基于Landsat 8 和隨機森林的青海門源天然草地地上生物量遙感估算

趙翊含，侯蒙京，馮琦勝，高宏元，梁天剛*，賀金生，2，錢大文

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

草地生物量是反映草地生長狀況和生態環境評估的一項重要指標［1］，準確地估測草地產量一直是畜牧業科研工作者持續探索的目標［2］。在高寒地區，由于過度放牧等因素引起的草地退化、水土流失等環境惡化問題日益嚴重；并且生態系統受干擾后植物生長和恢復速度緩慢［3］。為了確保高寒地區畜牧業的可持續發展，對于準確估算高寒草地地上生物量十分必要［4］。

傳統地面草地調查的方法耗費成本太高［5］，刈割法對草地的破壞性大，不利于草地的再生與可持續。基于遙感技術對草地資源監測省時、省力，能快速且客觀地對草地生長現狀做出評估，適合大尺度估測草地地上生物量（above-ground biomass，AGB）。自20 世紀80 年代，有很多國內外學者利用多源遙感數據開展AGB 估算研究。如姚興成等［6］利用中分辨率成像光譜儀（moderate-resolution imaging spectroradiometer，MODIS）和地面實測資料，以云南省為例建立草地AGB 遙感估算模型，用草地植被群落特征優化了模型精度；于惠等［7］構建古浪縣生物量估算模型，結果表明基于MODIS 短波紅外波段的歸一化耕作指數（normalized difference tillage index，NDTI）與草地AGB 的相關性最好；馮琦勝等［8］基于青藏高原地面實測樣本數據與MODIS-NDVI 構建生物量風干重估算模型，研究表明MODIS-NDVI 的指數函數模型是最優模型。盡管利用MODIS 遙感數據可以估測草地生物量，但是中分辨率影像存在空間分辨率較粗（大于250 m）和混合像元等問題，致使草地AGB 模型的精度不高，不適合縣域尺度高精度精細化的遙感監測。

Landsat-8 OLI（operational land imager，OLI）于2013 年發射，空間分辨率達30 m，與MODIS 相比其空間分辨率有實質性的提升，并且包含了9 個波段，其中包括一個15 m 的全色波段，適合對草地AGB 進行更精細的遙感估算。Landsat-8 OLI 改進了波段數量、光譜范圍及輻射分辨率等方面，在性能上相比之前幾代Landsat 衛星擁有明顯改善，對草地AGB 估測的準確性有提升作用［9］。李斌等［10］利用Landsat 8 遙感影像對不同處理階段［數字量化值（digital number，DN）、輻射定標、大氣糾正］的歸一化植被指數（normalized difference vegetation index，NDVI）進行對比研究，表明應該采用經過大氣校正后的值作為植被覆蓋的定量指標；Li 等［11］利用Landsat-8 增強型植被指數（enhanced vegetation index，EVI）建立了與植物干物質含量的回歸模型，評價分析了整個青藏高原草甸和草原群落地上生物量的變化特點，并比較了植物性狀。烏如汗［12］以內蒙古正藍旗為試驗區，利用Landsat-8 OLI 遙感影像和地面實測數據，確定了NDVI 植被指數的一元二次回歸模型作為草地地上生物量反演模型。以上研究表明Landsat-8 可以實現對草地AGB 的有效估測。然而，現有研究多采用傳統參數化回歸方法構建模型，反演誤差相對較大，因此在提高草地AGB 監測精度上仍然有必要進一步開展深入研究。近年來，機器學習方法在植被遙感監測方面取得了良好的效果［13］。它通常基于完整的光譜集進行建模，能夠充分利用光譜信息，具有非線性、預測準確率高、泛化能力強的特點［14］，因此在遙感領域基于非參數化的機器學習算法估測AGB 逐漸受到重視［13］。

基于以上因素考慮，本研究利用Landsat-8 OLI 遙感數據結合青海省門源縣的地面實測數據，將計算出的植被指數作為自變量，構建傳統回歸模型和隨機森林模型，對比分析模型精度，確定遙感反演生物量的最佳模型，以期為當地草地資源可持續利用和草畜科學管理提供理論依據。

1 材料與方法

1.1 研究區概況

研究區位于青海省海北藏族自治州門源回族自治縣，地處青藏高原腹地，地理位置為北緯37.5°-38.0°，東經102.0°-102.5°（圖1）。門源回族自治縣地勢西北高，東南低，平均海拔2800 m 以上，屬于高原大陸性氣候，夏季涼爽短暫，冬季寒冷漫長，年平均氣溫-1.7 ℃，年平均降水500 mm。該地區天然草地類型主要為高寒草甸類，主要以莎草科牧草嵩草（Kobresia myosuroides）、黑褐苔草（Carex atrofusca）和禾本科牧草為優勢種［15］。門源縣是河湟地區和河西走廊重要的水源涵養區和補給地，基于農牧資源豐富的區位優勢，積極貫徹“立草為業、畜牧強縣”的發展方針，大力發展高寒地區現代高效畜牧業及很多特色種養業。

1.2 數據來源

1.2.1實測數據獲取 本研究共獲得202 個地面實測調查樣本數據，包括青海海北高寒草地生態系統國家野外科學觀測研究站（簡稱海北站）多年連續的觀測數據和全縣范圍內開展的草地AGB 外業調查數據（圖1）。觀測站數據獲取時間為2019 年5-8 月、2020 年7 月和2021 年7-8 月，共得到178 個樣本數據，全縣范圍內AGB 外業調查共獲得24 個樣點數據，調查時間為2021 年7 月1 日-7 月8 日。在對草地生長期內的AGB 進行采集時，每個樣點選取3 個隨機樣方，樣方大小為50 cm×50 cm，齊地剪取地上生物量，并記錄各樣方測量時間、經緯度、海拔、草地蓋度、草層高度、鮮重，將樣本裝袋編號帶回實驗室進行65 ℃烘干48 h 后稱重。最終，取3 個樣方的平均值作為樣本生物量，總計有202 個觀測樣本。

圖1 研究區海拔及采樣點空間分布Fig.1 Spatial distribution of altitude and sampling points in the study area

1.2.2遙感數據獲取 遙感影像數據選用Landsat-8 OLI 數據，是美國陸地衛星計劃的第8 顆衛星，傳感器類型為陸地成像儀（OLI），時間分辨率為16 d，空間分辨率為30 m。下載的陸地產品級別為二級產品（Collection 2 Level-2），該產品已經經過輻射定標、大氣校正和幾何精校正。為了使影像數據的成像時間與地面實測數據的調查時間差異較小，下載影像時選擇影像完全覆蓋樣地，成像時間與采樣時間最接近的影像。選取2019-2021 年的7-8 月覆蓋門源縣采樣點及全境、云量覆蓋<30%的Landsat-8 OLI 影像，總計7 景（表1）。研究區實測草地生物量樣本數據有202 個，剔除受天氣狀況影響的樣本，實際參與模型構建的樣本數為199 個。

表1 所用Landsat-8 OLI 影像信息Table 1 Landsat-8 OLI image information used

1.3 模型構建及精度驗證

1.3.1單變量回歸模型建立 利用ArcGIS 10.2 軟件結合樣地經緯度坐標，導出Landsat-8 OLI 波段2~7，并利用Extraction 工具提取與199 個地面采樣點對應的像元值。根據公式計算29 種植被指數（表2）。分別構建29種植被指數與AGB 的單因素回歸模型，將植被指數作為自變量，AGB 作為因變量，二者之間建立線性、對數、指數、乘冪4 種擬合模型；采用10-fold 交叉驗證的方式對模型進行精度評價。

表2 植被指數計算公式Table 2 Calculation formula of vegetation indexes

1.3.2變量篩選 為了篩選出合適的建模因子，本研究在RStudio 中利用“glmnet”程序包實現基于最小絕對收縮篩選方法（the least absolute shrinkage and selection operator，LASSO）的變量選擇。這種方法可以將變量的數量顯著減少，從而實現降維，被廣泛應用于高維數據中選擇敏感變量［16］。LASSO 變量篩選方法以模型系數的絕對值函數作為懲罰項施加在普通最小二乘法的損失函數上，以此來壓縮模型系數，使一些解釋因變量效果不好的系數變小甚至壓縮至0 從而移出模型，因此LASSO 方法可以提供一種稀疏解，能夠同時進行變量篩選和參數估計［17-18］。

1.3.3隨機森林模型構建 隨機森林算法（random forest，RF）是一種新型機器學習算法，是運用多棵決策樹訓練樣本并集成預測的一種非參數機器學習算法，是從原始樣本中，利用bootstrap 重抽樣技術隨機抽取數據后構造多個樣本，然后對每個重抽樣樣本采用節點的隨機分裂技術構造N 棵決策樹［19］。在每棵樹生長過程中，從全部特征變量中隨機抽選mtry 個進行內部節點劃分；最后，將N 棵決策樹的預測結果集合起來，采用投票的方式決定新樣本的類別［20］。ntree 參數值上限一般設置為1000，大量研究已證明該值對許多RF 程序有效［21］。利用“random forest”程序包實現隨機森林算法，需要定義2 個參數：決策分類樹的數目（ntree）和節點分割的特征數目（mtry），參數影響著模型的精度，本研究采用的參數決策樹數量ntree 是500、特征個數mtry 是4。

1.3.4精度驗證 為了比較分析傳統回歸模型和RF 模型對AGB 的反演效果，這兩種模型均采用十折交叉驗證的方法［16］，將所用的樣本數據隨機分成10 份，其中9 份作為訓練集，1 份作為驗證集，分別計算均方根誤差（root mean square error，RMSE）、決定系數（coefficient of determination，R2）評價AGB 估測模型的精度［22］。計算公式如下：

式中：Yi表示樣本實測值，Yi表示模型反演值，yi表示樣本生物量實測平均值，n表示樣本數量；RMSE 越小，表示擬合精度越好；R2越接近1，擬合程度越好，參考價值越高［1］。

2 結果與分析

2.1 生物量樣本數據分析

2019-2021 年門源縣用于建模的199 個實測樣本的草地生物量統計結果見表3。從該表可以看出，海北觀測站 生物量最大值為4142 kg·hm-2，最小值為187.2 kg·hm-2，平均值 為2061 kg·hm-2，變異系 數（coefficient of variation，CV）為0.47。外業調查數據覆蓋門源縣主要草地類型，生物量最大值為4666.67 kg·hm-2，最小值為637.07 kg·hm-2，平均值為1824.01 kg·hm-2，變異系數（CV）為0.56。

表3 觀測樣本生物量匯總Table 3 Biomass statistics of observed samples

2.2 單因素植被指數回歸模型

續圖2 單因素植被指數回歸模型預測值與實測值關系Continued Fig.2 Relationship between predicted value and measured value of single-factor vegetation indexes regression model（n=199）

從結果來看，NDVI 與AGB 的對數模型精度達到最優，驗證集R2達0.50，RMSE 為702.89 kg·hm-2。單因素植被指數構建的模型的驗證集R2介于0.37~0.50，RMSE 介于702.89~792.95 kg·hm-2。所有植被指數最優模型 統 計 分 析 的 結 果 表 明，NDVI、EVI、RVI、GNDVI、DVI、RDVI、SAVI、MSAVI、OSAVI、PVI、TVI、TVI2、MSR、GSAVI、GOSAVI、TSAVI、RGBVI、VDVI、GRVI、GCI、RBNDVI、BNDVI 最優擬合模型為對數，其他植被指數的最優擬合模型為線性。圖2 為所有植被指數最優模型預測值和實測值的關系。驗證集R2較高的植被指數前10 個分別是NDVI、RBNDVI、TVI、GNDVI、MSR、TSAVI、GCI、GRVI、RVI、panNDVI。

圖2 單因素植被指數回歸模型預測值與實測值關系Fig.2 Relationship between predicted value and measured value of single-factor vegetation indexes regression model（n=199）

2.3 隨機森林模型分析

本研究對比了RF 模型變量篩選前后的模型精度。以29 種植被指數建立的模型，其驗證集R2為0.61，RMSE為621.14 kg·hm-2（圖3a）。利用LASSO 變量篩選方法，采用十折交叉驗證的方法確定模型均方誤差最小時λ0=0.25，并取距離λ0一個標準差位置的λ1（1.22）作為最終的懲罰系數（圖4）。從29 種植被指數中選出11 種植被指數，其重要性排序分別是MNLI、MSAVI、GRNDVI、MVI、SAVI、GOSAVI、GRVI、RVI、RGBVI、GSAVI、GCI（圖5）。通過變量篩選所構建的AGB 估測模型為最優，驗證集R2為0.62，RMSE 為621.95 kg·hm-2（圖3b）。由此可見，在保證精度的基礎上LASSO 變量篩選對草地AGB 反演模型可以有效進行降維和簡化。以上2 種RF 模型的精度均高于單因素回歸模型，其中最優RF 模型的R2比最優單因素模型提高0.12，RMSE 降低了80.95 kg·hm-2。

圖3 隨機森林模型預測值與實測值關系Fig.3 Relationship between predicted value and measured value of RF model

圖4 LASSO 篩選結果Fig.4 LASSO results

圖5 篩選出的植被指數的重要性Fig.5 The importance of screened vegetation indexes

2.4 門源縣AGB 空間分布

經過變量篩選的隨機森林草地AGB 估算模型最優，利用該模型對研究區2019-2021 年的AGB 進行了反演（圖6）。可以看出，草地AGB 高值主要集中在西北部，東南部相對較低，總體呈中部高，四周低的趨勢。2019 年草 地AGB 值 主 要 在0~800 kg·hm-2和1500~2000 kg·hm-2；2020 年 主 要 在0~800 kg·hm-2和1000~1500 kg·hm-2；2021 年主要分布在1000~1500 kg·hm-2和1500~2000 kg·hm-2。2019-2021 年全縣天然草地總產草量介于4.2827 萬~8.9776 萬t，平均單產介于1063.49~1484.82 kg·hm-2。門源縣草地類型以山地草甸、溫性草原和高寒草甸三類為主，其中分布最廣泛的是高寒草甸，2019-2021 年產草量介于4.0825 萬~5.6653 萬t，平均AGB為1060.38~1471.94 kg·hm-2。 山 地 草 甸2019-2021 年 產 草 量 介 于973.7687~1571.6790 t，平 均AGB 介 于1036.81~1637.43 kg·hm-2；溫 性 草 原2019-2021 年產草量介于746.6281~1112.7140 t，平均AGB 介于1198.72~1786.63 kg·hm-2（圖7）。

圖6 門源縣2019（a）、2020（b）和2021（c）年天然AGB 空間分布Fig.6 Spatial distribution of AGB of natural grassland in 2019（a），2020（b）and 2021（c）in Menyuan County

圖7 門源縣3 種主要草地類型的平均地上生物量統計Fig. 7 Average aboveground biomass of three main grassland types in Menyuan County

3 討論

現階段國內外研究基于植被指數構建生物量估算模型的方法較為常見，其數值可以代表植被活力，比單波段更具有靈敏性［23］。NDVI 是目前在遙感領域應用最廣泛的一種植被指數，可以捕獲植被整個生長季的動態，可用于植物生物量的預測，在各類植被指數中，NDVI 可以較好反映AGB，對植被生長的相關信息敏感，如水分含量和植被覆蓋等［24］。從本研究對29 種植被指數單因素建模的研究結果來看，NDVI 用于草地AGB 的估測模型在傳統統計回歸模型中精度最優，這與許多學者［25-27］的研究結果一致。也有很多學者研究認為SAVI 比NDVI 更適合估算草地地上生物量，因為SAVI 適應植被密度變化及消除土壤影響的能力較強［28-29］。但是，本研究區分布范圍最大的草地類型是高寒草甸，其植被覆蓋度較高且均一，所以土壤背景對植被估測的影響不如其他學者的研究區域明顯［4］，并且這是單因素估算模型，其反演結果不如NDVI 與AGB 估測效果，這與楊鵬萬等［4］的研究結果類似。

除NDVI 以外，單因 素 回 歸 模 型R2較高的前10 個植被指數依次分 別 是NDVI、RBNDVI、TVI、GNDVI、MSR、TSAVI、GCI、GRVI、RVI、panNDVI。經過LASSO 變量篩選后構建研究區草地AGB 反演模型的植被指數有RVI、SAVI、MSAVI、MNLI、GSAVI、GOSAVI、RGBVI、GRVI、GCI、MVI、GRNDVI。由此可見，傳統回歸結果與RF 建模篩選出的變量有一定差異。值得一提的是，在29 種單因素植被指數回歸模型中，NDVI 與AGB 的相關性最高，但是在隨機森林模型構建篩選變量的過程中NDVI 沒有被篩選出來，篩選出來的相近植被指數是GRNDVI。這是因為GRNDVI 是基于NDVI 基礎上改進的植被指數，GRNDVI 與NDVI 具有較好的線性相關關系，利用這兩種植被指數估算葉綠素可以達到相似的精度［30］。GRNDVI、RVI、GRVI、GCI 在變量篩選中被篩選出來，同時也是單因素回歸模型中與生物量擬合精度較好的幾種植被指數。但在RF 模型中其他被篩選出的植被指數（SAVI、MSAVI、MNLI、GSAVI、GOSAVI、RGBVI、MVI），其單因素模型的精度并不高。多種植被指數性能具有交互影響和互補作用可能是導致這種現象的重要原因。兩種模型采用的是完全不同的變量篩選方法，單因素模型是植被指數各自分別與AGB 的相關性，不能涵蓋多種因素的綜合影響，基于單因素構建的回歸模型的誤差偏大；而在RF 模型中是所有篩選出來的植被指數共同影響并決定了模型的精度，明顯優于傳統回歸模型。

比較分析最優RF 模型涉及的多種植被指數，可以看出：1）基于近紅外和紅光波段組合計算的植被指數，是反演草地AGB 的重要變量。本研究由LASSO 篩選出的11 個植被指數作為建模因子構建了最優RF 模型，其中，有4 個植被指數是基于Landsat-8 近紅外和紅光波段計算的，分別為RVI、SAVI、MSAVI、MNLI。已有研究表明，利用紅邊和近紅外光譜計算的植被指數可以最大限度地減少大氣和水分吸收等因素的影響［31］，研究區具有高寒陰濕的氣候特點，這可能是這類植被指數對高寒地區草地AGB 有顯著影響的關鍵因素之一。2）由綠波段參與構建的植被指數對高寒草地AGB 也具有重要影響。大部分植被指數多基于近紅外和紅光波段反射率組合計算，然而近年來研究發現以近紅外光和綠光波段構建的植被指數有時更加敏感［32-33］。本研究篩選出的GSAVI、GOSAVI、RGBVI、GRVI、GCI、GRNDVI 這幾種植被指數的計算均有綠波段，這表明基于綠波段構建的植被指數對高寒地區的草地生物量的精準估算也具有顯著影響。3）MVI 和MNLI 這兩個植被指數對干旱區的植被更為敏感，這可能與本研究區地處青藏高原腹地，屬半干旱地區密切相關。首先短波紅外波段參與構建的植被指數對特定區域的草地植被的AGB 更為敏感。MVI 植被指數是基于紅光、近紅外和短波紅外波段計算的植被指數，于惠等［7］研究表明短波紅外波段對干旱區草地植被更為敏感；此外，MNLI 是非線性植被指數的改進指數，它與葉綠素含量具有較高的相關性，在半干旱環境下，MNLI 可以將植被指數與地表生物物理參數的非線性化關系線性化［34］。故MVI 和MNLI 被篩選出來，對處于半干旱區的植被生物量估算也具有一定的作用。綜上所述，本研究通過變量篩選后，用選出的多種類型的植被指數建模，這不僅簡化了模型，也從葉綠素含量、消除水分影響、植被密度等方面綜合反映了植被特征，多種植被指數與草地AGB 之間具有顯著的相關性，并且不同植被指數之間也具有一定的互補性，從而可以綜合地反映草地生物量的狀況。

相較于傳統回歸分析方法，機器學習算法更適用于較復雜的運算，可以更好地進行變量篩選和組合，很大程度地提升草地AGB 估測模型的精度。隨機森林模型可以組合不同含義的變量特征，且有效解決“過飽和”和共線性的問題［35］。有很多學者利用隨機森林模型進行了農作物、森林、竹林的生物量估算和茶園提取［36-37］，均取得了不錯的研究成果。但是，RF 模型也有其局限性，特別是在構建回歸分類樹的方式上，會低估超出訓練集范圍內的高生物量值［32］。此外，機器學習基于數據驅動，通常需要大量的數據，若研究的數據量較小，可能對模型準確性造成一定的影響［38］。

Landsat-8 OLI 遙感數據因其相較于低分辨率衛星具有更高的分辨率，相較于之前的Landsat 衛星具有更多的波段和覆蓋范圍，有避免大氣吸收特征的干擾，可用于海岸帶觀測和云檢測，近紅外波段、短波紅外波段與MODIS 對應波段接近等優勢，在我國草地資源的監測方面有較廣泛的應用［22］。本研究通過機器學習算法和Landsat-8 OLI 數據提升了傳統草地生物量估算模型的精度，但是還存在一定的局限性。Landsat-8 OLI 的空間分辨率是30 m，比起MODIS 的影像空間分辨率（250~1000 m）高很多，但其時間分辨率是16 d。另外，研究區地處青藏高原，受云量和當地天氣影響較大，衛星過境重訪周期長，不能很好地獲得高質量的影像，可能會對模型精度和AGB 動態監測造成一定的影響。在未來研究中，可以獲取更高時間分辨率的Landsat 數據，進一步提高模型的精度。另外，草地生物量受地理位置、氣候、土壤、地形等多種因素的影響［39］。其中氣候因素包括光照、氣溫和降水等；土壤因素包括土壤營養元素、土壤結構和肥力等。并且還受到草地類型與分布、物種多樣性等生物因素，以及放牧、圍欄封育、輪牧等管理因素的共同影響［40］。未來可以考慮土壤、地形、氣象因素等變量參與模型的構建，優化現有RF 模型。

4 結論

本研究基于Landsat-8 OLI 遙感數據，分別構建并比較了29 個植被指數與草地AGB 的單因素回歸模型和RF模型，主要得出以下結論：1）綜合并確立了29 個國內外廣泛應用的植被指數與高寒草地AGB 的線性最優關系，表明多種植被指數與草地AGB 之間均具有顯著的相關性。2）多種植被指數之間具有一定的互補性，機器學習算法可以更好地進行變量篩選和組合，因此基于機器學習的多種植被指數的組合應用在很大程度上可以提升草地AGB 模型的反演精度，從而更加精準地反映草地生物量的時空變化狀況。3）得到了2019-2021 年門源縣30 m分辨率的草地AGB 分布。空間分布特征為西北部較高，東南部相對較低；大體呈中部高，四周低的狀況，其中高寒草甸分布最多。