基于多時相Landsat數(shù)據(jù)的汶川地震后王家?guī)r滑坡植被恢復估計

葉席， 匡鴻海

西南大學 地理科學學院，重慶 400715

近年來，全球因地震引發(fā)的地質災害頻有發(fā)生[1-5]，這些災害對植被造成了巨大破壞，對當?shù)刈匀画h(huán)境和生態(tài)系統(tǒng)造成一定程度的影響．地震災區(qū)自然環(huán)境和生態(tài)系統(tǒng)的恢復是長期且復雜的過程，地表植被的生長狀況是評價這一過程的極為重要的指標[6]，因此對地震災區(qū)植被生長情況的監(jiān)測不僅能夠為該地區(qū)生態(tài)系統(tǒng)恢復的評估提供依據(jù)，還能預測其未來的發(fā)展趨勢．遙感數(shù)據(jù)因其觀測范圍廣、時效性強等特點在植被的動態(tài)監(jiān)測和定量評估過程中體現(xiàn)出了巨大的優(yōu)勢，是植被覆蓋恢復動態(tài)監(jiān)測與評價最有效的數(shù)據(jù)來源[7]．

學者們利用衛(wèi)星影像和航空數(shù)據(jù)對地震災區(qū)植被恢復做了大量研究，王飛龍等[8]利用多期Landsat影像數(shù)據(jù)分析了汶川地震震中附近震后崩滑體上的植被動態(tài)恢復變化，發(fā)現(xiàn)該地區(qū)經(jīng)過近9年的恢復，植被覆蓋度恢復到0.74，與震前相比差值為0.08，根據(jù)擬合模型預計，2022年植被覆蓋度能恢復到震前水平；李明威等[9]利用不同來源的遙感影像分析了北川縣泥石流流域崩滑體時空演變特征，流域內植被覆蓋度在2008 年“9.24”泥石流災害后呈穩(wěn)定恢復，到 2016年研究區(qū)植被覆蓋度已恢復至較高水平；李京忠等[10]利用中分辨率成像光譜儀—植被指數(shù)(MODIS-NDVI)對都江堰龍溪河流域植被恢復進行了定量評估，認為植被覆蓋度對地震損害的響應存在滯后現(xiàn)象；田穎穎[11]利用多其次滑坡編目數(shù)據(jù)、MODIS-NDVI數(shù)據(jù)、無人機航拍數(shù)據(jù)等分析了2015年尼泊爾地震后地震區(qū)植被演化特征，結果發(fā)現(xiàn)震區(qū)植被在地震發(fā)生后以減少為主，并在震后前兩年明顯增加后逐漸減少，截至2019年，震區(qū)內的植被和松散斜坡物質仍沒有達到穩(wěn)定狀態(tài)；趙旦等[12]利用機載高分影像對汶川地震影響較為嚴重區(qū)域的農田和森林進行監(jiān)測，發(fā)現(xiàn)震后5年農田的恢復較低，森林恢復情況總體較好；Yang等[13]基于SPOT數(shù)據(jù)利用灰色預測模型預測臺灣1999年“9.21”地震受損植被恢復情況，發(fā)現(xiàn)先鋒植被恢復時間約為2年，相對實際情況約有2個月的滯后期；Chuang等[14]基于SPOT數(shù)據(jù)利用馬爾科夫鏈模型對臺灣“9.21”地震受損植被恢復狀況做出了評估，結果表明草可以作為滑坡區(qū)植被恢復的參考指標和滑坡區(qū)植被恢復的先導物種進行單獨提??；Verdonen等[15]基于QuickBird和WorldView數(shù)據(jù)利用NDVI來監(jiān)測俄羅斯北部滑坡凍土的植被狀況，結果顯示，山體滑坡對苔原植被產(chǎn)生顯著而持久的影響，隨著氣候持續(xù)變暖，熱喀斯特滑坡及其對植被的影響可能會在西伯利亞西北部和其他北極地區(qū)變得越來越普遍；Yang等[16]利用多年滑坡調查數(shù)據(jù)和MODIS-NDVI數(shù)據(jù)對汶川地震震中附近受損植被恢復情況進行了研究，結果顯示，震后植被恢復可能促進了震后滑坡活動的減少，汶川地震對區(qū)域地震后滑坡頻率的影響可能在地震后20年內消失．

以上這些研究在地震后植被恢復監(jiān)測和分析方面做出了重要貢獻．然而，通過時間序列遙感數(shù)據(jù)分析受損植被恢復情況并對其恢復到地震前水平做出預估的相關研究目前相對較少．2008年5月12日，四川省汶川縣發(fā)生了MS8.0地震，震中位于四川省汶川縣映秀鎮(zhèn)(31.021°N，103.367°E)，震源深度約14 km[17]．汶川地震導致約56 000處同震滑坡發(fā)生，總面積約為811 km2[18]，災害造成相關地區(qū)植被損壞嚴重．汶川地震所引發(fā)的北川縣王家?guī)r的同震滑坡，導致北川縣大量的房屋被掩埋和人員傷亡[19]，自然環(huán)境也受到重創(chuàng)，因此對該地區(qū)的植被恢復情況進行研究，可以為地震后自然環(huán)境和生態(tài)系統(tǒng)恢復、區(qū)域規(guī)劃發(fā)展等提供借鑒．

1 研究區(qū)與研究方法

1.1 研究區(qū)概況

研究區(qū)域為四川省北川縣曲山鎮(zhèn)王家?guī)r及其周邊區(qū)域(圖1)，面積約9.67 km2．該區(qū)域屬于亞熱帶季風性濕潤氣候，四季分明，氣候溫和，區(qū)域內雨量充沛，年均降水量約為1 280 mm，降水集中于6～9月；研究區(qū)海拔為630～1 610 m，地貌類型以山地、河谷為主，該區(qū)域屬龍門山前山與后山交界地帶，緊鄰汶川地震的主發(fā)震斷裂映秀—北川斷裂的上盤，屬于較為典型的地質構造不穩(wěn)定區(qū)域[20]，主要由寒武系砂巖、頁巖和片巖組成[21]．

底圖來源于四川測繪地理信息局，審圖號：圖川審(2016)018號．圖1 研究區(qū)地理位置

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

由于Landsat影像數(shù)據(jù)具有時間序列長、易獲取、分辨率較高等特點，本研究采用Landsat影像為數(shù)據(jù)源．由于Landsat 7的ETM+傳感器的掃描線校正儀(SLC)于2003年5月31日發(fā)生故障，導致每景ETM+影像存在約22%的數(shù)據(jù)掃描空隙，產(chǎn)生條帶現(xiàn)象，使得該數(shù)據(jù)的數(shù)據(jù)質量與數(shù)據(jù)應用都受到嚴重的制約和影響[22]，因此本研究選取成像時間為2007/05/06，2008/07/11，2009/06/12，2010/03/27，2011/05/07獲取的軌道號為129/38的5期Landsat 5 TM影像，以及成像時間為2013/05/22，2014/07/12，2015/10/19，2016/04/28，2017/05/01，2018/06/05，2019/08/11，2020/08/13獲取的軌道號為129/38的8期Landsat 8 OLI影像為本研究的數(shù)據(jù)源(數(shù)據(jù)來源于https：//earthexplorer.usgs.gov/)．利用ENVI 5.3軟件對13個不同時相的遙感影像進行輻射定標、FLAASH大氣校正以及目標區(qū)域掩膜提取等預處理．經(jīng)預處理后的假彩色影像如圖2所示．

1.3 歸一化植被指數(shù)(NDVI)計算

歸一化植被指數(shù)(NDVI)是反映植被生長狀況的重要指標因子，是目前應用最為廣泛的指標因子之一[23]，常用于對區(qū)域植被生長狀況的定量評估．NDVI能減少地形地貌的干擾，克服紅光波段反射率特別小(接近0)的情況，能夠將綠色植被信息與土壤背景等信息進行區(qū)分，NDVI值介于-1～1之間，其中0和負值代表無植被覆蓋的表面[24]．NDVI數(shù)值越大代表植被的覆蓋狀況越好，植被的生物量越多[25]．NDVI通常被定義為近紅外波段與可見光紅光波段反射率之差與反射率之和的比值[26]，其表達式為：

(1)

其中：NIR為近紅外波段的反射率；R為可見光紅光波段的反射率．

圖2 研究區(qū)假彩色衛(wèi)星圖像

1.4 滑坡區(qū)域提取

由于地震前(2007/05/06)后(2008/07/11)2個時相遙感影像所得的NDVI數(shù)據(jù)存在一定的差異，本研究將研究區(qū)內地震前后時相的NDVI差值數(shù)據(jù)劃分為多個小區(qū)域，再分別調整每個小區(qū)域的閾值，提取出研究區(qū)內震后的滑坡分布范圍[27]．利用多時相進行變化檢測時容易將與滑坡光譜變化特征相似的區(qū)域識別為滑坡[28]，因此同時對研究區(qū)震后的Google Earth高分辨率影像進行滑坡的目視解譯，再將二者的結果相結合，得到最終的滑坡區(qū)域．

1.5 季節(jié)性調整

為保證研究精度，本研究應選取云量在30%以下的遙感影像，在同震滑坡所處時期的研究區(qū)內多期遙感影像的積云均較多，為滿足研究精度，需要選取各年份不同時段的遙感影像．本研究選取的遙感影像獲取時段處于3-10月之間．雖然所選影像避開了頭年11月至次年2月天氣較寒冷因而植被生長條件較差的時段，但由于研究區(qū)內的地貌和氣候條件使得該區(qū)域內植物生長的季候差異明顯，因此需要對所選取的遙感影像進行季節(jié)性調整處理．本研究采用Yang等[29]提出的通過簡單的季節(jié)性調整來適應多時相圖像的季節(jié)性變化的方法來修正選取的遙感影像之間的植被季候差異．該方法假設自然季節(jié)的影響在遙感影像上的空間分布是均勻的，以僅受季節(jié)性影響而未受同震滑坡影響的非滑坡區(qū)作為參考區(qū)域，從每幅遙感影像中提取非滑坡區(qū)域所有像元的NDVI并取其平均值作為背景值，計算每幅遙感影像背景值之間的差異，并以此作為季節(jié)性調整的偏移值[30]．利用所得的季節(jié)性調整偏移值對滑坡區(qū)域和研究區(qū)域的NDVI進行調整后，得到調整后的NDVI，進而對每幅遙感影像進行調整后的植被覆蓋恢復率(VRR)計算．

1.6 VRR計算

遙感數(shù)據(jù)是區(qū)域范圍內植被覆蓋度提取的最有效數(shù)據(jù)源[31]．植被覆蓋度(FVC)是反映植被生態(tài)功能的動態(tài)變量，可用于監(jiān)測環(huán)境的變化和評價受損植被的恢復狀況[29]．NDVI被廣泛應用于植被覆蓋度的計算，利用NDVI與光譜混合分析(NDVI-SMA)相結合的方法提取植被覆蓋度，其原理是假設給定像元由植被與非植被覆蓋的地表兩部分所構成，其光譜信息只由這兩個組分線性合成，它們各自的面積在像元中所占的比率即為各因子的權重，由它們的相對比例加權得到每個像元，其中植被所占像元的百分比即為該像元的FVC[32]．植被覆蓋度表達式為：

(2)

式中：NDVIveg為全植被覆蓋像元的NDVI值；NDVIsoil為裸土像元的NDVI值；FVC的取值范圍為0～1．因受噪聲、植被分布情況以及鄰近地物輻射等因素的影響，植被覆蓋度置信度的取值主要由研究區(qū)域遙感影像的實際情況來決定，在沒有大量實測數(shù)據(jù)作參考的情況下，以累積頻率為5%和95%作為置信區(qū)間[33]．本研究選擇累積頻率為5%的NDVI作為裸土像元的NDVI值，累積頻率為95%的NDVI作為全植被覆蓋像元的NDVI值來進行植被覆蓋度的計算．

植被覆蓋恢復率(VRR)可以定義為植被覆蓋度受損后恢復的速率，它是根據(jù)植被受損前后的差異來監(jiān)測和評估植被的恢復情況，以植被覆蓋度為基礎來進行計算[34]．植被覆蓋恢復率的表達式為：

(3)

式中：FVC0為同震滑坡前的植被覆蓋度；FVC1為同震滑坡后的植被覆蓋度；FVCi為同震滑坡后i(i=1，2，…，12)年植被覆蓋度．

2 結果

2.1 滑坡提取結果

首先根據(jù)地震前后時相NDVI差值的閾值提取出滑坡分布的范圍，同時結合震后Google Earth的高分辨率影像對研究區(qū)內滑坡進行目視解譯的成果，得到研究區(qū)內滑坡區(qū)域(圖3)，其中滑坡區(qū)域面積約為0.94 km2，非滑坡區(qū)域面積約為8.73 km2．

圖3 研究區(qū)滑坡提取結果

2.2 NDVI計算結果

本研究中的NDVI均通過逐像素計算得到，研究區(qū)NDVI提取結果如圖4所示．根據(jù)滑坡提取結果，本研究將研究區(qū)域劃分為非滑坡區(qū)域以及滑坡區(qū)域，并分別計算出研究區(qū)域、非滑坡區(qū)域以及滑坡區(qū)域的NDVI均值，其結果如表1所示．滑坡區(qū)域的NDVI在同震滑坡發(fā)生后立即急劇下降，并在隨后的年份出現(xiàn)波動，但是NDVI總體呈現(xiàn)上升趨勢，對比圖4，可以清晰看到滑坡區(qū)域植被恢復顯著．在地震發(fā)生前，滑坡區(qū)域的植被覆蓋良好，平均NDVI為0.732；由于研究區(qū)內建筑物、河流、坑塘、水壩等的存在，使得非滑坡區(qū)域的平均NDVI低于滑坡區(qū)域，非滑坡區(qū)域的平均NDVI為0.646；整個研究區(qū)平均NDVI為0.654，表明震前該區(qū)域植被條件較好，但地震以及滑坡發(fā)生后，滑坡區(qū)域的平均NDVI急劇下降至0.385．12年后的2020年8月13日，研究區(qū)、非滑坡區(qū)域、滑坡區(qū)域的平均NDVI分別上升到0.671，0.675，0.624．由于選取的遙感影像時段處于3月至10月之間，存在一定的物候差異，從2009年6月12日到2010年3月27日、2014年7月12日到2015年10月19日，研究區(qū)內的NDVI水平均呈現(xiàn)較大輻度的下降．因此通過一定的季節(jié)性調整方法對數(shù)據(jù)進行修正以降低植被的物候差異是十分必要的．

圖4 研究區(qū)NDVI提取結果

表1 研究區(qū)、非滑坡區(qū)以及滑坡區(qū)NDVI

3 討論與分析

3.1 VRR分析

本研究利用NDVI和經(jīng)過季節(jié)性調整后的NDVI計算出FVC，進而得到VRR和調整后的VRR，結果如表2所示．為了評估整個研究區(qū)和滑坡區(qū)域的VRR變化情況，對研究區(qū)(圖5a)和滑坡區(qū)域(圖5b)構建VRR回歸模型，并對其進行趨勢分析．

地震以及滑坡發(fā)生之后，研究區(qū)及滑坡區(qū)的VRR總體呈現(xiàn)上升趨勢，在震后早期，受損植被的恢復速率相對較快，然后逐漸趨于平緩．從整個研究區(qū)來看，2010年3月27日、2014年7月12日以及2015年10月19日其VRR存在明顯異常，其中2010年3月27日的VRR相對偏低，其值為0.051，而2014年7月12日和2015年10月19日的VRR相對偏高，其值分別為0.729和0.772．出現(xiàn)異常情況可能是由季節(jié)因素導致的植物物候差異、降水條件的改變，以及人類活動的影響等所致．從滑坡區(qū)域來看，震后VRR的波動范圍較小，總體呈現(xiàn)逐步上升的趨勢．決定系數(shù)R2和p值分別用來衡量季節(jié)調整VRR回歸模型的擬合程度以及顯著性水平，從研究區(qū)和滑坡區(qū)域的VRR回歸模型的p值來看，VRR回歸模型是顯著的．從決定系數(shù)R2來看，滑坡區(qū)域VRR回歸模型R2值為0.863，研究區(qū)VRR回歸模型R2值為0.688，滑坡區(qū)域的VRR回歸模型優(yōu)于整個研究區(qū)的VRR回歸模型．

表2 研究區(qū)和滑坡區(qū)季節(jié)性調整前后的VRR

3.2 調整后的VRR分析

震后研究區(qū)及滑坡區(qū)域經(jīng)季節(jié)性調整后的VRR值總體呈逐漸上升趨勢，研究區(qū)及滑坡區(qū)域調整后的VRR的波動范圍均較小，受損植被處于逐步恢復的狀態(tài)．由研究區(qū)和滑坡區(qū)域調整后的VRR回歸模型的p值可知，調整后的VRR回歸模型是顯著的．對于整個研究區(qū)(圖5a)，調整后的VRR回歸模型R2值為0.961，較原VRR模型的R2值(0.688)有顯著提高，提高了0.273．對于滑坡區(qū)域(圖5b)，調整后的VRR回歸模型R2值為0.957，較原VRR回歸模型R2值增加了0.094，存在明顯提升．滑坡區(qū)域調整后VRR回歸模型R2值為0.961，而研究區(qū)調整后VRR回歸模型R2值為0.863，從決定系數(shù)R2來看，研究區(qū)調整后的VRR回歸模型優(yōu)于滑坡區(qū)．此外，就整個研究區(qū)而言，2010年3月27日觀測到相對較低的VRR值為0.051，而2014年7月12日和2015年10月19日觀測到相對較高的VRR分別為0.729和0.772．經(jīng)季節(jié)性調整后，2010年3月27日相對較低的VRR向上調整為0.416，2014年7月12日和2015年10月19日相對較高的2個VRR分別向下調整為0.632和0.604．

圖5 研究區(qū)(a)和滑坡區(qū)(b)VRR回歸模型

在對研究區(qū)NDVI水平進行季節(jié)性調整后，對比調整前后所得的VRR數(shù)據(jù)，可以發(fā)現(xiàn)無論是對于整個研究區(qū)還是滑坡區(qū)域，經(jīng)季節(jié)性調整后，VRR的數(shù)值波動范圍均有所減小，修正了VRR的大幅異常波動情況．同時VRR模型在R2和p值上均優(yōu)于未作季節(jié)性調整的原始VRR模型，說明對NDVI進行季節(jié)性修正有助于后續(xù)的VRR計算以及未來植被恢復情況的預測，季節(jié)性調整能較好地估計植被恢復狀況．根據(jù)滑坡區(qū)域調整后的VRR模型，本研究預測該滑坡區(qū)域的受損植被大約需要26年才能完全恢復．

3.3 本研究方法對其他區(qū)域的適用性

本研究采用多時相Landsat遙感影像對研究區(qū)域植被的時空變化情況進行了監(jiān)測和分析：首先對研究區(qū)2007-2020年的NDVI進行了分析，再以植被覆蓋度為基礎建立VRR模型，以及經(jīng)季節(jié)性調整后的VRR模型，并進行相應的分析．研究區(qū)內受損植被震后早期的恢復速率相對較快，之后逐漸趨于平緩．根據(jù)滑坡區(qū)域經(jīng)季節(jié)性調整后的VRR模型數(shù)據(jù)，本研究預測該滑坡區(qū)域的受損植被大約需要26年才能完全恢復．本研究選取的研究區(qū)域為典型的山地環(huán)境同震滑坡區(qū)域，并且研究區(qū)內植被隨季節(jié)變化的物候差異明顯．本研究所采用的研究方法可以為山地環(huán)境同震滑坡引起的植被在時空范圍上的變化情況以及受損植被完全恢復的合理預估提供參考，可應用于其他類似區(qū)域的相應研究．

4 結論

本研究利用多時相Landsat遙感影像對汶川地震后北川縣王家?guī)r地區(qū)受到同震滑坡影響的植被的時空變化進行了動態(tài)監(jiān)測與分析，并基于2007-2020年的13景Landsat遙感影像，對研究區(qū)的NDVI進行了分析，建立了相應的VRR回歸模型，同時為消除季節(jié)變化對植被的影響，對研究數(shù)據(jù)進行了季節(jié)性調整．結果顯示，在地震后12年的植被演替過程中，滑坡區(qū)域的NDVI總體水平呈上升趨勢，但與地震發(fā)生前的植被條件相比還有一定差距；對于整個研究區(qū)和滑坡區(qū)域，季節(jié)性調整后的VRR的相關系數(shù)高于原始VRR系數(shù)，表明在植物生長季候明顯的區(qū)域，利用多時相遙感影像進行長期的植被觀測時需要保持相同的季節(jié)或對不同季節(jié)進行季節(jié)性調整．本研究根據(jù)滑坡區(qū)域調整后的VRR模型，估計該滑坡區(qū)域的受損植被大約需要26年才能完全恢復．本研究對于同震滑坡后山區(qū)植被恢復的長期跟蹤監(jiān)測與未來趨勢預測具有一定參考意義，可應用于其他類似區(qū)域的研究．

地震后同震滑坡體上的植被恢復是一個長期且復雜的過程，目前對于引起這類變化的控制因素的研究成果相對較少，本研究也僅基于長時間序列NDVI數(shù)據(jù)來統(tǒng)計分析震后滑坡體植被恢復特征及變化趨勢，進而評估震后滑坡體植被恢復完成的時間．例如，坡度是導致滑坡以及影響之后植被恢復的重要因素之一[35]，但本研究在分析同震滑坡后植被恢復的情況時，未結合坡度、海拔、坡向、降水、水系分布以及人類活動等影響植被生長和恢復的因子進行綜合評估，在今后的研究中，可以利用最新的衛(wèi)星圖像和不同系列的傳感器，在考慮季節(jié)變化和地形因子、降水等影響的基礎上，建立可靠的植被恢復回歸模型，對植被恢復進行持續(xù)監(jiān)測．