基于GEE云平臺的榆林市植被覆蓋度時空變化研究

白欣，賀炳彥，高婷，李雪，孫喬

（1.長安大學地質工程與測繪學院，西安 710000；2.中煤航測遙感集團有限公司，西安 710000；3.西安靈圖信息科技有限公司，西安 710000；4.航天建筑設計研究院有限公司，北京 100162）

植被作為生態大循環的重要參與者，在全球能量流動和物質循環中起著重要作用［1］。通常用植被覆蓋度（Fractional vegetation cover，FVC）表征地表植被覆蓋情況，一般將植被覆蓋度定義為植被冠層在地面上垂直投影面積與土地總面積的比值［2］。國內外學者常通過像元二分模型來估算植被覆蓋度，多使用歸一化植被指數來分解像元遙感信息，從而監測研究區的植被覆蓋變化。張舒婷等［3］采用Landsat數據，利用像元二分法估算植被覆蓋度，并監測了黃土高原丘陵溝壑區18年的植被覆蓋變化；張志強等［4］采用像元二分模型、Hurst指數分析了2000—2019年黃河流域的植被覆蓋度時空變化特征和發展趨勢。

榆林市位于黃土高原腹地，水土流失嚴重，是中國典型的生態脆弱區。針對黃土高原風蝕沙化嚴重、土壤鹽堿化等問題，從1999年開始，國家實施退耕還林還草措施，以期植被較快恢復，國內許多學者也開始關注榆林市退耕還林還草后的植被覆蓋狀況。楊波等［5］使用MODIS數據研究了榆林市自還林還草工程以來的植被變化和未來發展趨勢，發現18年來榆林市歸一化植被指數（Normalized difference vegetation index，NDVI）共增長了88.56%，未來有10%的增長潛力；馬鋒等［6］研究發現榆林市2000—2018年植被覆蓋度呈波動上升的趨勢，且南部黃土丘陵溝壑區的植被覆蓋度增速高于北部風沙草灘區。但少有學者使用高分辨率數據對榆林植被覆蓋進行長時序的監測，主要原因是植被覆蓋度計算過程中高分辨率數據存儲不便、計算量大、運行時間長等問題。

Google earth engine（GEE）平臺的出現可以有效解決上述問題。GEE是一個基于云計算的“遙感大數據”分析和呈現平臺，相較于單機地理信息處理平臺而言，GEE不僅有大量的影像數據可供用戶使用，還提供了高性能的集群服務器，實現海量影像數據在線可視化處理，大大提高了工作效率［7］。因此，本文基于GEE平臺，針對榆林市的生態環境問題和還林還草工程，采用Landsat surface reflectance（陸地衛星地表反射率數據）和像元二分模型反演了2000—2019年榆林市植被覆蓋度，分析植被覆蓋時空變化特征，同時分析了氣候因子（年平均氣溫、年降水量）和地形因子（高程、坡度）對植被覆蓋變化的影響，最后通過Hurst指數預測了未來研究區植被覆蓋的變化趨勢，為評估榆林市還林還草工程提供科學依據。

1 研究區與數據

1.1 研究區概況

榆林市地處陜西省最北部、陜北黃土高原與毛烏素沙地南緣交界處，東隔黃河與山西相望，南與陜西省延安市接壤，西鄰寧夏、甘肅，北連內蒙古，高程530～1 922 m（圖1），面積約43 578 km2。研究區屬溫帶大陸性季風氣候，年均溫度為10℃，年平均降水量為400 mm左右，氣象災害較多，立地條件較差［8］。以長城為界，北邊為風沙草灘區，南部為黃土丘陵溝壑區。土壤多以風沙土、黃綿土為主，水土流失和風蝕沙化嚴重，是退耕還林還草工程的典型區域。

圖1 榆林市位置

1.2 數據來源及預處理

本研究采用的遙感影像數據為Landsat surface reflectance，空間分辨率30 m、時間分辨率16 d，來源于GEE遙感計算云平臺。首先，在GEE平臺上調用“LANDSAT/LT05/C01/T1_SR”“LANDSAT/LE07/C01/T1_SR”“LANDSAT/LC08/C01/T1_SR”影像數據集，該產品已經進行了大氣校正［9，10］。利用Landsat質量評估波段（QA）檢測了像元的云、云影、雪等干擾因素，篩選出了3 042景低云量影像。依據GEE平臺上的CFMASK算法去云處理后，在線計算歸一化植被指數，采用最大值合成方法合成年度NDVI最大值影像，在此基礎上選用像元二分模型估算植被覆蓋度，最后采用qualityMosaic算法在線裁剪影像。

2000—2019年榆林市年平均氣溫和年降水量數據均來源于國家地球系統科學數據中心（http：//www.geodata.cn），空間分辨率為1 km。將該數據上傳至GEE平臺，用上述qualityMosaic算法在線裁剪，得到榆林市每年的年平均氣溫和年降水量分布數據。從地理空間數據云上獲取坡度和高程數據，空間分辨率為30 m。

2 研究方法

2.1 植被覆蓋度估算

植被指數通常用歸一化植被指數值來表示，用來衡量植被生長狀況和空間分布，是最常用的一種植被指數，其計算方法是近紅外波段的反射值（NIR）和紅外波段反射值（RED）之差與其之和的比值，計算公式如下：

采用像元二分模型估算植被覆蓋度，計算公式如下：

式中，NDVIsoil和NDVIveg分別為純裸土和純植被像元的NDVI值，它們的理論值應為0和1，但受地表濕度、土壤類型、土壤顏色等的影響，其實際值會隨空間和時間而變化［11］。因此，在年最大NDVI累積頻率表上取5%和95%作為NDVIsoil和NDVIveg，然后代入式（2）中，計算榆林市2000—2019年植被覆蓋度［12］。根據水利部2008年頒布施行的《土壤侵蝕分類分級標準》［13］，本研究將植被覆蓋度分為5個等級：低植被覆蓋度（0≤FVC＜30%）、中低植被覆蓋度（30%≤FVC＜45%）、中 植 被 覆 蓋 度（45%≤FVC＜60%）、中高植被覆蓋度（60%≤FVC＜75%）、高植被覆蓋度（75%≤FVC≤100%）。

2.2 植被覆蓋度變化趨勢分析

為了分析研究區植被覆蓋度在過去20年的變化趨勢及其時空特征，采用Sen+Mann-Kendall法分析FVC長時間序列變化趨勢。Sen斜率估計用于計算趨勢值，其計算公式如下［14］：

式中，xi、xj分別代表第i、第j年的植被覆蓋度；β是斜率，用于判斷時間序列的上升或下降，當β小于0時表示其呈現下降趨勢，相反表示上升。

Mann-Kendall屬于非參數檢驗方法，與其他參數檢驗的方法相比，樣本受異常值干擾小，不需要遵從一定的分布，更適合順序變量［15］。其檢驗統計量S的計算方法如下［16］

本文所用時間序列長度n值為20，故趨勢檢驗使用檢驗統計量Z，計算方法見式（6），取顯著水平α=0.05，Z1-α/2=Z0.975=1.96。

2.3 相關性分析

偏相關分析是當兩個變量與第三個變量相關時，控制第三個變量不變，只分析剩余兩個變量之間的相關性。選取年平均氣溫和年降水量作為植被覆蓋度的相關變量，分別計算其與植被覆蓋度的偏相關系數，公式如下［18］：

式中，rxy為變量x、y的相關系數，xi和yi是第i年的值，xˉ和yˉ分別是兩個變量的平均值。rxy*z為控制自變量Z不變時因變量x與自變量y的偏相關系數。顯著性檢驗采用t檢驗法，統計量t的計算公式如下：

式中，rxy*z和rxz分別為偏相關系數和相關系數，n為年數，m為自變量個數。

2.4 Hurst指數

Hurst指數是描述自相似性和長期依賴性的有效方法，并在水文、氣候、地質等領域廣泛運用［19］。本文采用Hurst指數來描述研究區FVC未來的變化趨勢，其具體計算過程參考嚴恩萍等［20］的方法。Hurst指數取值包括3種形式：當0.5≤Hurst時，表示時間序列具有長期持續性，未來的變化趨勢與過去相同；當0≤Hrust＜0.5時，說明時間序列具有反持續性，未來的變化趨勢與過去相反；當Hrust=0.5，表明時間序列是相互獨立的，無相關性。

3 結果與分析

3.1 植被覆蓋度空間分布特征

為了解榆林市植被覆蓋度空間分布情況，利用2000—2019年的最大植被覆蓋度，合成了榆林市20年平均植被覆蓋度分布圖，如圖2a所示。榆林市植被覆蓋度由東南向西北逐漸遞減，這與榆林市北部風沙草灘區、南部黃土高原丘陵溝壑區的地貌有關。低植被覆蓋度區域和中低植被覆蓋度區域主要分布在北部風沙草灘區，該區域地處毛烏素沙地南緣地帶，沙化風蝕嚴重，植被恢復困難。中植被覆蓋度區域、中高植被覆蓋度區域和高植被覆蓋度區域主要分布在南部黃土丘陵溝壑區，該地區空氣濕度、土壤肥力等自然環境較風沙區更好，加之退耕還林工程的實施，使得植被覆蓋更高。統計各等級植被覆蓋度所占面積比例，見圖2b。榆林市中植被覆蓋度所占面積比最大，占總面積的25.1%。

圖2 2000—2019年榆林市平均FVC空間分布

3.2 植被覆蓋度動態變化特征

圖3表示研究期內榆林市植被覆蓋度最大值的年際變化。榆林市植被覆蓋度整體呈波動增長的趨勢，2019年最大植被覆蓋度為52.3%，較2000年增長了11.2個百分點，年增長率為1.37%（P＜0.05）。最大植被覆蓋度在2001年為最低值，可能因為2001年退耕還林還草工程剛剛開始，耕地面積減少，而新種植的植被覆蓋較低［21］。

圖3 2000—2019年榆林市最大FVC年際變化

基于SEN+Mann-Kendall的植被覆蓋度趨勢分析結果表明（圖4），榆林市8.3%的地區植被覆蓋度呈嚴重退化趨勢，19.6%的地區植被覆蓋度呈輕微退化趨勢，9.1%的地區植被覆蓋度穩定不變，30.7%的地區植被覆蓋度呈輕微改善趨勢，32.3%的地區植被覆蓋度有明顯改善。植被覆蓋度呈改善的地區集中在黃土丘陵溝壑區的子州縣、吳堡縣、綏德縣、米脂縣、清澗縣，橫山區、相比于風沙草灘區，這些地區降水量更多，土壤條件更好，植被恢復最為明顯。呈退化趨勢的地區主要分布在定邊南部和靖邊南部、榆陽區西部和神木西部，屬于風沙區，說明這些地區沙地生態系統相對脆弱，植被生存條件不穩定［22］。

圖4 2000—2019年榆林市植被覆蓋度變化趨勢分布

3.3 植被覆蓋度變化與高程、坡度的分異特征

3.3.1 植被覆蓋度變化與高程的分異特征由于不同高程地區的氣溫、空氣濕度、土壤肥力等因素不同，因此對植被生長的影響程度也不同。從圖5a可以看出，榆林市平均植被覆蓋度隨高程增加呈先上升后下降的趨勢，平均植被覆蓋度的最大值出現在高程730～930 m的區域內，值為0.62，這里位于各川川谷，地勢平緩，土壤水分良好，適宜植被生長和植樹種草；1 130～1 330 m高程區域的平均植被覆蓋度最低，值為0.41。高程730～1 130 m范圍內的植被覆蓋度改善最為顯著，之后隨著高程增加，植被覆蓋改善效果逐漸下降，植被覆蓋退化漸為明顯。1 530～1 922 m高程范圍內的植被覆蓋退化區域面積占比最大，該區域溝壑縱橫，生態脆弱，水土流失嚴重，人工植樹造林難度大，因此植被退化較為明顯，可見大于1 530 m高程區域內的植被退化需引起進一步關注。3.3.2植被覆蓋度變化與坡度的分異特征由圖5b可知，榆林市植被覆蓋度隨坡度的增加呈緩慢上升的趨勢，平均植被覆蓋度在坡度為0°～5°的區域最低，為0.39，在坡度大于35°的區域最高，為0.59。一方面，坡度平緩的區域人類活動較為密集，坡度大的區域人類活動受限，植被生長受到的干擾減小，一定程度上導致植被覆蓋隨坡度增大而上升［23］；另一方面，退耕還林工程按照“先陡坡后緩坡”的原則實施，對于坡度大、水土流失嚴重的坡耕地優先造林，因此坡度大的地區植被覆蓋度越大，植被覆蓋度的改善面積所占比重也隨之增大。坡度為0°～5°的區域內植被覆蓋退化面積所占的比重最小，為31.8%，而改善的面積所占比重最大，為53.3%，可能是因為坡度平緩的地區植被覆蓋受人為影響干擾最大。

圖5 不同地形下植被覆蓋度及時序變化趨勢

3.4 植被覆蓋度與氣候因子相關分析

將榆林市2000—2019年植被覆蓋度與年平均氣溫和年降水量進行偏相關分析和顯著性檢驗，見圖6。可以看出，植被覆蓋度與降水的偏相關性更大，降水量更能決定植被的生長狀況。植被覆蓋度與年降水量呈顯著相關（P＜0.05）的像元占榆林市總像元數的11.5%，其中，呈顯著正相關的區域占9.0%，主要分布在中部的黃土丘陵溝壑區，該區域水土流失嚴重，植被生長對降水的依賴性較強。植被覆蓋度與年平均氣溫的偏相關關系中，有7.2%的區域通過了顯著性0.05的檢驗，呈顯著正相關（P＜0.05）的像元占榆林市總像元數的5.4%，主要分布在府谷縣、神木市、佳縣等東部地區；僅有1.8%的地區植被覆蓋度與年均氣溫呈顯著負相關，主要位于定邊縣南部、靖邊縣西南部，這些地方海拔較高，降水量較少，溫度升高導致植被表層土壤含水量減少，進而加劇植被根層干燥化，抑制植被的生長［24］。

圖6 植被覆蓋度與降水（a）和氣溫（b）的偏相關性

3.5 植被覆蓋度變化的持續性分析

利用2000—2019年榆林市的植被覆蓋度數據，計算Hurst指數，以預測植被覆蓋度的未來變化趨勢。榆林市Hurst指數均值為0.58，Hurst指數小于0.5的面積占19.1%，大于0.5的面積占80.9%，表明未來一段時期內榆林市植被覆蓋度變化趨勢與2000—2019年的變化趨勢保持一致的地區面積占比為80.9%。結合前面的FVC變化趨勢得到榆林市植被覆蓋度變化的持續性特征空間分布圖（圖7），一共得到6種變化情形：持續性嚴重退化、持續性輕微退化、持續性穩定不變、持續性輕微改善、持續性明顯改善、未來變化趨勢不確定。

圖7 基于Hurst指數的FVC變化趨勢分布

植被覆蓋度持續性輕微改善和明顯改善的區域面積分別占21.9%、31.9%，主要分布在府谷縣、佳縣、米脂縣、吳堡縣、綏德縣、清澗縣、子州縣、橫山區、神木市東南部和榆陽區東南部；持續性穩定不變的區域占9.1%，主要分布在毛烏素沙地南緣等風沙草灘區和靖邊縣最南端；持續性輕微退化和持續性嚴重退化區域面積分別占14.0%、8.0%，主要分布在定邊縣南部、靖邊縣南部以及神木市西北部；未來變化趨勢不確定的地區占15.1%，這部分地區需研究人員進一步關注和研究。同時，持續性退化的區域還需要研究人員重點關注其植被恢復受干擾的原因。

4 結論

1）榆林市植被覆蓋度由西北向東南遞增，北部風沙草灘區的植被覆蓋度低，南部黃土丘陵溝壑區的植被覆蓋度高；中植被覆蓋度所占面積最大，約占榆林市總面積的25.1%。

2）2000—2019年榆林市植被覆蓋度整體呈波動增長的趨勢，20年共增長了11.2個百分點，年增長率為1.37%；研究區內有63.0%的區域植被覆蓋得到了改善，其主要分布在土壤肥力、降水條件更好的黃土丘陵溝壑區；有27.9%的區域植被覆蓋呈退化趨勢。

3）榆林市植被覆蓋度受地形因子和氣候因子的影響，其高程、坡度分異特征明顯，植被覆蓋度隨高程增加呈先上升后下降的趨勢，隨坡度的增加呈緩慢上升的趨勢；在高程為730～1 130 m時植被改善區域面積占比最大，在坡度為0°～5°時植被退化區域面積占比最大；氣候因子中，降水對榆林市植被覆蓋的影響大于氣溫，植被覆蓋度與降水呈顯著相關的像元占比為11.5%，與氣溫呈顯著相關的像元占比為5.4%。

4）榆林市Hurst指數均值為0.58，有80.9%的植被覆蓋區域具有正向持續性，其中，持續改善的區域占53.8%，持續性穩定不變的區域占9.1%，持續退化的區域占22.0%，還有15.1%的區域未來變化趨勢不確定，需研究人員繼續關注。