任天晨,陳軍鋒,劉 楠,范 強
(太原理工大學水利科學與工程學院,山西 太原 030024)
植被作為生物圈生態可持續的重要評價因子,其生長變化不僅指示了氣候變化如何影響生物圈,同時也可以通過植被蒸騰作用、地表反照率和粗糙度影響陸地和大氣之間能量交換,對區域防止水土流失、調節氣候等有重要影響[1-3]。然而隨著氣候持續變暖和人類活動增強,對植被生長產生了顯著影響[2]。因此,全球已經實施了許多與植被恢復相關的項目來改善植被生長狀況[2-4]。黃河流域作為黃河水系從源頭到入海所影響的重要地理生態區域和我國灌溉農業發達地區,直接養活著1.07億人口,另外還有4 億人口依賴于黃河流域的水生活[4]。但黃河流域生態本底極為敏感脆弱,為了恢復黃河流域生態環境,自1999 年我國實施了一系列生態恢復工程并取得了階段性成果[4-6]。但有研究指出在生態工程實施過程中沒有綜合考慮區域氣候、人文等因子,目前高植被覆蓋區域的植被(主要為林地)已經達到區域可容納的最高閾值而呈退化趨勢[7-9],這可能會影響黃河流域生態可持續發展。因此,為了黃河流域植被恢復的長期成功,需要明確黃土高原植被覆蓋變化的轉折點及不同季節植被覆蓋變化的主要驅動因子,為黃河流域的生態恢復工程效應的評價提供一定參考。
歸一化植被指數(Normalized Difference Vegetation Index,NDVI)是目前最常用的反映大尺度地表植被狀況和變化的指數[8-11]。曾有研究發現20世紀80年代以來,全球植被生長季內平均NDVI和年最大NDVI均呈增加趨勢,植被覆蓋增長速率、幅度均以顯著增加趨勢為主,這一趨勢在植樹造林的黃河流域最為顯著[12,13]。其中,有研究發現植被覆蓋增加的主要驅動力是二氧化碳施肥效應,其次為碳沉降、氣候變化和人類活動,黃河流域植被變綠關鍵因素主要為植樹造林[14-16]。然而有研究也發現隨著黃土高原植被覆蓋增加,植被蒸騰作用會加強水循環過程,從而影響區域土壤水分、地表徑流和地表干燥化程度等,反過來抑制黃河流域植被生長[9,11]。如Fu 等[14]通過對黃土高原研究發現通過大規模的生態恢復工程減少了黃河泥沙量,生態環境也得到改善,但由于前期的植被恢復工程缺乏科學的指導,部分地區植被恢復不合理規劃導致土壤水分缺失、徑流減少和生態系統退化[9,14]。同時,有研究也發現新種植的植被使蒸發量增加,導致土壤水分減少從而造成土壤表層干旱化程度加?。?5-18],這將進一步影響黃土高原植被覆蓋的可持續發展。Xin 等[19]研究也發現1982-1999 年黃土高原西北部低覆蓋度地區植被覆蓋呈上升趨勢,而部分高覆蓋地區呈下降趨勢,這可能是該區域植被覆蓋已經達到該區域環境的承載能力從而趨向于飽和。其次,由于全球氣候變暖加劇了氣候系統的不穩定性,導致全球極端氣候事件出現的頻率和強度均呈現明顯的增加趨勢,這對黃河流域植被覆蓋的可持續生長產生了顯著影響[15,20]。如Jiang 等[21]發現在氣候變暖背景下極端干旱和極端濕潤交替出現對植被生產力的影響具有1 到5 年的記憶效應,彌補了過去的全球尺度研究強調極端干旱之后植被生產力下降得結論。因此,極端氣候事件對黃土高原植被生長造成的影響不容忽視。雖然目前針對人類活動、平均氣候和極端氣候等因素對植被生長產生的影響已有很多研究,但大部分研究是基于年尺度進行,很少有研究綜合分析不同季節植被生長對人類活動和氣候變化的敏感性程度。
基于此,基于GIMMS NDVI 和MODIS NDVI 數據探究了1982-2020 年間黃河流域不同季節植被覆蓋時空動態,進一步考慮到黃河流域部分地區水蝕和風蝕等影響嚴重并結合其他學者對黃河流域的研究(很多學者研究發現在黃河流域開展的生態恢復項目由于沒有考慮到植被蒸散發等作用,導致區域蒸散發增加,這使得地表干燥化加重。因此,本文探究了氣溫(日最高氣溫最小值、日最低氣溫最大值、日最低氣溫最小值、年均氣溫、日最低氣溫最大值)、降水(每年連續五天降雨量最大值、年均累計降水)、土地利用、風速、海拔、坡度對黃河流域植被生長的影響程度,旨在回答以下問題:①自1982 年以來黃河流域植被覆蓋時空動態如何?②不同季節NDVI 對極端氣候、平均氣候、人類活動和地形因子的響應程度如何,不同因子的相互作用對黃河流域植被生長影響如何?
黃河流域是我國第二大流域,西起巴顏喀拉山,東至渤海,北臨陰山,南到秦嶺,流域面積達7.52×105km2。地勢西高東低,從西到東橫跨青藏高原、內蒙古高原、黃土高原和黃淮海平原 4個地貌單元。地勢總體呈西高東低的態勢,西部河源地區平均海拔超過4 000 m,中部海拔介于1 000~2 000 m[1]。地貌復雜多樣,西部多山,常年積雪,中部為黃土高原地貌,水土流失嚴重,東部由黃河積平原組成。流域內不同地區氣候差異顯著,從西到南依次為半濕潤氣候、半干旱氣候和干旱氣候[3]。全年日照時間長、輻射強,氣溫日較差大,年較差小。降水量少而不均,濕度小,蒸發量大,多冰雹、沙暴天氣。土地利用類型主要為草地和農用地,其次為未利用地、林地和水域(圖1)。

圖1 黃河流域海拔和坡度空間分布圖Fig.1 Spatial Distribution of Altitude and Slope in the Yellow River Basin
(1)NDVI數據:本文使用了來自西部環境與生態科學數據下載中心(https:∕∕ecocast.arc.nasa.gov∕data∕)提供的1982-2015年15天8 km 的GIMMS NDVI3g數據集和NASA(https:∕∕lpdaacsvc.cr.usgs.gov∕appeears)提供的2001-2020 年的16 天250 m 的MOD13Q1 產品的NDVI數據集用于監測黃河流域植被覆蓋變化。為進一步去除NDVI噪聲,使用了Savitzky-Golay 濾波對GIMMS NDVI3g 和MODIS NDVI 數據進行時間序列重建,以平滑NDVI因云或其他因素導致的錯誤峰值[15]。并將MODIS NDVI 采樣為8 km 分辨率,然后基于最大值合成法(MVC)將兩種NDVI數據合成為月數據。最后用2001-2015 的MODIS NDVI數據對GIMMS NDVI重建后的數據進行建模,基于建模結果將GIMMS NDVI 數據延長至2020 年。并基于MVC 方法將12月、1月、2月的NDVI數據合成為冬季,將3月、4月、5月的NDVI數據合成為春季,將6月、7月、8月的NDVI數據夏季,將9月、10月、11月的NDVI數據合成為秋季。
(2)土地利用數據:1980、1990、2000、2005、2010、2015 和2020 年土地利用數據來源于中國科學院資源環境科學數據平臺,本文主要選取一級分類,包括草地、林地、耕地、建設用地、水域和未利用地。
(3)DEM 和坡度數據:90 m 的DEM 數據源于中國科學院資源環境科學數據平臺(http:∕∕www.resdc.cn)提供,為與本文NDVI數據分辨率保持一致,將90 m 的DEM 數據重采樣為8 km。坡度數據是基于DEM 數據利用ARCGIS 軟件提供的計算工具計算得到。
(3)氣候數據:本文氣候數據包括平均氣候數據和極端氣候數據。平均氣候[氣溫(TEM)、降水(PER)、風速(WIN)]來源于中國氣象數據網的月值數據(http:∕∕data.cma.cn∕data∕),時間跨度為1982-2020 年。極端氣候數據源于HadEX3 數據集(www.climdex.org),該數據集采用綜合觀測資源來量化晝夜溫度和降水變化,由29 個氣候極端指數組成,這些指數是通過世界溫度及降水站數據計算而來的,綜合反映了極端溫度和降水事件的頻率和強度,該數據被應用于極端天氣事件研究中[29,30]。根據黃河流域溫度和降水的實際情況,本文選擇了5個最能反映溫度日變化范圍[日最高氣溫最大值(TXX)、日最低氣溫最大值(TNX)、日最高氣溫最小值(TXN)和日最低氣溫最小值(TNN)]、長時降水指標[每年連續五天降雨量最大值(RX5DAY)]來分析黃河流域植被覆蓋對極端氣候指標的響應情況。將12 月、1 月、2 月定義為冬季,3 月、4 月、5 月定義為春季,6 月、7 月、8 月定義為夏季,9 月、10 月、11 月定義為秋季,最后利用專業氣象插值軟件ANUSPLINE 將季節性氣象數據插值為空間分辨率為8km的柵格數據集。
1.3.1 NDVI突變點檢驗——Mann-Kendall突變點檢驗
本文采用Mann-Kendall 方法對1982-2020 年的春、夏、秋和冬季的NDVI進行突變點檢驗。首先使NDVI時間序列構造一個X秩序列記為Sk,在時間序列為隨機的假設下[23],定義統計量:
其中,E(Sk)和Var(Sk)分別是Sk的均值和方差。UFk形成UF的特征曲線,通過可信度檢驗來判斷NDVI是否有明顯的變化趨勢。UFk為標準正態分布,在給定顯著性水平α上,取顯著性水平a=0.05,則U0.05 的臨界值為1.96 的絕對值。將UF、UB的曲線和0.05 的顯著水平在同一張圖上描述出來,UF>0,表示序列呈上升趨勢;反之,表明呈下降趨勢,大于或小于0.05 顯著的臨界線的,表示在95%置信水平上有顯著的上升或下降趨勢。NDVI突變開始的時間即UF與UB曲線出現相交且交點位于臨界線之間,則交點對應的時間為NDVI突變點[23]。
1.3.2 NDVI趨勢分析——Sen斜率法
采用Sen 斜率法趨勢分析計算1982-2020 年春、夏、秋和冬季的NDVI變化趨勢[8]。Sen 斜率法是一種非參數統計的計算方法,對測量誤差和離群數據不敏感[23]。具體公式如下:
式中:β為斜率,i和j代表年份,如果β>0,NDVI時間序列具有正趨勢;如果β<0,NDVI序列具有負趨勢。
1.3.3 NDVI對不同因子的敏感性——地理探測器
地理探測器模型由王勁峰等人[24]提出,是將自變量空間分布與潛在因素分布進行比較,適合用于測量空間分層非均質性程度的空間分析方法。本文利用地理探測器度量季節性NDVI對海拔、坡度、土地利用變化、氣候因子的敏感性,用q值解釋這個程度。公式如下:
式中:q是影響因子對NDVI時空變化的解釋力,h是不同驅動因子的分類或分區數據;L為影響因子的樣本數量,Nh和N分別是h和整個區域的單元數;δ2h和δ2是h和整個區域的方差。q值越大,代表該因子對NDVI影響程度越大。并使用交互作用探測識別不同因子之間的相互作用,并評估他們的組合效應以觀測任何一對因素共同作用是否會增加或降低對NDVI空間分布的解釋力。
兩種原始NDVI月均數據的趨勢和相關性極為顯著,兩者的相關系數(R2)為0.926 5,均方根誤差(RMSE)為0.039 5,但GIMMS NDVI 在NDVI極高值和極低值處明顯高于MODIS NDVI[圖2(a)]。從MODS NDVI 對GIMMS NDVI 重建后的月均NDVI來看,兩者的R2為0.980 3,RMSE為0.019 9,在NDVI極值處也較為接近[圖2(b)]。從像元尺度來看,兩種數據間的相關性也較好[R2= 0.917 1,RMSE=0.030 5,圖2(c)],利用MODIS NDVI 數據對GIMMS NDVI 重建后的數據的相關性顯著提高[R2= 0.950 5,RMSE=0.014 9,圖2(d)]。因此,可以使用MODIS NDVI數據對GIMMS NDVI數據進行延長。

圖2 黃河流域 GIMMS 3g NDVI 和 MODIS NDVI 數據的對比Fig.2 Comparison of GIMMS 3g NDVI and MODIS NDVI Data in the Yellow River Basin
圖3 為黃河流域春、夏、秋和冬季NDVI年際變化及其Mann-Kendall(M-K)檢驗曲線。從春季NDVI年際變化圖來看,1982-2020黃河流域春季NDVI值呈增加趨勢,線性變化傾斜率為0.002 7 ∕a,說明黃河流域近40 年春季植被覆蓋呈增加趨勢。其中1982 年和2019 年分別是該地區近40 年NDVI最小和最大的年份,其NDVI值分別為0.251 和0.331。由UF曲線可以看出,春季NDVI為持續上升趨勢,但在2000年之前上升趨勢呈鋸齒狀變化,2000 年之后為穩定上升趨勢。UF和UB曲線相交于2005年,說明2005年是春季NDVI突變的開始,也表明了黃河流域春季NDVI自2005 年之后NDVI上升速率顯著增加。秋季NDVI的UF曲線和春季相似,線性變化傾斜率為0.001 6 ∕a。2002年為秋季NDVI突變年,從曲線來看2002年之后秋季NDVI增加速率顯著大,其增加速率大于春季(UF變化值)。夏季NDVI值在近40年也呈增加趨勢,線性變化傾斜率為0.001 8 ∕a。其中1982 年和2018 年分別是該地區近40 年NDVI最小和最大的年份,其NDVI值分別為0.436 和0.531。由UF曲線可以看出,黃河流域夏季NDVI為持續上升趨勢。UF和UB曲線相交于2008 年,說明2008 年是夏季NDVI突變開始的年份,從UF曲線可以看出2008 年以前NDVI變化速率基本為穩定狀態,但在2008 年之后夏季NDVI迅速增加。冬季NDVI值線性變化傾斜率為0.0014∕a,為四季增加速率最小的季節。其中1985 年和2017年分別是近40年冬季NDVI最小和最大的年份,其NDVI值分別為0.182和0.266。UF和UB曲線相交于2007年,說明2007年是冬季NDVI突變的開始,且在2007年之后NDVI上升速率顯著增加。由UF曲線看出,冬季和夏季NDVI的UF曲線走勢較為一致,但冬季NDVI的UF曲線在2007 年的變化波動性較大,其變化值小于夏季。

圖3 黃河流域NDVI時序曲線和Mann-Kendall(M-K)突變點檢驗曲線圖Fig.3 NDVI time series curve and Mann Kendall (M-K) mutation point test curve in the Yellow River basin
總體而言,夏季NDVI增加速率最快,其次為秋季,再為春季和冬季。四季的NDVI突變年份均在2000 年之后,且突變年后NDVI增加速率均為上升趨勢,這也進一步說明黃河流域1999年之后實施的生態恢復工程對植被覆蓋度增加作用顯著。
從1982-2020年黃河流域NDVI不同季節的平均值發現(圖4):干旱區NDVI值最低,其次為半干旱區,濕潤區NDVI值最高,其次為半濕潤區,這表明NDVI的空間分布與水熱資源分布一致。春季NDVI最高值集中在半濕潤區和濕潤區,海拔高度基本在1 000 m 以上,NDVI 值大于0.4,植被覆蓋度較高。37.33%的NDVI值小于0.2,集中分布在干旱和半干旱地區,極端最小值分布在河套平原地區(NDVI<0.1)。但在半干旱的劉家峽地區的NDVI集中在0.3~0.4,高于周圍地區。內蒙古、寧夏、山西、陜西和甘肅北部地區的NDVI值基本在0.2 以下,其植被覆蓋度較低。而在其他地區NDVI值基本在0.2 以上,且植被覆蓋度較高,尤其在陜西省的南部的植被覆蓋有最高值,而在內蒙古的毛烏素沙漠的植被覆蓋度有最低值,基本在0.2 以下。究其原因可能是因為在山西、陜西、甘肅隴東及河南地區等地區林地分布范圍較大,土地荒漠化程度較黃土高原北部地區的內蒙古、寧夏、青海等地較小,植被綠化程度較高。夏季NDVI值較春季均有增加,尤其在黃河流域源頭地區,NDVI基本大于0.4。半干旱區的NDVI也從春季的小于0.2 增加到了大于0.2,但在干旱區的內蒙古的NDVI仍小于0.2,河套平原部分地區仍小于0.1。秋季NDVI與春季NDVI空間分布較為一致,但在黃河源區的NDVI值大于春季,集中在0.4 以上。冬季黃河流域除半濕潤和濕潤區的植被覆蓋度較高,57.530%的地區NDVI值在0.2以下,是四季植被覆蓋度最低月。

圖4 黃河流域NDVI空間分布和空間分布區間統計圖Fig.4 Spatial distribution and spatial distribution interval statistics of NDVI in the Yellow River basin
從40 年間黃河流域不同季節NDVI變化趨勢的空間分布、變化趨勢等級所占面積百分比(圖5)發現:近40年四季NDVI均為增加趨勢,尤其在陜西北部、山西西部及甘肅隴東地區增加速率較大,但不同季節的NDVI增加速率存在顯著的空間差異。春季NDVI僅有9.2%的區域為減小趨勢,集中在小浪底、河套平原和毛烏素沙地區。1.6%的區域NDVI增加速率大于0.005 ∕a,集中在黃河流域下游地區。其他地區的NDVI變化速率為干旱區<半干旱<濕潤地區,變化趨勢與水熱資源顯著相關。夏季NDVI的在黃土高原地區的變化速率集中在0.001 ∕a 以上,但在黃河源區的NDVI變化速率集中在-0.001 ∕a~0.001 ∕a,且大部分區域NDVI為減小趨勢。在小浪底、三門峽和壺口地區的植被NDVI也為減小趨勢。秋季8.1%的區域NDVI為減小趨勢,集中分布在三門峽地區,其他區域的變化情況與夏季較為一致。冬季24.9%的區域NDVI為減小趨勢,主要分布在干旱和半干旱地區,是四季中NDVI呈減小區域最多的季節,僅有16.4%的區域NDVI增加速率高于0.001 ∕a,其他區域NDVI基本為不變趨勢。

圖5 黃河流域NDVI變化趨勢圖Fig.5 NDVI Change Trend in the Yellow River Basin
從1980、1990、2000、2005、2010、2015 和2020 年土地利用分布情況(圖6、圖7)發現:1980-1900 年間草地和水域面積的減少主要轉換為耕地和建設用地,林地和未利用地基本沒有變化。1990-2000 年間草地面積有一定的增加(增加了2.649%),耕地、林地和水域面積分別減少了2.798%、0.045%和0.073%,建設用地面積仍為持續增長趨勢,增長了0.14%。2000-2005年間減少的耕地面積主要轉換為草地、建設用地和未利用地。這主要可能與1999開始的退耕還林(還草)政策有關,將坡耕地退耕還林還草、宜林荒山荒地造林,因此在2000 年之后四季植被覆蓋均出現轉折點。2005-2015 年間草地面積減少了1.72%,林地和耕地面積分別增加了0.42%和1.8%。2015-2020年間雖然林地和建設用地面積分別增加了0.15%和1.24%,但耕地和未利用地面積減少了2.90%和1.24%,使草地面積增加了0.87%。據統計數據顯示這與目前的農村人口減少和保護生態環境的意識增加,使植被覆蓋在2015-2021 年間有所改善[16]。有研究發現這是因為特定生態環境變化和人類活動導致部分草地轉為耕地、林地和其他土地[18],這主要是由于人口快速增長和三北防護林體系工程的實施,建設用地的發展,導致許多草地遭到破壞,并且長期過度放牧加上氣候變化和嚙齒動物破壞的影響,長江源區天然草地退化嚴重,將不可避免地導致長江源區荒漠化加劇,且研究發現黃河源區植被覆蓋度呈微弱的減小趨勢[14]。然而,隨著政府實施了一系列移民安置項目,導致黃河流域建設用地持續占用農田、草地和濕地。因此,小浪底等地區植被覆蓋為退化趨勢。

圖6 黃河流域土地利用變化圖Fig.6 Land use change map of different land use areas in the Yellow River basin

圖7 黃河流域不同土地利用面積百分比統計圖Fig.7 Percentage statistics of different land use areas in the Yellow River basin
海拔和坡度會影響區域水熱條件、植被類型等,因此基于參考文獻[22]將黃河流域海拔劃分為極低海拔區(<1 000 m)、低海拔區(1 000~2 000 m)、較低海拔區(2 000~3 000 m)、中低海拔區(3 000~4 000 m)和高海拔區(>4 000 m),將坡度分為平地(0°~2°)、平坡(2°~6°)、緩坡(6°~15°)、斜坡(15°~25°)和陡坡(>25°)。然后統計了黃河流域四季NDVI變化速率在不同分區的分布狀況,發現NDVI變化速率在不同高程梯度中存在顯著差異,總體表現為隨著海拔高度增加,NDVI增加速率減小。其中,春季NDVI變化速率在極低海拔區和低海拔區分布相對分散,均值分布在0.001 8 ∕a。海拔高于2 000 m 之后,春季NDVI變化速率集中在0.001 6~0.001 8 ∕a,分布相對集中。夏季NDVI變化速率在小于3 000 m 的地區均分布較春季分散,且隨著海拔高度增加其增加速率逐漸降低。秋季NDVI變化速率在海拔低于3 000 m 的時候下降較慢,在海拔為3 000~4 000 m 時,NDVI變化速率突然下降,之后隨著海拔高度增加NDVI變化速率呈穩定態勢。主要原因是因為海拔高于3 000 m 的區域人類活動較少,植被覆蓋度相對較高,但在海拔超多4 000 m 時,已遠越過林線,植被覆蓋度較低,且該區域的生態恢復工程實施困難,因此NDVI變化相對穩定,且增加速率最低。海拔小于2 000 m的區域NDVI增加速率較大,植被覆蓋改善效果明顯,主要是因為自1999 年之后在該區域大面積實施退耕還林還草等政策,使植被改善效果最為顯著。冬季NDVI增加速率在海拔小于1 000 m 時較高,均值為0.000 3 ∕a。在海拔大于1 000 m后,NDVI變化速率基本為0,且變化趨勢較為穩定[圖8(a)]。主要是因為黃河流域大部分地區冬季較為干旱,植被覆蓋度相對較低,因此在海拔高于1 000 m 的地區NDVI的增加速率基本0,植被覆蓋度改善不明顯[25]。

圖8 黃河流域NDVI變化速率隨海拔和坡度變化圖Fig.8 Variation of NDVI change rate with altitude and gradient in the Yellow River basin
從不同坡度區域來看,坡度對春季和冬季NDVI變化速率影響不大,但夏季和秋季NDVI變化速率隨著坡度增加而顯著減小,尤其在坡度大于6°之后夏季NDVI變化速率減小最明顯。主要原因是因為秋冬季人類活動相對較夏秋季少。其次,在坡度較緩區域的人類活動較為密集,坡度較大區域人類活動受到地形限制逐漸減少,植被受到的干擾也就越小,且坡度大于15°的區域多為林地分布區域,土地利用不易發生轉換,這就使隨著波段升高,NDVI的變化速率較為穩定且較小,而在坡度較小的區域的NDVI變化速率較大且分散[圖8(b)]。
通常認為人類活動(如土地利用變化)、地形(如坡度、海拔)和氣候因子(平均氣候和極端氣候)的共同或單獨作用導致植被覆蓋發生變化。雖然目前針對上述的單一因子對植被覆蓋的影響的研究很多,但定量描述不同因子對植被覆蓋度的貢獻大小和因子之間的相互作用對植被覆蓋度的影響的研究較少。因此,本文基于地理探測器的因子探測器和交互探測器分析了年均和季節性氣候、土地利用和地形對植被覆蓋的驅動機制。從因子探測器的結果可知(圖9),近40年不同因子對NDVI的決定力q值大小為:日最高氣溫最小值(TXN)的q決定力>日最低氣溫最大值(TNX)>日最低氣溫最小值(TNN)>年均氣溫(TEM)>土地利用(LUCC)>日最低氣溫最大值(TX)>每年連續五天降雨量最大值(RX5DAY)>風速(WIN)>海拔(DEM)>坡度(SLOPE)>年均累計降水(PER)。由此可見,溫度(平均溫度和極端溫度)和土地利用變化對NDVI的影響大于其他因子,年均降水量對植被NDVI的影響較小,其次為地形。

圖9 不同因子對年均NDVI變化的q決定力值分布圖Fig.9 Distribution of q driving force value of different factors on annual average NDVI change
通過不同因子交互探測結果我們發現(表1),任何兩種因子的交互作用都大于或等于單個因子對NDVI的影響。其中,RX5DAY∩DEM(0.72)、RX5DAY∩SLOPE(0.70)、TXN∩DEM(0.74)、TXN∩TEM(0.71)、TNN∩DEM(0.27)的交互作用對NDVI的影響最大,影響因子大于0.7。DEM∩SLOPE(0.27)的交互作用對NDVI的影響最小。極端降雨和地形、極端氣溫和地形的共同作用對NDVI的影響大于其他因子的交互作用對NDVI的影響。這主要是因為黃河流域的大面積的土壤為黃土層,黃土土質松散,垂直節理發育,遇水則容易溶解,因此在坡度較大的地區極容易造成泥石流、滑坡等自然災害,這對海拔較高的地區的林草地破壞極為嚴重(表1)[21]。其次,IPCC 第五次報告指出1998年之后全球地表溫度由快速增溫進入停滯期,這對地處干旱、半干旱和半濕潤的黃河流域植被生長來說是有利的,地表溫度增溫的停滯有利于減小土壤表層水分的蒸散發,從而使土壤含水量提高[8,11,22]。因此,土地利用變化和氣候等的綜合影響使NDVI發生了突變。

表1 黃河流域年間NDVI系數Tab.1 NDVI Interaction in the Yellow River Basin
從不同季節的驅動因子隨時間的變化特征來看(圖10),發現TEM和PER在1980-2015 年間對植被NDVI的q決定力呈下降趨勢,不同季節的極端氣候指標對NDVI的影響情況在2000年達到最低,2010年次之,其他年份的影響均較高。而LUCC和SLOPE的影響力呈增加趨勢,且LUCC的q值增加幅度較SLOPE大。以上原因主要是因為2000年之前黃河流域NDVI的變化很大程度上決定于氣候因子的變化,但從1999年之后開展了一系列工程(如退耕還林、營造林、草地治理、耕地治理和移民搬遷等)使LUCC發生了顯著變化,且其中的很多工程屬于生態恢復工程,這些工程的實施極大的改變了該地區植被覆蓋狀況。因此在2000 年之后的植被NDVI變化對LUCC的變化較為敏感。其次,夏季的SLOPE的q值大于春、秋和冬季,有研究發現主要是因為黃河流域的黃土高原大部分地區夏季容易形成及時雨、大暴雨等天氣,而由于黃土層的土壤性質,導致在SLOPE較大的地區發生滑坡、泥石流等自然災害,從而對植被覆蓋產生顯著的影響。以上結論與李曉麗等[25]研究結論一致。

圖10 不同因子對季節性NDVI變化的q決定力值分布圖Fig.10 Distribution of q driving force values of different factors on seasonal NDVI changes
基于GIMMS NDVI 數據分析了1982-2020 年間黃河流域植被覆蓋的季節性變化和土地利用、地形、極端氣候和平均氣候對季節性植被NDVI的貢獻度和不同因子的交互作用對NDVI的影響。主要結論如下:
(1)1982-2020 年黃河流域干旱區NDVI值最低,其次為半干旱區,濕潤區NDVI值最高,其次為半濕潤區,即NDVI的空間分布與水熱資源分布一致。
(2)1982-2020 年間黃河流域整體以增加趨勢為主。春季NDVI增加速率為0.002 7 ∕a,其中1982 年和2019 年分別是該地區近40年NDVI最小和最大的兩個年份,2005年NDVI變化的突變年。夏季NDVI增加速率0.001 8 ∕a,2008 年是NDVI突變年,且2008 年后夏季NDVI迅速增加。秋季NDVI增加速率為0.001 6 ∕a,2002 年為秋季NDVI突變年。冬季NDVI增加速率為0.001 4 ∕a,2007年是NDVI突變年。
(3)黃河流域NDVI隨著海拔高度增加,NDVI增加速率減小,且夏季的這種趨勢最明顯。坡度對春、冬季NDVI變化速率影響不大,但夏、秋季NDVI變化速率隨著坡度增加而顯著減小,尤其在坡度大于6°后夏季NDVI變化速率減小最明顯。
(4)近40 年不同因子對NDVI的q影響大小為:日最高氣溫的最小值(TXN)的q決定力>日最低氣溫最大值(TNX)>日最低氣溫最小值(TNN)>年均氣溫(TEM)>土地利用(LUCC)>日最低氣溫最大值(TXX)>每年連續五天降雨量最大值(RX5DAY)>風速(WIN)>海拔(DEM)>坡度(SLOPE)>年均累計降水(PER)。其中,任何兩種因子的交互作用都大于或等于單個因子對NDVI的影響。且TEM和PER在1980-2020 年對植被NDVI的決定力呈下降趨勢,而LUCC 和SLOPE的影響力呈增加趨勢,且LUCC的q值增加幅度較SLOPE大。極端溫度和降水對其影響沒有明顯的時間趨勢。