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

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

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

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

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

1 材料與方法

1.1 試驗設(shè)計

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

圖1 研究區(qū)位置圖Fig.1 Location of the study area

1.2 數(shù)據(jù)獲取與預(yù)處理

草地冠層光譜數(shù)據(jù)測量采用美國ASD公司的雙光譜儀系統(tǒng)，該系統(tǒng)由移動端和固定端兩臺ASD Field Spec 4光譜儀組成，其波段范圍為350～2 500 nm。雙光譜儀系統(tǒng)通過無線通訊，實時獲得反射率數(shù)據(jù)，并且能自動完成波長交叉校準，可消除不同時間太陽輻射變化帶來的誤差，因此在多變的天氣條件下也可以獲得良好的特征光譜圖[24]。測量草地冠層光譜要在晴朗少云的天氣下進行，測量時探頭垂直向下，距冠層高度大約1 m，每個小區(qū)隨機選取3個測量點，每個測量點測定9條光譜曲線，去掉異常光譜曲線，將其余光譜的平均值作為該小區(qū)的草地冠層反射率光譜。在小區(qū)的每個光譜測量位點放置樣方框，齊地面刈割樣方內(nèi)的植物地上部分，帶回實驗室在65℃烘箱中烘48 h至恒重，得到每個小區(qū)的草地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次外業(yè)調(diào)查，在2020年的7月和9月，由于B放牧地封鎖無法獲得實測數(shù)據(jù)，因此共獲得180個樣本數(shù)據(jù)，去掉無效或錯誤數(shù)據(jù)，共有179個樣本數(shù)據(jù)，具體的樣本分布和描述性統(tǒng)計如表1所示。

表1 AGB數(shù)據(jù)描述性統(tǒng)計結(jié)果Table 1 Descriptive statistical result of AGB data

1.3 特征變量

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

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

1.4 特征選擇與建模

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

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

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

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

1.5 模型評價

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

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

(1)

(2)

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

進行不同類型特征變量組合時，考慮到OR特征波段和FD特征波段都屬于單波段特征變量，因此可以將其作為整體，與光譜位置面積參數(shù)PA和植被指數(shù)VI任意組合，進而進行RFE特征選擇和RF模型構(gòu)建。表4結(jié)果顯示，在5種組合中，PA+VI組合的模型精度最高，R2達到了0.71，其對應(yīng)的RMSE也最小，為548.97 kg·ha-1；精度次之的是所有變量構(gòu)建的模型，擬合系數(shù)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)。從模型的建模特征及其重要性結(jié)果(圖5)來看，模型精度最高的PA+VI組合，是所有組合中特征數(shù)量最少的，僅有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估測模型結(jié)果Table 4 Results of RF estimation model based on different combination of class of features (n=179)

利用全集數(shù)據(jù)(n=179)，對不同類型特征變量組合RF模型的實測生物量與估測生物量進行相關(guān)性分析(圖6)，可以得出，在不同的特征變量組合中，基于光譜位置面積參數(shù)PA和植被指數(shù)VI的組合PA+VI的模型AGB估測效果最好，其線性決定系數(shù)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的決定系數(shù)R2為0.66。可見，不同類型特征變量組合模型的整體估測效果和其測試集10折交叉驗證的評價結(jié)果基本一致。

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

3 討論

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

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

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

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

4 結(jié)論

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