生態(tài)工程背景下黃土高原植被變化時空特征及其驅動力

修麗娜， 顏長珍, 錢大文， 幸贊品

(1.西北師范大學 地理與環(huán)境科學學院， 甘肅 蘭州 730070; 2.中國科學院 西北生態(tài)環(huán)境資源研究院，甘肅 蘭州 730000； 3.中國科學院 西北高原生物研究所， 青海 西寧 810008)

黃土高原作為我國歷史上開發(fā)時間較早、人工改造力度較大的區(qū)域之一，深受人類活動干擾的影響，尤其是過去幾十年來，毀林開荒、過度放牧等行為嚴重的破壞當地的生態(tài)環(huán)境[1-2]，再加上黃土高原地處半干旱半濕潤氣候帶，風力和水力侵蝕嚴重，導致了水土流失、沙漠化加劇等一系列的生態(tài)環(huán)境問題[3-4]。針對這些生態(tài)環(huán)境問題，1990年以來，我國開始實施退耕還林還草的一系列生態(tài)恢復工程，而黃土高原則是重點實施區(qū)域之一[5-6]。植被是陸地生態(tài)系統的重要組成部分，在陸地表面能量交換、水分循環(huán)和生物地球化學循環(huán)過程中起著至關重要的作用[7]。植被的時空變化會影響全球碳循環(huán)和氣候系統的穩(wěn)定性，監(jiān)測植被的時空變化趨勢對研究全球氣候變化背景下的陸地生態(tài)系統具有重要的意義[8]。一方面，植被變化是區(qū)域生態(tài)變化的有效指標，其時空變化特征可以作為生態(tài)恢復工程的評價依據，另一方面，植被改善能夠通過改變地表反射率和下墊面粗糙度以達到減緩徑流沖蝕、風力侵蝕和保持水土的作用，對區(qū)域生態(tài)環(huán)境起促進作用[9]。因此黃土高原地表植被的變化研究對于該地區(qū)環(huán)境保護和生態(tài)工程建設具有重要的指示意義。

歸一化植被指數(normal difference vegetation index, NDVI)是目前最常用表示植被生長狀況的指標，它與植被覆蓋度、生物量和光合作用強度等都密切相關[10-11]。很多學者利用NDVI對于黃土高原地區(qū)植被狀況的研究已經做了大量的工作，并取得了豐富的成果。高健健等[12]利用1981—2012年的NDVI數據遙感反演黃土高原的植被覆蓋度，并發(fā)現2001年以后該區(qū)生長季平均植被覆蓋度大于之前，呈現增加趨勢。李雙雙等[11]發(fā)現2000—2009年陜甘寧黃土高原地區(qū)NDVI增速遠快于三北防護林工程區(qū)，這一結果主要是有人類活動和氣候變化共同驅動，且人類活動作用明顯。肖強等[7]通過NDVI數據估算了黃土高原植被覆蓋度，其研究結果發(fā)現2000—2010年植被覆蓋度總體呈現東南高西北低，由東南向西北遞減的特征。然而，這些研究通常使用單一數據集來評估相對較短時期內或較粗分辨率的植被動態(tài)，導致缺乏對植被變化的長期監(jiān)測，使這些研究的代表性有限。此外，NDVI作為表征植被生長狀況的指標，以往研究僅從NDVI單方面對黃土高原地區(qū)植被進行評價，無法獲取全面的植被變化信息。植被生長受到環(huán)境驅動因素(例如氣候變化)和人類活動的影響，在已有研究中區(qū)分這些因素的方法主要分為3種： ①利用氣候與NDVI之間的相關性，并輔助相應的統計數據間接的區(qū)分氣候因素和人類活動對NDVI變化的影響[4,11]; ②使用水分利用效率變化模型提取人類影響，該方法只適合于植被變化主要受降水因素影響的干旱和半干旱區(qū)域[13]； ③利用回歸殘差法能較好的區(qū)分氣候因素和人類活動對植被的影響，目前應用較為廣泛[1,14]。但此方法只是單純的得出人類活動對植被NDVI的正或負面影響，對于人類和氣候共同影響的區(qū)域沒有明確得出。為此，本研究利用GIMMS NDVI和MODIS NDVI相關關系將MODIS NDVI時序延伸，獲得長時間序列較高分辨率的NDVI數據集，從植被質量角度，以生長季NDVI作為指標對1990—2015年黃土高原植被的時空演變規(guī)律進行分析，輔助土地利用數據從植被面積變化角度全面分析黃土高原植被時空變化特征，并從氣候變化和人類活動兩個方面分析本地區(qū)植被變化的驅動力，以期為生態(tài)環(huán)境保護相關政策的制定、區(qū)域可持續(xù)發(fā)展提供理論依據。

1 材料與方法

1.1 研究區(qū)概況

黃土高原位于我國中部偏北，黃河流域中部，地理位置為33°41′—41°07′N, 100°54′—114°07′E。黃土高原地跨山西、寧夏、陜西、甘肅、內蒙古、河南和青海省，總面積為6.24×105km2。區(qū)域內地勢西北高、東南低，千溝萬壑、地形破碎。氣候為大陸性季風氣候，冬季寒冷干燥多風沙，夏季炎熱多暴雨。植被由東南向西北可劃分為森林帶、森林草原帶、典型草原帶、荒漠草原和草原化荒漠帶。土壤類型主要有黃綿土、黑鈣土等[12,15-16]。

1.2 研究數據

1.2.1 NDVI數據來源及處理 本研究利用的遙感數據為NOAA/AVHRR和TERRA/AQUA MODIS的NDVI數據集，NOAA/AVHRR NDVI數據集來自戈達德航天中心(Goddard Space Flight Center， GSFC)全球監(jiān)測與模型研究(Global Inventor Modeling and Mapping Studies， GIMMS)數據中心。影像以32bit geotiff格式存放，采用Albers Conic投影方式，空間分辨率為8 km×8 km，時間分辨率為15 d，時間序列從1990年1月至2006年12月(可在ftp.glcf.umiacs.umd.edu免費下載)。MOD13A3是MODIS陸地產品組開發(fā)的月植被指數產品，包括月最大合成的NDVI和EVI植被指數，空間分辨率為1km，可以通過NASA對地觀測系統數據共享平臺(https:∥earthdata.nasa.gov/)下載。

已有研究[17]發(fā)現GIMMS和GIMMS 3g數據集在北半球存在較大差異，因此建議可以結合GIMMS和MODIS數據對該區(qū)域的植被變化進行分析。通過GIMMS NDVI和MODIS NDVI相關關系建立將MODIS NDVI時序延伸。首先將15 d GIMMS NDVI數據采用最大值合成法生成月NDVI數據集，然后將MODIS 1km月NDVI(MOD13A3)數據集按照GIMMS NDVI的8km格網分別計算每個格網的月NDVI均值和最小值。計算最小值主要是為排除非植被的區(qū)域做數據準備，求得的均值是與GIMMS NDVI數據做匹配。選取研究區(qū)內，2000—2006年7年間逐月最小值均>0的格網，分別計算7年每個月GIMMS NDVI同MODIS NDVI均值之間的線性相關系數[18-20]，結果詳見表1。選取最小值大于0的柵格主要是篩選出植被分布區(qū)域，由于1 km分辨率相較8 km的空間分辨率較精細，便于排除非植被區(qū)域NDVI值參與計算。然后依據此關系，使用1990—1999年GIMMS NDVI擴展MODIS NDVI的時間序列。為了反映植被NDVI的年際變化特征，本研究中采用生長季NDVI值(4—10月NDVI的平均值)來表征植被生長[8]。

表1 GIMMS NDVI與MODIS NDVI之間的相關關系

1.2.2 氣象數據來源及處理 氣象數據來源于中國氣象局提供的中國氣象科學數據共享服務網(http:∥cdc.cma.gov.cn/home.do)，包括1990—2015年黃土高原及其周邊地區(qū)共84個氣象臺站觀測的日平均氣溫和日降水量數據。本文利用ANUSPLIN 4.3軟件對整理好的月降水和溫度進行空間插值，空間分辨率為1 km。該方法已經廣泛應用于氣象因素的插值[21]，其精度要高于反距離權重(inverse distance weighted, IDW)和克里金插值(Kriging)[22]。

1.2.3 土地利用數據的解譯 本研究中用于土地利用解譯的影像包括沒有云覆蓋的1990，2000，2005，2010的TM影像和2015年的OLI影像，這些數據下載于美國地質調查局(United States Geological Survey，USGS)網站(http：∥glovis.usgs.gov/)。大多數選定的圖像是從6月到9月獲得的，丟失的圖像或質量差的圖像被替換為同期最近的月份或年份的圖像。在進行遙感解譯之前，影像已經通過幾何校正，輻射校正和大氣校正等預處理。本研究中使用的6種土地利用類型分別是：林地，草地，耕地，人工表面和其他類型，其中其他類型包括裸巖、戈壁、裸土、沙漠和鹽堿地[23]。通過eCognition 8.6軟件從影像中提取土地覆蓋信息，該軟件使用面向對象的分類方法[24-25]。

1.3 研究方法

1.3.1 變化趨勢分析 為了研究過去26 a來黃土高原植被的生長狀況及在空間上的差異，使用一元線性回歸對每個像元進行斜率分析[7]。當斜率為正時，表示隨時間變化植被指數升高；反之，斜率為負時，表示隨時間變化植被指數呈現下降趨勢。對所得的斜率進行顯著性水平檢驗，通過置信度檢驗的區(qū)域則認為植被變化的趨勢顯著。其計算公式為：

(1)

式中：n——模擬的時間序列長度(n=26)；i——時間序列的序號； NDVIi——第i年的NDVI值；Θslope——變化曲線的斜率。

1.3.2 土地利用類型凈變化率 年均土地利用類型凈變化率定量反映該土地利用類型面積變化速率，其對比較各土地利用類型變化的區(qū)域差異具有重要作用[26]。其計算公式為：

(2)

式中：K——某土地利用類型年均凈變化率；Ua，Ub——研究初期與末期某一種土地利用類型的面積；T——研究期時長，時段設定為年。

1.3.3 土地利用類型轉移矩陣 利用土地利用類型面積作為土地利用類型狀態(tài)轉移矩陣中的向量，分別統計1990—2000年和2000—2015年2期的轉移矩陣[27]，計算過程為：

(3)

式中：A——面積；i，j(i，j=1,2，…，n)——轉移前與轉移后的土地利用類型；Aij——土地利用類型從類型i變?yōu)閖的面積；n——轉移前后土地利用類型數。

1.3.4 多元回歸分析 氣候因素(溫度和降水)之間通常是相互作用的，并且對植被的影響不是對立的。本研究利用多元回歸分析方法來研究溫度、降水對NDVI的共同影響[28]。

(4)

利用F檢驗對上述多元線性回歸方程進行顯著性檢驗。

1.3.5 改進的殘差趨勢分析 殘差趨勢法由Evans和Geerken提出[13]，該方法利用每個像元全時間序列的年NDVI與降水數據建立響應關系，從而得到每個像元在全時間序列上實際NDVI與預測NDVI的殘差，從而判斷人類因素對植被生長狀況的影響[29]。曹鑫等[30]考慮到溫度對植被生長的影響，在已有模型基礎上加入了溫度因子。

ΔNDVI=NDVI真實值-NDVI模擬值=

(5)

為了研究更深入地研究氣候和人類活動對于植被變化的影響程度，田海靜[31]等在以上模型的基礎上，將年度NDVI的變化趨勢與氣候變化之間相關性不顯著的區(qū)域相對應。 ①由氣候變化引起的植被顯著恢復。 NDVI隨時間的變化趨勢顯著增加，NDVI與氣候因素顯著相關且NDVI殘差趨勢變化不顯著; ②由氣候變化引起的植被顯著退化。NDVI隨時間的變化趨勢顯著下降，NDVI與氣候因素顯著相關且NDVI殘差趨勢變化不顯著； ③由人類引起的植被顯著恢復。NDVI隨時間的變化趨勢顯著增加，NDVI與氣候因素不顯著相關且NDVI殘差趨勢變化顯著增加; ④由人類引起的植被顯著退化。NDVI隨時間的變化趨勢顯著下降，NDVI與氣候因素不顯著相關且NDVI殘差趨勢變化顯著下降; ⑤氣候因素和人類活動共同引起的植被顯著恢復。NDVI隨時間的變化趨勢顯著增加，NDVI與氣候因素顯著相關且NDVI殘差趨勢變化顯著增加； ⑥氣候因素和人類活動共同引起的植被顯著退化。NDVI隨時間的變化趨勢顯著下降，NDVI與氣候因素顯著相關且NDVI殘差趨勢變化顯著下降。

2 結果與分析

2.1 黃土高原地區(qū)植被質量時空變化

2.1.1 黃土高原地區(qū)植被NDVI年際變化特征 由圖1可以看出，1990—2015黃土高原地區(qū)植被的NDVI總體表現為上升趨勢，且變化趨勢較明顯的分為兩個時期：1990—1999年期間，黃土高原地區(qū)的NDVI值波動較大，在0.3304～0.367 0之間，有輕微的減少趨勢，特別是1998和1999年；2000年以后，該區(qū)域的NDVI表現為增加趨勢，范圍為0.286 4～0.372 3，最大值出現在2012年。

圖1 黃土高原1990-2015年生長季NDVI年際變化

2.1.2 黃土高原地區(qū)植被NDVI空間變化特征 為了便于分析1990—2015年黃土高原NDVI變化趨勢的空間分布規(guī)律，對這26 a以來生長季的NDVI逐像元進行一元線性回歸分析(圖2)，將所得相關系數進行0.05顯著性水平檢驗，所得變化趨勢分布圖(圖3)。黃土高原地區(qū)植被NDVI呈現增加趨勢的區(qū)域面積較大，主要分布在黃土高原的東部、中部和南部地區(qū)，包括陜西省的北部，甘肅和寧夏地區(qū)南部，內蒙古自治區(qū)的東部，以及山西省的西部，其他區(qū)域的植被也有不同程度的增加；NDVI降低的區(qū)域主要分布在寧夏和內蒙古的北部，陜西的南部，以及山西的部分區(qū)域。年均增長率最大的像元值達到0.020 0，年均負增長率最大的像元值達到-0.027 4(圖2)。

圖2 黃土高原1990-2015年生長季NDVI變化趨勢

從黃土高原NDVI顯著性變化趨勢分布圖(圖3)可以看出，1990—2015年黃土高原植被NDVI整體呈現的上升趨勢，顯著上升的區(qū)域主要集中在中部和東部地區(qū)。研究結果表明(表2)，1990—2015年期間，黃土高原地區(qū)53.45%的地區(qū)，植被NDVI呈上升趨勢，其中33.05%的地區(qū)植被為顯著上升。植被NDVI下降的區(qū)域主要分布在甘肅和寧夏地區(qū)的北部，陜西省的南部和山西省的東部等地，有46.55%的地區(qū)植被NDVI呈下降趨勢，其中24.61%的地區(qū)植被NDVI呈顯著下降。1990—1999年，黃土高原的植被NDVI總體呈現降低趨勢，上升區(qū)域面積占總面積的44.92%，其中顯著上升僅占6.45%；2000年以后，相較之前，黃土高原植被NDVI迅速上升，上升區(qū)域面積達到總面積的91.90%，顯著上升為65.78%。

圖3 黃土高原1990-2015年生長季NDVI顯著性變化分布

表2 各時間段黃土高原生長季NDVI變化類型像元個數及百分比

2.2 黃土高原地區(qū)植被面積變化特征

黃土高原地區(qū)植被面積總體表現為增加趨勢(表3)，林地主要分布在黃土高原的南部和東部地區(qū)，陜西和山西省境內，草地分布于黃土高原的大部分區(qū)域(附圖2)。

林地自1990年開始一直表現為持續(xù)增加，且2000年之后增加速度較快，2000—2005年林地的凈變化率達到0.31%，之后減慢，但是仍快于2000年以前，直到2010，林地面積開始有輕微減少。草地面積2000年以前總體表現為減少，凈變化率為-0.10%，之后表現為與林地相同的變化趨勢，2000—2005年面積的增加速度較快。

表3 黃土高原植被面積變化

1990—2000年增加的林地主要來自于耕地的轉入，面積為470.73 km2,其次為草地和其他，面積分別為68.46和68.86 km2；減少的草地主要轉化為耕地，面積達到4 104.05 km2, 其次為人工表面，面積為270.87 km2(表4)。隨著人口數量的增多，人工表面的迅速擴張，黃土高原地區(qū)的開荒行為持續(xù)增加，導致大量的草地被開墾為耕地，或退化為裸土等其他類型。

2000—2015年，增加的林地和草地主要來自耕地的轉入，面積分別為3 723.29和14 278.59 km2；此外還有大量的裸土等其他類型轉入草地，面積達到1 197.67 km2(表4)。盡管城鎮(zhèn)化進程持續(xù)加劇的行為會破壞植被，但是也導致農村農民對土地、山林的依賴性降低，大量農田撂荒，砍伐減小。加之隨著退耕還林(草)生態(tài)工程的實施，大面積的退耕還林(草)、植樹造林等一系列措施的實行都會導致黃土高原的植被面積持續(xù)增加。

表4 黃土高原不同時期土地利用類型轉移矩陣 km2

2.3 植被對氣候變化和人類活動的響應

通過改進的殘差趨勢分析得到，人類和氣候共同引起植被恢復的區(qū)域主要分布在黃土高原的中部，主要集中在陜西的北部，占總面積的9.41%；人類引起植被恢復的區(qū)域主要分布在內蒙古的東部和北部、甘肅和寧夏的南部，以及陜西和山西省的中部地區(qū)，占總面積的21.74%；氣候引起植被恢復的區(qū)域僅占總面積的1.91%，主要分布于青海省。人類和氣候共同引起的植被退化主要分布在內蒙古、山西和寧夏的北部，占總面積的4.71%；人類引起的植被退化面積占11.99%，主要分布于內蒙古的西北部、山西和陜西省的部分地區(qū)；氣候引起的植被退化主要分布于青海和甘肅省(圖4，表5)。

表5 黃土高原1990-2015年生長季NDVI變化驅動因素解析

圖4 黃土高原1990-2015年植被 顯著恢復和顯著退化區(qū)域分布

對黃土高原植被主要受氣候影響的區(qū)域進行進一步分析。氣候引起植被變化區(qū)域的溫度和降水有著相同的趨勢，都呈上升趨勢(圖5)。氣溫在過去26 a間升高幅度較大，且植被退化區(qū)域比恢復區(qū)域溫度上升較快。降水變化趨勢表現為輕微上升，且植被恢復區(qū)域較退化區(qū)域明顯。

圖5 1990-2015年氣候引起植被恢復和退化區(qū)域的年均氣溫和年降水變化趨勢

3 討論與結論

1990—2015年，黃土高原地區(qū)植被的NDVI總體表現為上升趨勢，這一趨勢與已有的研究結果一致[8,32]。具體來說，1990—1998年黃土高原的植被處于相對穩(wěn)定的時期，之后植被NDVI迅速下降，2002年植被NDVI進入迅速上升時期，這主要是由于隨著氣候干旱化趨勢發(fā)展，植被NDVI表現為減少且幅度小，1999—2001年降水明顯偏少(圖6)，導致植被NDVI減少，自2002年以來，隨著降水量的恢復和退耕還林還草政策的大規(guī)模實施，植被NDVI呈現明顯上升趨勢[4]。自2000年起，黃土高原植被面積表現為增加趨勢，且主要來自于耕地的轉入，這說明退耕還林(草)工程的實施，在黃土高原地區(qū)效果明顯，不僅植被面積持續(xù)增加，而且生長狀況也得到了改善。但值得注意的是，城鎮(zhèn)化進程的加劇，雖然會破壞植被，但是也導致農村農民對土地、山林的依賴性降低，大量農田撂荒，砍伐減小，進而對黃土高原植被的恢復起到積極作用。

圖6 黃土高原地區(qū)1990-2015年生長季 年均溫度和年降水變化圖

26 a間黃土高原植被NDVI整體呈現的上升趨勢，顯著上升的區(qū)域主要集中在中部和東部地區(qū)。結合人類殘差分析的結果得知，黃土高原植被恢復的原因主要是氣候和人類共同影響，且人類影響的程度較大，區(qū)域占總面積的21.74%，本研究中大量耕地轉入草地印證這一結果。結合李雙雙等[11]的研究可知，陜北地區(qū)2000—2009年退耕還林(草)面積為9.81×105hm2，榆林、延安退耕還林還草面積分別為3.09×105hm2，2.91×105hm2，也說明中東部地區(qū)的植被恢復主要得益于生態(tài)工程的實施。人類引起植被退化的區(qū)域主要是由于城鎮(zhèn)化擴張。羅娟等[33]對毛烏素沙地近10 a的土地利用動態(tài)變化進行了分析，發(fā)現居民地和工礦用地面積增加，且主要是由草地轉入；馬瞳宇等[34]選取黃土高原風蝕水蝕交錯區(qū)典型代表流域進行分析，1990—2010年期間，該地區(qū)煤礦開采活動活躍，影響流域內林草生長，并對侵蝕環(huán)境的人為惡化起著巨大的促進作用；卞鴻雁等[15]發(fā)現1980—2005年黃土高原南部地區(qū)建設用地凈增了1 238.29 km2。除了以上原因，對比黃土高原土地利用分布圖可以看出，部分人類引起植被退化區(qū)域是耕地分布的區(qū)域，由于研究中所使用的遙感數據分辨率有限，不可避免的會將一些區(qū)域耕地的NDVI值的變化帶入模型中，隨著退耕政策的實施，導致耕地的減少，近而影響到最后的模擬結果。氣候因素主導植被退化區(qū)域主要是由于氣溫升高較快，但降水增加不明顯，加速了相應區(qū)域氣候的暖干化，抑制了植被的生長。以上研究結果也側面的驗證了人類殘差分析在黃土高原的有效性。

盡管本研究中使用的NDVI相較其他長時間序列的NDVI產品具有更高的分辨率，但1 km的空間分辨率對于黃土高原地區(qū)而言仍然相對較低。在一些復雜的地形中，坡度和其他地形因素對植被的反射率有一定的影響，進而影響NDVI的觀測值。此外，黃土高原所處的地理位置，使得當某一地區(qū)的植被稀疏時，植被變化趨于退化，這些退化或暴露的斑點可能會影響遙感監(jiān)測的反射率，進而影響NDVI。

殘差趨勢分析法被廣泛用于分離氣候和人為因素對植被變化的影響[35]。殘差法以像素為分析單元，提出了一種充分考慮坡度，土壤和植被空間變異性對溫度—降水—植被NDVI的影響的分析方法。然而，這種方法也有一些缺點[36]，潛在的降水—溫度—NDVI關系的趨勢分析是基于相同的時間序列。由于遙感數據的限制，殘差趨勢法只能檢測到遙感數據時間序列內的退化。考慮到影響植被NDVI變化的因素還有很多，比如蒸散量、土壤水分等，這些因素都需要在未來的研究中仔細考慮，進而改進模型。