孟新月， 葛 靜， 侯蒙京， 馮琦勝， 金哲人， 高金龍， 梁天剛

(草地農(nóng)業(yè)生態(tài)系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室, 蘭州大學(xué)草地農(nóng)業(yè)科技學(xué)院， 甘肅 蘭州 730020)

植被蓋度是指地上植被占地表總面積的百分比，能夠很好地反映植被的生長(zhǎng)情況，也是生態(tài)系統(tǒng)變化的敏感指標(biāo)[1]。植被蓋度的大小能夠反映植被的茂密程度以及植被群落生長(zhǎng)動(dòng)態(tài)，同時(shí)也可以揭示區(qū)域生態(tài)環(huán)境的變化狀況，明確草原植被蓋度對(duì)評(píng)價(jià)草地資源、監(jiān)測(cè)草地退化狀況具有重要作用[2]。

地表實(shí)測(cè)法和遙感反演法是獲取植被蓋度的兩種廣泛使用的方法[3]。地表實(shí)測(cè)法以樣方尺度來測(cè)量，按其測(cè)量原理可分為目視估測(cè)法、采樣法及儀器法3類[4]。但是，傳統(tǒng)的地表實(shí)測(cè)法耗時(shí)長(zhǎng)、所需精力多、精度低，且只適用于較小的空間尺度[4-6]。近年來，“3S”(Remote Sensing,RS;Global Position System,GPS;Geographic Information System,GIS)技術(shù)迅猛發(fā)展，具有平臺(tái)廣闊、層次豐富、時(shí)相較長(zhǎng)等特點(diǎn)，在草地資源監(jiān)測(cè)中發(fā)揮著日益重要的作用[7]。吳見等[8]依據(jù)像元二分模型法收集了黃山市在各個(gè)年度的植物蓋度數(shù)據(jù)，并利用Landsat TM遙感技術(shù)影像，分析了黃山市植物蓋度空間變化特征，結(jié)果表明1988—2018年黃山市植被蓋度整體上呈現(xiàn)出增加趨勢(shì)；趙國忱等[9]利用Landsat-8數(shù)據(jù)和像元二分模型，并利用歸一化植被指數(shù)(Normalized difference vegetation index,NDVI)估算了遼寧北票市植被蓋度，結(jié)果表明2013—2015年間北票市東部為主要增長(zhǎng)區(qū)域，2015—2017年間，南部和北部有較大增長(zhǎng)。宋清潔等[10]通過地面實(shí)測(cè)數(shù)據(jù)結(jié)合中分辨率成像光譜儀(Moderate-resolution imaging spectroradiometer,MODIS)植被指數(shù)數(shù)據(jù)，分析了增強(qiáng)型植被指數(shù)(Enhanced vegetation index,EVI)和NDVI與草地蓋度的線性關(guān)系，篩選出甘南州草地蓋度最優(yōu)反演模型為基于EVI構(gòu)建的對(duì)數(shù)模型，其R2為0.707，RMSE為10.69%。但是，已有研究也存在草地蓋度反演誤差大、反演模型構(gòu)建方法有待改進(jìn)等方面的問題。雖然利用MODIS等遙感數(shù)據(jù)和像元二分模型方法可以快速模擬植被蓋度隨時(shí)間序列的動(dòng)態(tài)變化，但在估算蓋度時(shí)尚存在一定缺陷，如不能揭示植被覆蓋所體現(xiàn)的生態(tài)過程，并容易受純植被像元的參照光譜值的影響等[11]。由于單一的植被指數(shù)模型和多因素線性回歸模型精度較低、穩(wěn)定性差、存在局限性[12]，機(jī)器學(xué)習(xí)模型已逐步取代傳統(tǒng)算法成為構(gòu)建草地蓋度模型的主流方法。機(jī)器學(xué)習(xí)模型是基于數(shù)據(jù)驅(qū)動(dòng)的，其自動(dòng)檢索和解釋數(shù)據(jù)的方法比較靈活，可用于任何訓(xùn)練任務(wù)，準(zhǔn)確性更高。機(jī)器學(xué)習(xí)模型能夠估計(jì)多個(gè)變量間的復(fù)雜關(guān)系，與傳統(tǒng)模型相比更具穩(wěn)健性，能有效提高模型的預(yù)測(cè)精度[13-14]。為了改善傳統(tǒng)模型在實(shí)際應(yīng)用中存在的不足，一些學(xué)者使用機(jī)器學(xué)習(xí)模型在草地蓋度遙感反演方面進(jìn)行了初步探索，如陳黔等[15]使用Landsat數(shù)據(jù)與其他數(shù)據(jù)進(jìn)行北方沙地灌木覆蓋度的估算，研究表明利用機(jī)器學(xué)習(xí)算法可以很好的實(shí)現(xiàn)這一估算(R2=0.72，RMSE=13.73%)。Ge等[16]基于MODIS衛(wèi)星數(shù)據(jù)和野外實(shí)測(cè)蓋度數(shù)據(jù)，構(gòu)建了黃河源區(qū)4種遙感反演模型，結(jié)果表明支持向量機(jī)模型(Support vector machine,SVM)為草地蓋度反演的最優(yōu)模型，比單因素模型的R2提高0.08～0.14。然而，機(jī)器學(xué)習(xí)算法需要大量草地蓋度觀測(cè)數(shù)據(jù)和多種生態(tài)環(huán)境變量，存在變量篩選、模型調(diào)參等多個(gè)環(huán)節(jié)的大量分析，研究結(jié)果受諸多因素的影響，因此在天然草地蓋度最優(yōu)反演算法研究方面仍需要進(jìn)行深入探索。青藏高原海拔高，地貌豐富，面積廣大，草地資源異常豐富，由于自然和人為的雙重影響，青藏高原一些地區(qū)草地資源呈退化趨勢(shì)。我國已有多位學(xué)者針對(duì)青藏高原草地蓋度變化開展了大量研究。如邵偉等[17]利用統(tǒng)計(jì)資料分析表明西藏高原草地退化嚴(yán)重。高清竹等[18]基于遙感數(shù)據(jù)對(duì)藏北地區(qū)草地退化進(jìn)行遙感監(jiān)測(cè)，結(jié)果表明藏北地區(qū)1982—2004年草地退化較為嚴(yán)重，重度退化草地面積占草地總面積的8%。陸晴等[19]利用NDVI數(shù)據(jù)和地面實(shí)測(cè)數(shù)據(jù)，對(duì)青藏高原1982—2013年高寒草地覆蓋時(shí)空變化進(jìn)行了研究，結(jié)果表明青藏高原高寒草地生長(zhǎng)季NDVI表現(xiàn)為從東南到西北逐漸減少的趨勢(shì)。但前人研究范圍大多限于青藏高原部分區(qū)域，從整體上系統(tǒng)進(jìn)行青藏高原蓋度研究較少。

基于以上考慮，本研究利用2003—2018年草地植被生長(zhǎng)季野外實(shí)測(cè)數(shù)據(jù)，首先用最小絕對(duì)壓縮變量篩選方法(Least absolute shrinkage and selection operator，LASSO)篩選合適的變量，再分別建立單因素回歸模型和多因素機(jī)器學(xué)習(xí)模型，基于最優(yōu)模型進(jìn)一步分析青藏高原地區(qū)2001—2019年生長(zhǎng)季草地蓋度的時(shí)空變化特征。本研究旨在通過分析草地植被蓋度的時(shí)空變化特征，為青藏高原地區(qū)草地資源的長(zhǎng)期利用和動(dòng)態(tài)監(jiān)測(cè)提供科學(xué)依據(jù)[20]。

1 數(shù)據(jù)及方法

1.1 研究區(qū)概況

青藏高原(26°50′～39°19′N，78°25′～103°04′E)包括西藏自治區(qū)、青海省和四川省、甘肅省、新疆維吾爾自治區(qū)及云南省的部分地區(qū)，地域遼闊，平均海拔在4 000 m以上，有“世界屋脊”之稱[21]。氣候具有輻射強(qiáng)、晝夜溫差大等特點(diǎn)，年平均氣溫范圍為-6℃～3℃；降水多集中在6—9月，占全年降水量的70%左右，且東西部降水量差距較大，東南部降水最高可達(dá)1 000 mm，而西北部降水多為40～100 mm[22]。青藏高原草地總面積約為121 742 km2，占高原總面積的60%[23]。因獨(dú)特的氣候，高寒草甸、高寒草原、高寒荒漠為青藏高原主要的三種草地類型(圖1)[24]。

圖1 研究區(qū)草地類型及外業(yè)調(diào)查樣點(diǎn)空間分布Fig.1 Grassland types and spatial distribution of field investigation samples in the study area

1.2 樣地?cái)?shù)據(jù)獲取

地面實(shí)測(cè)數(shù)據(jù)的獲取時(shí)間為2003—2018年的7—9月，主要獲取了生長(zhǎng)季的草地植被蓋度數(shù)據(jù)，總計(jì)獲得4 376個(gè)樣點(diǎn)(圖1)。在研究區(qū)內(nèi)設(shè)定面積為100 m×100 m的樣地，選擇植被生長(zhǎng)狀況均一、地勢(shì)平坦、具有代表性的典型群落地段[25]。依據(jù)樣地內(nèi)草地植被的變化狀況，每個(gè)樣地選擇3～5個(gè)樣方采用目測(cè)法估計(jì)草地植被的蓋度。在采樣過程中，記錄通過GPS定位的每一樣方中心點(diǎn)的經(jīng)緯度、高程和優(yōu)勢(shì)種、草地植被蓋度、草層高度等指標(biāo)[26]。一個(gè)樣地的草地蓋度值由5個(gè)樣方觀測(cè)值的平均值代表。剔除明顯異常的數(shù)據(jù)，使數(shù)據(jù)采樣點(diǎn)信息一致，最終得到有效樣本記錄4 355條用于模型構(gòu)建(表1)。

表1 2003-2018年青藏高原地區(qū)天然草地蓋度統(tǒng)計(jì)結(jié)果Table 1 Statistical results of natural grassland coverage on the Tibetan Plateau from 2003 to 2018

1.3 MODIS數(shù)據(jù)獲取及處理

本研究下載了覆蓋整個(gè)研究區(qū)的2001—2019年5—9月的MCD43A4反射率產(chǎn)品，時(shí)間分辨率為每天，空間分辨率為500 m。本研究共利用7個(gè)波段，每個(gè)波段所對(duì)應(yīng)的波譜范圍依次是620～670 nm(Band1)，841～876 nm(Band2)，459～479 nm(Band3)，545～565 nm(Band4)，1 230～1 250 nm(Band5)，1 628～1 652 nm(Band6)和2 105～2 155 nm(Band7)。

利用MODIS數(shù)據(jù)重投影工具進(jìn)行數(shù)據(jù)格式轉(zhuǎn)換和投影定義，將數(shù)據(jù)處理為TIFF格式和WGS84投影。在ArcMap 10.2軟件中利用extract by mask工具裁剪出研究區(qū)范圍，用于后續(xù)分析，同時(shí)，采用最大值合成法(Maximum value composite,MVC)，得到了月最大反射率數(shù)據(jù)集[27]。

1.4 生態(tài)環(huán)境指標(biāo)及預(yù)處理

地形數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http://www.gscloud.cn)，下載航天飛機(jī)雷達(dá)地形任務(wù)(SRTM) V4.1 TIFF數(shù)據(jù)，將SRTM數(shù)字高程模型(DEM)數(shù)據(jù)重采樣至500 m空間分辨率，以和MCD43A4產(chǎn)品匹配。利用ArcMap 10.2軟件計(jì)算得到研究區(qū)經(jīng)度、緯度、海拔、坡向和坡度的空間分布數(shù)據(jù)[28]。

土壤數(shù)據(jù)來源于中國土壤特征數(shù)據(jù)集(http://globalchange.bnu.edu.cn/research/soil)，其中包括土壤深度0～30 cm的堿解氮(Alkali-hydrolysable N,AN)、有效鉀(Available K,AK)、有效磷(Available P,AP)、容重(Bulk density,BD)、砂土含量(Clay fraction,CL)、礫石含量(Gravel,GRAV)、PH值(PH)、孔隙度(Porosity,POR)、含泥量(Silt fraction,SI)、含砂量(Sand fraction,SA)、土壤有機(jī)質(zhì)(Soil organic matter,SOM)、全鉀(Total potassium,TK)、全氮(Total nitrogen,TN)及全磷(Total phosphorus,TP)。同時(shí)使用本研究區(qū)邊界對(duì)土壤數(shù)據(jù)進(jìn)行裁剪，以便進(jìn)一步分析。

溫度和降水?dāng)?shù)據(jù)來源于國家地球系統(tǒng)科學(xué)數(shù)據(jù)中心共享服務(wù)平臺(tái)(http://www.geodata.cn)。為了與遙感植被指數(shù)產(chǎn)品MCD43A4的時(shí)間分辨率相匹配，對(duì)應(yīng)該產(chǎn)品周期將日數(shù)據(jù)處理為月累積降水和月平均溫度。利用該數(shù)據(jù)使用插值軟件ANUSPLIN Version 4.3 進(jìn)行空間插值?？臻g插值數(shù)據(jù)選擇500 m空間分辨率和Albers投影。擬合降水表面時(shí)，將經(jīng)度、緯度和海拔作為三個(gè)獨(dú)立的樣條變量；擬合溫度表面時(shí)，高程作為獨(dú)立協(xié)變量。插值結(jié)果通過裁剪得到本研究區(qū)的氣候空間數(shù)據(jù)集，進(jìn)而提取采樣地對(duì)應(yīng)的降水和溫度值。

1.5 變量及篩選

本研究所涉及的變量包括4類，共55個(gè)。其中，遙感植被指數(shù)有31個(gè)(表2)，氣象指標(biāo)包括年降水量(prec)、年平均溫度(tmp)、大于0度年積溫(GDD)、1月至采樣月累積降水(per_samp)和1月至采樣月累積溫度(tem_samp)，地形指標(biāo)包括海拔(DEM)、坡位(TPI)、坡向(Aspect)、曲率(Curvature)和坡度(Slope)，土壤指標(biāo)包括AK，AN，AP，BD，CL，GRAV，pH，POR，SA，SI，SOM，TK，TN及TP。為了精簡(jiǎn)變量，使用最小絕對(duì)壓縮變量篩選(Least absolute shrinkage and selection operator，LASSO)方法來選擇敏感變量。這種方法可以顯著減少變量的數(shù)量，降低共線性，提高模型精度，還可以消除模型的冗余性，精簡(jiǎn)模型。LASSO回歸分析是基于最小二乘法的原理，利用模型中各變量的回歸系數(shù)構(gòu)造一個(gè)函數(shù)作為懲罰函數(shù)，并對(duì)各變量的系數(shù)進(jìn)行壓縮，直到模型的殘差平方和最小，去掉系數(shù)降為0的變量。因此，LASSO回歸分析方法被廣泛應(yīng)用于高維數(shù)據(jù)的變量篩選[29-30]。

表2 基于MODIS數(shù)據(jù)的草地植被指數(shù)變量Table 2 Grassland vegetation index variables based on MODIS data

1.6 模型構(gòu)建與分析

本研究采用單因素回歸、人工神經(jīng)網(wǎng)絡(luò)(Artificial neural network,ANN)、SVM和隨機(jī)森林(Random forest,RF)4種模型對(duì)草地蓋度進(jìn)行建模。根據(jù)帕累托原則[31-32]，首先將4 355個(gè)樣本數(shù)據(jù)以8∶2的比例分為兩部分，其中3 484個(gè)樣本用于建模，871個(gè)樣本用于空間模擬結(jié)果的精度驗(yàn)證，再將3 484個(gè)樣本數(shù)據(jù)以8∶2的比例分為兩部分，其中2 787個(gè)樣本作為訓(xùn)練集數(shù)據(jù)，697個(gè)樣本作為測(cè)試集數(shù)據(jù)(圖1)。

1.6.1傳統(tǒng)回歸模型 本項(xiàng)研究首先使用回歸分析法構(gòu)建草地蓋度模型，利用LASSO篩選出的變量作為自變量，分別與草地蓋度建立對(duì)數(shù)、線性、乘冪和指數(shù)4類回歸模型[27]。

1.6.2機(jī)器學(xué)習(xí)模型 本文共采用3種機(jī)器學(xué)習(xí)模型，分別為RF，SVM，ANN模型。RF模型是基于分類樹算法的一種算法。RF模型利用bootstrap進(jìn)行抽樣，抽出的樣本用于回歸樹的構(gòu)建。對(duì)訓(xùn)練樣本進(jìn)行連續(xù)篩選，得到最小殘差平方和，最后形成一棵完整的樹[33]。本研究在MATLAB 2019a中，通過Tree Bagger工具包構(gòu)建RF模型，主要測(cè)試決策樹的數(shù)目，在ntree中設(shè)置以100為間隔，從100開始增加到2 000，共迭代20次，當(dāng)十折交叉驗(yàn)證誤差最小時(shí)，確定為最優(yōu)ntree，用于最終模型的構(gòu)建。SVM是一種具有相關(guān)學(xué)習(xí)算法的監(jiān)督學(xué)習(xí)模型，由高維或無限維空間中的一組超平面構(gòu)成，可以用于分類、回歸和其他任務(wù)。本研究在MATLAB 2019a中構(gòu)建SVM回歸，選擇徑向基函數(shù)(RBF)作為核函數(shù)，采用遺傳算法(GA)智能尋找最佳cost (bestc)和gamma (bestg)參數(shù)。ANN是由輸入層、輸出層和一個(gè)或多個(gè)隱含層組成的多層網(wǎng)絡(luò)結(jié)構(gòu)[34]。本研究通過MATLAB軟件的神經(jīng)網(wǎng)絡(luò)工具箱完成BP神經(jīng)網(wǎng)絡(luò)建模及其驗(yàn)證，設(shè)置隱含層神經(jīng)元個(gè)數(shù)以10為間隔，從10增加到100，共迭代10次，取交叉驗(yàn)證誤差最小時(shí)的neuron個(gè)數(shù)為最優(yōu)neuron，構(gòu)建最終模型。

1.6.3模型評(píng)價(jià)指標(biāo) 模型的適應(yīng)度由訓(xùn)練數(shù)據(jù)集的決定系數(shù)(R2)來衡量，預(yù)測(cè)精度的評(píng)價(jià)指標(biāo)為相關(guān)系數(shù)(r)、均方根誤差(RMSE)和均方根誤差變異系數(shù)(CVRMSE)[35-36]。具體計(jì)算公式如下：

式中，E(yi)和yi分別為y的實(shí)測(cè)值和預(yù)測(cè)值，n等于樣本量(這里n=697)。RMSE越低，模型精度越高[37]。

1.7 草地蓋度空間分布與動(dòng)態(tài)變化

2001—2019年最大草地蓋度的時(shí)間變化采用坡度線性趨勢(shì)模型分析，計(jì)算各像元的植被蓋度變化[40-41]。當(dāng)slope>0時(shí)，草地蓋度呈增加趨勢(shì)；當(dāng)slope<0時(shí)，草地蓋度呈降低趨勢(shì)。計(jì)算公式為：

式中，n等于總年數(shù)(19)；i是2001年至2019年的年份序號(hào)(即1～19)；其中，coveri為最優(yōu)模型模擬的第i年年最大草地蓋度。

利用F檢驗(yàn)分析草地變化的顯著性[42]。如果F>F0.05(1,n-2)，在95%置信水平下變化顯著，此處F0.05(1,14)=4.60。研究區(qū)草地蓋度變化分為4類：顯著增加(slope>0,F>4.60)，增加(slope>0,F<4.60)，減少(slope<0,F<4.60)和顯著減少(slope<0,F>4.60)。F檢驗(yàn)公式為:

式中，n等于總年數(shù)19，R2為每個(gè)像元的草地蓋度與時(shí)間序列的相關(guān)系數(shù)r的平方；是1～19的平均值，此處等于10。coveri為最優(yōu)模型模擬的第i年年最大草地蓋度，其中，為2001—2019年平均最大草地蓋度。

利用RF最優(yōu)模型，反演出每年生長(zhǎng)季草地蓋度最大值。結(jié)合slope線性趨勢(shì)分析和F檢驗(yàn)分析2001—2019年草地蓋度的變化。在ArcMap 10.2中，利用Reclassify工具對(duì)19年草地蓋度增減情況進(jìn)行分級(jí)，分為增加、顯著增加、減少、顯著減少四級(jí)[25]。

2 結(jié)果與分析

2.1 變量篩選結(jié)果

通過LASSO篩選，共選出10個(gè)敏感變量用于天然草地植被蓋度模型的構(gòu)建與分析。其中，包括2個(gè)植被指數(shù)，分別為L(zhǎng)SWI，RVI；2個(gè)氣象指標(biāo)，分別為prec和tem_samp；1個(gè)地形參數(shù)，為Slope；5個(gè)土壤屬性指標(biāo)，分別為AN，AP，BD，POR，TP。由此可見，土壤、地形、氣象和植被指數(shù)這四類變量與草地蓋度均具有密切關(guān)系。

2.2 單因子模型分析

表3是利用LASSO篩選出的變量建立的回歸模型及精度分析結(jié)果。由表3可以看出，在所有草地蓋度單因子參數(shù)模型中，基于LSWI和RVI的單因子模型精度優(yōu)于其他模型，其R2均為0.52，RMSE介于15.56%～18.21%；其次為基于氣象指標(biāo)的單因子模型，R2介于0.02～0.39，RMSE介于17.77%～23.06%；而基于土壤、地形等因素的草地蓋度單因子參數(shù)模型精度均較低，R2介于0.1～0.26，RMSE介于19.35%～23.12%。

表3 草地蓋度與植被指數(shù)的單因子模型精度評(píng)價(jià)結(jié)果Table 3 Evaluation results of single-factor model accuracy of grassland coverage and vegetation index

LSWI和RVI與草地植被蓋度之間有較好的線性相關(guān)關(guān)系?；赗VI的線性模型的擬合決定系數(shù)(R2=0.52)高于指數(shù)模型(R2=0.50)、對(duì)數(shù)模型(R2=0.47)及乘冪模型(R2=0.42)，基于LSWI的線性模型的擬合決定系數(shù)(R2=0.52)高于指數(shù)模型(R2=0.47)，而基于RVI的均方根誤差低于LSWI。由此可見，基于RVI的線性模型為青藏高原草地植被蓋度的最優(yōu)單因子模型。

2.3 機(jī)器學(xué)習(xí)模型分析

從分析結(jié)果(表4)可以看出，3種機(jī)器學(xué)習(xí)模型精度均高于單因子最優(yōu)模型，R2提高了0.09～0.16，RMSE降低了1.52%～2.81%。RF模型的決定系數(shù)(R2=0.68)高于SVM模型(R2=0.66)和ANN模型(R2=0.61)，而RF模型的均方根誤差最低(RMSE=12.75%)，同時(shí)RF模型的CVRMSE為16.81%，代表其有較好的預(yù)測(cè)能力。研究結(jié)果說明，本研究建立的RF模型預(yù)測(cè)值與實(shí)地測(cè)量的草地蓋度值非常接近，可見RF模型優(yōu)于其他模型，為青藏高原草地植被蓋度的最優(yōu)估測(cè)模型。

2.4 草地蓋度空間動(dòng)態(tài)變化分析

圖2為利用RF模型反演的2001—2019年草地年最大蓋度的平均結(jié)果?？傮w上看，青藏高原草地植被覆蓋狀況較好，東部地區(qū)草地覆蓋度高，西部地區(qū)覆蓋度低。低蓋度(<40%)草地主要分布在西藏自治區(qū)的中部高海拔地區(qū)，其面積占比為16%；較低蓋度(40%～60%)草地多分布在青海省西部和西藏自治區(qū)的周邊，其面積占比為33%；較高蓋度(60%～80%)草地則主要分布在青藏高原中部地區(qū)，與中低蓋度草地交錯(cuò)分布，其面積占比為36%；高蓋度(>80%)草地多分布在四川省、甘肅省和青海省的東部地區(qū)，面積占比為15%。高和較高植被覆蓋區(qū)域面積占青藏高原草地總面積的15%和36%，共占青藏高原草地總面積的51%，說明青藏高原植被生長(zhǎng)季大部分區(qū)域的草地植被覆蓋率較高，且植被覆蓋狀況良好。

表4 青藏高原草地植被蓋度機(jī)器學(xué)習(xí)模型精度評(píng)價(jià)結(jié)果Table 4 Accuracy evaluation results of machine learning model for grassland vegetation coverage in Qinghai-Tibet Plateau

利用最大值合成法，通過對(duì)7—9月的最大草地蓋度趨勢(shì)分析的結(jié)果表明(圖3)，2001—2019年青藏高原大部分地區(qū)草地蓋度呈現(xiàn)增加趨勢(shì)，呈增加趨勢(shì)區(qū)域占比為55.4%，減少區(qū)域占比為44.6%。草地蓋度減少大部分發(fā)生于研究區(qū)中心地帶，少數(shù)分布于四川省部分地區(qū)；草地蓋度顯著減少的地區(qū)主要在西藏自治區(qū)的東南部地區(qū)和西藏自治區(qū)北部的少部分地區(qū)，面積占比為5.2%；青藏高原的東、西部地區(qū)多為草地蓋度增加區(qū)域；青海省的東北地區(qū)為蓋度顯著增加區(qū)域，另有少數(shù)分布在西藏自治區(qū)的西南部地區(qū)，面積占比為8.1%。

圖3 2001—2019年青藏高原草地蓋度變化趨勢(shì)Fig.3 Trends of grassland coverage in the Tibetan Plateau from 2001 to 2019

3 討論

本研究對(duì)比分析了單變量參數(shù)模型與非參數(shù)模型性能，確定了青藏高原地區(qū)最優(yōu)草地蓋度反演模型。單因素模型統(tǒng)計(jì)分析的結(jié)果表明，土壤、地形、氣象和植被指數(shù)這4類變量均與草地蓋度具有密切的關(guān)系。土壤、地形、氣象因素對(duì)草地蓋度有一定程度的影響，土壤屬性指標(biāo)與草地蓋度的R2介于0.01～0.26，地形參數(shù)R2為0.02，氣象因素R2達(dá)0.38。由此可見，單因素模型不能反映多種因素的綜合影響，因此在大范圍的草地蓋度分析時(shí)基于單因素構(gòu)建的回歸模型的誤差偏大。本研究采用了3種常見的機(jī)器學(xué)習(xí)模型。RF模型對(duì)參數(shù)的調(diào)整變化并不敏感，具有極佳的性能。經(jīng)研究發(fā)現(xiàn)，噪聲數(shù)據(jù)對(duì)RF模型影響較小，具有魯棒性[43-45]。與ANN，SVM模型相比，RF模型的訓(xùn)練集和測(cè)試集精度存在一定差異。SVM模型更適合處理小樣本數(shù)據(jù)，但當(dāng)樣本量過大時(shí)，SVM模型需耗費(fèi)更多運(yùn)算時(shí)間用于參數(shù)調(diào)整。參數(shù)調(diào)整對(duì)ANN模型的影響較大。最優(yōu)機(jī)器學(xué)習(xí)的RF模型自變量包括LSWI和RVI 2個(gè)植被指數(shù)，prec和tem_samp 2個(gè)氣象指標(biāo)，1個(gè)坡度地形參數(shù)，以及AN，AP，BD，POR，TP5個(gè)土壤屬性指標(biāo)。這些變量涵蓋了可代表草地生態(tài)環(huán)境的多種屬性，各指數(shù)之間具有優(yōu)勢(shì)互補(bǔ)的特點(diǎn)，因此利用LASSO變量篩選方法選出的最優(yōu)變量組合方式構(gòu)建的RF模型可以更好的反演整個(gè)青藏高原天然草地蓋度的變化。然而，基于RF的最優(yōu)模型仍然存在一定的局限性和不確定性。首先，青藏高原地區(qū)草地面積較大、地形復(fù)雜，大部分樣地分布在高原東部地區(qū)，而中西部地區(qū)由于受道路和海拔的限制，樣地?cái)?shù)量有限，而地面采樣點(diǎn)和遙感數(shù)據(jù)的時(shí)間匹配也存在很大不確定性，使得構(gòu)建的模型存在一定的誤差。同時(shí)，青藏高原氣象站分布不均勻，氣象資料的空間插值存在一定誤差[46]。其次，與單變量參數(shù)模型相比，RF模型算法更靈活，由基于高維數(shù)據(jù)訓(xùn)練的大樣本決策樹組成，具有較強(qiáng)的數(shù)據(jù)誤差容忍度，但通常需要大量的標(biāo)記和地面測(cè)量數(shù)據(jù)。因此，模型仍有一定的局限性和不確定性。

本研究分析了2001—2019年青藏高原草地蓋度空間分布格局與動(dòng)態(tài)變化特征，針對(duì)青藏高原草地蓋度變化，學(xué)者已展開了大量研究。如唐志光等[47]研究表明三江源地區(qū)草地蓋度空間分布呈現(xiàn)出東高西低的特點(diǎn)；于伯華等[48]的研究結(jié)果表明，雅魯藏布江流域草地蓋度總體上呈現(xiàn)出上升趨勢(shì)。丁明軍等[49]基于NDVI反演了1982—2009年青藏高原草地植被蓋度的年際變化，高原大部分地區(qū)草地蓋度呈現(xiàn)增加趨勢(shì)。本研究結(jié)果與前人在該地區(qū)所做研究結(jié)果基本一致，均得出該地區(qū)2001—2010年內(nèi)草地蓋度呈現(xiàn)整體增加的趨勢(shì)。馬琳雅等[50]利用MODIS植被指數(shù)數(shù)據(jù)，分析了甘南州2001—2011年草地蓋度空間變化特征，基于MODIS-EVI的對(duì)數(shù)模型為最優(yōu)模型，R2為0.47。本研究與之相比精度較高，這得益于我們不僅考慮了遙感植被指數(shù)，還充分考慮了土壤、地形、氣象等變量，能更好的反演草地蓋度的時(shí)空變化。草地蓋度不僅受氣候變化的影響，同時(shí)人為因素也會(huì)對(duì)草地植被蓋度產(chǎn)生很大的影響。2001—2019年青藏高原地區(qū)草地年平均最大蓋度空間分布整體上呈現(xiàn)自西向東、自北向南遞增的趨勢(shì)，并且增加區(qū)域大于減少區(qū)域。這與近年來國家加強(qiáng)對(duì)草地退化的重視有關(guān)。自1999年以來，國家先后實(shí)施了一系列以草原綜合治理和恢復(fù)為重點(diǎn)的工程措施，如退牧還草、草原生態(tài)保護(hù)補(bǔ)助獎(jiǎng)勵(lì)等措施[51]。因此，合理的人為干預(yù)對(duì)草地植被的恢復(fù)具有顯著的作用。

在本項(xiàng)研究基礎(chǔ)上，為提高模型精度，可通過科學(xué)手段改進(jìn)現(xiàn)有模型，如:合理擴(kuò)大觀測(cè)范圍、增加更多采樣點(diǎn)以提高衛(wèi)星影像的空間匹配性等。另外，傳統(tǒng)的地表實(shí)測(cè)法費(fèi)時(shí)、費(fèi)力、精度低，且不能反映較大空間尺度的植被結(jié)構(gòu)和空間變化等信息，近年來無人機(jī)(Unmanned aerial vehicle,UAV)技術(shù)普遍應(yīng)用于地表觀測(cè)，具有體積小、質(zhì)量輕、靈活性高、可在特殊地區(qū)探測(cè)等優(yōu)點(diǎn)[52-53]，因此利用基于UAV的大范圍多樣點(diǎn)的草地蓋度快速監(jiān)測(cè)方法，增加建模樣本的數(shù)量和時(shí)空代表性，是未來改進(jìn)研究區(qū)天然草地蓋度反演模型的重要研究方向。此外，青藏高原由于地形復(fù)雜，東西部氣溫、降水變化強(qiáng)烈，受氣候等諸多因素的影響，東西部的人口也存在巨大差距，氣象因子和土壤、地形因子也存在一定的空間差異性。因此，深入研究和分析氣象因素和人為因素對(duì)草地蓋度的影響也是今后研究的重要內(nèi)容。

4 結(jié)論

本研究基于2003—2018年青藏高原地區(qū)地形、土壤、氣候和植被指數(shù)數(shù)據(jù)，對(duì)比分析了青藏高原天然草地植被蓋度單因素回歸模型和非參數(shù)模型的精度，研究了近19年(2001—2019)青藏高原草地蓋度時(shí)空動(dòng)態(tài)變化。結(jié)果表明，在單因子草地蓋度遙感模型中，MODIS RVI與草地植被蓋度的相關(guān)性最好，基于RVI構(gòu)建的線性模型的R2為0.52，RMSE為15.56；在3種機(jī)器學(xué)習(xí)方法中，RF算法是青藏高原草地植被蓋度反演的最優(yōu)模型，R2達(dá)0.68，RMSE為12.75；RF模型的自變量包括2個(gè)植被指數(shù)、2個(gè)氣象指標(biāo)、1個(gè)地形參數(shù)及5個(gè)土壤屬性指標(biāo)，不僅有基于紅光和近紅外波段構(gòu)建的傳統(tǒng)植被指數(shù)RVI，也包括可以反映干旱、半干旱區(qū)域植被生長(zhǎng)狀況的LSWI，這些變量涵蓋了優(yōu)勢(shì)互補(bǔ)的植被指數(shù)，因此RF模型在反演類型多樣、地形、氣候和土壤異質(zhì)性強(qiáng)的青藏高原天然草地蓋度變化方面具有獨(dú)特的優(yōu)勢(shì)；2001—2019年間青藏高原地區(qū)草地植被蓋度整體呈現(xiàn)出由西向東、由北到南的增加趨勢(shì)，55.4%區(qū)域呈增加趨勢(shì)，44.6%區(qū)域呈減少趨勢(shì)。