殷崎棟,柳彩霞,田 野

1 陜西生態(tài)環(huán)境規(guī)劃設(shè)計院有限公司, 西安 710075 2 中國科學(xué)院空天信息創(chuàng)新研究院遙感科學(xué)國家重點實驗室, 北京 100101 3 中國科科學(xué)院城市環(huán)境研究所城市環(huán)境與健康重點實驗室, 廈門 361021

陸地植被是大氣、巖石圈、土壤圈、水圈和生物圈相互作用的產(chǎn)物。它是連接氣候變化、人類活動和生態(tài)系統(tǒng)的一個重要環(huán)節(jié)。氣候和人類活動是控制和影響植被空間分布及其變化的基本驅(qū)動力[1-2],辨別氣候和人類活動影響區(qū)域?qū)τ诶斫夂凸芾淼乇碇脖环浅Ｖ匾?特別是在干旱和半干旱地區(qū),因為在這些地區(qū),年降雨量相對較低,年降雨量波動性較大。在辨識人類活動造成的影響程度時,氣候的變異性很難將植被中的自然變化與直接人類活動引起的自然變化分開[3- 6]。鑒于土地退化或者綠化工程可能在幾十年內(nèi)逐步發(fā)生,對植被的觀測時間序列必須一致,因此遙感技術(shù)成為監(jiān)測并分析植被變化的一項可靠技術(shù)。遙感觀測的植被指數(shù),由于具備統(tǒng)一的觀測基準(zhǔn)和長時間序列特征,可以定量描述植被覆蓋和植被生產(chǎn)力的變化。利用歸一化差異植被指數(shù)(NDVI)和氣候因子的關(guān)系,可以研究人為和自然因素引起的植被變化[6- 8]。RESTREND算法是一種從NDVI趨勢中消除氣候影響的方法[9],基于氣候數(shù)據(jù)和植被指數(shù)進(jìn)行土地退化的檢測[3, 9]。早期的RESTREND算法通過計算年最大NDVI、生態(tài)系統(tǒng)生產(chǎn)力和降水之間的線性回歸來控制影響植被的氣候要素的變化[4, 9],然后從線性回歸關(guān)系中預(yù)測的NDVI與觀測NDVI之間的差異,得到NDVI殘差,之后對殘差進(jìn)行線性趨勢分析,該分析結(jié)果被認(rèn)為是人類活動引起的地表植被變化。然而,由于降水量和植被指數(shù)之間的線性關(guān)系存在不穩(wěn)定性。當(dāng)在時間序列的中間發(fā)生退化時,這種線性關(guān)系可能會發(fā)生改變,導(dǎo)致不可靠的結(jié)果[5]。因此,檢測植被變化是否引起植被與降水之間的線性關(guān)系(也稱為斷點)至關(guān)重要。遙感產(chǎn)品在時間序列上的斷點檢測已應(yīng)用于土地覆蓋制圖[10- 12]和森林管理[13-14]。已有研究表明,當(dāng)在植被與降水之間的關(guān)系中斷時檢測斷點時,可以改善植被變化檢測精度[3,5,15]。時間序列分割和殘差趨勢分析(Time Series Segmentation and Residual Trend, TSS-RESTREND)是由Burrell, Evans[6]提出,加入了斷點檢測的RESTREND算法,成功用于檢測澳大利亞全國的植被變化情況[6,8],能夠改進(jìn)對變化區(qū)域和變化方向的檢測。

陜西省地處中國大陸腹地,是水土流失、沙化等環(huán)境問題最嚴(yán)峻的地區(qū)之一。自1999年起,陜西省相繼啟動并實施了生態(tài)環(huán)境建設(shè)綜合治理工程、天然林資源保護(hù)工程、退耕還林工程、重點防護(hù)林工程、水土保持工程和天然草場恢復(fù)與建設(shè)工程等一批重點生態(tài)建設(shè)工程,使陜西省特別是陜北地區(qū)植被狀況有了明顯好轉(zhuǎn),抑制了陜北毛烏素沙漠南擴(kuò)的勢頭[16]。最新研究表明,近10多年來,陜西省部分區(qū)域的植被生長發(fā)生了重大變化。如欒金凱等[17]利用MODIS/Terra NDVI時間序列遙感影像(2000—2016年)和溫度、降水氣候數(shù)據(jù),分析了17年來陜西省榆林市的植被指數(shù)空間分布變化規(guī)律,并且利用復(fù)直線回歸法揭示了氣象因素對植被生長和演化的影響,榆林市一半以上的區(qū)域人類活動對植被生長起到了促進(jìn)作用。鄧晨暉等[18]基于MODIS/Terra NDVI時間序列遙感影像(2000—2015年)、DEM和氣象數(shù)據(jù),采取趨勢分析、多元回歸殘差法、偏最小二乘回歸法,反映了秦嶺地區(qū)植被覆蓋度變化特征,探究其對氣候變化和人類活動的雙重響應(yīng)機(jī)制。在區(qū)域尺度上,更早的研究集中在關(guān)中[19]、陜北[20- 24]等地。在全省尺度上,韓紅珠等[25]利用MODIS數(shù)據(jù)提取了2001—2016年間的植被物候時空特征,結(jié)果表明退耕還林工程效益顯著[26-27]。目前,還沒有相關(guān)研究報道利用公里空間分辨率數(shù)據(jù)從宏觀尺度分析陜西全省的植被綠度時空變化特征及其人類活動影響區(qū)域分析,人類活動對陜西全省植被的影響程度和范圍還缺乏全局了解。本文從全省尺度,細(xì)分為陜北、關(guān)中平原、陜南地區(qū),利用可靠的TSS-RESTREND算法解釋人類活動對陜西省NDVI(2001—2018年)的貢獻(xiàn)影響程度,算法剝離了氣候變化本身對NDVI的貢獻(xiàn)影響,對學(xué)者、管理者摸清陜西全省植被綠度的長時間序列變化具有重要意義。同時,2019年是陜西省全面實施退耕還林生態(tài)建設(shè)工程20周年,利用產(chǎn)品一致性較好的遙感時間序列產(chǎn)品,全面分析近20年來全省植被變化情況和客觀評價生態(tài)建設(shè)工程成效,可為后續(xù)生態(tài)環(huán)境管理提供可靠的本底數(shù)據(jù)。

1 研究區(qū)概況

圖1 研究區(qū)位置意圖Fig.1 Location map of the study area

陜西省地理范圍為105°29′—11°15′E, 31°42′—39°35′ N之間(圖1),地處東部濕潤地區(qū)和西部干旱區(qū)的交界地帶,省內(nèi)氣候類型多樣,各地的年均氣溫在7—6℃,屬于典型的大陸性季風(fēng)氣候區(qū)。省內(nèi)地形地貌類型豐富,地勢特點是南北高,中間低,由西向東傾斜,年降水量受山地地形影響顯著。自北向南可粗略劃分為三大自然區(qū)域,分別為陜北黃土高原區(qū)、關(guān)中平原區(qū)和陜南秦巴山地區(qū)。陜北黃土高原區(qū)包括榆林地區(qū)和延安地區(qū),海拔900—1600 m,地勢較高,地處干旱、半干旱地區(qū),年均氣溫7—11℃,年降水量約400—600 mm,地區(qū)分配不均,季節(jié)變化大。關(guān)中平原區(qū)由渭河干流以及兩岸支流日久沖積、沉積形成,主要包含西安市、銅川市、寶雞市、咸陽市和渭南市,地勢相對較低,海拔在300—800 m之間,中部較為平坦,年均氣溫11—13℃,年降水量約500—700 mm,為溫帶半濕潤氣候。陜南秦巴山地處于秦嶺和大巴山系,海拔200—1200 m,主要包括漢中市、安康市、商洛市,水熱條件相對較好,年均氣溫14—15℃,年降水量約700—900 mm。同時,陜西省地處黃河中游、長江上中游地區(qū),屬于國家生態(tài)環(huán)境建設(shè)重點區(qū)域。

2 數(shù)據(jù)和方法

2.1 數(shù)據(jù)來源與預(yù)處理

本文使用的NDVI數(shù)據(jù)是2001年1月到2018年12月的MODIS NDVI (MOD13A2)數(shù)據(jù),出自美國國家航空航天局NASA發(fā)布的EOS/MODIS數(shù)據(jù)產(chǎn)品,空間分辨率為1 km×1 km,時間分辨率為16 d。MODIS/NDVI產(chǎn)品經(jīng)過了水、云、重氣溶膠等預(yù)處理,保證了數(shù)據(jù)質(zhì)量,被廣泛應(yīng)用于區(qū)域植被覆蓋變化研究。本文數(shù)據(jù)處理在谷歌地球數(shù)據(jù)云平臺 (Google Earth Engine) 中處理,NDVI數(shù)據(jù)獲取ID為MODIS/006/MOD13A2,利用數(shù)據(jù)集的質(zhì)量控制字段(SummaryQA),選取所有最優(yōu)觀測值(SummaryQA等于0)的最大值作為月NDVI。

降雨數(shù)據(jù)為1999—2018年陜西省及其周邊省份33個氣象臺站的月數(shù)據(jù),來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn/)。利用薄板樣條空間插值算法[28-29],得到1999年1月至2018年12月的1 km空間分辨率的逐月降水柵格數(shù)據(jù)。

2.2 生長季NDVI和降水觀測時間序列統(tǒng)計

計算觀測異常值。以生長季(6、7、8月)降水異常值計算為例,利用生長季每個月的降雨插值圖,計算如下變量：

(1)計算每一年生長季降水均值,即3個月(6、7、8月)的平均值(mm)；

(2)計算2001—2018年生長季降水均值(mm)；

(3)用(1)除以(2)得到異常百分比(%)

同樣地,計算生長季NDVI異常百分比。為了規(guī)避觀測數(shù)據(jù)缺失帶來的誤差,我們沒有計算6、7、8月的平均值,而只用8月的NDVI來代表生長季的NDVI。

2.3 TSS-RESTREND方法

標(biāo)準(zhǔn)RESTREND是由Evans and Geerken[9]開發(fā)并經(jīng)Wessels, Prince[4]等人進(jìn)一步修改。算法基于像素的分析,用于揭示氣候要素和人類活動導(dǎo)致的植被生態(tài)系統(tǒng)的退化過程及其規(guī)律[3, 30]。利用最小二乘回歸(Ordinary Least Squares Regression, OLSR),計算植被凈初級生產(chǎn)力(Net Primary Productivity, NPP)與降水量的關(guān)系,即植被—降水關(guān)系(Vegetation Precipitation Relationship, VPR),其中NPP用最佳生長季節(jié)NDVI(NDVImax)表示[5-6],降水量用最優(yōu)累積降水量表示。根據(jù)OLSR模型得到的預(yù)測NDVImax與每次觀測NDVImax的差值稱為VPR—殘差[9]。

yi=β0+β1xi

(1)

式中,yi是VPR—殘差,xi是年份,β0是截距,β1是斜率。

一般情況下,認(rèn)為VPR—殘差的趨勢與降水無關(guān),是地表土地變化的表征[9]。利用標(biāo)準(zhǔn)RESTREND算法,必須滿足以下3個條件[6]：

1)VPR關(guān)系必須是正向(斜率>0)且顯著的。對顯著的推薦值為：在P<0.05水平上,R2>0.3。

2)VPR—殘差在時間序列上是單調(diào)函數(shù)[31]。

3)VPR與時間高度相關(guān),即VPR在時間序列上具有可比性,例如在研究區(qū)間內(nèi)生態(tài)系統(tǒng)沒有發(fā)生重大結(jié)構(gòu)性變化[5]。

在我們的研究中,最佳生長季NDVI發(fā)生在夏季,最優(yōu)累積則是通過逐像素計算累積周期(1—12個月)和偏移期(1—3個月)的組合與NDVImax的相關(guān)系數(shù),相關(guān)系數(shù)最高的降水累積量為最優(yōu)累積降水。標(biāo)準(zhǔn)RESTREND算法有一定的局限性,當(dāng)變化的速度和方向在時間序列中發(fā)生變化時,無法識別變化趨勢[32]。

TSS-RESTREND (Time series segmentation and residual trend) 能克服上述問題,最先被Burrell, Evans[6]提出并用于澳大利亞全境的植被變化分析,該方法的主要目的是首先利用BFAST(Breaks for Additive Season and Trend)算法[33-34]檢測出NDVI時間序列的斷點。不像之前的研究應(yīng)用[13, 35],在TSS-RESTREND中的BFAST算法不直接用于NDVI時間序列,而是用于NDVI—降水關(guān)系的殘差時間序列上,在此應(yīng)用中關(guān)閉BFAST的季節(jié)成分(Seasonal Component).

TSS-RESTREND首先在研究時間段內(nèi)利用BFAST算法找出NDVI時間序列的不一致年份。在每一個子時間段內(nèi),利用Chow檢驗準(zhǔn)則判定不一致年份對植被NDVI的影響程度,當(dāng)F檢驗在α=0.05水平時,BFAST檢測的不一致年份并沒有改變植被的結(jié)構(gòu),BFAST檢測結(jié)果被駁回,對檢測出來的時間段進(jìn)行前后合并。TSS-RESREND有四種情形：

(1)當(dāng)VPR (α=0.05) 和 VPR—殘差 (α=0.05)顯著相關(guān), 符合標(biāo)準(zhǔn)RESTREND算法。

(2)當(dāng)在VPR—殘差中檢測到不一致年份,Chow檢驗同樣要再用于VPR。當(dāng)VPR—殘差(α=0.05)顯著而VPR (α=0.05)不顯著時,采用分段RESTREND算法(segmented RESTREND):

yi=β0+β1xi+β2zi+β3xizi

(2)

式中,xi為年份,zi為啞變量(0 或者1).β0為截距,β1為斜率,β2為斷點處的抵消項,β3為斷點處的斜率。

(3)當(dāng)一個像元在VPR中有顯著斷點時,可能預(yù)示著生態(tài)系統(tǒng)的結(jié)構(gòu)性改變。在這種情況下,假定降雨對斷點兩側(cè)的影響是一致的假設(shè)就不合理,因此,時間序列NDVI需要分段分析。為了使兩側(cè)的數(shù)據(jù)有可比性,可以對降雨進(jìn)行歸一化得到一個標(biāo)準(zhǔn)項:

(3)

式中,zi為標(biāo)準(zhǔn)項,xi為觀測值,μ為時間序列均值,σ為標(biāo)準(zhǔn)差。利用這個標(biāo)準(zhǔn)項,可以建立如下多變量回歸關(guān)系：

yi=β0+β1xi+β2zi+β3xizii=2001,…,2018

(4)

式中,xi為式(3)計算的降水標(biāo)準(zhǔn)項,zi為啞變量(0 或者1)。β0為截距,β1為斜率,β2為斷點處的抵消項,β3為斷點處的斜率。對一個像素而言,像元綠度的變化需要包含殘差的變化。

表1 NDVI變化分類

根據(jù)線性擬合得到的F檢驗的顯著性水平和斜率，可以將NDVI變化趨勢分為增加、下降和基本不變(No Significant Change, NSC)；其中增加細(xì)分為劇烈增加(Dramatic Increase,I1)(α<0．01)， 明顯增加(Obvious Increase, I2)(0．01≤α<0．025)， 中等增加(Moderate Increase, I3) (0．025≤α<0．05)和輕微增加(Slight Increase, INC) (0．05≤α<0．1)。下降細(xì)分為劇烈下降(Dramatic Decrease, D1)(α<0．01)， 明顯下降(Obvious Decrease, D2)(0．01≤α<0．025)， 中等下降(Moderate Decrease, D3)(0．025≤α<0．05) 和輕微下降(Slight Decrease, DNC)(0．05≤α<0．1)

(4)不符合上述條件的像素被認(rèn)為不適宜用TSS-RESTREND算法。

根據(jù)上述回歸關(guān)系中的β1和顯著性水平α[7, 36],根據(jù)NDVI殘差與時間序列的線性擬合,我們利用F檢驗的α值將NDVI的變化趨勢分為9個等級(表1)。NDVI下降的趨勢等級分為4個等級,分別為D1 (α<0.01), D2 (0.01≤α<0.025), D3 (0.025≤α<0.05) 和 DNC (0.05≤α<0.1),NDVI增加的趨勢等級分為4個等級,分別為I1 (α<0.01), I2 (0.01≤α<0.025), I3 (0.025≤α<0.05)和 INC (0.05≤α<0.1)。NSC為顯著性水平不通過(α>0.1),表示NDVI無變化。

3 結(jié)果與討論

3.1 2001—2018年生長季NDVI時空分布與統(tǒng)計分析

圖2展示了陜西省生長季最大NDVI在2001—2018年間的空間分布。從圖中可以看出,NDVI大小隨緯度變化形成南北差異,隨地形起伏變化形成區(qū)域差異,總體NDVI的空間分布特征與陜西省不同的地形地貌空間(圖1)具有較好的一致性,NDVI時間序列南高北低,關(guān)中平原包含少量低值區(qū)。

時間維度上,全省NDVI平均值的年增長率為0.006/年,如圖3所示。由于2011年研究區(qū)域生長季受云影響較大,缺失值較多,因此統(tǒng)計時去除2011年的NDVI值。該增長率超過了三北防護(hù)林(0.0007/a)[37]、黃河上游(0.0023/a)[38]和三江源地區(qū)(0.0001/a)[39],這說明退耕還林工程對陜西省的影響非常顯著。截至2018年,相比18年來的平均值(圖4),77.29%的區(qū)域大于均值。其中,陜北的榆林、延安大于均值的區(qū)域較大,分別為97.52%和89.03%(表2),秦巴山區(qū)次之,為73.91%。特別是2012年之后, NDVI高值向北逐年推進(jìn)趨勢明顯(圖2)。

圖2 2001—2018年陜西省年最大歸一化植被指數(shù)Fig.2 Maximum NDVI(Normalized Difference Vegetation Index) in Shanxi in August from 2001 to 2018

表2 NDVI與18年均值(2001—2018年)比較

圖3 2001—2018年陜西省生長季年NDVI曲線Fig.3 Variations in annual maximum NDVI and abnormal (%) in Shaanxi in August from 2001 to 2018

圖4 2001—2018年陜西省生長季NDVI異常百分比Fig.4 Maximum NDVI abnormal (%) in Shaanxi in August from 2001 to 2018

3.2 人類活動造成的NDVI變化趨勢與分級

利用TSS-RESTREND算法計算得到NDVI變化量如圖5所示。2001—2018年間,陜西省植被綠度受劇烈的人類活動影響,NDVI呈增加態(tài)勢,NDVI增加的區(qū)域達(dá)71.77%,而陜北地區(qū)的增加量明顯大于關(guān)中平原區(qū)和陜南秦巴山地,其中陜北的榆林NDVI增加區(qū)域為72.11%,延安為86.44%,均超過了全省平均水平,這是由于陜西省北部NDVI增長與人工造林面積密切相關(guān),有關(guān)研究報道達(dá)96%以上,說明陜西省北部實施的退耕還林生態(tài)工程效果顯著,據(jù)統(tǒng)計,2000—2014年間,陜北地區(qū)山地和沙地造林面積達(dá)238.6×102km2[40]。

NDVI減少區(qū)域主要分布在城市周邊,如西安市、咸陽市、漢中市、安康市。如圖所示。關(guān)中平原NDVI減少量區(qū)域面積最大,占該區(qū)域的6.04%,秦巴山地次之(1.20%)。

圖6 2001—2018年陜西省綠度分級Fig.6 Greenness level determined by TSS-RESTREND in Shaanxi from 2001 to 2018 根據(jù)線性擬合得到的F 檢驗的顯著性水平和斜率，可以將NDVI 變化趨勢分為增加、下降和基本不變(No Significant Change, NSC)。其中增加細(xì)分為劇烈增加(Dramatic Increase,I1)(α<0．01)， 明顯增加(Obvious Increase, I2)(0．01≤α<0．025)， 中等增加(Moderate Increase, I3) (0．025≤α<0．05)和輕微增加(Slight Increase, INC) (0．05≤α<0．1)。下降細(xì)分為劇烈下降(Dramatic Decrease, D1)(α<0．01)， 明顯下降(Obvious Decrease, D2)(0．01≤α<0．025)， 中等下降(Moderate Decrease, D3)(0．025≤α<0．05) 和輕微下降(Slight Decrease, DNC)(0．05≤α<0．1)

(3)按照表1所示的綠度分級規(guī)則,可以將全省人類活動引起的NDVI變化量定性分為增長(INC,I1,I2,I3)、減少(D1,D2,D3)和無明顯變化(NSC)8類,總體上陜西全省呈變綠趨勢(圖6)。榆林市和延安市的變綠區(qū)域明顯多于關(guān)中平原和秦巴山地(圖6,圖7),在空間上呈連續(xù)分布,而關(guān)中平原和秦巴山地的變綠區(qū)域大部分呈點狀分布。表3所示的TSS-RESTREND檢測出的劇烈增長區(qū)域,延安為55.46%(11935.9km2),榆林為34.34% (13825.7 km2),而陜南為41.03% (14746.8 km2),這說明處于濕潤氣候區(qū)的陜南地區(qū)也有顯著變綠趨勢,與已有研究結(jié)果一致[18, 41]。在這里需要說明是,當(dāng)NDVI和降雨量關(guān)系不符合TSS-RESTREND算法適用條件時,該像素不做任何處理。因此,表3的面積統(tǒng)計未包含算法不適用像元。雖然陜西省位于東部濕潤和西部干旱過渡的半干旱地區(qū),降雨要素是影響植被生長變化的首要氣候要素,但是局部區(qū)域,如陜南秦巴山地,氣溫和水熱因素也是抑制植被生長的重要要素。雖然TSS-RESTREND算法可以在算法層面剔除算法不適用的像元,但是人類活動對陜西南部植被的影響結(jié)果還須進(jìn)一步論證。TSS-RESTREND算法內(nèi)置的BFAST時間序列分析方法可以檢測出VPR關(guān)系變化年份,從而可以擴(kuò)大RESTREND算法的應(yīng)用區(qū)域,并且可以反映NDVI變化帶來的生態(tài)累積效應(yīng)。圖7展示了自2001—2018年,TSS-RESTREND檢測的VPR斷點年份。整體上,斷點年份發(fā)生在2008年之后,這說明自2008年后,原有的植被NPP與降水的關(guān)系被打破,常規(guī)的RESTREND等趨勢分析算法將不適用。斷點年份發(fā)生在2008年之前的區(qū)域集中在關(guān)中平原西安市周邊以及秦巴山區(qū)的安康市南部。需要說明的是,本文研究區(qū)間雖然是在2001—2018年,但是檢測的斷點年份卻在2004—2016年,當(dāng)綠度變化發(fā)生在研究時間段的前3年或者后3年期間時,TSS-RESTREND算法無法檢測出斷點年份,這是由于當(dāng)生態(tài)系統(tǒng)結(jié)構(gòu)狀態(tài)的改變需要幾年的適應(yīng)周期,反映在NDVI上,顯示出一定的滯后性[36]。

表3 綠度分級統(tǒng)計

圖7 斷點年份檢測Fig.7 Break year detection in Shaanxi from 2001 to 2018

4 結(jié)論

本文分析了自退耕還林生態(tài)工程實施以來的2001—2018年,陜西省NDVI的時空變化規(guī)律,利用TSS-RESTREND算法剝離了氣候要素(降雨)對植被NDVI的影響,得到了人類活動對植被變化的影響程度和區(qū)域,主要結(jié)論如下：

(1) 2001—2018年間,陜西省NDVI呈顯著增加,全省平均增加速率為0.006/a;

(2) 相比18年來的平均值,77.29%的區(qū)域大于均值。其中,陜北的榆林、延安大于均值的區(qū)域較大,分別為97.52%和89.03%,秦巴山區(qū)次之,為73.91%。2012年之后,NDVI高值向北逐年推進(jìn)趨勢明顯。

(3)在只考慮降雨對植被的影響前提下,人類活動對植被的變化影響巨大,全省NDVI增加的區(qū)域達(dá)71.77%,而陜北地區(qū)的增加量明顯大于關(guān)中平原區(qū)和陜南秦巴山地,其中陜北的榆林NDVI增加區(qū)域為72.11%,延安為86.44%,均超過了全省平均水平。

(4)總體上陜西全省呈變綠趨勢。榆林市和延安市的變綠區(qū)域明顯多于關(guān)中平原和秦巴山地,延安和榆林的劇烈增長區(qū)域分別為55.46%和34.34%,而陜南為41.03%,說明處于濕潤氣候區(qū)的陜南地區(qū)也有顯著變綠趨勢。