基于斷點(diǎn)檢測的延安市富縣植被變化研究

季燁云，劉鵬舉，高 影，任 怡，余 濤，牛思圓

（1. 中國林業(yè)科學(xué)研究院 資源信息研究所，北京 100091；2. 中國林業(yè)科學(xué)研究院 國家林業(yè)和草原局林業(yè)遙感與信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100091；3. 國家林業(yè)和草原局 調(diào)查規(guī)劃設(shè)計院，北京 100029）

植被是生態(tài)系統(tǒng)的重要組成部分，對全球氣候變化和生態(tài)系統(tǒng)的結(jié)構(gòu)功能都有顯著影響[1-3]。研究地表植被覆蓋及其長期變化不僅可以揭示地表空間變化規(guī)律，而且有利于分析區(qū)域生態(tài)系統(tǒng)及環(huán)境變化情況[4]，對因地制宜開展植被的保護(hù)管理工作具有重要意義。黃土高原位于黃河中游地區(qū)，是世界上水土流失最為嚴(yán)重的區(qū)域之一。陜北是黃土高原的中心部分，水土流失嚴(yán)重、生態(tài)環(huán)境脆弱，但自三北防護(hù)林、天然林資源保護(hù)、退耕還林還草等工程實(shí)施以來，該區(qū)域生態(tài)環(huán)境得到了明顯改善。因此有必要開展工程實(shí)施效果階段性評估，便于為今后該區(qū)域的林業(yè)生態(tài)建設(shè)和保護(hù)提供理論依據(jù)。

時間序列的遙感數(shù)據(jù)可以動態(tài)監(jiān)測大面積的植被生長狀況[5-6]、發(fā)展變化趨勢特征[7-8]、擾動情況[9-10]等?，F(xiàn)有的大多數(shù)研究中，對歸一化植被指數(shù)(normalized difference vegetation index，NDVI)時間序列的分析多針對空間和時間總體變化展開[11]，如STOW等[12]提出用一元線性回歸方法來分析每個柵格像元多年的變化趨勢，以此來計算植被的綠度變化率(greenness rate of change，GRC)。當(dāng)植被遭遇砍伐、病蟲害、火災(zāi)等自然災(zāi)害或人類活動影響時，植被綠度變化趨勢可能會出現(xiàn)突然的轉(zhuǎn)變。然而，一元線性回歸方法只能表現(xiàn)植被變化的整體趨勢，不能監(jiān)測到以上趨勢發(fā)生突然變化的情況，具有很大的局限性[13-14]。因此，學(xué)者們提出了專用于處理遙感時序數(shù)據(jù)的斷點(diǎn)檢測方法，如LandTrendr(Landsat based detection of Trends in Disturbance and Recovery)方法[15]、BFAST(breaks for additive seasonal and trend)方法[16]等。由于上述方法受限于衛(wèi)星傳感器、數(shù)據(jù)類型和數(shù)據(jù)序列長度等方面的不足，JAMALI等[17]提出了DBEST(detecting breakpoints and estimating segments in trend)算法，該算法不僅對衛(wèi)星傳感器類型和數(shù)據(jù)序列長度沒有限制，而且支持周期性和非周期性植被指數(shù)時間序列，檢測效果較好且具有普遍的適用性[18-20]。

本研究以延安市富縣為例，利用2000—2019年Landsat影像構(gòu)建NDVI時間序列，運(yùn)用DBEST算法識別NDVI時間序列的斷點(diǎn)，分析斷點(diǎn)的時空分布，并研究植被變化多階段趨勢特征，為生態(tài)修復(fù)工程規(guī)劃和長期生態(tài)環(huán)境保護(hù)決策提供依據(jù)。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

富縣位于陜西省北部，延安市南部，35°44′06″～36°23′23″N，108°29′30″～109°42′54″E，總面積為4 182 km2，屬于渭北黃土高原丘陵溝壑地帶。該縣屬大陸性暖溫帶季風(fēng)氣候，四季分明，年平均氣溫為7～9 ℃，年平均降水量為500～600 mm。區(qū)域內(nèi)分布著豐富的森林資源，喬木樹種主要有櫟類Quercus spp.、樺木Betula spp.、刺槐Robinia pseudoacacia等；灌木主要有胡頹子Elaeagnus pungens、馬桑Coriaria nepalensis、白刺花Sophora davidii等。居民地、河流、鐵路和公路分布概況如圖1所示。

圖1 2019 年富縣基礎(chǔ)地理信息示意圖Figure 1 Basic geographic information of Fu County in 2019

1.2 數(shù)據(jù)來源及預(yù)處理

Landsat TM/ETM+/OLI多光譜遙感影像來源于美國地質(zhì)調(diào)查局提供的地表反射率(surface reflectance，SR)產(chǎn)品，即L2級別產(chǎn)品(https://earthexplorer.usgs.gov/)，空間分辨率為30 m，時間分辨率為16 d，光譜范圍包括可見光、近紅外、短波紅外和熱紅外波段[21]。由于2003年之后的Landsat 7 ETM+影像出現(xiàn)數(shù)據(jù)條帶丟失的問題，因此需先利用ENVI軟件中的Landsat_gapfill插件去除此類影像中的條帶[22]。然后利用Fmask云檢測算法對所有影像進(jìn)行云、云影識別提取[23]。最后利用紅色和近紅外波段計算NDVI，并基于最大值合成法[24](maximum value composite，MVC)合成了2000—2019年的NDVI時間序列數(shù)據(jù)。

本研究還收集了研究區(qū)2019年的林地類型，數(shù)據(jù)來源于國家森林資源智慧管理平臺中的“森林資源一張圖”數(shù)據(jù)，用于分析不同林地的植被變化情況。各類型面積如圖2所示：有林地2 664.77 km2，疏林地 82.80 km2，灌木林地 554.95 km2，未成林造林地 253.85 km2，其他用地 625.63 km2，包括人造地表、水體、耕地等。

圖2 2019 年富縣林地類型與分布示意圖Figure 2 Forest land types and distribution of Fu County in 2019

2 研究方法

2.1 斷點(diǎn)檢測方法

DBEST算法主要是由時間序列分解和趨勢分割兩部分組成[17]。時間序列分解主要是指植被指數(shù)時間序列數(shù)據(jù)的預(yù)處理過程，而趨勢分割是實(shí)現(xiàn)斷點(diǎn)檢測的過程。由于本研究對生長季NDVI數(shù)據(jù)采用最大值生成年度NDVI數(shù)據(jù)，因此時間序列分解部分未進(jìn)行以局部移動平均為基礎(chǔ)的季節(jié)趨勢分解(seasonal-trend decomposition procedure based on Loess，STL)，僅利用水平位移第 1 閾值、水平位移第2閾值、持續(xù)時長共3個參數(shù)檢測候選斷點(diǎn)。趨勢分割(即斷點(diǎn)檢測)包括3個過程：①找到時間序列的轉(zhuǎn)點(diǎn)(即峰/谷點(diǎn)、距離閾值點(diǎn)和水平位移點(diǎn))；②通過下一個轉(zhuǎn)點(diǎn)減去當(dāng)前轉(zhuǎn)點(diǎn)處的NDVI值來計算局部趨勢變化(trend local change，TLC)；③所有轉(zhuǎn)點(diǎn)按TLC降序重新排列，利用貝葉斯信息準(zhǔn)則(bayesian information criterion，BIC)確定斷點(diǎn)數(shù)量，若BIC準(zhǔn)則得到的斷點(diǎn)數(shù)量少于轉(zhuǎn)點(diǎn)數(shù)量，則最終的斷點(diǎn)數(shù)為BIC準(zhǔn)則得到的結(jié)果。

運(yùn)用DBEST算法檢測斷點(diǎn)，設(shè)置參數(shù)如下：水平位移第1閾值(為0.1，即當(dāng)時序NDVI數(shù)據(jù)變化超過0.1時認(rèn)為該點(diǎn)是候選水平位移點(diǎn)；水平位移第2閾值)為0.2，即當(dāng)候選水平位移點(diǎn)前后最小步長范圍內(nèi)的子序列均值變化超過0.2認(rèn)為候選水平位移點(diǎn)符合第2閾值；持續(xù)時長()為2 a，在該時段內(nèi)只能有1個水平位移點(diǎn)，以檢測NDVI短時間內(nèi)的平均變化；以上3個參數(shù)可確定時間序列中的水平位移點(diǎn)。變化級別()、泛化百分比()、斷點(diǎn)數(shù)目(m)共3個參數(shù)可根據(jù)需求三選一達(dá)到不同的目的。由于富縣植被覆蓋率較高，NDVI位于-0.2～0.2之間的變化在原始影像上的目視解譯不太明顯，因此本研究將設(shè)為0.2，可保留大于0.2的斷點(diǎn)，以便于檢測較為顯著的變化。算法輸出變化開始時間(即斷點(diǎn)位置)、變化結(jié)束時間、變化幅度和變化持續(xù)時長。

2.2 NDVI趨勢分析

2.2.1 總體趨勢 本研究利用一元線性回歸方法來擬合每個像元年內(nèi)NDVI最大值(INDV(max))的年際變化趨勢[24]：

其中：i為年序號；n為總年數(shù)；INDV(max)為年內(nèi)NDVI最大值；為第i年的INDV(max)；為NDVI年內(nèi)最大值的年際變化趨勢。將基于一元線性回歸方法得到的富縣2000—2019年NDVI時間序列的趨勢定義為總體趨勢。依據(jù)宋怡等[25]提出的分類標(biāo)準(zhǔn)將趨勢劃分為3種類型：退化(≤-0.001 0)、基本不變 (- 0.001 0＜SNDVI≤0.001 0)、改善 (SNDVI＞0.0010)。

2.2.2 多階段趨勢 ① 對 NDVI時間序列進(jìn)行斷點(diǎn)檢測，分析斷點(diǎn)及 NDVI變化幅度的時空分布。② 針對未檢測到斷點(diǎn)或NDVI變化幅度小于0.2的區(qū)域，研究其總體趨勢；針對斷點(diǎn)NDVI變化幅度大于0.2的區(qū)域，對其進(jìn)行多階段趨勢分析。多階段趨勢分析過程如下：以NDVI時間序列的所有斷點(diǎn)及對應(yīng)的變化結(jié)束時間對其進(jìn)行分段，分別進(jìn)行擬合，從而得到趨勢序列。斷點(diǎn)及對應(yīng)的變化結(jié)束時間生成的擬合段定義為變化段，除了變化段以外的擬合段定義為趨勢段。當(dāng)前階段趨勢定義為最末一段趨勢段。趨勢類型按2.2.1中標(biāo)準(zhǔn)分類。對于多趨勢段區(qū)域，研究其當(dāng)前階段趨勢，以分析植被的現(xiàn)狀。對于持續(xù)時長小于5 a的趨勢段，由于時間太短，不進(jìn)行趨勢分析。③ 根據(jù)林地類型分析林地總體趨勢和當(dāng)前階段趨勢的變化情況。

3 結(jié)果與分析

3.1 斷點(diǎn)檢測與分析

應(yīng)用DBEST算法對富縣NDVI時間序列數(shù)據(jù)進(jìn)行斷點(diǎn)檢測，時間序列中第1斷點(diǎn)NDVI變化幅度的空間分布與不同幅度面積占比如圖3所示。發(fā)生正向變化的區(qū)域占全縣面積的79.34%；負(fù)向變化的區(qū)域占20.66%，多位于道路兩側(cè)及居民點(diǎn)附近(圖1)。NDVI變化幅度為-0.2～0.2的區(qū)域占總面積66.47%，主要位于西部及東北地區(qū)，植被發(fā)生的變化較小；變化幅度絕對值超過0.2的區(qū)域占總面積的33.53%，位于富縣的東部及河流道路周圍區(qū)域，屬于NDVI劇烈變化的區(qū)域。

圖3 富縣歸一化植被指數(shù)變化幅度及其面積占比Figure 3 Change magnitude and area percent of NDVI in Fu County

圖4顯示：2000—2019年間富縣NDVI變化幅度大于0.2的區(qū)域斷點(diǎn)個數(shù)的空間分布及各類型的面積占比。其中，66.47%的區(qū)域斷點(diǎn)NDVI變化幅度小于0.2，主要分布在富縣西部和東北地區(qū)，說明NDVI變化幅度小，植被相對穩(wěn)定，因此不作進(jìn)一步分析。斷點(diǎn)NDVI變化幅度大于0.2的區(qū)域占總面積的33.53%，其中斷點(diǎn)數(shù)量為1或2個的區(qū)域占總面積的21.59%，位于東部地區(qū)；斷點(diǎn)數(shù)量多于4個的區(qū)域，僅占5.88%，集中分布于道路、河流沿線，植被變化頻繁，與人為活動有關(guān)。

圖4 富縣不同 NDVI變化幅度小于 0.2 與大于0.2的斷點(diǎn)數(shù)目及其面積占比Figure 4 Number of breakpoints and area percent of NDVI (＜0.2 and＞0.2) in Fu County

3.2 NDVI時間序列趨勢分析

3.2.1 總體趨勢 富縣 NDVI 總體趨勢有很強(qiáng)的空間差異。對趨勢值進(jìn)行分類，得到不同類別的分布及面積(圖5)。20 a間富縣超過97%的區(qū)域植被呈改善趨勢；退化區(qū)域占比不足2%，植被改善的面積遠(yuǎn)大于退化的面積。退化區(qū)域集中分布在富縣東部區(qū)域、道路及河流沿線，在西北部散落分布，其中東部部分人造地表退化面積較大；基本不變的區(qū)域面積最小，不足1%，分布較零散。

圖5 2000—2019 年富縣 NDVI 時間序列總體趨勢類別及其面積占比Figure 5 Overall trend level and area percent of NDVI time series in Fu County from 2000 to 2019

3.2.2 多 階 段 趨 勢 基 于 3.1 斷 點(diǎn) 檢 測 結(jié) 果 ，對2000—2019年富縣NDVI時間序列進(jìn)行多階段趨勢分析。由于66.47%的區(qū)域無斷點(diǎn)，因此該區(qū)域只分析總體趨勢；33.53%的區(qū)域具有1個及以上斷點(diǎn)，對該區(qū)域時間序列進(jìn)行分段，分析多階段趨勢中的當(dāng)前階段趨勢。當(dāng)前階段趨勢統(tǒng)計結(jié)果及分布(圖6A)表明：27.57%的區(qū)域當(dāng)前階段趨勢持續(xù)時長大于5 a，植被當(dāng)前階段改善的面積遠(yuǎn)大于退化的面積，改善面積占24.54%，分布于東部區(qū)域及道路河流沿線的周邊區(qū)域；退化面積僅占2.12%，基本不變面積不足退化面積的一半，分布較零散。當(dāng)前階段趨勢的開始時間(圖6B)均在2014年之前，時間分布跨度大；開始于2002年的區(qū)域面積最大，占12.76%，且均為改善趨勢，表明經(jīng)歷一次變化以后，區(qū)域植被一直呈現(xiàn)改善狀態(tài)；其次是開始于2014年的區(qū)域，占3.68%，主要位于道路、河流沿線；其余年份所占面積均不足2%。

圖6 2000—2019年富縣 NDVI時間序列當(dāng)前階段趨勢及其面積占比Figure 6 Current stage trend and area percent of NDVI time series in Fu County from 2000 to 2019

3.3 不同林地趨勢分析

如圖2所示：富縣植被退化的區(qū)域大部分分布在西北、西南、東北以及中部的一些地區(qū)，如今中部地區(qū)以及一些道路沿線已經(jīng)變?yōu)橛谰眯缘慕ㄔO(shè)用地；其他大部分均已改善，大多變?yōu)橛辛值睾褪枇值亍?019年富縣林地類型中，有林地占63.72%，疏林地占1.98%，灌木林地占13.27%，未成林造林地占6.07%，4種林地占富縣總面積的85.04%。如表1所示：有林地、疏林地、灌木林地、未成林造林地的改善面積占其總面積的比例分別為96.98%、94.32%、94.41%、92.97%，其他用地改善和基本不變的比例分別為86.66%和2.50%，退化的比例達(dá)10.84%，由于其他用地包括建設(shè)用地，說明大部分退化與城市和道路建設(shè)相關(guān)。整體來說，富縣生態(tài)建設(shè)工程取得了良好的效果，大部分林地處于改善的趨勢，特別是有林地和灌木林地。

表1 富縣林地類型趨勢情況統(tǒng)計Table 1 Statistics on trend of forest land types in Fu County

4 結(jié)論

本研究基于斷點(diǎn)檢測的多階段趨勢分析方法，研究了陜西省富縣近20 a植被的時空變化特征，得到以下結(jié)論：①第1斷點(diǎn)檢測結(jié)果表明：富縣發(fā)生正向變化的區(qū)域面積大，占79.34%。檢測結(jié)果中NDVI變化幅度小于0.2的區(qū)域占66.47%，主要分布在富縣西部和東北地區(qū)，植被相對穩(wěn)定；大于0.2的區(qū)域占33.53%，其中斷點(diǎn)數(shù)量為1或2個的區(qū)域比例為21.59%，位于東部地區(qū)；斷點(diǎn)數(shù)量多于4個的區(qū)域，僅占5.88%，集中分布于道路、河流沿線，變化頻繁，與人為活動有關(guān)。②植被的多階段趨勢表明：富縣共有27.57%的區(qū)域當(dāng)前階段趨勢持續(xù)時長大于5 a，當(dāng)前階段趨勢為改善的面積占24.54%，分布于東部區(qū)域及道路河流沿線的周邊區(qū)域；為退化的面積占2.12%，基本不變面積僅占0.91%，分布均較零散。當(dāng)前階段趨勢的開始時間均在2014年之前，時間分布跨度大，空間異質(zhì)性強(qiáng)，揭示了富縣植被變化的多樣性及復(fù)雜性。③近年來富縣有林地、疏林地、灌木林地、未成林造林地均是改善情況明顯，其中有林地和灌木林地中的改善比例最大。林地中呈改善態(tài)勢的面積占全縣總面積的81.83%，富縣生態(tài)建設(shè)工程取得了良好的效果。同時對于檢測到退化區(qū)域的林地，需要進(jìn)一步分析，采取有效措施提高質(zhì)量。

1999年以來，富縣積極響應(yīng)國家號召，實(shí)施退耕還林、荒山荒地造林、封山育林等措施。這些生態(tài)工程建設(shè)促進(jìn)了植被的恢復(fù)，然而在不同地區(qū)的效果是存在差異的。本研究獲得的NDVI的變化直接反映了植被恢復(fù)的效果，但DBEST算法參數(shù)設(shè)置對分析結(jié)果的影響、趨勢變化與影響因素的關(guān)系尚不清楚，在后續(xù)研究中，需從多尺度應(yīng)用DBEST算法，探究其應(yīng)用技術(shù)，并綜合考慮氣溫、降水、人為因素等對趨勢變化的影響。