劉 文, 莫興國*, 劉蘇峽
(1. 中國科學院地理科學與資源研究所陸地水循環及地表過程重點實驗室, 北京100101; 2. 中國科學院大學資環學院/中丹學院, 北京100049)
開展草地地上生物量(Aboveground biomass,AGB)的估算有助于理解草地生態系統的碳貯量和碳平衡。青藏高原草地是中國重要的畜牧業生產基地,也是當地牧民生產生活的重要保障。準確和及時地估算草地AGB能夠為理論載畜量的確定提供理論依據,有助于草地資源的利用與保護。
直接收獲法只適用于小范圍內AGB的監測。由于遙感資料覆蓋范圍廣且具有時空連續性,所以被用于區域尺度的AGB估算。借助遙感資料,AGB的估算方法可以分為兩種。一種是機理模型,例如CASA(Carnegie-Ames-Stanford approach)和GLO-PEM(Global Production Efficiency Model)模型[1-2]。另一種是數據驅動模型,包括回歸模型和機器學習方法,例如指數模型、人工神經網絡(Artificial neural networks,ANNs)和隨機森林算法[3-5]。該方法利用遙感植被指數、氣象、地形、植被覆蓋度、草地類型和生物量采樣點的經緯度等解釋變量與AGB實測值構建經驗關系。再基于此關系,利用這些解釋變量的網格數據模擬區域尺度的AGB。之前的研究在對生長季AGB的平均值或生長季逐月的AGB進行估算時,相應地采用生長季平均或月最大植被指數[4,6]。在對年AGB進行估算時,訓練階段采用實時觀測的植被指數;在區域尺度上開展模擬時,采用年最大植被指數[7]。機器學習方法的計算較機理模型簡便。使用該方法能夠直接估算年AGB,不需要對凈初級生產力(Net primary productivity,NPP)的季節變化進行模擬。機器學習方法也被用于草地NPP的模擬并且與機理模型的精度相當[8]。然而利用NPP估算AGB可能會受到根冠比的不確定性的影響[9]。輻射傳輸模型能夠克服植被指數的飽和效應導致的回歸模型的模擬值偏低的問題。但模型的準確性依賴于實地觀測的其他植被生物物理參數(例如,葉當量水厚度和比葉重的線性關系)[10]。機理模型的模擬精度取決于物理機制的假設的合理性,而使用機器學習方法能直接從數據中提取信息,不需要物理機制的假設。但機器學習方法受限于訓練資料的可獲得性和代表性。ANNs是眾多機器學習方法中的一種,擅于模擬自變量與因變量之間的非線性關系,已被成功應用于草地AGB的估算[4,11]。因此基于ANNs構建AGB估算模型具有實際意義。大多數研究表明ANNs的性能優于線性和非線性回歸[4,11-13]。例如,中國內蒙古錫林郭勒溫帶典型草原AGB模擬的研究表明,ANNs比線性回歸、對數、乘冪和指數模型的估算精度都高[14]。這可能是因為ANNs能夠以任意精度逼近任何連續的多變量函數[15],所以它能夠學習環境因子與AGB之間的復雜函數關系。
基于機器學習的模擬研究主要圍繞輸入變量的選取和校正探討改善模型性能的方法。例如,由于通過衛星遙感獲取的植被指數容易受到云層和氣溶膠的干擾,基于高光譜計算的歸一化差值植被指數(Normalized difference vegetation index,NDVI)可用于校正MODIS(Moderate-resolution Imaging Spectroradiometer)的NDVI,校正的NDVI被作為模型的輸入[7]。使用高光譜和光譜響應函數模擬衛星傳感器的多光譜反射率,能夠解決遙感影像與采樣時間或樣方尺寸不匹配的問題[16]。由于半干旱生態系統中植被覆蓋度較低(<50%),NDVI容易受到土壤背景信號的影響,增強型植被指數(Enhanced vegetation index,EVI)也被作為模型的輸入[17]。將遙感植被指數和合成孔徑雷達數據共同作為輸入能夠提高僅以其中一種數據作為輸入的模擬精度[18]。地面激光掃描(Terrestrial laser scanning,TLS)能夠準確提取植被的結構參數(例如,冠層的平均高度和覆蓋度),是AGB的有效解釋變量[19]。然而,由于TLS的觀測數據還未與近地面遙感(例如,無人機)和衛星遙感的數據建立關系,限制了區域尺度的模擬。太陽誘導葉綠素熒光(Solar-induced chlorophyll fluorescence,SIF)與植物實際的光合作用直接有關[20]。少數基于機器學習方法開展農作物產量模擬的研究表明,將SIF作為輸入能夠提高作物產量的模擬效率[21-22]。然而將SIF作為機器學習模型的輸入是否也有利于草地AGB的估算有待研究。
少數研究基于ANNs和隨機森林算法開展了青藏高原草地AGB的估算,然而機器學習與機理模型的模擬值的時空變化還未被比較。借助機理模型模擬的AGB時空變化能夠對ANNs的模擬結果進行驗證。單個機理模型的模擬結果具有不確定性,但多個模型的模擬結果可能提供了AGB時空變化的合理范圍。此外,AGB在不同草地類型和海拔地區中的時空變化還未被分析。之前的研究表明,AGB與氣溫(Ta)呈顯著(α=0.05)的正相關(R2:0.45~0.62),而與降雨(PRE)的相關性不顯著[6,23]。然而,AGB可能與前一年的PRE具有顯著的相關性[24]。除了PRE和Ta,大氣CO2濃度也是草地NPP增加的主要影響因子之一[25]。變暖導致青藏高原的飽和水汽壓差(VPD)顯著增加,VPD對草地生產力產生的消極影響的程度和范圍都在顯著增加[26]。但前一年的PRE和其他氣候因子對AGB變化的影響強弱還未知。此外,PRE或Ta與AGB的相關性包含了其他氣候因子的間接影響。每個氣候因子的直接影響的強弱有待比較。本研究基于遙感植被指數和其他輔助數據構建青藏高原草地AGB的機器學習模型;說明SIF的加入是否有利于模擬精度的提升;展示AGB在整個草地、不同草地類型和海拔區間中的時空變化;揭示ANNs與機理模型模擬值的時空變化的差異性;并識別AGB變化的主要氣候影響因子。ANNs的模擬資料為機理模型的驗證提供參考,為青藏高原理論載畜量的確定提供數據支撐。分析結果為草原管理提供理論依據,有助于促進畜牧業的可持續發展。
青藏高原是世界上海拔最高和面積最大的高原,平均海拔在4 000 m以上(圖1a)。年平均氣溫在-15~10℃。降雨主要集中在6—9月份。年降雨量呈東南地區高(>600 mm)和西北地區低(<200 mm)的分布特征。高寒草原約占高原面積的75%。草地類型(圖1b)主要包括[27]:溫帶禾草、雜類草草甸草原(TMS),溫帶叢生禾草典型草原(TAS),溫帶叢生矮禾草、矮半灌木荒漠草原(TDS),高寒禾草、苔草草原(AS),亞熱帶、熱帶草叢(STGC),溫帶禾草、雜類草草甸(TGM),溫帶禾草、雜類草鹽生草甸(TSM)和高寒嵩草、雜類草草甸(AM)。海拔來源于STRM(Shuttle Topography Radar Mission)V4.1版本的數據,分辨率為1 km(https://www.resdc.cn)。為了分析AGB隨海拔和草地類型的變化,將海拔和草地類型數據都重采樣為5 km。

圖1 青藏高原海拔、草地類型年降水和平均氣溫空間分布Fig.1 Spatial distribution of altitudes,grassland types, annual precipitation, and annual mean air temperature on the Tibetan Plateau
環境數據來源于遙感和再分析資料(表1)。ANNs的輸入為年最大EVI、年平均SIF,Ta、下行短波輻射(SRAD)、表層土壤水分(SM)和年PRE。其中,年最大EVI為16天EVI的年最大值。年平均SIF和SM分別為4天SIF和月平均SM的年平均。在訓練階段,為了與EVI數據的分辨率匹配,將其他輸入數據都重采樣為250 m分辨率。再根據采樣點經緯度和采樣年份提取了相應格點的環境數據。從采樣點周圍5 km范圍內選取具有代表性的EVI,使EVI與AGB達到較強的相關性。一共調整了82個采樣點的EVI。在預報階段,以5 km分辨率開展AGB的模擬,輸入數據都被重采樣為5 km分辨率。使用ArcGIS 10.5軟件以雙線性插值方法進行了重采樣。為了檢驗以SIF作為輸入變量對估算精度的影響,設置了2種輸入方案,即輸入變量不包括SIF和包括SIF,分別記作ANNs-S1和ANNs-S2。

表1 遙感和再分析數據信息Table 1 Remote sensing and reanalysis data information
AGB觀測資料包括從文獻[34-42]中和國家生態科學數據中心(海北站,www.cnern.org.cn)搜集的224個AGB實測值和本研究的8個實地觀測值,采樣點主要分布在西藏自治區、青海省和四川省,草地類型包括AM,AS和TSM。文獻和本研究中的AGB測定方法一致。選擇分布均勻的同一類型的草地作為樣地,面積為10 m×10 m~250 m×250 m;其中布置3~9個重復樣方,樣方大小為0.3 m×0.3 m~1 m×1 m,采用齊地刈割法收獲地上部分鮮草。將草樣放入烘箱以65℃烘干至恒重。采樣時間為一年中的7—9月份,該時段內的AGB基本達到了年最大值。在訓練階段,對同一格點內的多個觀測值進行了平均。
采用了MODIS、MuSyQ(Multisource Data Synergized Quantitative),CASA,GLO-PEM,VIP(Vegetation Interface Processes)等5套NPP資料(表2)。5套資料的分布時段與研究時段一致,并且空間分辨率較高。5個模型的計算原理不同。VIP模型基于植物的生物化學過程計算NPP,而其他4種模型都基于光能利用率計算NPP,但光能利用率的計算都存在差異。使用不同的NPP數據有利于對ANNs的結果進行充分地驗證。利用根冠比估算地上NPP,再將地上NPP除以系數0.45換算為干重,即AGB。草地類型AM,AS,TDS,TMS,TSM,TAS,TGM和STGC的根冠比分別為3.09,3.18,6.4,5.2,3.5,5.6,3.5和3.5[42-43]。為了與ANNs的模擬資料進行比較,將5套AGB資料都重采樣為5 km分辨率。5套資料與觀測值的相關性都較好(R2=0.65~0.72,P<0.01)。

表2 NPP資料信息Table 2 NPP data information
采用的ANNs結構類型為3層的前饋神經網絡,包括輸入層、隱含層和輸出層。訓練算法為基于LM(Levenberg-Marquardt)算法的反向傳播算法。將此種ANNs模型記作BP(Back propagation)-ANNs。采用試算的方法確定隱含層神經元的個數。同時保證參數的個數不多于訓練樣本的個數。隨機選取80%的樣本作為訓練集,剩余樣本作為驗證集。驗證集被用于評估訓練過程中和結束后的預測精度。當驗證誤差不再減小甚至開始增加時,訓練停止。采用該方法是為了避免過擬合。
基于十折交叉驗證[4]評估BP-ANNs的性能。將觀測資料隨機劃分為10組,每組樣本數量相等。每組樣本輪流作為測試集,其它9組樣本作為訓練資料。類似地,為了避免過擬合,從其他9組樣本中隨機選擇1組樣本作為驗證集,其他8組樣本用于訓練。將10組訓練和測試誤差的平均值作為BP-ANNs的性能指標。
AGB的影響因子主要考慮了前一年的PRE、當年的PRE,Ta,VPD,SRAD,風速(U)和CO2濃度。輻射是年代際AGB增加的主要影響因子之一[48]。U的增加有利于土壤蒸發,蒸發通過消耗土壤水分引起植被生長的水分脅迫,所以將U也作為AGB的影響因子。年平均CO2濃度來自莫納羅亞火山觀測站(https://gml.noaa.gov)。年平均U來自中國區域地面氣象要素驅動數據集(CMFD)。采用CMFD中的逐日比濕、氣溫和氣壓計算逐日VPD,再計算年平均VPD,計算方法參見[26]。為了減輕多重共線性對最小二乘估計的影響,基于嶺回歸[49]計算每個因子的標準化回歸系數。根據系數絕對值的大小識別每個因子對AGB變化影響的強弱。在AGB下降區域,輸入變量不包括CO2濃度。
年最大EVI與AGB的相關性較強,而相應格點上SIF與AGB的相關性較低(圖2)。這可能是由于SIF的原始分辨率較低,不能較好地反映采樣點的植被生長狀況。基于2種輸入方案的BP-ANNs在訓練階段的擬合值分別解釋了觀測值變化的92%和93%,均方根誤差(RMSE)分別為18.48和18.13 g·m-2(圖3a)。十折交叉驗證的結果表明,基于2種輸入方案的BP-ANNs在訓練和驗證階段的模擬值與實測值的決定系數(R2)在 0.90~0.92之間;在測試階段的R2分別為0.91和0.90,RMSE分別為23.34和24.52 g·m-2(表3,圖3b)。由于SIF的加入未改善AGB的估算精度,以下分析都基于ANNs-S1的模擬值。

圖2 AGB觀測值和年最大EVI(a)或年平均SIF(b)的散點圖Fig.2 Scatter plot of AGB observations and annual maximum EVI or annual mean SIF. 注:實線為對數據的最小二乘線性擬合;陰影為線性回歸估計值的置信區間,顯著性水平為α=0.05Note:The solid line is the least square linear fitting of the data. The shadow indicates the confidence interval for the linear regression estimates with a significance level at α=0.05

圖3 訓練期的擬合值與觀測值的散點圖(a)和基于十折交叉驗證的預測值與實測值的散點圖(b)Fig.3 Scatter plot of fitted values versus observations in the training period (a) and scatter plot of predictions versus observations based on 10-fold cross-validation (b)

表3 基于十折交叉驗證的BP-ANNs的模擬精度Table 3 Performance of BP-ANNs based on 10-fold cross-validation
草地AGB呈東南地區高、西南地區低的空間分布特點(圖4a),與5套AGB模擬資料的空間分布特征一致。西部和東部地區的AGB分別在50 g·m-2以下和100 g·m-2以上,與觀測值的空間分布一致。2001—2018年草地AGB的區域平均為(74.75±52.65) g·m-2。
草地AGB呈增加和減少趨勢的區域分別占草地面積的61%和39%。AGB呈顯著(α=0.05)增加趨勢的區域占草地面積的13%,主要分布在藏北和青海省東北部(圖4b);而呈顯著減少趨勢的區域僅占4%,零星地分布在西藏中部、青海省南部和四川省西北部。

圖4 2001-2018年草地AGB的多年平均(a)和年際變化(b)的空間分布Fig.4 Spatial distribution of annual mean (a) and interannual variations (b) of AGB of grassland during 2001-2018注:(a)中直方圖表示每個區間內AGB的頻率;圖(b)中“+”表示變化顯著(α=0.05)的區域Note:The inside histogram in (a) indicates the frequency of AGB in each interval;“+” in (b) mark the areas with significant (α=0.05) trends
STGC和TGM的AGB最高,其次為AM,TMS,TAS和AS的AGB(圖5a)。TDS和TSM的AGB最低。AM和TAS中的AGB的標準差最大,說明AGB的空間異質性最高;其他草地類型中AGB的空間異質性隨AGB的減少而減小。STGC,TAS,TDS和TSM的AGB都呈顯著增加趨勢(圖5b),其中STGC的AGB的增加速率最大,其次為TAS,TDS和TSM。

圖5 不同草地類型中多年平均AGB的區域平均(a)和區域平均AGB的年際變化(b)Fig.5 Regional averaged annual mean AGB and the interannual variations of regional averaged AGB in different grassland types注:(a)中的誤差線表示在每種草地類型中所有格點上的多年平均AGB的標準差;“*”和“**”表示趨勢的顯著性水平分別為α=0.05和α=0.01Note:Error bar in (a) represents the standard deviation of annual mean AGB of all pixels in each grassland type. “*” and “**” represent the significance level of the trend as α=0.05 and α=0.01,respectively
在海拔低于和高于3 800 m的區域中,AGB隨海拔的升高分別呈顯著(α=0.01)增加和減少趨勢(圖6a)。在海拔介于3 400~3 600 m的區域中,AGB的標準差最大,說明它的空間異質性最高;隨著海拔繼續升高,AGB的空間異質性逐漸降低。
在海拔位于4 600~4 800 m和高于5 600 m的地區中,AGB趨勢的平均值為負值;在其他區域中,AGB趨勢的平均值都為正值(圖6b)。在海拔高于5 000 m的區域中,AGB趨勢的空間異質性較低。AGB趨勢的平均值隨海拔的升高呈顯著地減小現象。

圖6 AGB和它的年際變化隨海拔的變化Fig.6 Changes in AGB and its interannual variation along with altitude注:(a)和(b)中的點分別表示每個海拔區間內多年平均AGB和AGB趨勢的區域平均;陰影表示在每個海拔區間中所有格點的多年平均AGB或AGB趨勢的標準差Note:oints in (a) and (b) represent the regional average of the annual mean AGB and AGB trends in each altitude interval,respectively. The shadow indicates the standard deviation of annual mean AGB or AGB trends of all pixels in each altitude interval
BP-ANNs的模擬值在5個機理模型模擬值的范圍內(56.12~160.98 g·m-2)。BP-ANNs的模擬值與實測值的相對誤差(RE)絕對值的平均值為57%,平均絕對誤差(MAE)的平均值為31.72 g·m-2,在機理模型的模擬誤差范圍內(RE:42%~77%;MAE:25.26~43.40 g·m-2)。在大部分草地區域中,BP-ANNs的模擬值低于VIP的模擬值而高于其他4套資料(圖7)。BP-ANNs與GLO-PEM模擬值的RE絕對值的區域平均最大,其次為MuSyQ,VIP,MODIS和CASA。BP-ANNs的模擬值與VIP,CASA和MODIS的模擬值的R2也較高,RMSE也較低。5套資料比BP-ANNs的模擬值都高(低)的區域占草地面積的1%(10%),相應地區的EVI的平均值也較高(低),MAE在54~108 g·m-2(25~84 g·m-2)之間。在三江源地區,BP-ANNs的模擬值(110.18 g·m-2)套資料的范圍內(43.66~166.39 g·m-2),與VIP和CASA的模擬值較接近,與GLO-PEM模擬值的相對誤差最大,為60%。

圖7 BP-ANNs的模擬值的多年平均分別與5套資料多年平均的相對誤差(RE)的概率密度曲線Fig.7 Probability density curves of relative error (RE) between the annual mean from BP-ANNs and the five datasets注:某個格點上的RE等于BP-ANNs的模擬值的多年平均與某一套資料的多年平均的差值除以BP-ANNs的模擬值的多年平均Note:RE in a pixel equals the difference between the annual mean from BP-ANNs and one of the datasets divided by the annual mean from BP-ANNs
VIP與MODIS模型的模擬值在藏北地區呈現顯著增加的趨勢;在東北部分地區,5套資料呈現顯著增加的趨勢。在BP-ANNs和5套資料都有值的區域中,BP-ANNs的模擬值呈顯著增加趨勢的面積比例(7%)與GLO-PEM、CASA和MuSyQ接近(6%,11%和10%),然而低于VIP和MODIS(38%和25%);呈顯著減少趨勢的面積比例(5%)與5套資料都接近(1%~5%)。BP-ANNs的模擬值呈顯著增加和減少的面積比例都與年最大EVI相似(顯著增加:9%,顯著減小:3%)。在區域平均水平上,BP-ANNs的模擬值比5套資料的年際波動性弱(圖8),但變化趨勢有相似之處。例如,2001—2010年BP-ANNs,MODIS,MuSyQ,GLO-PEM和CASA的模擬值都呈微弱的增加趨勢。

圖8 基于BP-ANNs和5個機理模型的2001—2018年青藏高原草地區域平均AGB的年際變化;CV表示變異系數Fig.8 Interannual variations of regional averaged AGB of grassland on the Tibetan Plateau during 2001—2018 based on BP-ANNs and the 5 mechanism models;CV represents coefficient of variation
在AGB增加的區域中,每個因子按重要性從高到低排序依次為CO2濃度、SRAD、前一年的PRE(PRE1),VPD,U,Ta和當年的PRE。在AGB減少的區域中,每個因子按影響從強到弱排序依次為SRAD,VPD,PRE1,U,當年的PRE和Ta。SRAD和當年的PRE在大部分草地(80%和60%)中的回歸系數為負值,說明2種因子都對AGB產生負效應,無實際意義。這可能是由于回歸系數受到了各因子之間多重共線性的影響。
基于隨機森林算法的1982—2010年的AGB為78.40 g·m-2[3],2000—2014年的AGB為77.12 g·m-2[6],這與BP-ANNs的模擬值接近。然而,Gao等[7]使用高光譜反射率計算的NDVI校正MODIS NDVI,基于隨機森林算法得到2000—2017年的AGB為59.63 g·m-2,比BP-ANNs和上述2項研究的模擬值都低。2003—2015年6—10月觀測的1 433個三江源中東部地區的AGB樣本的平均值為108.64 g·m-2[4],BP-ANNs的模擬值的誤差比5個機理模型都小。
在整個草地區域中或在BP-ANNs的模擬值比5套資料都低或都高的區域中,BP-ANNs的模擬值與EVI的R2(0.95~0.98)比5套資料與EVI的R2(0.11~0.75)都高。這可能是因為在訓練階段AGB觀測值與EVI的相關性較高。在整個草地區域中,BP-ANNs的模擬值與SIF的R2較低(0.71),但與MuSyQ和GLO-PEM的模擬值與SIF的R2(0.72和0.73)接近。在BP-ANNs的模擬值比5套資料都低或都高的區域中,BP-ANNs模擬值與SIF的R2在5套資料與SIF的R2的范圍之內。
與AGB實測值對應的EVI范圍為0.05~0.87。在BP-ANNs的模擬值比5套資料都低(高)的區域中,EVI在0.05~0.87范圍內的區域面積比例為99.59%(99.64%),說明與AGB實測值對應的EVI范圍能夠代表大部分區域中EVI的范圍。在5套資料比BP-ANNs的模擬值都高(低)且EVI處于0.05~0.87的區域中,每個機理模型模擬值的平均值比AGB觀測值的平均值都高(低)。由于BP-ANNs以AGB觀測值作為訓練資料,它在該區域中的模擬值的平均值與AGB觀測值的平均值最接近。因此,在該區域中,5套資料也比BP-ANNs模擬值的平均值都高(低)。總之,5套資料比BP-ANNs的模擬值都高(低)是由于5套資料的平均值比BP-ANNs訓練資料的平均值都高(低)。例如,青海、四川和甘肅省交界處的EVI高于0.5,VIP和CASA模型在大部分地區的模擬值大于250 g·m-2。然而與高于0.5的EVI對應的AGB觀測值的平均值為195 g·m-2。將這些地區中的EVI輸入BP-ANNs,獲得的模擬值將小于VIP和CASA模型的模擬值。因為AGB觀測資料有限,它們的平均值是否具有代表性可能存在不確定性。
CO2濃度對AGB的影響較強,且產生正效應。之前的研究也表明,CO2濃度對草地總初級生產力產生正效應,比Ta、U和輻射的影響都強[49]。U對AGB變化產生負效應,并且U對AGB的影響比Ta強(表4)。在中國半干旱和干旱地區,U對NDVI也產生負效應,并且NDVI對U的敏感性高于對Ta的敏感性[50]。

表4 在AGB增加和減少區域中,每個因子的標準化回歸系數的區域平均Table 4 The regional average of standardized regression coefficients of each factor in the area with increased and decreased AGB,respectively
在AGB上升和下降區域,氣候因子分別解釋了AGB變化的(59±18)%和(51±17)%。還有一部分未解釋的變化可能與土壤養分(例如,總碳、總氮、總磷等)、土壤酸堿度和人類活動的干擾(例如,放牧、道路建設等)有關[51-52]。此外,氣候變量的年總量和年平均值不能反映它們季節變化的年際差異。
基于AGB觀測值,使用BP-ANNs構建了遙感植被指數和其他環境因子與AGB的關系,展示了不同草地類型和海拔地區中AGB的時空變化。揭示了ANNs與機理模型模擬值的時空變化的差異。2001—2018年青藏高原草地AGB的區域平均為(74.75±52.65) g·m-2,呈微弱的增加趨勢。AGB增加的草地面積比例較高,為61%。AGB從草叢到草甸再到草原依次減少。AGB隨海拔的升高先增加后減少,3 400~3 800 m海拔區間內的AGB最高。AGB的增加趨勢隨海拔的升高而減弱。使用BP-ANNs能獲得與基于生物化學過程或光能利用率的機理模型的模擬精度相當的模擬值。大氣CO2濃度、VPD、前一年的PRE對AGB變化的影響最強,其次為U和Ta。