1999-2018年黃河流域NDVI時空變化及驅動力分析

付含培, 王讓虎, 王曉軍

(1.山西大學 環境與資源學院, 太原 030006; 2.山西大學 黃土高原研究所, 太原 030006)

植被是聯結大氣圈、土壤圈、水圈和生物圈的重要紐帶，對于監測生態環境的變化具有重要指示作用[1]，因此植被覆蓋的時空變化特征及驅動因子分析在區域生態環境變化研究中具有重要意義[2-3]。歸一化植被指數(NDVI)[4]作為表征植被生長狀況的重要指示因子，與葉面積指數、生物量存在明顯的相關性，可以對地表植被生長狀況進行快速評估[5-8]。基于MODIS，SPOT/VEGETATION等衛星遙感影像得到的NDVI數據具有易獲取、分辨率高的特點，是長時序植被覆蓋研究的重要數據來源。目前國內外眾多學者已利用不同時段的NDVI數據采用不同方法對區域尺度植被覆蓋的時空變化規律及對氣候變化的響應進行了深入研究[9-10]。

黃河流域是我國第二大流域，生態環境較為脆弱，近年來在氣候變化和人類活動等因素作用下，流域植被覆蓋狀況發生了轉變[11]。目前，黃河流域植被覆蓋變化監測及驅動機制研究取得了一定的進展。張亞玲等[12]采用1998—2012年SPOT-NDVI數據發現黃河流域NDVI平均值與緯度變化存在一定的負相關性，并隨著時間增長呈現出增加的趨勢。袁麗華等[13]利用2000—2010年的MOD13Q1數據發現，黃河流域NDVI呈現出北部較低，而東南部、西部較高的特征。趙安周等[14]基于MODIS-NDVI數據發現在2000—2014年期間，黃河流域NDVI保持了增長的趨勢，并且達到了6.93%/10 a的增速。劉啟興等[15]利用2000—2016年的MODIS-NDVI數據，綜合采用多種方法(肯德爾法檢驗法等)探討了氣候因子對于黃河流域植被覆蓋特征的變化，研究結果顯示該區域NDVI在2000—2016年保持了較高的穩定性；其中大部分區域的植被NDVI保持增長趨勢，其比例達到了70.4%，而年增長率基本保持在0～0.004范圍內。郭帥等[16]分析1982—2015年GIMMS NDVI數據及相關氣象數據，得出黃河流域NDVI與年平均氣溫、年降水量以正相關為主，其中顯著正相關分別為22.39%，21.99%。葉培龍等[17]通過分析1980—2018年黃河上游各氣象站點的數據，發現黃河上游呈現暖濕化趨勢，其有利于黃河上游植被生長。但鮮有學者對影響黃河流域植被覆蓋的自然和人文因子展開系統性分析和研究，且缺少相關驅動力分析。地理探測器是探測空間分異性，以及揭示其背后驅動力的一組統計學方法[18]，已有學者將該統計方法應用于植被NDVI驅動機制的研究[19-21]。因此，本研究基于SPOT NDVI遙感數據，采用Sen趨勢分析、Manna-Kendall檢驗、Hurst指數、空間轉移矩陣分析1999—2018年黃河流域NDVI時空變化特征，并采用地理探測器模型對黃河流域NDVI空間分異特征與驅動因子進行探測分析，為黃河流域植被恢復與生態建設提供科學依據。

1 研究區概況

黃河流域位于95°50′29″—119°6′26″E，32°6′53″—41°48′18″N，流域面積為79.46萬km2，地跨9省。海拔0～6 255 m，由西向東依次為西部青藏高原區，多為山地，海拔較高；中部地區以黃土高原區為主，黃土地貌，地表破碎，水土流失嚴重；東部黃淮海平原區，主要由黃河沖擊而成，海拔較低，整體地勢西高東低。黃河流域西部為大陸性干旱氣候，中部為季風性半干旱氣候，東部屬海洋性半濕潤氣候。黃河流域年均溫9℃，流域大部分地區年降水量在200～650 mm，西南部和東南部多于650 mm，深居內陸的西北部部分地區則少于150 mm。除氣候及地形等方面的因素外，該區域存在多種類型的地貌，同樣促進了不同植被的生長。

2 數據與方法

2.1 數據來源與預處理

在本次研究中采用的植被NDVI數據來自于1 km植被指數(NDVI)空間分布數據集，并通過最大值合成法形成；氣象氣候數據來自于氣象背景數據集，基于全國1 915個站點的氣象數據插值生成500 m×500 m的空間分布數據集；地貌、土壤數據分別來自于中國1∶100萬地貌類型數據、土壤類型空間分布數據；90 m DEM來自于SRTM V4.1(Shuttle Radar Topography Mission，SRTM)數據，通過拼接重采樣生成1 km柵格數據，并利用ArcGIS進一步獲取坡度和坡向數據；土地利用類型數據主要來自于中國土地利用現狀數據集，分辨率1 km；人口密度、GDP密度數據分別來源于中國人口空間分布公里網格數據集、GDP空間分布公里網格數據集，均為1 km網格。所有的數據均通過中國科學院資源環境科學數據中心(http:∥www.resdc.cn/)獲取，并通過ArcGIS軟件利用黃河流域范圍矢量數據進行裁剪、投影變換等處理。

2.2 研究方法

2.2.1 趨勢分析 采用Sen氏趨勢分析與M-K檢驗結合，分析像元尺度上NDVI變化趨勢。Sen氏趨勢分析對數據誤差有較強的抵抗能力，且分析結果較為科學可信[22]。其計算公式如下：

當斜率S大于0時，表示呈現增長趨勢，小于0則呈現下降趨勢。

Manna-Kendall檢驗(M-K檢驗)可判斷趨勢的顯著性，且能夠排除少數異常值的干擾。因此，本文將二者相結合，用來判斷黃河流域內NDVI在像元尺度的變化趨勢及其顯著性。

2.2.2 Hurst指數 Hurst指數多用于描述自然界中長時間序列的自然現象[23]。具體公式見相關文獻[13]。若通過上述公式證明時間序列中存在Hurst現象，可利用最小二乘法擬合得到H值。

NDVI序列的完全隨機性和持續性由H值決定。當滿足條件0

2.2.3 地理探測器 目前在空間分異性研究中廣泛采用了地理探測器模型[18]。在本研究中即利用該工具分析黃河流域NDVI變化的驅動因子單個及其交互關系。已有研究表明，氣候因子對NDVI有顯著影響[15]。黃河流域面積廣闊，地形復雜多樣。西部高原山地，植被覆蓋較為原始，受人類影響較小，西北部內蒙古高原及北部河套平原，歷史上為農業灌溉區，受人類活動影響較大。因此，在考慮系統性、典型性等因素下，本研究選取氣候、地表及人類活動三大類因素，14個因子，探測其對黃河流域NDVI的影響。在ArcGIS 10.2中，按照1 km×1 km網格，隨機生成7 520個采樣點，結合空間位置關系對采樣點NDVI數據進行關聯。在本次研究中共選取14個影響因子(表1)。本文將平均降水量、濕潤指數、干燥度、積溫(≥10℃)、積溫(≥0℃)、年均溫分為6個等級，將高程分為7個等級，將坡度、地貌、土壤、GDP、人口、土地利用類型分為9個等級，將坡向分為10級，劃分過程中采用自然斷點法(圖1)。

表1 黃河流域NDVI影響因子

地理探測器具體包括如下4種類型：

(1) 分異及因子探測。其可以對NDVI的空間分異性進行探測，在此基礎上可以進一步分析不同因子產生的影響。具體公式如下：

(2) 交互作用探測。基于交互作用探測可以對人類活動、氣象以及地表等因素之間的交互作用進行分析，可以評價多因子的交互作用是否能夠影響到對于NDVI的解釋能力。

(3) 風險區探測。據此可以評價不同子區域屬性均值的差異性，用于鑒別各區域植被覆蓋情況，用t統計量來檢驗[20]。具體公式如下：

(4) 生態探測。生態探測需要借助F統計量，衡量兩因子X1和X2對NDVI空間分布的影響是否有顯著差異。

式中：L1，L2分別為X1，X2分層的數量；Nx1，Nx2分別為二者的樣本數目；SSWx1，SSWx2對應著二者的分層、層內方差和[18]。

3 結果與分析

3.1 植被NDVI空間分布格局

基于1999—2018年的NDVI數據，結合黃河流域實際情況，本文利用等間距法將植被狀況劃分為5種類型(表2)。由植被覆蓋空間格局圖(圖2)可見，黃河流域NDVI分布具有明顯的空間分異特征，總體表現為南高北低，與緯度表現出一定的負相關性，即緯度越低，則NDVI值越高。其中，較高及以上植被覆蓋度占黃河流域面積的47.13%，主要分布于流域的東部和南部地區，包括山西西部呂梁山脈、陜南及關中地區；中植被覆蓋度占26.80%，多位于中部地區，呈東北—西南走向帶狀分布；較低植被覆蓋度占21.41%，集中在流域的北部；低植被覆蓋度面積僅占4.65%，分布在河套平原兩側。黃河流域西南部主要為山地，植被以森林、草地以及草本濕地為主，有較高的NDVI值；東南部屬于海洋性半濕潤氣候，水熱較充足，且林地和農作物分布較為廣泛，所以NDVI值較高。北、中部為內蒙古高原、河套平原，生態環境脆弱，氣溫高、降水少，植被破壞嚴重，現有植被多為低矮的農作物，因此NDVI值較低。

圖1 黃河流域因子空間分布

表2 植被覆蓋類型

圖2 1999-2018植被覆蓋空間格局

3.2 植被NDVI時間變化趨勢

3.2.1 區域尺度NDVI的時間變化趨勢 由區域尺度植被NDVI的時間變化趨勢圖(圖3)可見，1999—2018年黃河流域整體呈增加態勢，平均年際變化率為0.007。20 a間，年均NDVI最高值出現在2012年，為0.44，最低值出現在2000年，為0.34。就季節來看，夏季NDVI值最高，秋季、春季其次，冬季最低。黃河流域春、夏、秋季NDVI變化趨勢與年際變化趨勢基本一致，呈波動上升趨勢，平均年際變化率分別為0.005，0.007，0.007，而冬季NDVI的年際變化相對平穩，年際變化率只有0.000 07，這可能因為當地植被主要以落葉植被為主，冬季植被大都凋落，增長不夠明顯。從月份NDVI變化趨勢來看，最低值出現在冬季的2月份，為0.16；NDVI最高值出現在夏季的8月份，為0.53；2—8月NVDI值上升明顯，8月—次年2月下降明顯，呈波峰和波谷狀。

圖3 1999-2018年區域尺度NDVI時間變化趨勢

3.2.2 像元尺度NDVI的時間變化趨勢 根據Sen趨勢分析結果將NDVI變化趨勢分為7級(圖4)，根據M-K檢驗結果(0.05置信水平)將NDVI變化顯著性分為6級(圖5)。由變化趨勢結果可見，黃河流域20 a間NDVI時間變化趨勢以極顯著緩慢增長為主，整體變化趨勢向好。其中，極顯著緩慢增長區域所占面積最大，占整個流域面積的53.12%；不顯著變化的區域其次，占29.98%；極顯著較快增長的區域占15.47%；快速增長的占1%；而極顯著較快退化和極顯著緩慢退化合計共占黃河流域面積的1.43%。極顯著較快增長和極顯著快速增長的區域主要位于黃河流域的東部，黃河兩岸，地理位置上屬山西西部呂梁山脈和陜西北部黃土高原，這與該區域持續退耕還林，人工在造林有密不可分的關系，一定程度上減少了黃土高原的水土流失，改善了黃河生態環境。基本不變的區域則位于黃河流域北部，河套平原周邊地區以及黃河下游地區，該區域城鎮、耕地面積較大，因而基本沒有變化。

圖4 1999-2018年NDVI趨勢變化

圖5 1999-2018年NDVI變化顯著性

3.2.3 NDVI變化的可持續性分析 從黃河流域NDVI的Hurst指數(圖6)可以看出，H值低于0.3的區域零星分布在流域的中部及南部地區，僅占流域面積的1.46%；H值大于0.7的區域位于中部和北部地區，占流域面積的0.67%；0.5～0.6的區域占28.78%，0.6～0.7占6.32%，主要分布于山西西部和陜西中部；0.3～0.4的區域占流域的18.13%，0.4～0.5的區域占黃河流域的45.15%，其反退化占比較高，主要分布于甘肅隴中黃土高原、寧夏地區和陜西、山西黃土高原區，這些區域本身為黃土高原或干旱沙地，植被覆蓋差，多為季節性降水，加之人類活動頻繁，其植被生長持續性不強。0.5～0.6的區域占流域面積的28.78%，0.6～0.7占6.32%。

本文將Hurst指數圖與趨勢分析圖在ArcGIS中疊加，得到黃河流域NDVI未來趨勢變化圖(圖7)。結果表明，持續快速退化和持續較快退化面積極小，僅占0.01%，幾乎可以忽略不計；持續較緩慢退化的面積占0.14%，大多分布于城鎮及周邊，是城區擴張造成的；持續基本不變的面積占30.19%；波動變化占整個流域面積的44.86%，是占比最大的區域；緩慢退化區和緩慢增長區域是波動變化的主要發生地；持續緩慢增長占17.67%，與Hurst指數0.5～0.6分布基本一致；持續較快增長占6.66%，持續快速增長面積僅占0.47%，二者占比較小，多分布于持續緩慢增長區內。

圖6 1999-2018年NDVI Hurst指數

圖7 未來變化趨勢

3.3 植被NDVI空間分布變化

在ArcGIS中通過疊加分析得到了1999—2018年植被覆蓋空間轉移矩陣表(表3)，結果表明，1999—2018年期間高植被覆蓋度和較高植被覆蓋度大幅增加，中植被覆蓋度增加明顯，主要是由低一級植被覆蓋度轉化而來。演化趨勢為：較低植被覆蓋度→中植被覆蓋度，中植被覆蓋度→較高植被覆蓋度，較高植被覆蓋度→高植被覆蓋度。20 a來，植被覆蓋度保持穩定的區域占31.6%，持續增加的占65.99%，減少的占2.41%。總體來看，黃河流域植被增長明顯，生態恢復成效顯著，生態環境整體向好。

表3 黃河流域植被覆蓋空間轉移矩陣 %

3.4 NDVI空間分布的地理探測

3.4.1 不同驅動因子對NDVI空間分異的驅動分析 用地理探測器對NDVI的空間分異及變化趨勢進行歸因，因子探測結果見表4，不同因子對黃河流域NDVI空間分布的解釋能力如下：年平均降水量>濕潤指數>干燥度>土壤類型>年均溫>土地利用>積溫(≥0℃)>積溫(≥10℃)>高程>地貌>坡度>人口>GDP>坡向，這也反映了地處內陸半濕潤、半干旱、干旱交匯區的黃河流域，水分為主的氣候因素是影響NDVI的主導因素。其中，年平均降水量、濕潤指數、干燥度、土壤類型的解釋力均超過30%，分別為59.6%，46.1%，34.9%，33.9%，是主要影響因子；年均溫、土地利用、積溫(≥0℃)、積溫(≥10℃)、高程、地貌解釋力大于20%，分別為27.9%，26.9%，26.2%，25.7%，24.2%，20.7%，是次要影響因素；坡度、人口、GDP、坡向解釋力均低于20%，影響力較低。

表4 因子探測結果

3.4.2 交互作用探測分析 由表5可知，交互作用解釋力最強的是平均降水量和土地利用，它們雙因子交互q值達到了0.704，這說明了土地利用的變化影響降水量,并以此作為因變量解釋黃河流域NDVI較強。解釋力最差的是坡向與GDP的交互，這可能是由于本身坡向及GDP對NDVI解釋力較差，在這樣兩個因子交互下，解釋力不足。同時，降水量與大多數因子相結合都能產生較高的q值，最低值是平均降水量與坡向的組合，卻也達到了0.638。坡向、GDP、人口密度與多數因子交互解釋力較低。

表5 交互作用探測結果

由表6可知，絕大多數因子交互呈現相互增強或非線性增強，表現為兩因子q值相加大于單因子q值。其中以年平均降水量為代表的主導因子表現更為明顯，X4∩X12(0.704)>X4∩X11(0.698)>X4∩X10(0.688)>X4∩X7(0.687)>X4∩X6(0.669)>X4∩X1(0.665)>X4∩X2(0.654)>X4∩X8(0.640)>X4∩X14(0.632)>X4∩X5(0.627)>X4∩X13(0.624)>X4∩X3(0.615)>X4∩X9(0.609)，結果表明，年平均降水量與土地利用類型、土壤類型、地貌、高程、年均溫等因子的相互作用呈現相互增強的關系。X9∩X4(0.609)>X9∩X5(0.479)>X9∩X3(0.368)>X9∩X11(0.359)>X9∩X6(0.295)>X9∩X12(0.290)>X9∩X1(0.279)>X9∩X2(0.276)>X9∩X7(0.254)>X9∩X10(0.224)>X9∩X8(0.128)>X9∩X14(0.114)>X9∩X13(0.047)，且X9∩Xn均大于X9+Xn，可見坡向與年平均降水量、濕潤指數、干燥度、土壤類型等因子的相互作用效應呈現非線性增強的關系。

研究結果表明，僅有兩對因子交互對NDVI的影響是相互獨立的，分別是年平均降水量與土地利用類型交互及高程與人口密度的交互(表6)，X4∩X13=0.624=X4+X13，X7∩X14=0.339=X7+X14。

3.4.3 生態探測分析 由表7可知，氣候要素中各因子大多與其他氣候要素因子存在顯著性關系，年平均降水量尤為突出，而與地表因素和人類活動因素極少存在顯著性關系。比如積溫(≥0℃)與干燥度、年平均降水量、濕潤指數和土地利用類型存在顯著性關系，年平均降水量與積溫(≥0℃)、積溫(≥10℃)和干燥度存在顯著性關系。在地表因素中，大多數因子均與土壤和土地利用類型之間存在顯著性差異，比如高程與土壤類型和土地利用類型存在顯著性關系；坡向與地貌、土壤類型、土地利用類型和人口密度有顯著性差異；地貌與坡度、坡向、土壤類型和土地利用類型存在顯著性差異。由此可見，NDVI受三大類要素共同影響，氣候因素仍為主導。其中，年平均降水量、土壤和土地利用類型對NDVI影響最大，GDP對NDVI影響最小。

表6 影響植被NDVI主導因子交互作用分析

表7 生態探測結果

3.4.4 影響因子的適宜范圍或類型 基于風險探測，我們能夠得出各類影響因子中，植被NDVI最適宜的范圍(表8)。氣象因子方面，年平均降水量在658～1 039 mm，年均溫-2.6～1.8℃時，NDVI均值最高，可見水熱條件較好的地區植被覆蓋度更好。高程在3 216～4 005 m，坡度在15.29°～27.56°時植被覆蓋度最高，與地貌的大起伏中山、大起伏高山，土地利用的林地、草地相對應，說明在山高坡陡的地區，植被覆蓋度更高。草甸土和砂姜黑土廣泛分布于山地和山前洼地中，從側面反映出海拔對植被覆蓋的影響。

表8 不同影響因子的適宜類型或范圍

總體來看，由風險探測我們可知，氣象因素、地表因素是影響黃河流域NDVI的主要因素，而人類活動對NDVI的影響主要體現在土地利用方式的改變上，退耕還林與植樹造林工程發揮了重要作用。

4 討 論

黃河流域地跨我國三級階梯，地貌類型復雜多樣，自然與人為因素共同作用，對流域內植被覆蓋產生重要影響。本研究表明黃河流域自實施退耕還林、生態修復工程實施以來，流域NDVI變化趨勢以極顯著緩慢增長為主，植被覆蓋空間格局整體向好，而氣象因素、地表因素與人類活動因素共同影響了植被分異，其中氣象因素依舊為主導因素，平均降水量與濕潤指數解釋力均大于45%，占主導地位。由此可見降水是植被活動的主要驅動力，這與前人的研究結果一致[24-25]。土壤類型是植被NDVI變化的重要因子[20]，在以降水為主要限制因素的地區，土壤類型對植被生長和雨水再利用效率影響顯著，如淋溶土、半淋溶土地區的NDVI值較高；而漠土與干旱土NDVI值較低，不適宜植被生長；初育土廣泛分布于陜北黃土高原區，該區是黃河流域生態修復與水土保持的重點區域，因此NDVI增速較快。人類活動中，土地利用方式對NDVI影響較大，自1999年退耕還林還草工程實施以來，黃河流域NDVI持續增長，但增速較慢，這可能與城市化和過度放牧有關。

氣候因素與地表因素、人類活動因素相交互，會顯著提高對NDVI的解釋力。本文研究發現，降水量與土地利用解釋力最強，q值為0.704，降水量與其余各因子交互q值均大于0.6，濕潤指數與其余各因子交互q值均大于0.47，也反映了水分的主導作用。這與相關研究[26]認為的社會經濟因子與自然因子交互，會顯著提高其對植被NDVI空間分布的影響力一致。

像元尺度上的區域NDVI時空變化及其驅動因素分析一直是地理學和生態學研究的重點，本文利用SPOT NDVI數據可觀揭示了黃河流域20 a間植被覆蓋度變化的時空特征，本研究表明黃河流域NDVI變化的可持續性不強，主要以波動變化為主，與其他已有研究的相同時期對比，基本趨勢是一致的[13]。已有研究表明，SPOT系列衛星NDVI數據在光譜波段設計、空間精度和幾何糾正等方面都比其他數據更有優勢，因此其應用也較為廣泛[27]。因此，本研究是可信的。但是NDVI存在飽和現象，特別在夏季或植被覆蓋度高的林區尤為明顯。EVI能很好地反映植被生長，即使在生物含量較高的月份也能反映出植被生長的季節性特征[28]，因而EVI更適用于長序列季節性研究。本文大多基于年際NDVI數據展開分析，因此受飽和現象影響較小。

地理探測器是探測空間分異性，以及揭示其背后驅動力的一組統計學方法[18]。地理探測器在各領域應用廣泛，特別是在對NDVI的探測上。本文采用地理探測器對黃河流域NDVI空間分異特征展開探測，結論與同流域研究基本一致[29]。但是，在選取因子指標和對各因子數據進行離散化處理時，存在一定的主觀性。

5 結 論

(1) 1999—2018年黃河流域植被覆蓋整體較好，多年平均NDVI值為0.39，呈現南高北低的分布格局，其中較高及以上植被覆蓋度區域占流域面積的47.13%。從NDVI空間變化特征來看，1999—2018年期間，高植被覆蓋度和較高植被覆蓋度大幅增加，中植被覆蓋度增加明顯，主要是由低一級植被覆蓋度轉化而來。

(2) 20 a來黃河流域NDVI整體呈增加態勢，年增率為0.007；季節變化特征明顯，最低值出現在冬季的2月份，NDVI值為0.16，最高值出現在夏季的8月份，NDVI值為0.53。從NDVI像元尺度的時間變化趨勢來看，1999—2018年黃河流域NDVI變化趨勢以極顯著緩慢增長為主，占整個流域面積的53.12%，整體變化趨勢向好，但可持續性不強；極顯著較快增長和快速增長的區域占15.47%，主要位于山西西部呂梁山脈和陜西北部黃土高原，這與當地生態恢復工程緊密相關。需要注意的是，未來持續波動區域所占面積比例達到44.86%，其中以緩慢退化和緩慢增長為主。

(3) 黃河流域的NDVI空間差異主要由年平均降水量、濕潤指數、干燥度、土壤類型決定，其因子解釋力均超過30%，這表明氣候因素仍然是影響黃河流域NDVI的主導因素。絕大多數驅動因子間交互呈現相互增強或非線性增強，交互解釋力最強的是年平均降水量和土地利用，q值為0.704。只有年平均降水量與土地利用類型交互及高程與人口密度的交互是相對獨立的。