基于隨機森林的高寒草地地上生物量高光譜估算

高宏元, 侯蒙京, 葛 靜, 包旭瑩, 李元春, 劉 潔, 馮琦勝, 梁天剛*, 賀金生,2, 錢大文

(1. 蘭州大學草地農業科技學院, 蘭州大學草地農業生態系統國家重點實驗室, 蘭州大學農業農村部草牧業創新重點實驗室, 蘭州大學草地農業教育部工程研究中心, 甘肅 蘭州 730020； 2. 北京大學城市與環境學院, 北京 100871；3. 中國科學院西北高原生物研究所, 青海 西寧 810008)

草地地上生物量(Above ground biomass,AGB)是指某一時刻單位面積內的草地植物在地上部分的有機物的總量[1]，它不僅是衡量草地群落生長狀況與生產力水平的關鍵參數，還是表征群落能量流動和物質循環的重要指標之一[2]。在高寒地區，草地是家畜的主要食物來源與能量供給[3]，草地AGB變化還能反映草地的退化與土壤的侵蝕程度[4]。因此，測定高寒草地AGB十分有必要。傳統的草地AGB測定通常采用實地收獲法，耗時費力，時效性較差，且對草地破壞性大。

高光譜遙感能在不直接接觸目標物體的同時，獲得其豐富的光譜信息[5]。對草地高光譜數據進行挖掘和分析，可以獲得植被的物理化學組分和生理生態情況等指標，這使實時監測草地植被成為了可能，如紀童等[6]用高光譜數據實現了草坪草葉綠素含量的監測，高金龍等[7]利用高光譜數據構建了高寒天然草地氮、磷養分的估測模型，王磊[8]用高光譜反演了草地的葉面積指數，韓萬強等[9]用高光譜數據識別了3種草地主要植物等。對于草地AGB的高光譜監測，目前也有不少相關研究。胥慧等[10]的研究表明基于光譜紅谷吸收深度D和光譜綠峰反射高度H的高光譜特征參數(D-H)/(D+H)和草地生物量有較高的相關性；夏浪等[11]使用NDVI反演生物量并得到較高的模型精度；馬維維[12]用高光譜衛星數據提取NDVI等5個植被指數用于反演草地生物量。這些研究大多直接以植被指數為建模特征，少有考慮原始光譜特征及其轉換參數對模型的影響。張凱等[13]篩選了甘南地區草地冠層的6個光譜特征變量并分別進行AGB的線性、對數等模型的構建；安海波等[14]用不同植被指數對內蒙古天然和人工草地進行了AGB指數、對數等的非線性回歸建模。上述研究構建的AGB反演模型大多為單因素或多因素參數統計模型，更注重數據的空間分布，這在研究樣本較少且草地類型均一時是有效的，當樣本較多且研究區情況復雜時模型的預測能力存在不穩定的問題[15]。因此，在進行草地AGB的高光譜研究時，建模特征的選擇和模型算法是影響AGB估算模型的關鍵因素，在模型簡化和穩定性方面有重要意義。

目前，用于高光譜研究的算法較多，連續投影算法(Successive projections algorithm,SPA)由于其良好的光譜冗余信息消除能力，常用于連續光譜數據的選擇，楊晨波等人[16]研究表明SPA可以極大減小原始光譜數據的維度，從而簡化模型；遞歸特征消除算法(Recursive feature elimination,RFE)是一種用于篩選最優特征子集的貪心算法，一般與機器學習結合使用，在降低數據維度的同時可找到精度最高的模型[17-18]；隨機森林(Random forest,RF)是一種決策樹集成算法，具有不易過擬合和普適性廣的優點[19]，在植被生理生態等指標的研究中，常用于估測模型的構建[20-22]。綜上，本研究基于SPA和RFE特征選擇方法和RF模型，開展高寒地區草地AGB的高光譜研究，以期對高寒地區的放牧強度、草畜平衡和生態環境實時監測提供科學依據。

1 材料與方法

1.1 試驗設計

本研究的試驗區位于青海海北高寒草地生態國家野外科學觀測研究站(簡稱海北站，37°37′ N，101°19′ E)附近(圖1)，坐落于青藏高原東北隅，草地類型屬于典型的高寒草地，年平均氣溫-1.7℃，年降水量范圍為426～860 mm，植被類型是以金露梅(Potentillafruticosa)為建群種的高寒灌叢草甸和以嵩草屬(Kobresia)植物為建群種的高寒嵩草草甸[23]。試驗區A，B是兩塊典型的高寒草地，每塊試驗地有15個小區(每個小區面積0.2 ha)，設有禁牧(CK)、輕度放牧(Light，L)、中度放牧(Medium，M)和重度放牧(Heavy，H)4個放牧梯度，對應的放牧強度分別為0頭·ha-1(CK)、0.5 頭·ha-1(L)、1.0 頭·ha-1(M)和2.0 頭·ha-1(H)。放牧家畜為牦牛，放牧方式為輪牧，時間為每年7—9月。放牧設置主要模擬了天然草地的復雜情況，能反映本研究構建的AGB模型的普遍適用性。

圖1 研究區位置圖Fig.1 Location of the study area

1.2 數據獲取與預處理

草地冠層光譜數據測量采用美國ASD公司的雙光譜儀系統，該系統由移動端和固定端兩臺ASD Field Spec 4光譜儀組成，其波段范圍為350～2 500 nm。雙光譜儀系統通過無線通訊，實時獲得反射率數據，并且能自動完成波長交叉校準，可消除不同時間太陽輻射變化帶來的誤差，因此在多變的天氣條件下也可以獲得良好的特征光譜圖[24]。測量草地冠層光譜要在晴朗少云的天氣下進行，測量時探頭垂直向下，距冠層高度大約1 m，每個小區隨機選取3個測量點，每個測量點測定9條光譜曲線，去掉異常光譜曲線，將其余光譜的平均值作為該小區的草地冠層反射率光譜。在小區的每個光譜測量位點放置樣方框，齊地面刈割樣方內的植物地上部分，帶回實驗室在65℃烘箱中烘48 h至恒重，得到每個小區的草地AGB。

由于光譜曲線有一定的噪聲，用Savitzky-Golay平滑法[25]進行平滑去噪處理，得到原始光譜曲線，進而得到一階微分光譜。此外，考慮到水分、氧氣和儀器敏感性波動等因素的影響，剔除了1 301～1 450 nm，1 801～2 050 nm和2 301～2 500 nm波段的光譜曲線[26]。本研究在2019年5月、6月、7月、8月、9月和2020年7月、9月總計開展了7次外業調查，在2020年的7月和9月，由于B放牧地封鎖無法獲得實測數據，因此共獲得180個樣本數據，去掉無效或錯誤數據，共有179個樣本數據，具體的樣本分布和描述性統計如表1所示。

表1 AGB數據描述性統計結果Table 1 Descriptive statistical result of AGB data

1.3 特征變量

本研究的特征變量有原始光譜(Original spectrum,OR)、一階微分光譜FD(First derivative spectrum,FD)、光譜位置面積參數(Spectral parameters of spectral position and area,PA)和植被指數(Vegetation indices,VI)4類(表2)，這些不同類別的特征變量，將用于后續的特征選擇和模型構建。

表2 特征變量及其定義Table 2 Variables and definition

1.4 特征選擇與建模

本研究先用連續投影算法SPA對原始光譜和一階導光譜數據進行特征波段反射率的提取，再用遞歸特征消除算法RFE對特征波段反射率和其他特征變量進行特征選擇，最后用隨機森林RF算法構建草地AGB的反演模型，以上算法均由Python編程語言實現。

SPA是一種使矢量空間共線性最小化的前向變量選擇算法，它的優勢在于可提取全波段的幾個特征波段，能夠消除光譜矩陣中冗余的信息[28]，因此本研究用SPA提取OR和FD的光譜特征波段，從而降低數據維度，使模型更簡單高效。

RFE是一種尋求最優特征子集的貪心算法[29]，基本思想是構建底層模型進行初始特征集的訓練，并給每個特征賦予權重，然后去掉權重最小的特征，將其他特征組成新的特征子集，再進行訓練，遞歸重復此過程直至達到最終所需的特征數目。在本研究中，RFE構建RF底層模型時選取默認參數，選取訓練結果決定系數R2最大時對應的特征子集作為構建模型的特征組合。

RF是多棵決策樹構成的集成模型，模型的最終輸出結果由森林中的每一棵決策樹共同決定，在RF中以每棵決策樹輸出的均值為最終結果。RF算法的具體過程[30]如下：在原始訓練特征集中用Bootstrap抽樣方法獲得n個特征子集；對每個特征子集選擇m個特征，并對每個訓練特征子集構建決策樹，得到n個決策樹模型，建立起隨機森林；計算每棵決策樹的結果，將n棵決策樹輸出結果的均值作為最終結果。此外，RF算法提供了特征重要性(總和為1)的接口，方便比較建模特征的重要性。經過參數優選和多次訓練，確定RF的主要參數決策樹數量為1 000。

1.5 模型評價

為了減少訓練樣本劃分偶然性帶來的結果誤差，本研究采用10折交叉驗證確保模型的穩定性。在10折交叉驗證中，試驗數據樣本被劃分為10份，輪流將其中9份作為訓練集，1份作為測試集，將測試集結果取平均值作為模型最終的評價結果。

本研究采用決定系數(Coefficient of determination,R2)和均方根誤差(Root mean square error,RMSE)評價AGB估測模型的精度，計算公式如下：

(1)

(2)

2 結果與分析

2.1 原始光譜和一階微分的SPA特征波段提取

圖2 SPA特征波段提取及波段分布Fig.2 SPA feature bands extraction and bands distribution

2.2 基于不同類型特征變量的AGB模型分析

對OR特征波段、FD特征波段、光譜位置面積參數PA和植被指數VI(表2)4類特征變量分別先進行RFE特征選擇，將RFE選擇出的特征子集作為RF的建模特征，再經過10折交叉驗證，計算得到訓練集和測試集的RMSE和R2(表3)，并利用RF自帶的特征重要性接口分析建模特征的重要性(圖3)。

圖3 不同RF模型的特征重要性Fig.3 The importance of features in different RF models

表3 不同類別特征變量的RF估測模型結果Table 3 Results of RF estimation model based on different class of features (n=179)

從模型精度來看，基于植被指數VI構建的RF模型精度最高(R2=0.70，RMSE=557.87 kg·ha-1)，對應7個建模特征為NDBleaf，OSAVI，TCARI，MTCI，PRI，SIPI，VARIg；光譜位置面積參數PA和一階光譜特征波段FD的模型精度次之，分別為R2=0.64(RMSE=596.42 kg·ha-1)和R2=0.57(RMSE=685.96 kg·ha-1)；精度最差的是原始光譜OR特征波段構建的RF模型，測試集R2僅為0.52(RMSE=700.06 kg·ha-1)。從RF模型構建的特征來看，經過RFE特征選擇后，4種模型的特征數量都有所降低，尤其是VI構建的RF模型，特征數量從26減少到7，在較大降低數據維度的同時，模型精度也最高，其中由R800和R450構建的結構不敏感色素指數SIPI的重要性最高(0.47)，其次是由R531和R570構建的光化學反射系數PRI(0.25)，這兩個特征為模型貢獻了72%的重要性。

對全集數據(n=179)進行實測生物量與估測生物量的相關性分析，結果表明，在4類特征變量中，基于植被指數VI的RF地上生物量反演模型的估測效果最好，其線性決定系數R2達到0.72，其次是光譜位置面積參數PA(R2=0.68)，再次是一階光譜特征波段FD(R2=0.59)(圖4)，模型估測效果最差的是原始光譜特征波段OR，其實測AGB和估測AGB的決定系數R2為0.56，這和測試集的10折交叉驗證評價結果相一致。

圖4 不同RF估測模型的實測AGB與估測AGB的相關性Fig.4 Correlation between measured AGB and predicted AGB in different RF estimation models

2.3 基于不同類型特征變量組合的AGB模型分析

進行不同類型特征變量組合時，考慮到OR特征波段和FD特征波段都屬于單波段特征變量，因此可以將其作為整體，與光譜位置面積參數PA和植被指數VI任意組合，進而進行RFE特征選擇和RF模型構建。表4結果顯示，在5種組合中，PA+VI組合的模型精度最高，R2達到了0.71，其對應的RMSE也最小，為548.97 kg·ha-1；精度次之的是所有變量構建的模型，擬合系數R2=0.70(RMSE=562.93 kg·ha-1)；OR+FD組合與OR+FD+PA組合的模型精度居中，R2分別為0.69和0.67(RMSE分別為579.00和591.17 kg·ha-1)；OR+FD+VI的模型精度最差(R2=0.64，RMSE=620.05 kg·ha-1)。從模型的建模特征及其重要性結果(圖5)來看，模型精度最高的PA+VI組合，是所有組合中特征數量最少的，僅有5個建模特征H,Dr,VI2,PRI和SIPI，其中特征重要性最高的是VI2和PRI，分別為0.37和0.30，二者為整個模型貢獻了67%的重要性；模型精度次之的全部變量組合，有6個建模特征，其中PRI和VI1的重要性最高，分別為0.30和0.26；OR+FD組合有7個建模特征，最重要的特征是R350，重要性為0.24；OR+FD+PA組合和OR+FD+VI組合的建模特征最多，均為10個，重要性最高的特征分別是VI2和PRI。

表4 不同類別特征變量組合的RF估測模型結果Table 4 Results of RF estimation model based on different combination of class of features (n=179)

利用全集數據(n=179)，對不同類型特征變量組合RF模型的實測生物量與估測生物量進行相關性分析(圖6)，可以得出，在不同的特征變量組合中，基于光譜位置面積參數PA和植被指數VI的組合PA+VI的模型AGB估測效果最好，其線性決定系數R2達到0.73，均方根誤差RMSE為476.93 kg·ha-1；估測效果次之的是OR+FD+PA+VI所有變量組合(R2=0.72)；OR+FD組合和OR+FD+PA組合的AGB模型估測效果一致(R2=0.68)；OR+FD+VI組合的RF模型的AGB估測效果最差，其實測AGB和估測AGB的決定系數R2為0.66。可見，不同類型特征變量組合模型的整體估測效果和其測試集10折交叉驗證的評價結果基本一致。

圖6 不同RF估測模型的實測AGB與估測AGB的相關性Fig.6 Correlation between measured AGB and predicted AGB in different RF estimation

3 討論

與傳統的統計方法相比，機器學習模型適用于較大規模數據和復雜場景，更關心模型的可用性與預測能力[15]。本文用RF構建的最優高寒草地AGB估測模型(PA+VI組合)的測試集精度為R2=0.71，全集的預測精度為R2=0.73，在不破壞草地的同時能夠快速測定AGB。

特征選擇是高光譜數據模型構建前的一個重要步驟，主要起數據降維和模型簡化的作用。本研究先后用SPA和RFE算法進行草地AGB建模變量的篩選，其中SPA算法主要用于連續光譜的特征波段提取。本研究表明，SPA將OR和FD的波段數量從1 551個分別減少到9個和19個(圖2)，可極大地降低數據的維度，但在特征波段提取時，SPA并未選擇RMSE最小時的波段數量，而是綜合考量波段數量和RMSE對模型的影響，選擇了適宜數量的特征波段，這和王承克[31]、吳迪[32]等的研究結論相一致。此外，OR的SPA特征波段主要集中在740，950，1 100和1 300 nm等附近，其中740 nm屬于紅邊范圍，950，1 100和1 300 nm屬于近紅外范圍。有研究表明紅邊包含了可以表征生物量、葉綠素等參量的光譜信息[33]，而920～1 120 nm之間的近紅外波段很有可能與葉片水分和干物質的吸收有關[34]，1 300 nm波段則和楊晨波等[16]對冬小麥的地上干生物量的敏感波段的提取結果相近。

不同類型特征變量的RF估測模型結果表明，PA和VI的建模精度好于OR，FD特征波段的建模精度，這可能與PA和VI是由多個波段組合計算的指數，而OR，FD特征波段是單波段有關，在一定程度上說明了PA和VI等波段組合指數比OR等單波段包含的光譜信息豐富。在不同類型特征變量組合的RF估測模型結果中，精度最高的是PA+VI組合，其僅有5個建模特征，是所有組合中建模特征數量最少的組合，其中特征VI2和PRI為整個模型提供了67%的重要性。VI2是紅谷吸收深度D和綠峰反射高度H的歸一化指數(表2)，能較好地反映草地生物量等參量，如胥慧等[10]的研究表明(D—H)/(D+H)和草地生物量的相關系數達到0.660。PRI指光化學反射系數，是R531和R570的歸一化指數[35]，能反映植物光合作用過程中光能利用效率，可用于研究植被生產力[36]，道日娜[37]研究發現，重度放牧時草地生物量最小，此時的PRI也最小，因此PRI也能較好的反映草地的AGB。

然而，本研究也存在一定的局限性。例如，機器學習是個典型的數據驅動的問題，通常需要大量的數據[15]，而本研究的樣本數為179條，這對機器學習而言數據量較小，可能對模型的精度造成一定的影響。放牧干擾行為也可能影響到模型最終結果：蘇日古嘎[38]研究發現，與未放牧草地相比，放牧會改變草地植被的群落物種多樣性、優勢度和植物種群的空間異質性，且不同放牧家畜對草地植被的影響程度不同；還有研究表明，由于長期圍欄放牧家畜的踐踏和選擇性取食，導致植被環境斑塊化[39]，以上行為都會改變草地植被冠層的結構組成，使獲取的高光譜反射率呈現不同的特征，進而對RF模型的精度造成一定的影響。此外，由于ASD光譜儀測定的是目標物“點”的光譜數據，是350～2 500 nm全光譜段(間隔1 nm)的觀測結果，由大量窄波段組成，所選最優模型變量也由其中部分特征波段組成，但目前可用的高時空分辨率的光學遙感衛星通常只有寬波段的觀測結果，后續將進一步研究最優模型變量的特征波段和遙感影像波段的對應關系，進而實現大面積草地AGB的模型構建和反演。無人機高光譜成像技術具有獲取數據快、操作靈活等特點，現被越來越多地應用到農作物等的監測，如陶惠林等[40]借助無人機高光譜遙感平臺，獲取了冬小麥各生育期的無人機影像，提取了植被指數和紅邊特征，最終構建了冬小麥產量模型。與ASD地物光譜儀相比，無人機高光譜成像技術的覆蓋面積更廣，因而可以構建較大尺度的生物量估算模型，這對高寒草地AGB的遙感估測具有參考意義。

4 結論

SPA可有效降低高光譜數據維度，使原始光譜OR和一階微分光譜FD的波段數量分別減少了99.4%和98.8%。在4種特征變量分別構建的草地AGB估測模型中，基于VI的RF模型精度最高(測試集R2=0.70，RMSE=557.87 kg·ha-1)，建模特征為NDBleaf，OSAVI，TCARI，MTCI，PRI，SIPI，VARIg，其中特征SIPI和PRI共為模型貢獻了72%的重要性，實測AGB和估測AGB的線性決定系數R2達到0.72。不同類型特征變量組合構建的草地AGB估測模型中，PA+VI組合的RF模型精度最高(R2=0.71，RMSE=548.97 kg·ha-1)，建模特征為H，Dr，VI2，PRI和SIPI，其中特征VI2和PRI共為模型貢獻了67%的重要性，實測AGB和估測AGB的線性決定系數R2達到0.73。整體而言，相對OR，FD等特征變量，PA，VI的特征變量組合構建的草地AGB的RF模型精度較高，建模效果更好。