運用Sentinel-2遙感影像數據估測森林蓄積量1)

李坤 吳達勝 方陸明 劉建超

(浙江省林業智能監測與信息技術研究重點實驗室(浙江農林大學)，杭州，311300) (浙江農林大學暨陽學院)

森林是陸地生態系統的重要組成部分，在全球碳循環中起著至關重要的作用，全面快速地掌握森林信息對國民經濟和綠色可持續發展有著積極作用。森林蓄積量是衡量森林質量的一個重要因子，也是反映森林資源水平的一個基本指標，精準估測森林蓄積量有助于監測森林質量，為林業資源管理決策提供科學依據。傳統的森林資源調查不僅費時費力，而且難以反映森林的動態發展情況。因此，林業工作者及研究人員一直致力于探索如何利用飛速發展的遙感技術對蓄積量進行估測[1-2]。

目前，遙感技術應用于森林蓄積量估測主要體現在三方面：一是選取光學遙感影像，對其進行主成分分析后選取主成分并提取紋理特征，探討不同分辨率不同窗口大小提取到的紋理特征與森林蓄積量的相關性；二是結合不同遙感數據源(光學遙感、雷達遙感)構建估測模型；三是使用不同建模方法，從回歸模型發展至人工神經網絡、k最近鄰算法等經典機器學習模型[3]。劉明艷等[4]利用主成分分析法對老禿頂子自然保護區進行蓄積量估測，模型預估精度達到92.18%；楊明星等[5]利用多窗口遙感影像紋理特征對蓄積量進行建模評估，使用隨機森林模型對蓄積量進行建模評估，精度達到75.46%，估測效果較好；蔣馥根等[6]基于方差優化K近鄰法對蓄積量進行反演，該模型決定系數(R2)達到0.69，相對均方根誤差(RRMSE)為32.04%。機器學習模型為非參數模型，相較于傳統的經驗模型，不僅能夠克服特征因子之間的共線性問題，而且具有更好的泛化能力，能夠對類別特征進行處理，更適用于處理數據量大、特征繁多的林業數據。

國內外學者通常運用當地的遙感影像和外業調查數據對當前區域的蓄積量進行預測或者反演[7-9]，很少將同一模型應用于不同區域，檢驗模型的普適性；在選取特征建模時，由于傳統經驗模型不能直接處理類別特征，需要根據類別特征將數據劃分為不同的類別，再根據不同類別分別建立不同的模型[10-11]，這樣的模型建立繁瑣，且適用性不強。

本研究從非參數模型泛化性能較強的特點出發[6，12]，結合機器學習算法處理類別特征的能力[13]，提出最小絕對收縮和選擇算子(Lasso)特征選擇方法的堆疊(Stacking)融合模型。根據Sentinel-2遙感影像、二類調查數據與數字高程模型(DEM)數據，首先選擇K最近鄰算法(K-NN)、梯度提升迭代決策樹(GBDT)、極端梯度提升(XGBoost)、梯度增強集成分類器(CatBoost)4種單一模型進行建模估測，然后探討類別特征的加入對模型指標的影響，最后建立堆疊(Stacking)融合模型。采用十折交叉驗證法檢驗模型性能指標，比較融合模型與單一模型的精度，驗證不同模型的普適性，探討類別特征對模型優化的影響，為森林蓄積量估測中特征選取以及模型構建提供新的思路。

1 研究區概況

淳安縣位于浙江省西部，地處北緯29°11′～30°2′，東經118°20′～119°20′，陸域面積4 417.48 km2，屬中亞熱帶季風氣候。全縣森林覆蓋率達65.00%，生物種類繁多、資源豐富，年平均溫度為17 ℃，年平均降水量1 430 mm，地形以丘陵山地為主，海拔50～1 500 m，土壤以紅壤、黃壤、巖性土、水稻土為主，其中紅壤為面積最大的土屬。

臨海市位于浙江省東南沿海，地處北緯28°40′～29°4′，東經120°49′～121°41′，陸地總面積2 203 km2，屬亞熱帶季風氣候，年平均氣溫17 ℃，年平均降水量1 550 mm，四季分明，適合植被生長，境內森林覆蓋率為64.20%。地形以丘陵和山地為主，山地面積占總面積的70.7%，平原面積占22.8%，呈“七山一水兩分田”地貌。

2 研究方法

2.1 遙感影像數據和地形數據

本研究使用Sentinel-2遙感衛星影像，成像時間為2017年11月27日，共12景，根據研究區行政區邊界裁剪拼接后的遙感影像如圖1(a)和圖1(c)所示。Sentinel-2號衛星(Sentinel-2A、Sentinel-2B)攜帶一枚多光譜成像儀(MSI)，可覆蓋13個光譜波段(見表1)。具有光學遙感衛星常見的可見光、近紅外和短波紅外波段，還包含3個應用在高光譜數據中的紅邊波段，具有不同的空間分辨率。Sentinel-2號數據是在紅邊范圍內的3個波段對監測植被的輻射強度非常有效[14-15]。

表1 Sentinel-2遙感影像波段信息

本研究中的地形數據來源于2017年研究區的二類調查數據，以及30 m分辨率的數字高程模型(DEM)數據。根據研究區行政邊界拼接裁剪后的數字高程模型數據如圖1(b)和圖1(d)所示。

2.2 特征因子提取

由于遙感影像、數字高程模型數據、二類調查數據坐標系不同，在數據整合前需進行坐標系轉換，本研究統一使用CGCS_2000_3_Degree_GK_CM_120E坐標系[16]。

光譜特征因子：Sentinel-2共包含13個波段，B1波段與B10波段與本研究無關，故選取剩余11個單波段因子。基于單波段因子進行波段運算，提取12個常見植被指數，即土壤調節植被指數、比值植被指數、非線性指數、歸一化水指數、歸一化植被指數、改進的歸一化植被指數、歸一化紅外指數、歸一化差異綠度指數、增強型植被指數、差值植被指數、類胡蘿卜素反射指數、綠光葉綠素指數。根據Sentinel-2特有的紅邊波段提取7個植被指數：比值紅邊植被指數、歸一化紅邊植被指數1、歸一化紅邊植被指數2(、改進的歸一化紅邊植被指數、紅邊葉綠素指數、非線性紅邊指數、反轉紅邊葉綠素指數[17]。植被指數的具體計算公式見表2。

紋理特征因子：本研究采用灰度共生矩陣法提取遙感影像的紋理特征[19]。以B2、B3、B4和B8等4個10 m分辨率波段組合后進行主成分分析，使用ENVI軟件提取第一主成分的8個紋理特征因子：均值、方差、協同性、對比度、相異性、熵、角二階矩、相關性。

表2 植被指數

地面數據特征因子：從二調數據中獲取每個小班的坡向、坡位、土層厚度、腐殖質厚度、植被覆蓋度、群落、優勢樹種、年齡、郁閉度等9個自變量因子，其中坡向、腐殖質厚度、群落、優勢樹種為類別特征因子。從數字高程模型數據中提取海拔、坡度和坡向等3個自變量因子。

本研究共提取50個特征因子(包含四個類別特征因子)。使用ArcMap軟件以森林資源細班為單元整合所有自變量因子，獲得淳安縣數據133 734條、臨海市數據59 636條。刪除字段缺失數據，再根據3倍標準差原則去除部分數據，共獲得淳安縣研究數據40 216條、臨海市研究數據19 840條。

2.3 基本模型

K最近鄰算法(K-NN算法)：K最近鄰算法是比較成熟的非參數方法中最典型的方法之一。該算法通過數據所有的特征計算新數據與數據集中數據點的距離(該距離通常采用歐式距離以及曼哈頓距離)，按照距離遞增次序進行排序，返回K個點的加權值作為預測值。該算法沒有過多的理論推導與步驟，實現較為方便。

圖1 研究區衛星影像圖和數字高程模型數據

梯度提升迭代決策樹(GBDT)：梯度提升迭代決策樹是一種集成學習算法，由決策樹和梯度提升組成，GBDT采用基函數的線性組合以及不斷減小訓練過程產生的殘差(預測值與真實值的差值)來完成數據回歸。其核心是計算上一輪迭代的負梯度，減少殘差的梯度方向上建立新的決策樹，通過不斷迭代改進來得到更準確的預測結果。該算法很好地降低了多種因素之間的信息冗余程度。

極端梯度提升算法(XGBoost算法)：極端梯度提升算法是一種基于決策樹的集成機器學習算法。XGBoost算法的基本思想是首先建立一個基分類/回歸器，逐步加入新的分類/回歸器，每加入一個分類/回歸器后，再次計算其目標函數的值，保證在迭代過程中，目標函數值逐步下降，以此不斷提升模型的表達效果。該算法有較強的泛化性能，能降低過擬合，適用性較強，在速度和效率上都優于其他梯度提升算法。

梯度增強集成分類器(CatBoost)：傳統的梯度增強對樣本梯度計算依賴于樣本自身，噪聲點會帶來預測偏移，最終導致過擬合，為解決這一問題，CatBoost首先對整個數據集進行若干次排序，去除第i條數據，針對前i-1條數據，分別計算損失函數以及梯度，并建立殘差樹，最后將殘差樹累加到原模型上，有效避免了預測偏移，減少過擬合。而且CatBoost算法能自動處理類別特征，并根據特征的內在聯系將原有類別特征進行組合，豐富特征維度，提升預測結果的精確性。

堆疊法(Stacking)模型融合：Stacking是一種集成學習方法，基于K折交叉驗證的思想，通過元分類器或元回歸器來融合多個分類或回歸模型。首先基于完整的訓練集對基礎模型進行訓練，然后將每個基礎模型的預測值作為元模型的訓練集，用新的訓練集對元模型進行訓練。Stacking通常會使用不同的機器學習算法來保證基礎模型的多樣性，并通過元模型以最優方法對初步回歸結果進行整合，可以充分利用基礎模型的預測結果，而且可以發現糾正基礎模型中的部分偏差，使得準確率得到提高，同時降低過擬合的風險。本研究使用CatBoost和GBDT作為基礎模型，XGBoost作為元模型。

圖2 Stacking模型融合示意圖

Lasso特征選擇：最小絕對收縮和選擇算子(Lasso)是一種基于L1范式的特征選擇方法，屬于嵌入式特征選擇。Lasso通過加入L1懲罰項，將一些不重要特征的系數壓縮為0，即剔除該特征，從而達到稀疏化和特征選擇的目的。與其他變量選擇方法比較，Lasso方法能夠較準確地選擇出重要特征，而且還具備特征選擇的穩定性。

類別特征編碼：本研究共選取四個類別特征。由于K-NN、GBDT和XGBoost不直接支持類別特征，因此需要對類別特征進行編碼。根據類別特征屬性，本研究采用兩種類別編碼方式。

標簽編碼：使用腐殖質厚度、群落為定序變量，采用標簽編碼，將類別特征依照特質以數字劃分。將腐殖質厚度的{“薄”，“中”，“厚”}、群落的{“簡單”，“較完整”，“完整”}轉換為{0，1，2}。

獨熱編碼：使用坡向取值為{“陰坡”，“陽坡”，“無坡向”}、優勢樹種取值為{“杉木”，“馬尾松”，“闊葉混交林”等8種優勢樹種}為定類變量，采用獨熱編碼，其方法是使用N位狀態寄存器來對N個狀態進行編碼，每個狀態都有獨立的寄存器位，在任意時候，其中只有一位有效。坡向編碼結果如圖3所示。

圖3 獨熱編碼圖

2.4 模型構建

因本文研究數據涉及到的自變量因子有46個(不包括4個類別特征)，為了找出對森林蓄積量具有更大的影響作用的因子，進一步提升模型的估測精度及其泛化能力。使用Lasso特征選擇方法，從繁多的特征因子中選取部分重要特征，建立K-NN、GBDT、XGBoost和CatBoost等4種單一模型，并使用控制變量法探索類別特征對模型指標的影響，然后使用Stacking方法對單一模型進行融合，分析比較融合模型與單一模型的優劣。本文選用的Lasso特征選擇以及機器學習模型均通過網格搜索(GridSearchCV)進行參數調優。

建模的具體步驟如下：(1)將所有樣本按照7∶3的比例隨機劃分訓練集與測試集；(2)通過網格搜索對4個單一模型進行反復調參，建立最優參數估測模型；(3)使用最優參數模型對測試集進行估測，得到單一模型預測結果；(4)檢驗類別特征對模型性能指標的影響，分析模型的普適性；(5)使用Stacking方法融合多個單一模型，并與單一模型估測結果進行對比。

2.5 模型精度評價

本研究選用十折交叉驗證的方式來評估模型的精度。采用的性能指標有決定系數(R2)、均方誤差(MSE)、平均絕對誤差(MAE)、平均百分比誤差(MAPE)，計算公式如下：

由于使用K-NN算法前需要對數據進行歸一化，而樹模型并不需要對數據進行歸一化，從而導致數據量綱不統一。使用不同量綱下的MAE和MSE不能體現模型之間的差異，故使用均方誤差與因變量平均值的比值(相對均方誤差)替代性能指標的均方誤差，平均絕對誤差與因變量平均值的比值(相對平均絕對誤差)替代性能指標的平均絕對誤差。

3 結果與分析

3.1 數值特征選擇

本文共提取數值特征46個，類別特征4個。為了減少信息冗余，降低共線性，使用Lasso方法對46個數值特征進行特征選擇。經參數調優后，參數阿爾法取值0.003，系數閾值0.05，并選取系數值大于0.05的自變量因子作為建模特征因子(包括郁閉度、比值紅邊植被指數、改進的歸一化植被指數、年齡、方差、相關性、相異性和熵等特征因子)。

3.2 類別特征對單一模型的影響

為比較類別特征(坡向、腐殖質厚度、群落和優勢樹種)對于模型性能指標的影響，使用淳安縣數據，在未加入類別特征的情況下，經篩選后的8個數值特征進行建模和運用4個類別特征及篩選后的8個數值特征進行建模，模型性能指標見表3。

由表3可知，類別特征的加入對模型各項性能指標的提升有著顯著影響。在加入類別指標后，GBDT模型的R2提高了0.12，CatBoost模型的R2提高了0.16；GBDT模型的相對均方誤差從55.55%降低到32.19%，CatBoost模型的相對均方誤差從54.0%降低到26.7%；GBDT模型的相對平均絕對誤差從28.32%降低到20.17%，K-NN模型的相對平均絕對誤差從34.02%降低到24.47%；K-NN模型的平均百分比誤差從36.74%降低到33.75%，CatBoost模型的平均百分比誤差從36.07%降低到22.57%。因此，加入類別特征后，4種模型的各項性能指標均有顯著改善，CatBoost模型總體表現最好，各項指標優于其他3種模型。加入類別特征后的4種模型估測結果散點圖如圖4所示。

表3 淳安縣森林蓄積量估測模型的性能指標

由圖4可知，CatBoost、XGBoost、GBDT、K-NN等4種模型的R2逐漸下降，預測點與真實值之間差值越來越大，觀察點在回歸直線附近愈發稀疏，模型的可解釋性降低。

3.3 單一模型的普適性

為探究模型跨空間普適能力，并進一步驗證類別特征的加入對模型性能指標的影響是否具有普適性。將4種模型應用于臨海市森林蓄積量估測，估測結果性能指標如表4所示。該結果進一步驗證了類別特征對模型的優化作用具有一定的普適性，且XGBoost和CatBoost模型優于K-NN以及GBDT模型。

圖4 4種模型對淳安縣森林蓄積量預測的擬合

表4 臨海市森林蓄積量估測模型的性能指標

雖然4種模型的R2相較淳安縣的森林蓄積量估測結果有所下降，但是相對均方誤差，相對平均絕對誤差以及平均百分比誤差等3項指標均好于淳安的估測結果，說明模型具有一定的跨空間普適性，在不同地區均能取得較好表現。加入類別特征后的蓄積量估測散點圖如圖5所示。

3.4 融合模型

為驗證融合模型與單模型的優劣，基于前述篩選的類別特征及數值特征，以CatBoost和GBDT為基礎模型，XGBoost為元模型，進行Stacking模型融合，與相同情況下單模型中表現最好的CatBoost進行對比(見表5)。由表5可知，模型融合后，除相對均方誤差外，其余性能指標較CatBoost均有所改善，說明融合模型的綜合效果優于單一模型。

圖5 4種模型對臨海市森林蓄積量預測的擬合

4 結論與討論

遙感影像包含了更廣域的、動態的地面覆蓋物信息，為森林蓄積量估測提供了范圍更廣、實時性更強的數據來源，結合DEM、森林資源二類調查數據，為森林資源動態監測提供更為有效的支持。本研究使用2017年淳安縣及臨海市的Sentinel-2遙感影像、二類調查以及DEM數據，經過數據預處理以及特征選擇，建立了單模型和多模型融合的蓄積量預測模型。通過控制變量法探索類別特征對于模型指標的影響，證明類別特征的加入能有效提高模型性能；實驗中，使用臨海市森林資源數，進一步據驗證模型的間普適性以及類別特征對模型精度的影響；Stacking融合模型與單模型相比，綜合指性能標更佳，估測精度更高。

在所選取的特征因子中，獲取數值特征46個，Lasso特征選擇方法有效減少了特征數量，使建模所用數值特征縮減到8個，加快了模型訓練速度、提高了模型泛化能力，并使模型具有更好的解釋性；類別特征的加入明顯優化了模型，相比于常見的根據類別特征將數據劃分為不同類別，建立模型方法更加簡便易用；單模型中，CatBoost綜合表現最佳，但融合模型Stacking的綜合性能指標更優于CatBoost模型；Sentinel-2遙感影像特有的三個紅邊波段對于森林蓄積量的估測結果有較大的影響，與前人研究結果一致[19-20]。

但受遙感影像精度及與二類調查數據獲取時間差異等資源的限制，本研究仍存在一些不足，需要進行優化改進：研究所選取的研究區域空間跨度較大，但是仍未超出省級界限，可選擇空間跨度更大的地區作為研究區域，進一步分析模型的普適性；使用了Sentinel-2的單波段與基于部分波段提取的紋理特征因子，但是缺少對紋理特征因子的深入探究，可使用不同窗口大小、不同步長、不同波段組合來提取紋理特征，探究不同紋理特征對模型精度的影響[21]；本研究使用遙感影像分辨率為10 m，若能獲取更高分辨率的衛星遙感影像和DEM數據進行研究，則可以進一步提高預測精度；除光學遙感外，雷達遙感影像也常被用來作為蓄積量估測的數據源，雷達遙感影像可穿透森林冠層，獲取植被垂直結構參數，可應用于蓄積量預測[19]，進一步探討雷達影像對蓄積量估測模型精度的影響。