基于時序突變檢測的植被空間變化特征識別方法
——以海河北部山區(qū)為例

海 月,楊廣斌,李若男,鄭 華

1 貴州師范大學(xué)地理與環(huán)境科學(xué)學(xué)院, 貴陽 550025 2 貴州省山地資源與環(huán)境遙感應(yīng)用重點實驗室, 貴陽 550001 3 中國科學(xué)院生態(tài)環(huán)境研究中心城市與區(qū)域生態(tài)國家重點實驗室, 北京 100085

植被是覆蓋地表植物群落的總稱,是生態(tài)系統(tǒng)的重要組成部分。植被既受大氣、土壤、水分等環(huán)境要素的綜合影響,又能反作用于環(huán)境[1-2]。研究植被長期變化不僅能夠揭示生態(tài)系的演變過程并指示生態(tài)系統(tǒng)健康狀況；更有利于人類對植被進行合理的保護、利用及管理[3-4]。隨著遙感技術(shù)的發(fā)展、數(shù)據(jù)獲取的便捷及分析手段的提升,基于植被長時序趨勢特征的研究越來越多,趨勢分析成為常用的分析方法。歸一化植被指數(shù)(Normalized Vegetation Index,NDVI)能夠反映植被的生長狀態(tài)、覆蓋情況、生物量及環(huán)境要素影響[5-8],并作為植被質(zhì)量評估指標應(yīng)用十分廣泛。長時序NDVI分析的優(yōu)勢不僅在于能有效的反應(yīng)植被變化的情況,還在于不同數(shù)據(jù)源能夠混合使用,使得大空間尺度的植被質(zhì)量評估成為可能[9-11]。

在全球氣候變化的大背景下,利用長時序NDVI的空間特征評估氣候變化帶來的影響是關(guān)注的重點。許多研究已從不同空間尺度上(全球、區(qū)域、流域等)研究了氣候變化對植被及土地利用變化的影響[12-17]。從研究結(jié)果看,NDVI空間特征與氣候變化密切相關(guān),氣候的時空異質(zhì)性對NDVI的變化趨勢有十分顯著的影響[18-23]。此外,長時序NDVI分析還被廣泛用于分析土地利用演變,包括生態(tài)恢復(fù)工程帶來的植被變化[24-27],以及城市擴張帶來的植被退化等[28-33]。從時間尺度上看,除年際間長時序分析外[34],部分研究還分析了NDVI在季節(jié)和年內(nèi)的變化趨勢[35],并在小空間尺度上,結(jié)合不同物種的物候信息揭示氣候和人類活動對植被生長變化的影響[36-37]。

現(xiàn)有的大多數(shù)研究中,對NDVI的趨勢分析多針對時間和空間總體變化展開[34,38],無法描述精細變化特征,不能反映長時序數(shù)據(jù)中的波動特征,從而無法滿足植被保護及恢復(fù)差異化決策和管理需求,使得長時序NDVI分析結(jié)果受到質(zhì)疑[34, 39-40]?；诖?本研究建立了一種新的空間特征識別方法,利用2000—2018年MODIS NDVI產(chǎn)品數(shù)據(jù)集,將時序突變檢測與趨勢分析相結(jié)合,基于柵格尺度識別NDVI空間突變拐點,分析了突變拐點前后NDVI變化趨勢與常規(guī)全時序變化趨勢的異同,分析了兩種趨勢分析的空間特征差異,并探討了趨勢變化的影響因素。

1 研究區(qū)域

海河流域北部(東經(jīng)115°30′至119°45′,北緯39°10′至42°40′)由灤河及潮白河流域組成,北接內(nèi)蒙古壩上高原,南連京津冀城市群。區(qū)域內(nèi)年均溫1℃至11℃,年降水量400—700 mm,屬暖溫帶半干旱-半濕潤季風氣候。灤河及潮白河流域面積89870 km2,其中上游山區(qū)面積占73.9%(66376 km2),跨河北、內(nèi)蒙古、北京、天津、遼寧五省(市)。區(qū)域內(nèi)西北部以草地覆蓋為主,中部及東部以林灌覆蓋為主,其中林灌草的覆蓋面積占78.1%(圖1)。研究區(qū)域內(nèi)分布著眾多針對森林、草地、地質(zhì)、生物等國家、省、縣級保護區(qū),總面積4093 km2,占區(qū)域面積的6.2%。從行政區(qū)域來看,海河流域北部山區(qū)主要涵蓋承德市大部分區(qū)域,以及北京、唐山北部,張家口西部,秦皇島東部,是京津冀都市圈重要的生態(tài)涵養(yǎng)區(qū)；此外,潮白河及灤河還是首都水源地的源頭和首都應(yīng)急水源地所在地。區(qū)域內(nèi)的植被具有十分重要的生態(tài)屏障作用,其空間變化特征不僅關(guān)系到本地生態(tài)系統(tǒng)健康,同時還會對下游都市圈的生態(tài)安全產(chǎn)生影響。

圖1 研究區(qū)域土地利用及氣象站點分布示意圖Fig.1 Distribution of land use and meteorological stations

2 研究方法

為了揭示植被在不同時段內(nèi)的空變化特征,從而精確評估變化趨勢,為明確區(qū)域人類活動對植被的影響以及差異化的植被恢復(fù)策略提供支持。本研究構(gòu)建了基于柵格尺度的NDVI時序突變檢測方法,采用Pettit突變分析方法識別NDVI空間突變的拐點,采用最小二乘擬合計算突變拐點前后NDVI的變化趨勢,分析各時段趨勢變化空間特征,并與常規(guī)全時序NDVI趨勢分析及時空變化特征對比,提出NDVI空間變化特征差異,并分析差異的影響因素(圖2)。

圖2 基于時序突變檢測的NDVI空間特征識別框架Fig.2 Framework of NDVI spatial feature recognition based on temporal mutation detection

2.1 數(shù)據(jù)來源與處理方法

NDVI數(shù)據(jù)來自Google earth engine中MODIS13Q1數(shù)據(jù)集,空間分辨率為250 m,時間分辨率為16 d。本研究中使用數(shù)據(jù)的時間段為2000年至2018年。在提取數(shù)據(jù)集中的各期NDVI后,使用Savitzky-Golay濾波方法平滑數(shù)據(jù)并排除異常值,識別研究區(qū)域多年NDVI季節(jié)變化特征,采用每年4—9月生長季節(jié)均值作為NDVI的表征。

氣象數(shù)據(jù)來自資源環(huán)境數(shù)據(jù)云平臺,采用2000—2018年研究區(qū)域及周邊59個國家級氣象站的年均值,利用ANUSPLIN插值軟件對年均氣溫和降水進行插值,空間分辨率為250 m,進而提取海河北部山區(qū)氣象因子的空間分布。

2.2 突變分析及突變拐點識別

突變分析用于識別基于柵格尺度的NDVI在2000—2018年間的突變時間,采用Pettit突變分析方法。該方法不僅能夠識別長時間序列中突變點的位置及數(shù)量,還能夠判斷突變點的統(tǒng)計意義是否顯著,方法如下：

(1)

(2)

K=max|Ut,n|,t=1,2,…,n

(3)

(4)

式中,Ut,n為定義的統(tǒng)計量；n為時間序列長度；xj,xi為檢驗樣本中的隨機變量；K為突變年份；P為顯著性水平檢驗值；在給定顯著性水平α情況下,若P≤α,則認為該突變點滿足顯著性水平。

確定各柵格的突變年份后,統(tǒng)計各年突變柵格的數(shù)量,得到累積突變柵格比例,將累計曲線拐點對應(yīng)年份作為時序分析的分界點。

2.3 趨勢分析及分級

NDVI趨勢分析采用最小二乘擬合計算植被變化趨勢的斜率,計算方法如下：

(5)

式中,SNDVI為變化趨勢,SNDVI>0表明植被變化呈現(xiàn)增加趨勢,SNDVI<0表明植被變化呈現(xiàn)降低趨勢,SNDVI=0表明植被變化無明顯變化趨勢；n為時間序列長度；i為研究年份序號；NDVIi為第i年的生長季平均NDVI值。從顯著性水平來看,選取0.05顯著性水平將NDVI變化分為3個等級(表1)。

表1 NDVI變化趨勢等級

3 結(jié)果

3.1 NDVI突變分析及空間變化

圖3 年度突變柵格及累計突變比例Fig.3 Annual mutation grids and cumulative percentage

基于柵格尺度的突變空間分析結(jié)果表明,顯著突變(P<0.05)年份的空間分布呈現(xiàn)出較為明顯的空間異質(zhì)性：初期(2004—2008年)突變區(qū)域較小,主要集中在區(qū)域東南部；中期(2009—2011年)突變范圍急劇擴大,由東南向西北覆蓋絕大部分區(qū)域；后期(2012—2014年)突變區(qū)域減少且向西北部集中。從面積比例上看,NDVI存在顯著突變出現(xiàn)的時間為2004—2014年,其中突變集中的時間為2009—2011年(圖3)。從突變累積曲線看,2011年顯著突變的柵格出現(xiàn)拐點。至2011年,研究區(qū)域內(nèi)累積共有69.1%的面積發(fā)生了突變,占突變面積的88.7%,基于以上結(jié)果,本研究以2011年作為分界點,分析2000—2011年和2011—2018年NDVI的變化趨勢,并將其與2000—2018年的全時序趨勢分析結(jié)果進行對比。

3.2 NDVI變化趨勢及其與氣象因子相關(guān)性分析

圖4為海河北部山區(qū)不同時間段NDVI的變化趨勢。對比全時序和突變拐點前后的結(jié)果可以看出,全時序變化趨勢中,區(qū)域整體呈現(xiàn)改善趨勢(87.2%),退化區(qū)域集中中部以及東南部(1.7%)。突變拐點前的趨勢變化中,呈現(xiàn)改善的區(qū)域分布在中部和東部(占54.3%),退化區(qū)(占1.6%,)主要分布在西北部草原區(qū)(正藍旗、克什克騰旗、多倫縣),中部承德市周邊(雙灤區(qū)、雙橋區(qū)),東南部唐山北部(遵化市、遷西縣、遷安市、灤平縣)。突變拐點后,改善區(qū)域主要分布在北部和西部(占52.7%),退化區(qū)(占1.2%)主要分布在東部寬城縣、青龍縣、以及遷西縣和遷安市北部。

圖4 海河北部山區(qū)NDVI變化趨勢Fig.4 Classification of NDVI variation trend in the mountain area of Northern Haihe Basin

對比不同時間段NDVI變化趨勢的重疊率可以看出(表2),突變拐點前后顯著改善和基本不變的區(qū)域與全時序分析的重疊率較高(均大于50%),顯著退化區(qū)的重疊率較低,特別是2011—2018年間,顯著退化區(qū)域與全時序分析的重疊區(qū)僅為15.8%。這表明基于時序突變識別的變化趨勢與全時序分析的識別結(jié)果存在空間差異,特別是顯著退化區(qū)的識別結(jié)果差異最大,而退化區(qū)正是植被恢復(fù)并需要人類正向干預(yù)的重點區(qū)域。

表2 突變前后NDVI變化趨勢占全時序分析比例

圖5 海河北部山區(qū)NDVI與氣象因子相關(guān)性Fig.5 Correlation of NDVI to meteorological factors in the mountain area of Northern Haihe Basin

對比全時序和時序突變拐點前后時間段,NDVI與氣象因子相關(guān)性也存在明顯的空間異質(zhì)性(圖5)。氣溫與NDVI總體上呈現(xiàn)負相關(guān),其中突變拐點前呈現(xiàn)負相關(guān)的面積最大,占77.9%,突變后負相關(guān)面積降至50.6%,集中在西部和東部；相比之下,全時序分析的負相關(guān)面積比例為72.3%,顯著負相關(guān)占比(11.7%)大于突變拐點前后(1.7%和3.4%)。三種時段中,降水與NDVI的正相關(guān)面積比例均在75%以上；突變后,呈現(xiàn)負相關(guān)的面積有所增加,占總面積的24.46%,主要分布在西部赤城、豐寧縣和東南部遵化、青龍、寬城等地。結(jié)合NDVI變化趨勢來看,不同時段內(nèi),氣象因子對NDVI顯著改善的作用并不大,但卻是突變后中東部局部呈現(xiàn)退化的原因之一。

3.3 典型退化區(qū)時序突變與傳統(tǒng)檢測比較

從識別結(jié)果來看,退化主要集中在三個區(qū)域(圖6),分別為西北部草原區(qū)(包括多倫縣、圍場縣、以及正藍旗和克什克騰旗部分區(qū)域),中部林灌區(qū)(包括隆化縣、灤平縣、承德縣、雙灤區(qū)和雙橋區(qū)),和東南林農(nóng)區(qū)(包括寬城縣、興隆縣、遵化市、遷西縣、遷安市、以及灤縣部分區(qū)域)。

圖6 典型退化區(qū)空間分布Fig.6 Spatial distribution of typical degradation area

從面積上看(表3),三個典型退化區(qū)的面積占全時序分析結(jié)果的67.0%,其中中部林灌區(qū)的占比最大,為43.1%,其次為東南林農(nóng)區(qū)(15.1%)和西北草原區(qū)(8.8%)。對比兩種趨勢分析結(jié)果,退化區(qū)的分布存在明顯的空間異質(zhì)性。突變拐點前后顯著退化區(qū)的比例與全時序突變趨勢的重疊比例均低于50%,特別是2011—2018年,平均重疊比例僅為16.3%。結(jié)合區(qū)域氣候變化及經(jīng)濟發(fā)展政策來看,氣候變化帶來的退化影響在2000—2011年主要影響西北部草原區(qū),在2011—2018年主要影響東南林農(nóng)區(qū)；城鎮(zhèn)擴張導(dǎo)致的退化主要是在2011—2018年間影響中部林灌區(qū)和東南部林農(nóng)區(qū)。

表3 典型退化區(qū)面積比例

4 討論

從面積上看,基于時序突變分析與常規(guī)全時序分析在海河北部山區(qū)NDVI趨勢空間特征上存在明顯差異。常規(guī)全時序趨勢分析方法中,顯著改善的區(qū)域占80%以上,退化區(qū)域占1.7%；基于時序突變分析方法中,突變拐點前后時段內(nèi),改善區(qū)面積均下降了30%左右,基本不變的區(qū)域增加了45%左右,退化區(qū)占1.2%—1.5%左右。從空間分布上看,2000—2018年區(qū)域NDVI整體呈現(xiàn)顯著上升的趨勢,植被處于改善狀態(tài),研究結(jié)果與海河及我國北方2000年以后NDVI的變化趨勢一致。何航等[38]對中國北方NDVI變化特征的研究中表明,在氣候變化背景下,1982—2015年間華北地區(qū)植被程良性改善狀態(tài),顯著增加面積大于40%。王永財?shù)萚41]針對海河流域1998—2011年不同土地覆被類型的NDVI變化趨勢進行了分析,顯示各類土地覆被的NDVI均呈現(xiàn)上升趨勢,其中森林、農(nóng)田、濕地、草地的上升趨勢明顯。基于時序突變趨勢分析結(jié)果表明,雖然改善區(qū)仍占主導(dǎo)(>50%),但面積顯著減少,突變拐點前后空間格局差異明顯。突變前,改善區(qū)主要分布在中部和南部,突變拐點后主要分布在北部和西部；退化區(qū)則由突變前的西北和東南部變?yōu)橥蛔兒蟮臇|部地區(qū)。對比前人的研究成果,基于時序突變的分析方法準確的識別了突變拐點前錫林郭勒東部草原草甸區(qū)的退化[42],以及突變拐點后唐山北部林農(nóng)區(qū)的退化[43]。

NDVI的變化受到多種因素的共同影響,主要包括氣候變化和人類活動兩方面。在全球暖濕化的總體趨勢下,氣象因子對NDVI的影響一直是研究的熱點。結(jié)合區(qū)域氣象因子變化,海河北部山區(qū)NDVI變化趨勢與降水的時空變化密切相關(guān)。突變拐點前,降水下降區(qū)分布在西北部草原區(qū),突變拐點后則轉(zhuǎn)移至東部林農(nóng)區(qū)。從NDVI與氣溫的相關(guān)性變化趨勢看,呈現(xiàn)負相關(guān)的區(qū)域主要分布在承德周邊的林灌覆蓋區(qū),表明植被覆蓋對局地氣溫起到了一定調(diào)節(jié)作用。相比突變前,突變后NDVI與氣溫呈負相關(guān)的區(qū)域減少(由77.9%變?yōu)?0.6%),但呈現(xiàn)顯著負相關(guān)的區(qū)域增加(由1.7%變?yōu)?.4%),表明雖然植被對氣候調(diào)節(jié)作用的面積有所減少,但強度依然明顯。從相關(guān)性分析來看,氣象因子對NDVI的變化的空間異質(zhì)性十分明顯,呈現(xiàn)促進作用的區(qū)域集中在承德西部、北京北部和張家口赤城縣,主要源于該區(qū)域在2011年后降水及氣溫呈現(xiàn)整體增加趨勢,有利于植被的生長。呈現(xiàn)負面作用的區(qū)域集中在東部(遷西、遵化、遷安、青龍和寬城),其原因可能主要由于東部在2011年后降水呈現(xiàn)下降趨勢,其中遷西、遷安和青龍等地呈現(xiàn)出顯著下降趨勢。從不同植被類型的面積來看,森林仍然是受影響最大的區(qū)域,其次為灌叢和耕地；從植被類型的比例看,55%以上的森林和灌叢、47.0%的耕地及27.7%的草地受到了氣象因子的正向影響,僅有4.5%以上的森林灌叢、9.9%的耕地和10.7%的草地受到負向影響。相比氣象因子,人類活動是海河北部山區(qū)NDVI整體呈上升趨勢的主要因素,這主要是源于近幾十年來在我國北方開展的生態(tài)造林及森林保育工程。從研究區(qū)域整體來看,基于京津風沙源和退耕還林還草等工程,2000年至2015年間,區(qū)域內(nèi)林灌由41868 km2增加至42086 km2,草地由9598 km2增加至9,752km2。統(tǒng)計表明,承德市自2002年京津風沙源治理工程建設(shè)以來,共造林61.3萬hm2,截至2015年,承德市的森林面積占京津冀三省市林地面積的36%[44]。張家口市轄赤城縣在2001—2013年間通過飛播或人工種植共造林4.0萬hm2[45]。隨著工程后期植被恢復(fù)增速的降低以及區(qū)域發(fā)展等因素,人類活動對NDVI帶來的負面影響也逐步體現(xiàn),特別是2011年“首都經(jīng)濟圈”概念提出后,唐山、承德等城市周邊的城市化進程也直接導(dǎo)致了2011—2018年中部承德周邊和東南部唐山-秦皇島周邊NDVI的退化。綜合來看,海河北部山區(qū)的植被整體改善主要是由于人類活動的影響,而突變點前后局部植被退化的空間異質(zhì)性則是氣象因素和人類活動共同作用的結(jié)果。

從研究結(jié)果看,常規(guī)全時序分析能夠識別總體趨勢變化,但卻無法識別由于數(shù)據(jù)區(qū)間波動帶來局部時空特征,無法準確評估生態(tài)系統(tǒng)的退化狀況和識別退化區(qū)域,從而不能準確判別退化的影響因素,無法為生態(tài)工程效果評估和后續(xù)政策的調(diào)整提供科學(xué)依據(jù)?；跁r序突變分析方法在識別NDVI分時段趨勢的基礎(chǔ)上,能夠揭示植被空間變化的動態(tài)趨勢,為明確氣象因素和人類活動對植被的影響(特別是負面影響),并為制定差異化修復(fù)策略提供了有效的決策信息。

5 結(jié)論

(1)結(jié)合歸一化植被指數(shù)NDVI,建立了基于時序突變檢測的植被空間變化特征識別方法,并與常規(guī)長時序趨勢檢測方法及結(jié)果對比；結(jié)果表明,基于時序突變檢測的方法能夠識別NDVI精細的空間變化特征,從而為后續(xù)修復(fù)提供完整有效的信息。

(2)以海河北部山區(qū)為例,通過時序突變分析表明海河北部山區(qū)NDVI在2011年存在突變拐點；突變拐點前后NDVI變化趨勢均以改善為主,但退化區(qū)域呈現(xiàn)出空異質(zhì)性,突變拐點前的退化區(qū)主要位于西北草原區(qū),突變拐點后的退化區(qū)主要位于中東部林農(nóng)區(qū)。

(3)分析表明,區(qū)域內(nèi)植被整體改善主要源于人類活動的影響(生態(tài)工程等),氣候因素的作用有限；但區(qū)域內(nèi)突變拐點前后植被局部退化的空間變化特征則受到氣候因素和人類活動共同影響。

(4)基于時序突變分析方法與常規(guī)長時序分析方法的評估結(jié)果存在空間差異,特別是退化區(qū)域的識別結(jié)果差異顯著；對比前人研究成果,基于時序突變分析結(jié)果能夠更為準確的識別區(qū)域生態(tài)系統(tǒng)的變化及影響因素。