基于多时相Landsat数据的汶川地震后王家岩滑坡植被恢复估计

叶席， 匡鸿海

西南大学 地理科学学院，重庆 400715

近年来，全球因地震引发的地质灾害频有发生[1-5]，这些灾害对植被造成了巨大破坏，对当地自然环境和生态系统造成一定程度的影响．地震灾区自然环境和生态系统的恢复是长期且复杂的过程，地表植被的生长状况是评价这一过程的极为重要的指标[6]，因此对地震灾区植被生长情况的监测不仅能够为该地区生态系统恢复的评估提供依据，还能预测其未来的发展趋势．遥感数据因其观测范围广、时效性强等特点在植被的动态监测和定量评估过程中体现出了巨大的优势，是植被覆盖恢复动态监测与评价最有效的数据来源[7]．

学者们利用卫星影像和航空数据对地震灾区植被恢复做了大量研究，王飞龙等[8]利用多期Landsat影像数据分析了汶川地震震中附近震后崩滑体上的植被动态恢复变化，发现该地区经过近9年的恢复，植被覆盖度恢复到0.74，与震前相比差值为0.08，根据拟合模型预计，2022年植被覆盖度能恢复到震前水平；李明威等[9]利用不同来源的遥感影像分析了北川县泥石流流域崩滑体时空演变特征，流域内植被覆盖度在2008 年“9.24”泥石流灾害后呈稳定恢复，到 2016年研究区植被覆盖度已恢复至较高水平；李京忠等[10]利用中分辨率成像光谱仪—植被指数(MODIS-NDVI)对都江堰龙溪河流域植被恢复进行了定量评估，认为植被覆盖度对地震损害的响应存在滞后现象；田颖颖[11]利用多其次滑坡编目数据、MODIS-NDVI数据、无人机航拍数据等分析了2015年尼泊尔地震后地震区植被演化特征，结果发现震区植被在地震发生后以减少为主，并在震后前两年明显增加后逐渐减少，截至2019年，震区内的植被和松散斜坡物质仍没有达到稳定状态；赵旦等[12]利用机载高分影像对汶川地震影响较为严重区域的农田和森林进行监测，发现震后5年农田的恢复较低，森林恢复情况总体较好；Yang等[13]基于SPOT数据利用灰色预测模型预测台湾1999年“9.21”地震受损植被恢复情况，发现先锋植被恢复时间约为2年，相对实际情况约有2个月的滞后期；Chuang等[14]基于SPOT数据利用马尔科夫链模型对台湾“9.21”地震受损植被恢复状况做出了评估，结果表明草可以作为滑坡区植被恢复的参考指标和滑坡区植被恢复的先导物种进行单独提取；Verdonen等[15]基于QuickBird和WorldView数据利用NDVI来监测俄罗斯北部滑坡冻土的植被状况，结果显示，山体滑坡对苔原植被产生显著而持久的影响，随着气候持续变暖，热喀斯特滑坡及其对植被的影响可能会在西伯利亚西北部和其他北极地区变得越来越普遍；Yang等[16]利用多年滑坡调查数据和MODIS-NDVI数据对汶川地震震中附近受损植被恢复情况进行了研究，结果显示，震后植被恢复可能促进了震后滑坡活动的减少，汶川地震对区域地震后滑坡频率的影响可能在地震后20年内消失．

以上这些研究在地震后植被恢复监测和分析方面做出了重要贡献．然而，通过时间序列遥感数据分析受损植被恢复情况并对其恢复到地震前水平做出预估的相关研究目前相对较少．2008年5月12日，四川省汶川县发生了MS8.0地震，震中位于四川省汶川县映秀镇(31.021°N，103.367°E)，震源深度约14 km[17]．汶川地震导致约56 000处同震滑坡发生，总面积约为811 km2[18]，灾害造成相关地区植被损坏严重．汶川地震所引发的北川县王家岩的同震滑坡，导致北川县大量的房屋被掩埋和人员伤亡[19]，自然环境也受到重创，因此对该地区的植被恢复情况进行研究，可以为地震后自然环境和生态系统恢复、区域规划发展等提供借鉴．

1 研究区与研究方法

1.1 研究区概况

研究区域为四川省北川县曲山镇王家岩及其周边区域(图1)，面积约9.67 km2．该区域属于亚热带季风性湿润气候，四季分明，气候温和，区域内雨量充沛，年均降水量约为1 280 mm，降水集中于6～9月；研究区海拔为630～1 610 m，地貌类型以山地、河谷为主，该区域属龙门山前山与后山交界地带，紧邻汶川地震的主发震断裂映秀—北川断裂的上盘，属于较为典型的地质构造不稳定区域[20]，主要由寒武系砂岩、页岩和片岩组成[21]．

底图来源于四川测绘地理信息局，审图号：图川审(2016)018号．图1 研究区地理位置

1.2 数据来源及数据预处理

由于Landsat影像数据具有时间序列长、易获取、分辨率较高等特点，本研究采用Landsat影像为数据源．由于Landsat 7的ETM+传感器的扫描线校正仪(SLC)于2003年5月31日发生故障，导致每景ETM+影像存在约22%的数据扫描空隙，产生条带现象，使得该数据的数据质量与数据应用都受到严重的制约和影响[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影像为本研究的数据源(数据来源于https：//earthexplorer.usgs.gov/)．利用ENVI 5.3软件对13个不同时相的遥感影像进行辐射定标、FLAASH大气校正以及目标区域掩膜提取等预处理．经预处理后的假彩色影像如图2所示．

1.3 归一化植被指数(NDVI)计算

归一化植被指数(NDVI)是反映植被生长状况的重要指标因子，是目前应用最为广泛的指标因子之一[23]，常用于对区域植被生长状况的定量评估．NDVI能减少地形地貌的干扰，克服红光波段反射率特别小(接近0)的情况，能够将绿色植被信息与土壤背景等信息进行区分，NDVI值介于-1～1之间，其中0和负值代表无植被覆盖的表面[24]．NDVI数值越大代表植被的覆盖状况越好，植被的生物量越多[25]．NDVI通常被定义为近红外波段与可见光红光波段反射率之差与反射率之和的比值[26]，其表达式为：

(1)

其中：NIR为近红外波段的反射率；R为可见光红光波段的反射率．

图2 研究区假彩色卫星图像

1.4 滑坡区域提取

由于地震前(2007/05/06)后(2008/07/11)2个时相遥感影像所得的NDVI数据存在一定的差异，本研究将研究区内地震前后时相的NDVI差值数据划分为多个小区域，再分别调整每个小区域的阈值，提取出研究区内震后的滑坡分布范围[27]．利用多时相进行变化检测时容易将与滑坡光谱变化特征相似的区域识别为滑坡[28]，因此同时对研究区震后的Google Earth高分辨率影像进行滑坡的目视解译，再将二者的结果相结合，得到最终的滑坡区域．

1.5 季节性调整

为保证研究精度，本研究应选取云量在30%以下的遥感影像，在同震滑坡所处时期的研究区内多期遥感影像的积云均较多，为满足研究精度，需要选取各年份不同时段的遥感影像．本研究选取的遥感影像获取时段处于3-10月之间．虽然所选影像避开了头年11月至次年2月天气较寒冷因而植被生长条件较差的时段，但由于研究区内的地貌和气候条件使得该区域内植物生长的季候差异明显，因此需要对所选取的遥感影像进行季节性调整处理．本研究采用Yang等[29]提出的通过简单的季节性调整来适应多时相图像的季节性变化的方法来修正选取的遥感影像之间的植被季候差异．该方法假设自然季节的影响在遥感影像上的空间分布是均匀的，以仅受季节性影响而未受同震滑坡影响的非滑坡区作为参考区域，从每幅遥感影像中提取非滑坡区域所有像元的NDVI并取其平均值作为背景值，计算每幅遥感影像背景值之间的差异，并以此作为季节性调整的偏移值[30]．利用所得的季节性调整偏移值对滑坡区域和研究区域的NDVI进行调整后，得到调整后的NDVI，进而对每幅遥感影像进行调整后的植被覆盖恢复率(VRR)计算．

1.6 VRR计算

遥感数据是区域范围内植被覆盖度提取的最有效数据源[31]．植被覆盖度(FVC)是反映植被生态功能的动态变量，可用于监测环境的变化和评价受损植被的恢复状况[29]．NDVI被广泛应用于植被覆盖度的计算，利用NDVI与光谱混合分析(NDVI-SMA)相结合的方法提取植被覆盖度，其原理是假设给定像元由植被与非植被覆盖的地表两部分所构成，其光谱信息只由这两个组分线性合成，它们各自的面积在像元中所占的比率即为各因子的权重，由它们的相对比例加权得到每个像元，其中植被所占像元的百分比即为该像元的FVC[32]．植被覆盖度表达式为：

(2)

式中：NDVIveg为全植被覆盖像元的NDVI值；NDVIsoil为裸土像元的NDVI值；FVC的取值范围为0～1．因受噪声、植被分布情况以及邻近地物辐射等因素的影响，植被覆盖度置信度的取值主要由研究区域遥感影像的实际情况来决定，在没有大量实测数据作参考的情况下，以累积频率为5%和95%作为置信区间[33]．本研究选择累积频率为5%的NDVI作为裸土像元的NDVI值，累积频率为95%的NDVI作为全植被覆盖像元的NDVI值来进行植被覆盖度的计算．

植被覆盖恢复率(VRR)可以定义为植被覆盖度受损后恢复的速率，它是根据植被受损前后的差异来监测和评估植被的恢复情况，以植被覆盖度为基础来进行计算[34]．植被覆盖恢复率的表达式为：

(3)

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

2 结果

2.1 滑坡提取结果

首先根据地震前后时相NDVI差值的阈值提取出滑坡分布的范围，同时结合震后Google Earth的高分辨率影像对研究区内滑坡进行目视解译的成果，得到研究区内滑坡区域(图3)，其中滑坡区域面积约为0.94 km2，非滑坡区域面积约为8.73 km2．

图3 研究区滑坡提取结果

2.2 NDVI计算结果

本研究中的NDVI均通过逐像素计算得到，研究区NDVI提取结果如图4所示．根据滑坡提取结果，本研究将研究区域划分为非滑坡区域以及滑坡区域，并分别计算出研究区域、非滑坡区域以及滑坡区域的NDVI均值，其结果如表1所示．滑坡区域的NDVI在同震滑坡发生后立即急剧下降，并在随后的年份出现波动，但是NDVI总体呈现上升趋势，对比图4，可以清晰看到滑坡区域植被恢复显著．在地震发生前，滑坡区域的植被覆盖良好，平均NDVI为0.732；由于研究区内建筑物、河流、坑塘、水坝等的存在，使得非滑坡区域的平均NDVI低于滑坡区域，非滑坡区域的平均NDVI为0.646；整个研究区平均NDVI为0.654，表明震前该区域植被条件较好，但地震以及滑坡发生后，滑坡区域的平均NDVI急剧下降至0.385．12年后的2020年8月13日，研究区、非滑坡区域、滑坡区域的平均NDVI分别上升到0.671，0.675，0.624．由于选取的遥感影像时段处于3月至10月之间，存在一定的物候差异，从2009年6月12日到2010年3月27日、2014年7月12日到2015年10月19日，研究区内的NDVI水平均呈现较大辐度的下降．因此通过一定的季节性调整方法对数据进行修正以降低植被的物候差异是十分必要的．

图4 研究区NDVI提取结果

表1 研究区、非滑坡区以及滑坡区NDVI

3 讨论与分析

3.1 VRR分析

本研究利用NDVI和经过季节性调整后的NDVI计算出FVC，进而得到VRR和调整后的VRR，结果如表2所示．为了评估整个研究区和滑坡区域的VRR变化情况，对研究区(图5a)和滑坡区域(图5b)构建VRR回归模型，并对其进行趋势分析．

地震以及滑坡发生之后，研究区及滑坡区的VRR总体呈现上升趋势，在震后早期，受损植被的恢复速率相对较快，然后逐渐趋于平缓．从整个研究区来看，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．出现异常情况可能是由季节因素导致的植物物候差异、降水条件的改变，以及人类活动的影响等所致．从滑坡区域来看，震后VRR的波动范围较小，总体呈现逐步上升的趋势．决定系数R2和p值分别用来衡量季节调整VRR回归模型的拟合程度以及显著性水平，从研究区和滑坡区域的VRR回归模型的p值来看，VRR回归模型是显著的．从决定系数R2来看，滑坡区域VRR回归模型R2值为0.863，研究区VRR回归模型R2值为0.688，滑坡区域的VRR回归模型优于整个研究区的VRR回归模型．

表2 研究区和滑坡区季节性调整前后的VRR

3.2 调整后的VRR分析

震后研究区及滑坡区域经季节性调整后的VRR值总体呈逐渐上升趋势，研究区及滑坡区域调整后的VRR的波动范围均较小，受损植被处于逐步恢复的状态．由研究区和滑坡区域调整后的VRR回归模型的p值可知，调整后的VRR回归模型是显著的．对于整个研究区(图5a)，调整后的VRR回归模型R2值为0.961，较原VRR模型的R2值(0.688)有显著提高，提高了0.273．对于滑坡区域(图5b)，调整后的VRR回归模型R2值为0.957，较原VRR回归模型R2值增加了0.094，存在明显提升．滑坡区域调整后VRR回归模型R2值为0.961，而研究区调整后VRR回归模型R2值为0.863，从决定系数R2来看，研究区调整后的VRR回归模型优于滑坡区．此外，就整个研究区而言，2010年3月27日观测到相对较低的VRR值为0.051，而2014年7月12日和2015年10月19日观测到相对较高的VRR分别为0.729和0.772．经季节性调整后，2010年3月27日相对较低的VRR向上调整为0.416，2014年7月12日和2015年10月19日相对较高的2个VRR分别向下调整为0.632和0.604．

图5 研究区(a)和滑坡区(b)VRR回归模型

在对研究区NDVI水平进行季节性调整后，对比调整前后所得的VRR数据，可以发现无论是对于整个研究区还是滑坡区域，经季节性调整后，VRR的数值波动范围均有所减小，修正了VRR的大幅异常波动情况．同时VRR模型在R2和p值上均优于未作季节性调整的原始VRR模型，说明对NDVI进行季节性修正有助于后续的VRR计算以及未来植被恢复情况的预测，季节性调整能较好地估计植被恢复状况．根据滑坡区域调整后的VRR模型，本研究预测该滑坡区域的受损植被大约需要26年才能完全恢复．

3.3 本研究方法对其他区域的适用性

本研究采用多时相Landsat遥感影像对研究区域植被的时空变化情况进行了监测和分析：首先对研究区2007-2020年的NDVI进行了分析，再以植被覆盖度为基础建立VRR模型，以及经季节性调整后的VRR模型，并进行相应的分析．研究区内受损植被震后早期的恢复速率相对较快，之后逐渐趋于平缓．根据滑坡区域经季节性调整后的VRR模型数据，本研究预测该滑坡区域的受损植被大约需要26年才能完全恢复．本研究选取的研究区域为典型的山地环境同震滑坡区域，并且研究区内植被随季节变化的物候差异明显．本研究所采用的研究方法可以为山地环境同震滑坡引起的植被在时空范围上的变化情况以及受损植被完全恢复的合理预估提供参考，可应用于其他类似区域的相应研究．

4 结论

本研究利用多时相Landsat遥感影像对汶川地震后北川县王家岩地区受到同震滑坡影响的植被的时空变化进行了动态监测与分析，并基于2007-2020年的13景Landsat遥感影像，对研究区的NDVI进行了分析，建立了相应的VRR回归模型，同时为消除季节变化对植被的影响，对研究数据进行了季节性调整．结果显示，在地震后12年的植被演替过程中，滑坡区域的NDVI总体水平呈上升趋势，但与地震发生前的植被条件相比还有一定差距；对于整个研究区和滑坡区域，季节性调整后的VRR的相关系数高于原始VRR系数，表明在植物生长季候明显的区域，利用多时相遥感影像进行长期的植被观测时需要保持相同的季节或对不同季节进行季节性调整．本研究根据滑坡区域调整后的VRR模型，估计该滑坡区域的受损植被大约需要26年才能完全恢复．本研究对于同震滑坡后山区植被恢复的长期跟踪监测与未来趋势预测具有一定参考意义，可应用于其他类似区域的研究．

地震后同震滑坡体上的植被恢复是一个长期且复杂的过程，目前对于引起这类变化的控制因素的研究成果相对较少，本研究也仅基于长时间序列NDVI数据来统计分析震后滑坡体植被恢复特征及变化趋势，进而评估震后滑坡体植被恢复完成的时间．例如，坡度是导致滑坡以及影响之后植被恢复的重要因素之一[35]，但本研究在分析同震滑坡后植被恢复的情况时，未结合坡度、海拔、坡向、降水、水系分布以及人类活动等影响植被生长和恢复的因子进行综合评估，在今后的研究中，可以利用最新的卫星图像和不同系列的传感器，在考虑季节变化和地形因子、降水等影响的基础上，建立可靠的植被恢复回归模型，对植被恢复进行持续监测．