基于長時間序列NDVI的黃土高原延河流域及其溝壑區植被覆蓋變化分析

賈云飛, 李云飛, 范天程, 曾 競, 趙建林

(長安大學 地質工程與測繪學院, 西安 710054)

黃土高原位于我國中部偏北部，屬半干旱大陸性季風氣候區，氣候變化敏感，生態環境脆弱[1]，是中國水土流失的嚴重區域之一。嚴重的水土流失會影響社會經濟的發展，高泥沙輸移會威脅黃河下游地區[2]。相關研究表明黃河近90%的輸沙量來源于黃土高原[3]。

為了控制水土流失和高輸沙量，過去50 a中國在黃土高原實施了一系列水土保持措施：其中包括坡面區域的大量梯田建設和生態修復項目，以及溝壑區域大量淤地壩的建設[4]。其中最為顯著的生態修復項目是自1999年開展的“退耕還林(草)”。相關研究表明自“退耕還林”政策實施以來，黃土高原植被覆蓋度顯著提高，植被狀況明顯好轉；郭敏杰等[5]研究表明1982—2006年黃土高原植被覆蓋度呈現穩定的增加趨勢，增速為0.075%/a；黑哲[6]研究表明2000—2019年超過90%的區域植被NDVI呈增加趨勢，其中71.56%的地區呈顯著增加；郭永強等[7]研究表明實施“退耕還林”(2000—2015)后，植被覆蓋率年均增幅為0.59%，年植被覆蓋率從1987年的41.78%提高到2015年的53.23%。胡春宏等[8]研究表明新中國成立后，經過治理，入黃沙量由1919—1959年16億t/a銳減至2000—2018年約2.5億t/a。黃土高原具有典型的丘陵溝壑地貌，而侵蝕過程和泥沙主要來源于溝壑區與坡面。其中坡面以面蝕、細溝侵蝕以及耕作侵蝕為主，而溝壑區主要以溝蝕和重力侵蝕為主，研究表明在黃土高原流域中，溝壑間侵蝕產沙大于坡面產沙。Zhao等[9]研究表明當溝壑區的密度大于30%時，溝壑區對流域泥沙的貢獻大于75%。以往的研究表明坡面的植被覆蓋情況有了明顯的改善，但坡面的植被覆蓋情況改善是否會誘發坡面間的溝壑區植被覆蓋增加，目前還尚未有研究表明。

本文以黃土高原退耕還林重點區域延河流域為研究對象，基于Google Earth(GE)平臺采用系統樣本法和人工交匯法提取該流域的溝壑區與坡面地貌；基于長時間序列NDVI數據，采用趨勢分析法、Pettitt突變點分析、殘差分析法以及地學統計方法，在分析延河流域過去20 a(2000—2019年)整體的植被覆蓋變化情況以及降雨量和人類活動對NDVI變化的影響的基礎上，分析黃土高原坡面和溝壑區的植被覆蓋變化情況，探索坡面與溝壑區植被變化的差異性及相關性。

1 研究區域概況

本文的研究區域是延河流域，延河是黃河的一級支流，全長286.9 km。延河流域(108°45′—110°28′E，36°23′—37°17′N)位于黃土高原中部，流域面積7 725 km2，地勢西北高、東南低。延河流域屬于黃土丘陵溝壑區，地表破碎，溝壑密度在2.1～4.6 km/km2，水土流失嚴重。流域的植被區劃屬于森林草原地帶，主要是人工種植的次生植被和干旱草本植物。該流域屬暖溫帶半濕潤地區，年平均降水量約為520 mm，降水主要集中在6—9月份，約占全年降雨量的75%。

2 數據來源與研究方法

2.1 數據來源

2.1.1 NDVI數據 本文使用歸一化植被指數(Normalized Difference Vegetation Index，NDVI)來分析植被覆蓋變化。NDVI與植被的初級生產力、葉覆蓋、生物量具有很好的相關性，已被廣泛用于量化植被的生長趨勢和過程[10]。本文使用的NDVI數據為MODIS 16 d合成的NDVI產品MOD13Q1，空間分辨率為250 m，時間跨度為2000—2019年，由Google Earth Engine平臺提供，使用最大值合成法合成影像，得到延河流域各年的NDVI最大值柵格數據集，并將NDVI像元初始值轉化到-1～1[11]。

2.1.2 降雨數據 為分析降雨對延河流域NDVI植被覆蓋變化的貢獻，本文采用1 km分辨率的降雨數據集，該數據來源于中國科學院資源環境科學與數據中心的中國1980年以來逐年年降水量空間插值數據集，該數據集是基于全國2 400多個氣象站點日觀測數據，通過整理、計算和空間插值處理生成。本文使用了該數據集2000—2019年的降雨數據。

2.1.3 溝壑區樣本提取與柵格化 為了區分溝壑區和坡面植被覆蓋變化情況，本文首先在延河流域提取溝壑區樣本，提取過程如下：首先，在延河流域系統均勻地分布272個3 km×3 km正方形樣本區域，該樣本覆蓋整個延河流域。然后基于Google Earth平臺，人工勾繪出每個樣本區內的溝壑區域，并假設樣本區內除溝壑區域外就是非溝壑區域即坡面區域。在此基礎上，對提取的矢量溝壑區和坡面進行柵格化。柵格化過程如下：首先創建與NDVI像元一致的緩沖區，基于緩沖區對溝壑區矢量與坡面區矢量進行分割，使得分割后的每一個區域能夠正好覆蓋所對應的像元，然后計算得到分割后的溝壑區和坡面區內像元對應的矢量面積，再結合像元面積得到每個像元的溝壑區面積占比與坡面面積占比。具體流程見圖1。

注：A為272個樣本區中的某一樣本區；B為基于Google Earth人工勾繪的樣本區內溝壑區域；C，D為溝壑區矢量和坡面區矢量分割結果；E，F為溝壑區和坡面區像元面積占比。

由于本文采用的NDVI數據的空間分辨率為250 m，故柵格化后的溝壑區和坡面的空間分辨率與NDVI像元一致，柵格化后的像元值為對應柵格的溝壑區面積占比(GD)，計算公式如下：

(1)

式中：GD為每個像元值，即溝壑區面積占比；Sgully為每個像元中人工勾繪溝壑的面積(m2)。

坡面面積占比(SD)計算公式如下：

(2)

式中：SD為每個像元值，即坡面面積占比；Sslope為每個像元中非溝壑區即坡面的面積(m2)。

2.2 研究方法

2.2.1 延河流域溝壑區和坡面NDVI 為分析延河流域溝壑區和坡面NDVI變化趨勢，本文分別計算了272個樣本區2000—2019年的樣本區(NDVImean)，坡面(NDVIslope)和溝壑區(NDVIgully)的平均NDVI年最大值。其中，樣本區平均NDVI等于樣本區內各像元的平均值。溝壑區NDVI采用面積加權法計算，計算公式如下：

(3)

式中：i=1,2,…,272；NDVIgully(i)為第i個樣本區的溝壑區平均NDVI最大值；NDVIij為第i個樣本區第j個像元的年NDVI最大值；GDij為第i個樣本區第j個像元對應的溝壑區面積占比。

坡面NDVI同樣采用面積加權法計算，計算公式如下：

(4)

式中：i=1,2,…,272；NDVIslope(i)為第i個樣本區的坡面平均NDVI最大值；NDVIij為第i個樣本區第j個像元的年NDVI最大值；SDij為第i個樣本區第j個像元對應的坡面面積占比。

在此基礎上，本文采用ANOVA分析方法，分析延河流域溝壑區NDVI(NDVIgully)和坡面區NDVI(NDVIslope)是否具有顯著差異(以p<0.05檢驗其顯著性)。以及采用線性回歸方程和R2分析樣本區中兩者相關性即：

NDVIgully=a+b×NDVIslope

(5)

式中：a和b為線性回歸系數，如果b>1則表明樣本區中溝壑區NDVI大于坡面NDVI，反之溝壑區NDVI小于坡面NDVI。

2.2.2 趨勢分析法 本文采用一元線性趨勢分析法分析整個延河流域20 a的NDVI變化趨勢。基于像元尺度擬合NDVI變化，并整合每個像元的時變特征和區域空間變化，以反映研究區域NDVI的20 a時空格局演變[12]。一元線性趨勢線斜率的計算公式為：

(6)

式中：n為研究數據的時間長度，本文為20 a；NDVIj為整個延河流域第j年的NDVI最大值；k為像元回歸方程趨勢線斜率，若k>0，表示研究時段內NDVI變化呈增加趨勢，反之，則呈減少趨勢[13]。為了更好地判斷整個研究區域植被覆蓋的動態變化趨勢，根據k值的大小，可定義7個NDVI趨勢變化水平[14]：重度退化、中度退化、輕度退化、基本不變、輕度改善、中度改善以及明顯改善(表1)。

表1 NDVI趨勢變化水平

2.2.3 突變點分析 為了檢驗延河流域過去20 a植被覆蓋變化是否連續，本文采用Pettitt[15]突變檢驗方法檢驗延河植被覆蓋變化是否存在顯著突變點(p<0.05)。Pettitt突變檢驗是一種非參數檢驗法，其假設樣本容量中存在突變點t，采用Mann-Whitney的統計量Ut,n檢驗樣本序列突變點t前后兩個子樣本，x1,…,xt和xt+1,…,xn二者累積分布是否存在顯著差異。

(7)

式中：xt，xi分別為第t年和第i年的NDVI相鄰時間段斜率樣本；n為數據系列長度；Ut,n為將根據第一個樣本序列超過第二個樣本序列次數的統計組成新的序列。

Pettitt法原假設H0為序列不存在突變點。若t時刻滿足：

(8)

則t點處為突變點。

同時計算統計量：

(9)

若概率p≤0.05，則認為檢測出的突變點在統計意義上是顯著的。

2.2.4 NDVI變化驅動力分析 由于延河流域過去20 a人類活動頻繁，特別是退耕還林項目的實施，同時降雨也會對區域植被覆蓋產生影響。因此，本文采用由Evans等[16]提出的殘差分析方法研究降雨量和人類活動對延河流域NDVI變化的影響及相對貢獻。該方法的主要步驟為：首先基于2000—2019年NDVI和降雨量時間序列數據，以NDVI為因變量，以降雨量為自變量，建立一元線性回歸模型，計算模型中的指標參數(斜率和截距)；基于降雨量與回歸模型的參數，計算得到NDVI的預測值(NDVIPRE)，表示降雨量對NDVI變化的影響；最后計算NDVI觀測值與NDVIPRE的差值，即殘差(NDVIHA)，表示人類活動對NDVI變化的影響。具體計算公式如下：

NDVIPRE=c+d×PRE

(10)

NDVIHA=NDVIOBS-NDVIPRE

(11)

式中：NDVIPRE為回歸模型的NDVI預測值，即降雨量對NDVI變化的貢獻；PRE為年降雨量(mm/a)；NDVIHA為殘差，即人類活動對NDVI變化的貢獻；NDVIOBS為基于NDVI柵格數據集的觀測值；c，d為模型參數。

3 結果與分析

3.1 延河流域NDVI時空變化特征

(1) 2000—2019年延河流域年平均NDVI整體呈現增加趨勢，在研究時段內，延河流域年平均NDVI由2000年的0.415增加到2018年的0.750，其增長率達到44.60%，且年平均趨勢率為1.30%(p<0.001)，相關系數R2=0.544，表明2000—2019年延河流域的植被覆蓋狀況改善較為明顯。基于Pettitt′s突變點分析表明延河流域NDVI突變年份在2008年(p<0.001)，因此將延河流域NDVI變化整體趨勢分為兩個階段：2000—2008年和2009—2019年。2000—2008年延河流域年平均NDVI增加顯著，平均趨勢率為2.00%(p<0.001)，相關系數R2=0.343；2009—2019年延河流域年平均NDVI增加趨勢減緩，平均趨勢率為0.70%(p<0.001)，相關系數R2=0.105(圖2)。

(2) 2000—2019年延河流域NDVI變化趨勢具有一定的空間異質性，延河流域大面積植被覆蓋恢復良好，年平均NDVI呈增加和減小的區域分別占98.90%，0.80%。其中呈現明顯改善(k>0.009)的區域面積占比約為81.27%，約6 314.20 km2；其次是中度改善的區域面積占比約為14.44%，約1 121.84 km2，主要分布在西北部、西南部及東南部；重度退化(k<-0.009)的面積占比約為0.31%，約24.42 km2(表2)，主要分布在中部地區，該地區為延安市區。

圖2 2000-2019年延河流域NDVI年際變化及分段變化

通過對比可以發現2000—2008年延河流域植被覆蓋狀況增加趨勢更加明顯，整個延河流域植被覆蓋情況都有所改善，退化的區域面積可以忽略不計；而2009—2019年植被覆蓋狀況增加趨勢明顯減緩，且具有很大的空間異質性；除西北部、中部以及東南部地區均出現不同水平的退化趨勢外，其他地區的增加趨勢也有所減緩(圖3)。2000—2008年延河流域年平均NDVI呈增加和減小趨勢的區域面積分別占總面積的98.57%，0.90%。其中呈明顯改善(k>0.009)的區域面積占比約為88.85%，約6 903.16 km2；其次是中度改善的區域面積占比約為7.41%，約575.86 km2；退化的面積占比都在0.50%以下。2009—2019年延河流域年平均NDVI變化趨勢具有很大的空間異質性；年平均NDVI呈增加和減小趨勢的區域面積分別占總面積的87.13%，7.33%。其中明顯改善的區域面積占比約為33.70%，約2 617.78 km2；其次是中度改善的區域面積占比約為34.01%，約2 642.71 km2；再次是輕度改善的區域面積占比約為19.42%，約1 508.60 km2；基本不變的區域面積約占5.54%，約430.51 km2；輕度退化的區域面積約占4.65%，約361.20 km2(表2)。

圖3 不同年份延河流域NDVI變化趨勢空間分布對比

表2 2000-2019及分段時期延河流域趨勢變化分布

3.2 NDVI變化驅動力分析

延河流域NDVI變化驅動力分析表明2000—2019年延河流域降雨量對NDVI變化的影響一直保持穩定，平均趨勢率為0.037%/a(p<0.05)，且相關性較弱(R2=0.001)；而2000—2019年人類活動對NDVI變化的影響呈現明顯增加的趨勢，平均趨勢率為1.29%/a(p<0.001)，相關性較強(R2=0.755)。這在一定程度上說明人類活動對延河流域植被覆蓋的改善產生了積極影響。人類活動對延河流域NDVI變化的影響也呈現分段特征，2000—2008年人類活動對延河流域NDVI變化的貢獻明顯大于2009—2019年(圖4)。

3.3 延河流域溝壑區與坡面植被覆蓋變化特征

基于272個樣本區內平均溝壑區NDVI和坡面NDVI分析表明，隨著坡面NDVI的增加，溝壑區NDVI也呈現明顯增加的趨勢，溝壑區與坡面的線性趨勢率為1.02(p<0.001)，即溝壑區植被覆蓋要略大于坡面植被覆蓋，兩者相關系數R2=0.97(圖5)。

然而，兩者相關性存在年際差異，且與延河流域整體NDVI變化趨勢存在關聯。如圖6所示，基于每年樣本區中溝壑區NDVI和坡面NDVI線性相關分析表明，2000—2008年，各樣本區中坡面NDVI與溝壑區NDVI的線性趨勢率均≥1，且相關系數R2均>0.90，即當整體植被覆蓋呈現顯著增加時，溝壑區植被覆蓋要好于坡面區域植被覆蓋，且兩者相關性較強；而2009—2019年，各樣本區中坡面NDVI與溝壑區NDVI的線性趨勢率均<1，即溝壑區植被覆蓋小于坡面區域植被覆蓋，且相關系數R2除了2010年、2011年、2015年外，均小于0.90。

注：NDVIP，NDVIH分別為降水、人類活動對NDVI變化的貢獻。

圖5 坡面NDVI與溝壑區NDVI線性相關

非獨立樣本t檢驗的結果表明溝壑區NDVI與坡面NDVI存在明顯差異(p<0.05)。2000—2019年溝壑區NDVI與坡面NDVI差異的均值為0.012(p<0.001)；2000—2008年溝壑區NDVI與坡面NDVI差異的均值為0.009(p<0.001)；2009—2019年溝壑區NDVI與坡面NDVI差異的均值為0.014(p<0.001)。

4 討論與結論

4.1 討 論

本文研究表明，2000—2019年延河流域整體植被覆蓋改善較快，但具有一定的空間異質性，中部地區城市建設導致植被覆蓋減少，其他地區都植被覆蓋都呈現增加趨勢。延河流域植被覆蓋變化具有分段特征，2000—2008年植被覆蓋恢復速度較快，2009—2019年植被覆蓋恢復速度有所減緩，恢復趨勢趨于平穩。退耕還林還草分為1999年起實施的前一輪退耕還林還草和2014年起實施的新一輪退耕還林還草。1999—2001年是退耕還林還草的試點示范階段，2002年全國全面啟動退耕還林還草工程，2008年起開始鞏固各工程省區的退耕還林成果，為擴大退耕還林、退牧還草規模，2014年開始了新一輪退耕還林還草[17-18]。與此同時降雨和人類活動對NDVI變化歸因分析表明，降雨量對延河流域NDVI變化的貢獻保持穩定，而人類活動對延河流域NDVI的貢獻有著明顯增加的趨勢。1999年來開展的如“退耕還林”等一系列的生態工程是2000—2019延河流域NDVI明顯增加的主要原因。

自退耕還林項目實施以來，延河流域丘陵溝壑區植被恢復態勢較為顯著，坡面植被覆蓋狀況也得到了明顯改善[19]。基于人工提取的大量溝壑區和坡面樣本，本文初步分析了退耕還林項目對延河流域溝壑區和坡面區植被覆蓋影響。初步結果表明樣本區中溝壑區與坡面區的植被覆蓋具有同步性，兩者相關性較強，這表明當坡面植被覆蓋增加時，也間接增加了溝壑區植被覆蓋。這與前人研究一致，相關研究表明溝壑區由于棄耕時間長且人為擾動較小，形成了比較穩定的植物群落，因此在坡面植被覆蓋改善時，溝壑區的植被覆蓋也隨之改善，且溝壑植被覆蓋優于坡面植被覆蓋。然而，坡面上的植被恢復與坡度、水土保持政策密切相關，近年來也得到了極大的改善[20]。本文研究表明，延河流域坡面NDVI與溝壑區NDVI具有較強的相關性，且同樣具有分段特征，坡面NDVI增加的趨勢越快，溝壑區NDVI增加的趨勢也越快，且相關性較強。

圖6 2000-2019各年坡面NDVI與溝壑區NDVI的線性關系

4.2 結 論

(1) 2000—2019年延河流域植被NDVI增加顯著，平均趨勢率為1.30%/a(p<0.001)，其中，中部地區是城市用地NDVI呈減小趨勢；其他地區NDVI都呈現增加趨勢，北部相比南部NDVI增加較明顯。2000—2008年延河流域植被NDVI增加顯著，平均趨勢率為2.00%/a(p<0.001)，且整個延河流域植被覆蓋都呈現改善狀況；2009—2019年植被NDVI增加趨勢減緩，平均趨勢率為0.70%/a(p<0.001)，且西北部、中部和東南部植被覆蓋呈現輕度退化趨勢。

(2) 降雨和人類活動對延河流域NDVI變化的貢獻具有明顯的差異性，2000—2019年降雨對NDVI變化的貢獻保持穩定，且相關系數R2=0.001，相關性較弱；而人類活動對NDVI變化的貢獻較大，人類活動對延河流域NDVI變化的貢獻的平均趨勢率為1.30%/a(p<0.001)，相關系數R2=0.755，存在較強的相關性。

(3) 2000—2019年延河流域溝壑區域NDVI隨著坡面NDVI的增加而增加，年平均趨勢率為1.02(p<0.001)，且相關系數R2=0.967，相關性較強。坡面NDVI與溝壑區NDVI的線性關系也呈現分段趨勢，2000—2008各年坡面NDVI與坡面NDVI的平均趨勢率均>1，且相關性較強，2009—2019各年的平均趨勢率均<1，且相關性較2008年以前有所減弱。