黃土高原植被的時空變化及其驅動力分析研究

董 鐿 尹冬勤 李 淵 嚴泰來 王紅說,2*

(1.中國農業大學 土地科學與技術學院，北京 100083；2.農業農村部農業災害遙感重點實驗室，北京 100083)

植被作為構成地球生態系統的主要部分，是連接大氣、水體、土壤的自然紐帶，在全球物質循環和能量流動中扮演重要角色[1-2]。在全球氣候變化和城市化的背景下，監測植被的時空變化，評估氣候變化和人類活動對植被變化的影響程度，對確定合理的生態工程布局和適應性管理對策具有重要的實用價值，已經成為生態系統可持續發展的熱點問題[3]。遙感衛星在比較短的時間可以進行大規模的監測，得到的數據具有信息量大、區域廣的特點[4]。遙感等技術手段在植被監測中得到了廣泛的運用。歸一化植被指數NDVI定義為遙感影像的近紅外波段(NIR)和紅波段(R)反射率的比值[5-6]，是作為植被生長狀況的最有效指導因子之一。利用長時間序列的NDVI影像數據，可以定量的監測植被的生長變化[7]。

黃土高原水土流失嚴重，是我國的生態脆弱區和水土保持重點區域[8]。為了改善黃土高原日益惡化的生態環境，我國開始實施退耕還林、還草的政策，樹立和踐行“綠水青山就是金山銀山”的基本理念，堅持生態優先、綠色發展[9]。作為世界上最大的黃土區域，黃土高原植被覆蓋是維系我國生態平衡不可缺少的重要生態屏障。因此，研究黃土高原地區的植被變化特征，探討氣候變化和人類活動對植被的影響變得十分重要和必要。

植被的變化是一個長期的過程，需要長期的觀察才能得到其時空變化趨勢科學的結論[10]。已有研究表明植被的生長受到環境驅動因素的影響[11]。Ichii等[12]利用全球NDVI和氣候數據，分析植被NDVI與氣候之間的關系，研究發現：在北半球中、高緯度地區NDVI的增加與溫度升高有關，而南半球半干旱地區NDVI的減少則是由于降水減少所引起；降水量與植被生長之間的關系具有空間異質性，在干旱地區，由于降水量小，蒸發量大，植被NDVI與降水呈正相關，而在潮濕地區，植被NDVI與降水之間呈負相關[13]。除氣候變化外，人類活動也是植被變化的主要驅動力。隨著人口的增長和經濟的發展，人類活動對植被的影響越來越大[14]。全球的植被變化主要由人類活動造成，其中中國和印度在全球的植被變化中占據主導地位[15]。自20世紀80年代以來，中國在黃土高原上實施了一系列控制水土流失的政策，這些政策在改善植被里發揮了重要作用，1999—2010年，黃土高原累計造林面積已達1.89×107hm2[16]。

已有研究表明植被變化與氣候變化、人類活動有著密切的關系，但是關于評估和度量人類活動和氣候變化對黃土高原植被影響程度的研究卻很少。因此，本研究擬利用黃土高原1982—2015年長時間序列GIMMS NDVI 3g數據及同時期氣象資料，輔助土地利用數據和人口空間數據，并結合趨勢分析、Mann-Kendall檢驗、Hurst指數等方法，研究黃土高原的植被變化特征及未來變化趨勢，以期更加全面的揭示黃土高原地區近年植被的時空變化，利用殘差分析定量的分析氣候因子和人類活動對植被變化，為生態環境改善、植被恢復以及區域經濟發展提供科學依據。

1 研究區概況與數據獲取

1.1 研究區概況

黃土高原地區位于北緯33°41′～41°16′，東經107°54′～114°33′。位于我國中北部，總面積64.87×104km2(圖1)。其行政區劃上包括寧夏回族自治區、內蒙古自治區、青海省、甘肅省、陜西省、山西省和河南省。黃土高原是內陸地區，是溫帶大陸性季風氣候，地勢由西北向東南傾斜。黃土高原地區地質地貌復雜，氣候的多變性和人類的密集活動，使得黃土高原的水土流失和荒漠化十分嚴重。

1.2 數據來源及預處理

(1)選取GIMMS—NDVI時序數據為美國國家航空航天局戈達德航天中心全球監測與模型組發布的15 d合成的最大值NDVI數據集(https:∥ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/)，其空間分辨率為0.083°×0.083°，時間跨度為1982—2015年，該數據經過加工校正，使得數據精度變得更高、誤差變得更小，適用于植被覆蓋變化的長期監測。采用國際上通用的最大合成法(MVC)對GIMMS NDVI進行處理得到月數據，最后再由月數據平均得到年數據。

(2)氣象數據。氣象數據來源于國家氣象局氣象數據共享網(http:∥data.cma.cn/site/index.html)，選取1982—2015年黃土高原及周邊共101個氣象臺站的月降水和月平均溫度數據。氣象數據處理基于ANUSPLIN 4.3軟件，該軟件基于薄盤光滑樣條函數理論，能夠很好的實現空間插值，可結合DEM數據對氣象數據進行插值。

(3)土地利用數據。該數據來自中國科學院地理科學與資源研究所數據共享中心 (http:∥www.resdc.cn/)。空間分辨率為1 km×1 km，該數據主要基于30 m×30 m分辨率的Landsat TM/ETM遙感影像，人工目視解譯后生成，包括耕地、林地、草地、水域、居民地和未利用土地6 個一級類型以及25 個二級類型。1995—2010年全國人口空間分布公里網格數據 (http:∥www.resdc.cn/)也是中國科學院地理科學與資源研究所數據共享中心提供，是根據全國人口普查數據創建的。將土地利用數據和人口空間分布重采樣和投影變換，以匹配相對應的NDVI數據。

圖1 研究區地理位置及氣象站點分布Fig.1 Study area and meteorological station location

2 研究方法

2.1 趨勢線分析法

用一元線性回歸計算出黃土高原地區1982—2015年的變化趨勢，其計算公式如下：

(1)

式中：Slope為像元NDVI回歸方程的斜率；i代表年序號；n為年跨度；MNDVI,i為第i年最大的NDVI值。Slope>0,表示隨時間變化NDVI值呈上升趨勢，且Slope值越大，表示上升趨勢越明顯；反之，若Slope<0，表示隨時間變化NDVI值呈下降趨勢。

2.2 Theil-Sen median趨勢分析與Mann-Kendall檢驗

除了采用線性回歸方法外，本文還采用Sen方法分析植被變化趨勢的空間分布特征，并通過Mann-Kendall方法進行檢驗。該方法能夠適用于數據缺失的資料，能夠剔除異常值的干擾，對離散數據和測量誤差有較強的規避能力，其分析公式如下：

(2)

式中：β代表NDVI的變化趨勢，β>0代表NDVI呈現上升的趨勢，β<0表示NDVI呈現下降的趨勢；xj、xi為序列數據。

Mann-Kendall方法是一種非參數統計檢驗的方法，在長時間序列的趨勢檢驗和分析中得到了廣泛的應用，其統計檢驗方法如下：

(3)

(4)

(5)

式中：S為檢驗統計量；Z為標準化后的檢驗統計量；xj、xi為時間序列數據；n為序列樣本數；當n≥8時，S近似為正態分布，其方差計算公式如下：

(6)

標準化后的Z為標準正態分布，若|Z|>Z1-a/2，則表示存在顯著的趨勢變化。Z1-a/2為標準正態函數分布表在置信度水平α下所對應的值。

2.3 Hurst指數分析

Hurst指數可預測分析時間序列預測未來發展趨勢[17]，在水文、氣象等多個領域得到廣泛應用[18]。本研究采用重標極差(R/S)分析方法計算Hurst指數。其計算方法步驟如下：

1)對于時間序列[(NDVI(t))],t=1,2,…,n，定義其均值序列：

(7)

2)計算累計離差公式為：

(8)

3)計算極差序列公式為：

(9)

4)計算標準差序列公式為：

(10)

5)計算Hurst指數公式為：

(11)

Hurst指數值范圍0～1：若0

2.4 氣象因子與NDVI的相關分析

氣候變化是影響該研究區NDVI變化的主要因素之一。因此，通過逐像元計算年最大NDVI與年均氣溫和降雨之間的相關系數，來表征氣候因子與NDVI的相關程度。其計算公式為：

(12)

(13)

式中：Rxy,z為自變量z固定后因變量x與自變量y的偏相關系數。偏相關系數的顯著性檢驗，一般采用T檢驗法，其計算公式如下：

(14)

式中：tp為偏相關系數顯著性檢驗的統計量；Rxy,z為偏相關系數，n為樣本量，m為自變量的數目。

氣候因素之間通常是相互作用的，并且對植被的影響不是對立的，利用多元回歸分析方法研究溫度、降水氣象要素對NDVI的共同影響[19]，并且利用F檢驗對上述多元回歸方程進行顯著性檢驗。

復相關分析的表達式如下：

(15)

式中：x是因變量，y、z為自變量，為復相關系數。

復相關系數的顯著性檢驗采用F檢驗法，其統計量計算公式為：

(16)

式中：Fm為檢驗顯著性的統計量，Fx,yz為復相關系數，n為樣本數，k代表自變量數。

2.5 殘差分析

殘差分析是由Evans和 Geerken[20]提出，殘差值是指實際觀測值與預測值之間的差異值。本文通過剔除年NDVI變化中氣象因子因素的影響，得到植被生長變化中的人為因素。進行NDVI值與氣象因子做回歸分析，得到預測值。在不考慮其他因素的情況下，計算出NDVI的預測值和真實值之間的差值，以此推出人為因素對植被覆蓋變化的影響。其表達式如下所示：

y=a+bx1+cx2

(17)

ε=NDVIobs-NDVIpre

(18)

式中:y為因變量NDVI的值,x1x2為自變量，即降水和氣溫2個氣候因子；a和b分別被稱為截距和斜率；ε為NDVI中人類活動所貢獻的部分，即NDVI殘差值,NDVIobs代表影像實際值，NDVIpre代表預測值。當ε>0時，表示人類產生正面影響；當ε<0,表示人類活動產生負面影響；當ε=0，表示人類活動微弱，基本受到氣候因素的影響。然后對年殘差進行趨勢分析，如果殘差趨勢變化明顯，則NDVI受到人類活動的影響。

3 結果與分析

3.1 黃土高原NDVI時空變化特征

從NDVI的年際變化特征進行分析黃土高原1982—2015年植被指數年際間的變化趨勢，結果見圖2。由圖2可見黃土高原地區植被的NDVI總體表現為波動上升趨勢，且變化趨勢較為明顯的分為2 個時期。1982—1999年黃土高原地區植被變化出現較為緩慢的增長趨勢，而自從開始了退耕還林政策后，2000—2015年開始較為快速的增長，在2012的上升最為明顯。

圖2 黃土高原1982—2015年NDVI時間變化趨勢Fig.2 Time trend of NDVI in the Loess Plateau from 1982 to 2015

為了研究黃土高原地區NDVI的空間分布格局，計算了1982—2015年的平均值，并繪制了空間分布圖，由圖3中可見，黃土高原植被格局呈現由東南向西北遞減的態勢，高值區主要分布在黃土高原東南部地區，該地區的植被類型主要分布為針葉林、落葉闊葉林、灌叢，植被長勢較好。低值區主要分布在西北部區域，這些地區多為荒漠地區、牧區，植被較為稀疏。自從21世紀以來，植被覆蓋率逐漸提高，退耕還林之后NDVI均值由0.478增加到0.51。

圖3 3期黃土高原年均歸一化植被指數空間分布Fig.3 Spatial distribution of annual average NDVI in the Loess Plateau in three periods

3.2 黃土高原NDVI變化趨勢分析

為了研究NDVI的變化趨勢，采用一元線性回歸分析法，通過95%的顯著性檢驗，逐像元模擬黃土高原NDVI植被變化特征(圖4)。由圖4(b)可見，1982—2015年以來，黃土高原大部分地區植被呈現增加趨勢，植被呈現增加的區域占總面積的87.90%，有71.20%的區域通過了(P<0.05)的顯著性檢驗，植被減少大都出現在省會區域或者重工業區域，比如西寧，西安，包頭等，由于經濟的快速發展，人口的大量集中，使生態環境遭到破壞。而陜北高原、黃土高原的中部、東南部，NDVI出現顯著的增加趨勢，表明該地區的生態環境得到明顯的改善。

根據Sen趨勢分析和Mann-Kendall檢驗結合得到變化趨勢的分布圖圖4(a)，黃土高原的大部分區域在近三十年表現出較強的增加趨勢。植被呈現增加的區域占總面積的87.50%，經過Mann-Kendall檢驗(P<0.05)后，有66.32%區域通過了假設檢驗，植被減少的區域大多出現在經濟較發達區域，NDVI的減少與這些地區的快速城市化發展有著密切的聯系。

圖4 黃土高原NDVI Sen’s趨勢變化(a)與Slope趨勢變化(b)Fig.4 NDVI Sen’s (a) and Slope (b) trend changs in the Loess Plateau

3.3 黃土高原NDVI趨勢預測

基于Hurst指數和研究區的實際情況，對黃土高原NDVI的可持續進行統計分析發現大部分區域是反持續發展的。根據圖5分析黃土高原NDVI的Hurst指數平均值是0.416，呈現持續性和反持續性的區域分別占13.70%和86.30%，表明黃土高原地區反向特征要強于同向特征。從空間分布來看，NDVI呈現持續性發展的地區大多集中在中部地區，未來的植被將可能改善；而反持續的主要分布在北部丘陵溝壑區，未來植被可能會出現退化趨勢，需要注意。

圖5 黃土高原NDVI可持續性空間分布Fig.5 Spatial distribution of NDVI persistent in the Loess Plateau

3.4 植被變化與氣候因子相關分析

黃土高原地區植被NDVI和降水和平均溫度的相關分析見圖6。由圖6可知，NDVI與降水的相關系數在-0.53～0.77，與溫度的相關系數是-0.76～0.70。經過分析，NDVI與降水的相關系數的平均值為0.21，呈現正相關的區域占總面積的84.51%，大多分布在干旱區、半干旱區所在的城市，15.49%的區域呈現負相關，分布在東南部半濕潤地區，只有30.82%通過了(P<0.05)顯著性檢驗。NDVI與平均溫度的相關系數平均值0.14，呈現正相關的區域占總面積的79.94%，負相關的區域約占20.06%，分布在半干旱區的包頭市和呼和浩特市；半濕潤區東南部的西安市和渭南市；半濕潤區中東部的山西中部太原市、呂梁市、延安市等，只有28.76%通過了(P<0.05)顯著性檢驗。因此相比于黃土高原植被與溫度的關系，黃土高原的降水對植被變化的影響更加明顯。

圖6 黃土高原NDVI與年降水量、年平均氣溫間的相關系數(a)、(c)與及其顯著性(b)、(d)Fig.6 Correlation between NDVI and annual precipitation (a) and annual average temperature (c) and their significances (b) and (d) in the Loess Plateau

進一步采用偏相關分析來分析降水和氣溫對植被的影響。據圖7可知NDVI與降水的偏相關系數在-0.51～0.80，與平均溫度的偏相關系數-0.76～0.71。經過分析，NDVI與降水的偏相關系數的平均值是0.26，呈現正相關的區域占總面積的88.56%，分布在干旱區、半干旱區的絕大部分城市，11.44%的區域呈現負相關，分布在東南部，總共有37.74%通過了(P<0.05)顯著性檢驗；NDVI與平均溫度的偏相關系數的平均值是0.21，呈現正相關的區域占總面積的81.45%，負相關的區域占18.55%，大約有35.6%通過了(P<0.05)顯著性檢驗。

采用偏相關系數分析降水和植被的單獨作用后，進一步采用復相關系數分析降水和氣溫對植被的綜合影響。NDVI與降水、平均氣溫的復相關系數在0.005～0.80，平均值是0.40。有48.95%通過了顯著性檢驗，復相關系數比偏相關系數、相關系數更高表明氣候因子對植被生長的復合影響力更大。

圖7 黃土高原NDVI與年降水量、年平均氣溫間的偏相關系數(a)、(c)與及其顯著性(b)、(d)Fig.7 Partial correlation between NDVI and annual precipitation (a) and annual average temperature (c) and their significances (b) and (d) in the Loess Plateau

圖8 黃土高原NDVI與降水量、年平均氣溫的復相關系數(a)與顯著性(b)Fig.8 Multiple correlation between NDVI and annual precipitation and annual average temperature in the Loess Plateau (a) and its significance (b)

3.5 人類生產活動對植被的影響

除了氣候變化因素外，人類生產活動(人工生態恢復工程、人口密度)也會影響植被的生長變化。

從圖9可以看出，黃土高原地區殘差值呈上升趨勢，且殘差值由負轉正，說明人類活動對黃土高原的植被影響由消極轉變成積極。從2000年開始，趨勢明顯加快，說明退耕還林有著改善植被作用。由圖10可知，黃土高原大部分地區的NDVI殘差呈現增長的趨勢，增長的區域達到88.94%，殘差降低的區域為11.06%。NDVI殘差呈現減少的區域主要出現在一些經濟發展較快的區域，如西安、銀川、西寧、呼和浩特等周邊地區。這些地區經濟發展較快，人口比較集中化，以及城市化快速發展，從而導致該經濟發達地區的NDVI值減少。上述地區是我國重要工業基地，在發展的過程中，不免會破壞當地的生態環境。NDVI殘差增長的區域基本是黃河灌溉區、陜北高原等地區。NDVI殘差的增加說明體現人類活動對植被變化起到了促進作用，例如陜北高原、鄂爾多斯高原是明顯的增加，是退耕還林、退耕還草的重點區域。說明人類活動對植被變化起著很重要的作用，為生態治理提供了科學依據。

圖9 NDVI殘差變化趨勢時間分布Fig.9 Trend of NDVI residual time in the Loess Plateau

圖10 NDVI殘差空間變化趨勢Fig.10 Spatial distribution of NDVI residuals in the Loess Plateau

為了更深入研究氣候變化和人類活動對植被的影響程度，通過疊加NDVI的趨勢圖、年最大NDVI與年降水量、年均氣溫的復相關系數和殘差趨勢圖可以看出，人類活動和氣候因子引起的植被恢復約占總面積的11.32%，人類活動引起的植被恢復約占總面積的32.74%，氣候因子引起的植被恢復約占總面積的4.56%。而引起植被退化的面積占總面積非常小，說明整體氣候與人類活動都是促進植被恢復，與氣候因子相比，人類活動對植被恢復的影響更加明顯，說明研究區的退耕還林、還草政策實施較好，人為因素起了作用。

4 討 論

4.1 土地利用的變化

人類活動引起的土地利用也是植被變化的主要因子，從土地利用數據表明，1980—2010年以來，耕地、草地、水體和未利用土地面積減少，而林地和建設用地面積增加，說明一系列的生態工程活動在黃土高原地區起到了良好的作用。耕地、草地、水體和未利用土地面積分別1980年的211 791、268 249、10 456和46 605 km2減少至2010年的210 234、264 759、9 566和45 787 km2；林地和建設用地面積分別從1980年的97 901和14 172 km2增加至2010年的101 538和17 285 km2。黃土高原土地利用類型具有明顯的分階段變化特征，以2000年為分界點，1980—2000年主要變化特征為草地面積減少，耕地、建設用地增加，隨著人口的增加，城市化的發展，大量耕地轉變為建設用地，大量草地被開墾成耕地，草地面積減少；2000—2010年主要變化特征是耕地和草地面積都減少，林地面積增加。增加的林地主要來自于耕地和草地，隨著退耕還林、還草等生態工程的實施，有大量的耕地轉變為草地與林地，使得黃土高原植被面積持續增加。

表1 1980—2010年黃土高原土地利用轉移矩陣Table 1 The conversion of land use in the Loess Plateau during 1980-2010 km2

4.2 人口的規模與變化

人口的遷移也會影響著植被，根據人口普查數據，1995—2010年，黃土高原地區的總人口呈現逐年增長的趨勢，1995年黃土高原的總人口數為 9 162.77 萬人,2000年總人口10 052.43萬人,2000年后人口增速變緩,2010年總人口達到 10 771.93 萬人[21]。由于自然地理和經濟發展的限制，人口分布表現出很大的空間異質性，人口集中區域基本集中省會等經濟發達地區。根據圖12和NDVI變化趨勢圖疊加分析可知，可以表明人們傾向于生活在植被狀況良好的區域，一定的人口可以維持植被的穩定生長，但是人口稠密也會對植被產生破壞作用，大規模的城市化建設，使得以前的耕地、草地、林地等轉化為建設用地。人口的遷移還會給土地利用、植被帶來變化。對于經濟發達地區，大量人口遷入，會使建設用地不斷增加，因此在經濟發展較快的地區人類活動對植被NDVI起到了抑制的作用。

圖12 1995—2010年黃土高原人口密度(a)及1995—2010年黃土高原人口密度變化幅度圖(b)，人/km2Fig.12 Spatial distributions of population density (a) and its changes (b) in the Loess Plateau from 1995 to 2010

5 結 論

總體上，黃土高原地區生態環境變化趨勢向好，而變化原因相對于自然因素，人類活動是黃土高原地區生態環境變化的主要驅動因素，國家采取的一系列措施有效地明顯改善了該地區的生態環境，而城鎮化也帶來了一些負面影響。以上結果說明，基于氣象數據以及遙感NDVI數據采用的統計方法，對于黃土高原地區生態環境變化及其驅動因子的分析是有效的，總體結果符合實際情況。

由于黃土高原特殊復雜地形地貌的影響，數據、特別是遙感數據帶有較大的不確定性，部分地區可能誤差較大。因此，以上研究結果只是在總體上反映黃土高原的生態環境的狀況。研究的深入還要有賴于更準確的數據以及詳盡的實地調查。