基于NDVI变化检测的滑坡遥感精细识别

基于NDVI变化检测的滑坡遥感精细识别

郭擎1,朱丽娅1,2,李安1,顾铃燕1,2

（1.中国科学院空天信息创新研究院；

2.南京大学金陵学院）

基金资助：国家自然科学基金面上项目“多源多时相遥感图像光谱特征鲁棒性融合研究”。

作者简介：郭擎(-)，女，河南驻马店人，博士，研究员，主要从事遥感信息提取与滑坡灾害监测研究。E?mail：guoqing

1引言

全球范围的地震断裂带现阶段已进入比较活跃的时期，而且随着经济的高速发展，人们对资源的不合理开采及对环境的破坏，使得自然灾害频繁发生，特别是滑坡灾害，它不但可以由强降雨和过度挖掘开采引发，还作为地震主要的次生灾害给人类的生存环境带来巨大破坏，造成不可挽回的经济损失。因此如何在滑坡发生后组织滑坡信息的有效提取，是进行灾后救援的重要组成部分。

灾后航空影像全野外调查识别滑坡很及时，可根据实际需要制定针对性的航线进行成像，但也存在着不足之处，航空影像的图幅范围较小，获取成本较高，由于航空影像一般由飞机拍摄，畸变较大。传统的滑坡遥感解译—目视解译，主要根据地物的色彩、色调、纹理和阴影等所形成的几何形态特征进行解译，解译时需要有经验的专业判读人员。这两种滑坡信息提取方法费时费力，大量的遥感信息不能得到充分的发掘和提取，难以满足滑坡灾害调查尤其是严峻的滑坡灾害综合防治的需要。因此，基于遥感影像的光谱、几何、纹理等特征的滑坡灾害信息识别新方法的探索开始出现。

基于单时相影像的滑坡识别方法，利用一幅某一时相的遥感数据进行滑坡识别，如万保峰等首先对研究区影像进行预处理和特征提取，然后采用灰度共生矩阵纹理分析方法对重庆武隆县的滑坡进行识别，识别效果明显，可大致确定滑坡形态，为滑坡危害性评估提供依据；朱娇等提出一种基于归一化植被指数（NormalizedDifferenceVegetationindex，NDVI）、土壤亮度指数、第一主成分信息的快速、高效的滑坡信息提取思路，分别对面向像元和基于对象的方法进行了验证。该类方法充分利用了遥感影像的光谱、几何、纹理等信息。但单时相的滑坡识别方法容易将与滑坡光谱特征相似的地物（建筑、裸土和道路等）视为滑坡，特别地，在存在较多岩石、裸地、道路等的山区，滑坡识别结果的错分率、误分率较高。

多时相的滑坡识别，利用多幅不同时间的遥感影像进行滑坡识别，如李松等以多时相遥感影像为数据源，结合比率变换和纹理分析进行变化检测来识别滑坡灾害信息，突破了影像光谱特征的局限；胡德勇等通过多期遥感图像的多尺度、多层次分割，建立研究区的滑坡识别知识规则，对马来西亚热带雨林中发生的滑坡进行识别；Behling则在分析滑坡发生的相关地表变化与其他土地覆盖变化之间的NDVI时序变化轨迹的基础上进行滑坡识别；虎振兴等以多源影像作为数据源，在Behling等方法的基础上，研究了滑坡发生后植被恢复速率的量化方法，构建了利用NDVI时序变化信息识别历史滑坡的方法，对我国云南鲁甸地区进行了滑坡识别。唐世超等则是利用多时相雷达图像对震后滑坡进行识别。随着近几年现代遥感技术的发展，时间序列遥感数据逐渐丰富，利用这种数据进行滑坡识别的方法也开始涌现。

为了提高遥感滑坡应急的时效性，并有效降低滑坡的误分率、漏分率，在分析当前滑坡灾害识别方法的基础上，以多时相遥感影像的变化检测为基础，结合面向对象的几何规则进行滑坡识别，利用滑坡前后NDVI差异较大的特点检测滑坡，排除光谱特征与滑坡相似的地物，并结合遥感影像的几何特征，做到高效提取滑坡的同时，提高滑坡遥感识别结果的准确性和可靠性。

2研究区与数据

2.1研究区概况

九寨沟位于四川省阿坝藏族羌族自治州九寨沟县境内，地处青藏高原东缘向四川盆地过渡地带，是一条纵深50km的山沟谷地，总面积.74km2。该区域内地势南高北低，谷深切，高差悬殊，区域北缘九寨沟口的海拔仅m，中部峰岭的海拔均m以上，南缘达m以上。九寨沟地质背景复杂，碳酸盐岩分布广泛，褶皱断裂发育，新构造运动强烈，地壳抬升幅度大，造就了多种地貌。

北京时间年8月8日21时19分46秒，四川省北部阿坝州九寨沟县发生7.0级地震，震中位于33.20°N，.82°E。九寨沟核心景区西部5km处比芒村，截至年8月13日20时，地震造成25人死亡，人受伤，6人失联，人受灾，间房屋不同程度受损（其中倒塌76间）。

此次地震是发生在走滑型断裂上的走滑型地震，震中及邻区发育有近南北向的岷江断裂、北西走向的东昆仑断裂东段，两条断裂分别为巴颜喀拉块体的东边界和北边界。发震区域构造复杂，位于南北地震带的中断，是巴颜喀拉块体边界断裂持续活动的结果。本次地震中强烈的震动和显著的地壳形变导致了大量同震滑坡的发生，且滑坡主要以中小型浅层滑坡和崩塌为主。

结合九寨沟震后滑坡的分布和滑坡快速识别方法的特点，确定九寨沟县的上四寨村—九寨天堂段为研究区，如图1所示。该地区发生多处滑坡，且有云量较少的遥感影像覆盖。

图1研究区概况

2.2数据及其预处理

研究中收集到覆盖该研究区的两个时相的高分辨率图像，分别是高分一号（GF-1）8m分辨率多光谱相机PMS2在年8月1日的影像（震前）以及年8月16日的影像（震后）。将多光谱影像按辐射定标、大气校正、几何校正的步骤完成预处理，同时以震后的遥感影像作为基准影像，对震前影像进行配准，以尽量减少误差，最终裁剪出预先设定的研究区。

图2影像预处理流程

3基于变化检测的滑坡精细识别

3.1滑坡预选区确定

本研究滑坡识别主要用到的是基于遥感影像变化检测的方法，该方法利用多时相获取的覆盖同一地表区域的遥感影像以及其他辅助数据来分析地表的变化。首先，在遥感影像经过预处理的基础上，选择滑坡发生前后的两幅NDVI影像（即年8月的九寨沟地区滑前高分影像得到归一化植被指数图与8月的滑后归一化植被指数图），进行基于归一化植被指数的变化检测，并设定一个扰动变化参数I：

I=1+(1+Npre)-(1+Npost)(1+Npre)+(1+Npost)(1)

其中：Npre为滑坡发生前影像的NDVI；Npost为滑坡发生后影像的NDVI。

对同一研究区，利用研究区双时相影像进行NDVI变化检测，由式（1）确定的植被扰动变化参数I初步判断非滑坡要素。地震发生前滑坡未发生，归一化植被指数大；滑坡发生后，滑坡区域的归一化植被指数则明显降低。因此，归一化植被指数的变化越大，扰动变化参数I也随之变大，利用该特点可以有效排除非滑坡要素：由于两幅影像数据都来源于8月夏季的九寨沟，因此，没有发生滑坡且覆盖研究区大部分的植被可首先被排除；同理，NDVI在滑前滑后并未发生变化的建筑及裸地等也被排除。

对发生变化的区域进行后续滑坡识别，依据Behling方法中的植被扰动分类规则（阈值由统计分析的结果决定）进行分类，分类依据如表1所示。

依据表1的扰动分类规则，按照植被扰动参数和归一化植被指数的值将每个研究区像素分为9类，对9类像素多次组合实验，特别是选择植被扰动参数较大的像素，并结合目视解译的结果，选出像素组合分布中能精确表达滑坡实际分布的部分，将其称为滑坡预选区。

3.2滑坡的精细提取

确定滑坡预选区后，影响滑坡提取精度的主要因素是一些误判，如部分道路被误判为滑坡等，为了排除这些因素的影响进一步提高滑坡提取精度，结合面向对象的几何特征进行分类，进一步进行滑坡提取。

对滑后研究区影像进行分割，与多尺度分割不同，滑坡的进一步提取只需把尺度确定在高效提取出对提取精度造成影响的地物因素即可，如道路。因此对研究区多次尝试使用不同的尺度参数进行分割，直到确定一个尺度参数，在该尺度参数下能比较精确地得到目标地物分布范围。并在该分割尺度基础上，结合地物的几何特征，通过找出对滑坡提取精度造成影响的目标地物，完成地物的排除。本研究所用几何特征是对象的圆度和伸长率。道路对象的圆度一般较低，圆度越小，提取的道路越准确。道路对象的伸长率一般较高，伸长率越小，提取的地物越不包含道路。

对象的圆度Roundness正比于多边形的面积与多边形最长直径平方的比率。最长直径指的是多边形外接矩形框的长轴长度。圆的圆度为1，正方形的圆度为4/π，圆度的计算公式为：

Roundness=4*Areaπ*Maxaxislen2（2）

其中：Area为对象的多边形的总面积减去多边形中洞的面积；Maxaxislen为对象的长轴长，由多边形的外接矩形确定。

对象的伸长率Elongation，表明多边形长轴和短轴的比率。长轴和短轴来自于多边形的外接矩形框。正方形的伸长率为1，矩形的伸长率大于1。伸长率的计算公式为：

Elongation=MaxaxislenMinaxislen（3）

其中：Minaxislen为对象的短轴长，由多边形的外接矩形确定。

最后将地物排除效果与滑坡预选区的结果互相结合与叠加，完成最终的滑坡提取。

4实验结果与分析

4.1实验设置

为了验证本研究方法的有效性，对方法进行绝对验证和相对验证。将本研究方法结果与人工目视解释结果对比进行绝对验证，将本研究方法结果与同样面向对象和规则的单时相滑坡识别方法进行对比进行相对验证。面向对象和规则的单时相滑坡识别方法主要包括影像分割、滑坡分类规则的建立以及滑坡信息提取。通过影像分割使同质像元组成大小不同的对象，利用对象的光谱特征、纹理特征和几何特征等进行分类，充分利用遥感影像信息，有效克服像元层次分类的不足。根据研究区的地形地貌特征和滑坡特征，建立面向对象的滑坡识别规则。

经过多次试验尝试，根据文献确定规则为0.02NDVI0.4，0.07蓝光波段反射率0.3，NDVI纹理均值0.，1伸长率2.7，0.19圆度0.65时，滑坡识别效果最佳。

研究的基于变化检测和几何规则的多时相滑坡识别方法，依据3.1节尤其是最后一段描述的方法，在经过多次实验后，选择植被扰动分类中的5、6、9等3类像素组合作为滑坡候选区，此时对滑坡识别精度造成影响的是部分道路区域，而圆度和伸长率能不同程度地对道路进行排除，从实验结果看，对于本文中的研究区，圆度对道路的排除效果更好。因此进一步建立0.19圆度0.65几何规则对道路进行排除，在以上条件下，能较精确地确定滑坡的分布区域。

为了定量分析与比较滑坡识别的精度，将目视解译的结果作为正确分类的标准，计算混淆矩阵，对最终的滑坡识别结果进行分类精度评价。

4.2实验结果分析

单时相对比方法和本文多时相变化检测方法最终提取的滑坡结果以及目视解释的滑坡识别结果如图3所示。从图3所提取出的总体滑坡分布可以看出，面向对象和规则的单时相滑坡识别方法与基于变化检测和几何规则的多时相方法相比，滑坡漏分与误分的现象比较严重，光谱特征与滑坡相似的地物没有被排除，如上四寨村（见图3中绿框）和图3中右上角大多数建筑不能通过地物识别规则排除，部分裸地也被误判为滑坡，直接导致滑坡提取精度不高。单时相滑坡识别方法中每条规则的确立比较复杂，规则阈值确立过程中的不确定因素高，同时由于不同地区地形地貌因子不同，这种单时相滑坡识别方法的规则建立每次都需要因地制宜地发生改变，较为费时费力。

图3滑坡识别结果

而本研究提出的方法能很好地解决这一问题，该方法通过归一化植被指数的变化，识别出研究区内植被发生较大变化的像素并快速建立滑坡候选区，排除大部分与滑坡光谱特征相似的地物（如建筑、裸地等），因为这些地物的归一化植被指数在滑坡前后没有发生较大变化，根据植被扰动参数可判断为非滑坡。这种通过归一化植被指数变化确定滑坡的方法，不需要考虑过多的地形地貌因子的限制，滑坡识别时效性高，可以适用于大部分地区。并且本研究方法滑坡识别的错分与误分主要发生在部分道路上，通过进一步建立识别规则可以有效降低。

局部道路和建筑的排除效果细节示例分别如图4和图5所示。从图4可以看出，多时相变化检测方法利用圆度几何规则能对道路进行有效地排除。而图4中的红色道路在单时相方法图3（c）中都被误判为滑坡。结合图5中建筑的分布，将图3（b）与图3（c）同区域的建筑排除效果进行对比，可以发现基于变化检测的多时相方法对建筑的排除效果较好，利用建筑区的归一化植被指数在滑坡前后没有发生较大变化的特点，能够有效提升滑坡的识别精度。

图4部分道路的排除效果

除了上面对结果的视觉定性分析，还进行了客观定量评价。滑坡提取的精度评价结果如表2所示，可以看到，用基于变化检测和几何规则相结合的多时相方法进行滑坡提取的精度较高，和面向对象的单时相方法相比，前者制图精度高达88.80%，用户精度高达81.19%，都大幅超过面向对象单时相法的精度，漏分误差及错分误差分别下降23.22%和11.72%。综合上面的绝对验证和相对验证，视觉定性分析和客观定量评价的结果可知，本研究方法可为有效组织滑坡灾后救援与重建工作提供可靠依据。

图5部分建筑的排除效果

5结语

研究提出了基于多时相遥感影像归一化植被指数变化检测的滑坡遥感精细识别方法，并结合地物的几何信息，有效降低传统滑坡识别方法中的错分误差和漏分误差，提升滑坡遥感识别结果的准确性和可靠性。以九寨沟滑坡为例，识别的制图精度和用户精度分别达到88.80%和81.19%，比单时相滑坡识别方法的漏分误差和错分误差分别降低了23.22%和11.72%。滑坡的遥感识别方法解决传统目视解译识别方法耗费大量人力物力以及过度依赖解译专家经验的问题，为有效组织灾后救援及灾后重建等提供了宝贵的时间。

与面向对象和规则的单时相遥感影像滑坡检测方法相比，该方法更好地排除了与滑坡光谱特征相似的地物，比如绝大部分建筑、道路以及裸地等，有效地反映出了低植被覆盖区域滑坡引起的NDVI的变化。该方法的识别步骤更为简洁，对遥感数据要求更低，效率更高，并充分利用了遥感影像丰富的几何信息，在一定程度上保证了滑坡识别结果的可靠性。

然而如果缺失同一季节的遥感影像，该方法可能会将森林、田野等的自然变化误判为滑坡预选区，需进一步通过获得多源遥感影像或量化植被生长速率等方法来解决。