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

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

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

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

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

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

1 研究區(qū)域概況

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

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源

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

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

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

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

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

(1)

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

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

(2)

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

2.2 研究方法

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

(3)

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

坡面NDVI同樣采用面積加權(quán)法計算，計算公式如下：

(4)

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

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

NDVIgully=a+b×NDVIslope

(5)

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

2.2.2 趨勢分析法 本文采用一元線性趨勢分析法分析整個延河流域20 a的NDVI變化趨勢?；谙裨叨葦M合NDVI變化，并整合每個像元的時變特征和區(qū)域空間變化，以反映研究區(qū)域NDVI的20 a時空格局演變[12]。一元線性趨勢線斜率的計算公式為：

(6)

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

表1 NDVI趨勢變化水平

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

(7)

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

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

(8)

則t點處為突變點。

同時計算統(tǒng)計量：

(9)

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

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

NDVIPRE=c+d×PRE

(10)

NDVIHA=NDVIOBS-NDVIPRE

(11)

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

3 結(jié)果與分析

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

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

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

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

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

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

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

3.2 NDVI變化驅(qū)動力分析

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

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

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

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

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

圖5 坡面NDVI與溝壑區(qū)NDVI線性相關(guān)

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

4 討論與結(jié)論

4.1 討 論

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

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

圖6 2000-2019各年坡面NDVI與溝壑區(qū)NDVI的線性關(guān)系

4.2 結(jié) 論

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

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

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