基于MODIS數據與機器學習的青藏高原草地地上生物量研究

金哲人，馮琦勝，王瑞涇，梁天剛

（蘭州大學草地農業生態系統國家重點實驗室，蘭州大學農業農村部草牧業創新重點實驗室，蘭州大學草地農業教育部工程研究中心，蘭州大學草地農業科技學院，甘肅蘭州 730020）

我國幅員遼闊，擁有約960萬km2的陸地面積，其中有400萬km2的草原，占陸地面積的40%；同時，我國包含有被稱為“世界屋脊”的青藏高原，占陸地面積的1/4，其中約70%是草原［1-3］。草地資源不僅是陸地上綠植中所占面積最大的可再生資源，而且對畜牧業的發展也起到重要的基礎作用［4］。草地地上生物量（aboveground biomass，AGB）作為評估草地生態環境的重要載體可監測草地退化程度。草地地上生物量在評估區域碳循環中作用重大，是用來判斷草地能否可持續健康發展的重要指標［5］。研究青藏高原草地地上生物量對當地草地資源的監測與良性發展具有積極影響。

運用遙感技術分析草地資源擁有比野外勘測更加高效和便捷的優勢，可以對大范圍的草地資源進行監測。國內外利用衛星影像對草地資源的研究已有很多，很多研究采用實測數據與植被指數結合進行草地生物量反演，對草地生物量變化進行評估。趙慧芳等［6］研究發現基于MOD13Q1數據的青海省草地地上生物量與植被指數NDVI具有從東南部向西北部逐級遞減的空間分布特征。喬蕻強等［7］利用MODIS-NDVI和MODIS-EVI通過植被指數對肅南縣草原生產能力進行估算，建立研究區草地估產模型。杜玉娥等［8］結合實測數據并利用MODIS植被指數NDVI和EVI資料，估算了2002-2009年青藏高原草地生物量變化。馬青青等［9］結合瑪曲2016年地面實測數據和MODIS影像資料，得出NDVI和地上生物量在6-9月相關性較高，其中9月的擬合水平（R2=0.5102）最高。陸蔭等［10］基于實測數據和MOD13Q1植被指數產品，反演了2000-2019年甘南州草地地上生物量時空分布特征，結果表明MODIS EVI植被指數適合于當地草地生物量變化研究，模型決定系數R2=0.52149，RMSE=527.9 kg·hm-2。姚興成等［11］利用MOD13Q1數據，建立對云南省草地AGB遙感估算模型，之后引入高度和蓋度采用優化的估算模型分析其空間分布格局，模型的估算精度由原來的35.0%提升到43.7%。Gao等［12］利用MOD13Q1數據和隨機森林模型與多因子（NDVI、緯度、經度、海拔、草地類型）結合，估算了2000-2017年青藏高原高寒草地地上生物量和蓋度，并提高了估算精度。Zeng等［13］利用隨機森林算法結合野外觀測數據、MOD13A2遙感植被指數（NDVI、EVI）、氣象數據和地形數據，對2000-2014年青藏高原草地地上生物量進行估算，所建立的隨機森林模型在草地地上生物量估算中表現良好，能解釋觀測數據86%的變化。Zhou等［14］采用支持向量機（support vector machine，SVM）、隨機森林（random forest，RF）和高精度地表建模（high accuracy surface modeling，HASM）3種機器學習模型對2001-2019年中國三江源地區草地AGB進行了模擬，HASM模型比RF和SVM模型表現更為優秀（R2=0.8459＞0.72＞0.5858；RMSE=29＜41＜56）。Chen等［15］使用多層感知器神經網絡的機器學習方法，結合現場觀測和氣候數據，從Sentinel-2影像中估算奶牛場牧草地上生物量，大約60%的生物量變化得到了解釋。通過以往研究發現，目前所選用的遙感數據多為月數據，實測數據與遙感影像時間匹配跨度大，會存在實測數據與遙感數據不完全匹配的問題；并且構建草地地上生物量的模型精度較低。

基于上述原因，本研究基于MODIS逐日反射率數據構建的多種植被指數，比較分析了不同機器學習方法反演青藏高原草地地上生物量的適宜性，評估了2000-2020年青藏高原草地地上生物量的變化，研究結果有助于對草地資源變化狀況進行了解，為青藏高原地區草地生物量估算提供科學依據，服務于對青藏高原草地地上生物量的相關研究。

1 材料與方法

1.1 研究區域

青藏高原（25°59′26″-40°1′6″N，67°40′37″-104°40′43″E）位于我國西部，東至祁連山東端、橫斷山等山脈東緣，南抵喜馬拉雅山等山脈南麓，西起興都庫什山脈和帕米爾高原西緣，北至帕米爾高原北緣、西昆侖山和祁連山山脈北麓，平均海拔約4320 m，總面積為3.08×106km2，其中我國境內的青藏高原面積約2.58×106km2［16］。以青藏高原地區為研究區域（圖1），整個青藏高原地區天然草地面積約為1.5×108hm2，在高原氣候的特殊影響與作用下，整個地區的草原類型以高寒草甸和高寒草原為主［17］。青藏高原不僅影響我國的氣候與生態環境，甚至對世界的環境變化也具有重要作用，其中的草地資源尤其突出表現在維持該地區生態系統的穩定性方面，發揮了高原生態屏障的保護功能［18］。青藏高原的高寒草甸、高寒草原等大面積的草地資源作為高原本身的一道生態屏障，對防控當地土地沙化和水土流失具有積極的影響［19］。

1.2 數據采集與預處理生物量數據集

本研究中所用的草地外業調查數據來自2003-2020年5-9月青藏高原的實地觀測調查。樣地設置在地勢平坦、草地植物群落均一性較好的地段，每塊樣地隨機選擇3～5個大小為0.5 m×0.5 m的樣方。通過GPS測定樣方的經緯度和海拔，并記錄樣地所隸屬的行政區、樣方內的物種名稱、蓋度、高度等植物群落信息。地上生物量采用齊地面刈割的方法，現場稱其鮮重，之后帶回實驗室在70℃烘干至恒重稱其干重。在本研究中，共采用4823條采樣點的實測數據來評估該地區的草地生物量（圖1），其中有1338條采樣點的實測數據是使用Xia等［20］估算青藏高原草地生物量和周轉時間的研究中所使用的數據。

氣象數據來源于國家科技基礎條件平臺-國家地球系統科學數據中心（http：//www.geodata.cn）提供的2000-2020年逐月氣象差值數據集。該氣象數據集包括逐月平均氣溫（℃）、逐月最低氣溫（℃）、逐月最高氣溫（℃）和逐月降水量（mm）數據，空間分辨率為0.0083333°（約1 km）。數據是據CRU（氣候研究中心Climatic Research Unit）發布的全球0.5°氣候數據以及WorldClim發布的全球高分辨率氣候數據，通過Delta空間降尺度方案在中國地區降尺度生成。并用496個獨立氣象觀測點數據進行驗證，驗證結果可信［21］。

表1 本研究中使用的植被指數Table 1 The vegetation index used in this study

1.3 機器學習方法

本研究主要運用R語言的caret包（classification and regression training）進行機器學習模型的構建，該包提供的工具有數據分割、預處理、特征選擇、使用重采樣的模型調整、變量重要性評估等。本研究使用了caret包提供的機器學習方法包括：基于樣條函數的廣義加性模型（generalized additive model using splines，bam）、貝葉斯規整化神經網絡（bayesian regularized neural networks，brnn）、提升樹（boosted tree，bstTree）、彈性網絡（elasticnet，enet）、基于徑向基函數核的高斯過程（gaussian process with radial basis function kernel，gaussprRadial）、廣義線性模型網（generalized linear model net，glmnet）、k-臨近算法（k-nearest neighbors，kknn）、基于正向選擇的線性回歸（linear regression with forward selection，leapForward）、單 調 多 層 感 知 機 神 經 網 絡（monotone multi-layer perceptron neural network，monmlp）、非負最小二乘（non-negative least squares，nnls）、嶺回歸（ridge regression，ridge）、基于徑向基函數核的支持向量機（support vector machines with radial basis function kernel，svmRadialSigma）、稀疏偏最小二乘（sparse partial least squares，spls）、條件推斷隨機森林（conditional inference random forest，cforest）、平行隨機森林（parallel random forest，parRF）、分位數隨機森林（quantile random forest，qrf）、基于ranger的隨機森林（random forest，ranger）、基于Rborist的隨機森林（random forest，Rborist）、基于randomForest的隨機森林（random forest，rf）、基于randomForest和RRF的正則化隨機森林（regularized random forest，RRF）、基于RRF的正則化隨機森林（regularized random forest，RRFglobal）、極端的梯度提升1（extreme gradient boosting 1，xgbDART）和極端的梯度提升2（extreme gradient boosting 2，xgbTree）共23種，探索適合特殊地區青藏高原的機器學習方法。各個方法具體描述見caret包的說明文檔（https：//topepo.github.io/caret/index.html）。

1.4 變量的選擇

本研究在選取變量的時候考慮了兩個方面，一方面是植被指數，一方面是溫度降水，本研究認為溫度降水對植被敏感，所以本研究用了溫度和降水數據，通過查閱相關文獻使用了大家常用的反映植被生長狀況的植被指數。

將2000-2020年的青藏高原各個相關數據集整合，變量共包含MCD43A4產品的7個波段（Band1、Band2、Band3、Band4、Band5、Band6、Band7）、10種 植 被 指 數（NDVI、EVI、EVI2、DVI、RVI、SAVI、MSAVI、OSAVI、SATVI、NDPI）、1-12月的每月降水量（prec_01、prec_02、prec_03、prec_04、prec_05、prec_06、prec_07、prec_08、prec_09、prec_10、prec_11、prec_12）、年降水量（prec_all）、1-12月的每月最高溫度（tmax_01、tmax_02、tmax_03、tmax_04、tmax_05、tmax_06、tmax_07、tmax_08、tmax_09、tmax_10、tmax_11、tmax_12）、1-12月的每月最低溫度（tmin_01、tmin_02、tmin_03、tmin_04、tmin_05、tmin_06、tmin_07、tmin_08、tmin_09、tmin_10、tmin_11、tmin_12）、

1-12月的每月平均溫度（tp_01、tp_02、tp_03、tp_04、tp_05、tp_06、tp_07、tp_08、tp_09、tp_10、tp_11、tp_12）和年平均溫度（tp_year），一共67個變量。

本研究采用R語言的CAST包中的前向特征選擇方法（the forward feature selection，ffs）進行變量篩選，此篩選方法先基于全變量中任意2個變量訓練模型，選出精度最高的模型所用變量，之后逐個添加特征變量，并選取模型精度提高最大的特征變量，當模型精度不再提高時停止變量篩選，最終得到最優特征變量組合。

1.5 草地地上生物量變化趨勢分析

采用Theil-Sen median趨勢分析、Mann-Kendall以及Hurst指數方法，研究青藏高原覆蓋區域草地生物量的空間分布特征、時間變化趨勢、變化趨勢特征和可持續性特征［32-33］。Theil-Sen median是一種非參數統計的趨勢計算方法，效果穩健并可以減少數據異常值的影響［32］。Mann-Kendall作為一種非參數統計檢驗方法，優點是不需要樣本服從特定的分布，可以直接檢驗時間序列的單調變化趨勢［34-36］。Theil-Sen median趨勢分析和Mann-Kendall趨勢檢驗結合可判斷長時間序列數據趨勢，能夠有效分析青藏高原草地生物量時間序列的變化趨勢［37］。Hurst指數是定量描述時間序列信息長期依賴性的有效方法［37］。Theil-Sen median趨勢分析、Mann-Kendall以及Hurst指數（即H值）方法的計算公式詳見文獻［32］。將結果劃分為輕微降低、顯著降低、穩定不變、輕微增加和顯著增加5種類別。H值取值越接近1，持續性越強；H=0.5，說明不存在長期相關性；H越接近0，反持續性越強。

2 結果與分析

2.1 采樣點數據分析

按草地類型統計青藏高原逐年實地采樣數據的結果表明（表2）：在采樣數據中，風干重由高到低排序，前3位依次是山地草甸類、沼澤類和高寒草甸類，分別達到2001.48、1812.50和1533.51 kg·hm-2。風干重較低的3位是高寒荒漠草原、溫性荒漠和高寒荒漠，分別為406.76、393.85和304.06 kg·hm-2，同時高寒荒漠的標準偏差也最低，只有89.62 kg·hm-2。山地草甸的蓋度達到87.13%，高寒草甸的蓋度達到83.15%，沼澤的蓋度達到75.50%，是所有草地類型中蓋度較高的前3位。蓋度最低的3位是溫性草甸草原、高寒荒漠草原和溫性荒漠，分別為35.62%、29.16%和25.24%。各草地類型的高度大多集中在4～15 cm，低地草甸平均高度最高，達到了26.12 cm，高寒草甸草原平均高度最低，只有4.74 cm。

(3) 電網輸配電效益下降問題。分布式電源使電網輸配電效益下降幾乎是不可避免的,主要表現為:分布式發電影響電網的售電量,從而影響電網的效益;增大公共電網網損率,降低公共電網的設備利用率。

表2 不同草地類型的調查樣本數據基本情況Table 2 Survey sample data and basic statistics of different grassland type

2.2 變量相關分析

利用ffs算法篩選變量的結果如圖2所示。經過篩選得出變量“prec_05”、“prec_06”、“tp_12”、“NDPI”、“prec_04”、“tmax_01”、“prec_08”、“prec_12”與生物量關系較大，其中篩選出的8個變量中有5個與月降水量有關，它們各自之間的相關系數關系如圖3所示。結果表明：與生物量相關性最大的變量是“NDPI”，相關系數值為0.589，其次是“prec_05”和“prec_04”，即4、5月的降水量，與生物量的相關性達到了0.400以上。

2.3 基于機器學習的最優模型

運用R語言中caret包，將67個變量全部帶入所選用機器學習模型中運算與生物量的關系，結果如圖4所示。全變量與生物量決定系數最高的模型是rf，R2達到了0.6360。決定系數較低的后3位模型是glmnet、nnls和leapForward，R2都在0.46以下。按照R2由高到低排序，前8位模型依次為rf、RRF、RRFglobal、Rborist、parRF、ranger、qrf、cforest，它們都屬于隨機森林類方法，R2都在0.61以上。其次R2較高的模型是xgbTree、xgbDART、kknn、svmRadialSigma和bstTree，這些模型的R2也都達到了0.57以上。說明對于研究全變量和生物量的關系，隨機森林rf模型為最優模型。

基于篩選出來的8個變量構建機器學習模型，根據圖4可知，在所用的機器學習的方法中，R2由高到低前5位依 次是Rborist、RRF、ranger、RRFglobal和parRF，R2分 別 為0.6484、0.6477、0.6474、0.6473和0.6460。R2在0.40以下的模型有spls、glmnet、enet、ridge、leapForward和nnls模型，說明這些方法運行效果欠佳，并不適合研究篩選后的變量與生物量的關系。結果表明Rborist、RRF、ranger、RRFglobal、parRF、rf、qrf和kknn這8種方法基于篩選變量的模型比全變量模型的R2值高，說明減少變量可以消除數據冗余，從而提高模型精度。由此本研究選取出最優的機器學習模型為Rborist。篩選前后各機器學習方法參數取值如表3所示。

2.4 青藏高原2000-2020年生物量變化分析

基于選擇出的最優機器學習模型，本研究模擬了青藏高原2000-2020年草地地上生物量（圖5）。結果表明，青藏高原的草地地上生物量由西北向東南方向逐步增加，東南部的草地地上生物量普遍在900 kg·hm-2以上，要明顯好于西北部。

采用Theil-Sen median趨勢分析和Mann-Kendall方法，可有效反映出青藏高原地區2000-2020年的草地生物量的時空變化趨勢，如圖6所示。研究結果表明，青藏高原草地生物量輕微降低區域占草地總面積的9.83%、顯著降低區域占比1.33%、穩定不變區域占比64.33%、輕微增加區域占比18.24%，顯著增加區域占比6.26%。草地顯著降低的面積最少，穩定不變的面積最大，并且從圖中可看出草地生物量降低的區域主要集中在青藏高原東部和中部，增加的區域分布在青藏高原東北部。

將Theil-Sen median趨勢分析和Mann-Kendall檢驗結果與Hurst指數相結合，得到青藏高原草地地上生物量變化趨勢與持續性的耦合信息，如圖7所示。研究結果表明，青藏高原草地地上生物量變化趨勢不穩定的區域占比61.38%，持續性輕微降低的區域占比4.67%，持續性明顯降低的區域占比1.19%，持續性穩定不變的區域占比19.71%，持續性輕微增加的區域占比7.87%，持續性明顯增加的區域占比5.18%。總的來說，青藏高原2000-2020年草地地上生物量增加的區域要大于降低的區域，總體向好發展。

?

將草地類型按青藏高原行政區劃統計（表4），結果表明各省份不同草地類型變化趨勢不確定的區域面積占比較大，其中甘肅省、青海省和四川省持續性增加的草地區域面積都達到了20%以上，四川省持續性降低的草地區域面積超過了10%。山地草甸類在各省份均持續性輕微降低，其中西藏持續性輕微降低的區域面積占比最大，達到15.57%。甘肅省持續性明顯增加的草地類型前3位依次是溫性荒漠草原類、高寒草原類和溫性荒漠類。西藏各草地類型生物量均有所降低，其中溫性草甸草原類持續性輕微降低的區域面積達到20%以上，持續性增加總面積比率與退化的比例相持平。整體上各個省份大部分區域的各草地類型在持續性增加。

表4 2000-2020年青藏高原各類草地變化特征分析Table 4 The change characteristic analysis of different grassland from 2000-2020 on the Tibetan Plateau（%）

續表Continued Table

3 討論

3.1 機器學習模型的優劣

本研究基于地面實測數據和MCD43A4產品構建的植被指數，使用ffs篩選得到變量“prec_05”、“prec_06”、“tp_12”、“NDPI”、“prec_04”、“tmax_01”、“prec_08”、“prec_12”和生物量的關系較大，運用最優模型Rborist，對2000-2020年青藏高原地區的草地生物量進行反演，并對青藏高原2000-2020年的草地地上生物量時空變化趨勢進行分析。Meng等［38］基于遙感數據和機器學習算法對2011-2016年甘南地區草地地上生物量建模，用隨機森林算法構建了最優估計模型，該模型可解釋牧草生長季AGB的89.41%的變化。Zeng等［13］利用隨機森林算法對2000-2014年青藏高原草地生物量反演，該研究在RF建模過程中選擇MOD13A2遙感植被指數（NDVI和EVI）、氣象變量（年平均溫度和年平均降水量）和地形變量（海拔、坡度和坡向），所建立的隨機森林模型在草地地上生物量估計中表現良好，能解釋觀測數據86%的變化。Lin等［39］利用優化算法和多維特征改進高寒草地植被覆蓋度估算，其中隨機森林算法（R2：0.861，RMSE：9.5%）的植被覆蓋度反演效果最好。Tang等［40］利用機器學習方法對2001-2020年黃河源地區反演草地生物量，研究結果表明隨機森林模型優于其他3種模型。本研究主要是利用多種機器學習方法對草地生物量進行反演。同樣發現隨機森林方法效果最佳，R2由高到低排序前10位中有8位是隨機森林的方法，前5位都是隨機森林的方法。隨機森林優點是可以準確描述生物物理參數與復雜環境因素之間的復雜非線性關系，并且易于實現，因為它的參數簡單，預測變量不需要預處理，且具有較好的泛化性和魯棒性［41］。本研究所選用的最優模型是隨機森林方法中基于Rborist包的隨機森林（random forest，Rborist）。Rborist包能夠實現隨機森林算法，并且擁有較高性能。這個包是Arborist項目的衍生品，Arborist項目是一個針對各種決策樹方法的多語言項目。Arborist在實現隨機森林的基礎上具有可擴展和高性能特點，支持回歸和分類決策樹模型，能夠在訓練和測試中利用多核并行運算，并且具有集群友好性。

3.2 研究結果的一致性

植被指數對地表植被狀況敏感，張雅等［42］基于Landsat 8數據提取常用的植被指數，利用統計分析方法建立紫泥泉牧場陰坡和陽坡的草原生物量遙感估算模型，并進行了生物量空間反演與驗證，其中NDVI與生物量的相關性最高。葛靜等［43］基于ADC（農業數字相機agricultural digital camera）和MODIS遙感數據對黃河源區的高寒草地地上生物量進行監測研究，結合MODIS NDVI構建草地地上生物量反演模型。Wang等［31］提出的NDPI植被指數表現優異。Gan等［44］利用MODIS MOD09A1產品提取的植被指數結合多種方法監測2009-2013年黃淮冬小麥（Triticum aestivum）返青期（green-up dates，GUDs），結果表明NDPI與地面數據的一致性優于其他植被指數。Su等［45］利用Sentinel-2衛星遙感影像，通過對植被光譜反射率與降塵吸收率（the amount of dust absorption，ADA）的相關性分析，建立降塵反演模型，獲取常綠植被的降塵分布，結果表明NDPI更適合建立反演模型。An等［46］基于MODIS MOD09A1和不同算法分別監測青藏高原植被生長季的開始和結束時間，研究結果指出基于EVI和NDPI的物候指標時空穩定性更好。本研究選用的10種植被指數與草地生物量之間的相關系數如圖8所示，可以看出NDPI與草地生物量關系更密切。原因在于NDPI使用加權的red-SWIR組合來代替NDVI中的紅波段，克服了土壤背景異質性的不利影響；同時NDPI通過添加SWIR波段來響應葉片水分含量，由于水分在草本植被總重量中所占比例較高，因此NDPI能夠捕捉葉片水分含量的變化，考慮了草地葉片水分變化對生物量的影響［47］。本研究結果與Xu等［47］認為NDPI在草地地上生物量估算中具有優越性相一致。

本研究經過篩選得出8個變量中有5個與月降水量有關，說明降水量對青藏高原草地生物量影響很大。卓嘎等［48］對2000-2016年青藏高原植被覆蓋時空變化及其對氣候因子的響應的研究結果表明降水量對植被覆蓋的影響大于平均氣溫對植被覆蓋的影響。Zhang等［49］基于MODIS MOD09A1產品和氣象數據，監測了青藏高原的綠度變化趨勢及其驅動因素，結果表明青藏高原東北部植被對降水相對敏感，而高原南部植被對溫度相對敏感。Zeng等［13］利用隨機森林算法估算青藏高原草地地上生物量，研究指出在整個研究區草地AGB與氣溫、降水呈顯著正相關。草地AGB與年平均降水量的相關性為0.54（P＜0.05），遠高于年平均溫度（R=0.38，P＜0.05）。Shen等［50］利用NDVI和氣象數據研究生長季前溫度和降水對青藏高原中東部草地春季物候的影響，研究結果表明溫度比降水對青藏高原草地返青期空間格局的影響更大。韓炳宏等［51］基于2000-2018年MODIS NDVI數據及氣象數據，探討了青藏高原地區NDVI的時空變化趨勢及其與氣溫、降水的關系，研究發現青藏高原地區NDVI與同期氣溫和降水均顯著相關，但與氣溫的關系更密切。綜上，目前對影響青藏高原地區草地生物量的氣象因素，有學者認為主要是氣溫，有學者認為是降水。出現差異的原因很可能是青藏高原不同地區對氣溫和降水敏感度不一樣，而且研究的年份、使用的數據和研究的方法不同都可能導致結果不同［49］。

3.3 研究展望

本研究證明了NDPI在青藏高原地區反演草地生物量有出色的表現，并且證實了Rborist方法對青藏高原草地生物量建模的有效性。今后將繼續深入研究NDPI的適用性，并進一步分析溫度和降水與草地生物量的關系。在未來的工作中，考慮把土壤性質等其他環境因素納入建模過程，以提高預測性能，并且最優模型的選取是否適合運用在其他地方還有待進一步驗證。

4 結論

本研究基于機器學習與MODIS數據研究青藏高原草地地上生物量的情況，得出以下結論：

1）構建的機器學習模型中，Rborist模型精度最高，基于篩選后變量的R2達到0.6484。“prec_05”、“prec_06”、“tp_12”、“NDPI”、“prec_04”、“tmax_01”、“prec_08”、“prec_12”這8個變量與生物量相關。

2）青藏高原東南部的生物量要高于西北部，呈現由東南向西北遞減趨勢。

3）2000-2020年間青藏高原草地生物量穩步增長，整體向好發展。青藏高原61.38%的草地變化趨勢不具有可持續性，4.67%的草地持續性輕微惡化，持續性明顯惡化的區域占比1.19%，呈穩定或恢復趨勢的區域占比32.76%。