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

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

1 貴州師范大學地理與環境科學學院, 貴陽 550025 2 貴州省山地資源與環境遙感應用重點實驗室, 貴陽 550001 3 中國科學院生態環境研究中心城市與區域生態國家重點實驗室, 北京 100085

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

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

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

1 研究區域

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

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

2 研究方法

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

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

2.1 數據來源與處理方法

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

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

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

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

(1)

(2)

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

(3)

(4)

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

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

2.3 趨勢分析及分級

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

(5)

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

表1 NDVI變化趨勢等級

3 結果

3.1 NDVI突變分析及空間變化

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

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

3.2 NDVI變化趨勢及其與氣象因子相關性分析

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

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

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

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

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

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

3.3 典型退化區時序突變與傳統檢測比較

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

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

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

表3 典型退化區面積比例

4 討論

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

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

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

5 結論

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

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

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

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