Interpretation method for regional co-seismic collapses based on multi-feature fusion of optical remote sensing
-
摘要: 同震崩塌滑坡的解译及定位是震区灾后恢复工作中需要重点解决的问题，如何在灾害快速、自动解译的基础上，不断提高解译精度，是目前同震崩塌滑坡解译的研究热点之一，也是促使地质灾害早期识别向智能化、科学化发展的必要前提。文章在团队前期所提出的遥感影像局部阈值二值化方法的基础上，针对同震崩塌滑坡解译结果假阳率偏高的问题，分析了假阳性地物的光学和几何特点，提出了融合目标区域光学影像灰度特征、区域坡度信息、NDVI指数特征及解译地物主轴特征的同震崩塌滑坡多特征融合解译方法。为验证所提出模型的准确性，以2014年云南鲁甸地震龙头山镇为研究区，利用震后获取的高分一号（GF-1）卫星影像数据及数字高程模型对该同震崩塌滑坡进行了解译识别，结果表明，文中提出的方法准确解译出了同震崩塌滑坡区域，并有效去除了假阳性地物干扰，提高了解译精度。Abstract: Interpretation of co-seismic collapse landslides is a key problem that needs to be solved in the post-disaster recovery work in earthquake areas. The issue regarding continuously improvement of interpretation accuracy for rapid and automatic interpretation of disasters is currently a hot topic, which is also a prerequisite to promote the development of early recognition of geological disasters towards intelligence and scientific. Based on the local threshold binarization method of remote sensing image proposed by the team in the early stage, this paper analyzes the optical and geometric characteristics of false positive features and proposes a fusion for the high false positive rate of the interpreted results of co-seismic mountain collapse. The multi-feature fusion interpretation method of the co-seismic mountain collapse with the gray feature of the optical image of the target area, the regional slope information, the NDVI feature and the interpretation of the main axis feature of the ground feature. In order to verify the accuracy of the proposed model, based on the 2014 Ludian earthquake in Yunnan, a case study was carried out in the Longtoushan town area. The Gaofen-1 (GF-1) satellite image data obtained after the earthquake and the digital elevation model were used for the earthquake in this area. The interpretation and recognition of the collapse of the cracked mountain shows that the method proposed in this paper accurately interprets the collapsed area of the cracked mountain body, effectively removes the false positive ground object interference, and improves the accuracy of interpretation.
-
0. 引 言
我国是世界上受地震灾害影响最为严重的国家之一。尤其是青藏高原、四川盆地、云贵高原等第一、第二级阶梯艰险山区，地震不仅能够在震时造成严重的破坏，同时也引发大量的同震崩塌滑坡，造成严重的破坏[1]。以2014年8月3日云南鲁甸地区6.5级地震为例，同震崩塌滑坡灾害所造成的死亡和失踪人数约占总人数的30%[2-4]。此外，同震崩塌滑坡形成的堰塞坝和松散堆积物也是导致溃坝洪水、次生泥石流等链生灾害的主要诱因[5-6]，使得主震后一段时间内地质灾害发灾趋势显著增强。因此，如何准确、快速地对同震崩塌滑坡进行识别和定位是灾后恢复与建设所需优先解决的问题之一。
在遥感技术广泛运用于地质灾害领域前，同震崩塌滑坡的调查通常以人工踏勘为主。虽然识别结果详细准确，但效率较低，且震后区域地形陡峭、人力难至、危险性较高，因此覆盖面和适用性受限[7-8]。随着遥感技术逐渐成熟，星载和机载的光学遥感影像等数据获取成本显著降低，基于光学遥感影像提取区域内同震崩塌滑坡、泥石流等灾点信息的技术已成为主要手段。例如，汶川地震之后，遥感手段就广泛用于同震地质灾害的解译中，为同震地质灾害的快速编目提供了有力支撑[9]。总体而言，现有的识别方法可以分为三类：（1）变化检测方法。这种方法需要有同一位置的两期或多期遥感数据，找出因同震崩塌滑坡引起的变化区域，数据源可以是二维光学影像，也可以是三维地形数据[10-11]。变化检测对新生崩塌体的应用效果较好，但是需要具有时间序列的遥感数据，通常此类数据较难获取；（2）机器学习方法。常用的有贝叶斯、逻辑回归、支持向量机和人工神经网络等算法[12-13]。这类方法在数据准备部分需要提取所使用数据的各类相关特征，然后使用各种分类器进行同震崩塌滑坡识别，其自动化程度相较于前述方法更高。但机器学习方法的特征选择和超参数调试工作量较大；（3）特征阈值方法。这种方法多使用基于像素的方法或基于对象的多尺度分割方法，对同震崩塌滑坡区域的光谱、纹理、地貌或地形等特征进行统计，设置一种或者多种阈值进行崩塌体识别[14-16]。其判断精度高，工作量相对较小，但是由于其判断标准为特定区域特定特征的统计值，目前，无论是规范还是已有研究都未形成统一标准。
目前基于光学遥感影像解译同震崩塌滑坡已逐渐发展成为一种较为可行的解决方案。我们在前期研究[17-19]中，针对传统二值化分割方法[20-21]在遥感影像解译中所存在的受高灰度地物干扰等问题，提出了基于蒙特卡洛迭代来对全色波段影像进行随机分块，并对局部影像采用大津（Otsu）方法进行二值化分割，从而区分震裂物源前景与其他地物背景。该方法能够克服目前基于全局阈值二值化分割方法在区域影像光照不均、山体阴影遮挡或复杂地物背景等因素对解译结果的干扰。但是，该方法依然存在误报率（或假阳率）较高的问题，即一些相似灰度的地物被误判为同震崩塌滑坡。我们在研究[19]中发现，以2008年汶川地震北川地区的解译结果为例，误判为同震崩塌滑坡的地物通常是一些裸露的河漫滩、开垦的斜坡等与同震崩塌滑坡具有相似的灰度特征的地物，尤其是在地形DEM模型分辨率较低时，无法利用我们所提出的光学影像与地形坡度的融合方法进行区分。因此，如何进一步降低解译结果中的假阳性问题，从而提高同震崩塌滑坡解译精度是目前需要研究与解决的关键问题。
围绕上述问题，本研究针对复杂艰险山区地质灾害区域遥感影像特点和已有解译方法中产生的假阳性地物几何特征进行了分析，进而提出了一种遥感影像多特征策略融合的同震崩塌滑坡解译方法。该方法在我们已提出的局部阈值二值化识别分析的方法基础[18]上，进一步引入归一化植被指数（NDVI）特征[22-25]和主轴特征[26-27]，通过融合地形坡度特征、主轴轴向特征、NDVI特征，来排除灰度特性与同震崩塌滑坡相似的平缓区域对解译结果的干扰、河岸边裸露的河漫滩等沿水系流向的长条形地物的干扰以及植被区对识别结果的干扰，从而实现高精度同震崩塌滑坡的检测解译识别。最后，以2014年云南鲁甸地震龙头山镇区域为研究区，利用震后获取的高分一号（GF-1）卫星影像数据及数字高程模型对该区域同震崩塌滑坡进行了解译识别，从而验证本文所提出的方法的准确性和可行性。
1. 方法与模型
本研究首先利用局部阈值二值化方法对影像同震崩塌滑坡进行初步解译，再针对识别结果出现的假阳性地物，通过坡度信息、主轴特征因子和NDVI指数进行二次分析与判断。通过目标区域的坡度筛选，排除灰度特性与同震崩塌滑坡相似的平缓区域对解译结果的干扰；辅以目标区域的主轴分析，排除长条形地物对解译结果的干扰；最后结合目标区域的NDVI指数特征，避免将高植被覆盖区域误检为同震崩塌滑坡，进而提高同震崩塌滑坡解译精度。方法具体操作流程见图1。
1.1 局部阈值二值化
同震崩塌滑坡识别的本质是对光学影像中潜在的物源区域进行分割操作。针对同一幅影像上不同部分具有不同灰度的现象，利用单一的全局阈值无法对影像细节做出较好地分割。因此，本研究采用局部阈值二值化方法对影像进行初步分割，来进行同震崩塌滑坡识别。
局部阈值二值化方法首先利用蒙特卡洛迭代来对影像进行随机分块，对分块后的局部影像采用经典大津（Otsu）二值化方法来区分识别目标的前景与背景。根据蒙特卡洛迭代结果的置信度来判断像元是否属于崩塌体。基于全色影像与地形数据融合的崩塌体识别方法具体操作思路如下：
（1）对
$H \times W$ 大小的影像进行蒙特卡洛随机分块，共设置N步蒙特卡洛计算步，H和W分别指代影像随机分块的高度和宽度。（2）通过每次计算步获得的蒙特卡洛随机分块图像(大小为
$H \times W$ )。（3）对上述图像块进行经典大津（Otsu）二值化阈值处理，获得每个像元在第
$k$ 次蒙特卡洛计算步的二值化分割结果${f_k}\left( {x,y} \right)$ ，其中$k < N$ 。（4）对每次蒙特卡洛计算步获得的像元分割结果进行二值化概率分布统计。
$$ p\left( {x,y} \right) = \frac{1}{N}\mathop \sum \limits_{k = 1}^N {f_k}\left( {x,y} \right) $$ (1) （5）设置置信度
${T_P}$ 为概率阈值，大于${T_P}$ 为前景目标，小于${T_P}$ 则为背景杂波。我们在研究[18] 中发现，当${T_P} = 0.80$ 具有较好的效果，即假设100次蒙特卡洛迭代计算中，某像元有80次以上被判识为目标前景，则该像元在最终结果中被判为同震崩塌滑坡。基于局部阈值二值化的全色影像同震崩塌滑坡识别方法的具体思路如图2所示。
![]()
图 2 全色影像局部阈值二值化初步检测流程图Figure 2. Flow chart of collapse identification for panchromatic image based on local threshold binarization1.2 多特征策略融合
同震崩塌滑坡的影像灰度特征与裸露的土壤极为相似，但裸露的土壤不一定是同震崩塌滑坡，如裸露的河漫滩、开垦的斜坡地、泥土道路等地物。这些地物很难有效地与同震崩塌滑坡进行区分，致使我们在研究[17]中提出的方法在同震崩塌滑坡解译结果上假阳率偏高，导致检测的准确性受到影响。
为了排除这些误判为同震崩塌滑坡的假阳性地物干扰，研究首先引入地形坡度策略对初步识别结果进行筛选。但是，受数字高程模型DEM精度的限制，地形坡度策略虽然能够去除少部分地形的干扰，但在河漫滩、坡脚等处仍然可能存在误判。因此，在地形坡度策略的基础上，研究引入了主轴特征对识别方法加以补充。最后，考虑到植被区对最终结果的影响，研究尝试结合以下特征策略（地形坡度、主轴轴向、NDVI指数）对初步的检测结果进行二次筛选，尝试进一步降低同震崩塌滑坡识别假阳率。
1.2.1 地形坡度特征
地形坡度表征每个像元到与其相邻的像元方向上的最大高程变化率。通常来说，坡度值越小，地势越平坦；坡度值越大，地势越陡峭。考虑到同震崩塌滑坡的发育特征，即通常发育在坡度较陡的山体，我们沿用已有研究[18]所建议的5°坡度作为阈值分界，对初步检测结果进行过滤，去掉坡度较缓的区域。
1.2.2 主轴轴向特征
同震崩塌滑坡的滑移方向多为由高到低的滑流，因此，其主轴轴向应当具有明显的方向性。同震崩塌滑坡受重力因素影响，其发生移动一般表现为从山坡到山谷的运动，多为从上到下的直线式滑行，通常轴向特征表现为近竖直向。沿山道路的主轴轴向则多以横向为主，为此通过目标体的主轴轴向可以有效区分同震崩塌滑坡与道路。此外，河岸边裸露的河漫滩的主轴方向也多为横向。因此，可以根据地物的主轴轴向特征降低裸露泥土道路、河漫滩对同震崩塌滑坡检测识别检测的影响。
目标区域主轴轴向的计算过程如图3所示。
设影像中的目标区域为R，一共有N个像素，
$\left( {x,y} \right)$ 处的像素值为$f\left( {x,y} \right)$ 。则目标的质心定义为目标像素坐标的均值$\left( {\bar x,\bar y} \right)$ ：$$ \bar x = \frac{1}{N}\mathop {\sum \sum }\limits_{\left( {x,y} \right) \in R} xf\left( {x,y} \right) $$ (2) $$ \bar y = \frac{1}{N}\mathop {\sum \sum }\limits_{\left( {x,y} \right) \in R} yf\left( {x,y} \right) $$ (3) 则中心矩
${m_{i,j}}$ 定义为：$$ {m_{i,j}} = \mathop {\sum \sum }\limits_{\left( {x,y} \right) \in R} {(x - \bar x)^i}{(y - \bar y)^j}f\left( {x,y} \right) $$ (4) 其中，i为x轴方向的阶数，j为y轴方向的阶数。
目标的惯性矩定义为：
$$ \begin{split} I\left( \theta \right) &= \mathop {\sum \sum }\limits_{\left( {x,y} \right) \in R} {\overline {AC} ^2}f\left( {x,y} \right) \\ &= \mathop {\sum \sum }\limits_{\left( {x,y} \right) \in R} {\left[ { {\left( {y - \bar y} \right)\cos \theta + } {\left( {x - \bar x} \right)\sin \theta } } \right]^2} \cdot f\left( {x,y} \right) \end{split} $$ (5) 其中，惯性主轴定义惯性矩I(θ)为最小的方向，对I(θ)求导：
$$ \begin{split} I'\left( \theta \right) &= \underset{\left(x,y\right)\in R}{\sum \sum }2\left[\left(y-\overline{y}\right)\mathrm{cos}\theta + \left(\overline{x}-x\right)\mathrm{sin}\theta \right]\cdot \left[ { - \left( {y - \bar y} \right)\sin \theta + \left( {\bar x - x} \right)\cos \theta } \right] \cdot f\left( {x,y} \right) \\ &= \mathop {\sum \sum }\limits_{\left( {x,y} \right) \in R} \left\{ {\left[ {{{\left( {x - \bar x} \right)}^2} - {{\left( {y - \bar y} \right)}^2}} \right]\sin 2\theta - } {2\left( {x - \bar x} \right)\left( {y - \bar y} \right)\cos 2\theta } \right\} \cdot f\left( {x,y} \right) = \left( {{m_{2,0}} - {m_{0,2}}} \right)\sin 2\theta - 2{m_{1,1}}\cos 2\theta \end{split} $$ (6) 令
$I'\left( \theta \right) = 0$ 可得：$$ {\rm{arctan}}2\theta = \frac{{2{m_{1,1}}}}{{{m_{2,0}} - {m_{0,2}}}} $$ (7) 可以发现，式（2）—（7）有两个解
$$ {\theta _1} = \frac{1}{2}{\rm{arctan}}\left[ {\frac{{2{m_{1,1}}}}{{{m_{2,0}} - {m_{0,2}}}}} \right] $$ (8) $$ {\theta _2} = \frac{1}{2}{\rm{arctan}}\left[ {\frac{{2{m_{1,1}}}}{{{m_{2,0}} - {m_{0,2}}}}} \right] + \frac{\pi }{2} $$ (9) 为了区别最大值和最小值，求I(θ)的二阶导数，得
$$ I''\left( \theta \right) = \left[ {2\left( {{m_{2,0}} - {m_{0,2}}} \right)\cos2\theta + } {4{m_{1,1}}\sin2\theta } \right] \cdot f\left( {x,y} \right) $$ (10) 将
${\theta _1}$ 和${\theta _2}$ 分别代入式（9）中进行计算，当$I''(\theta ) > 0$ 时，${\theta _i}$ 即为所求的惯性主轴方向，主轴轴向示意图如图4所示。利用主轴特征对解译结果进行修正的步骤中，关键在于如何得到主轴的位置以及确定主轴端点的坐标和高程。文中方法具体如下：首先计算得到目标区域的主轴轴向[28-30]，再对主轴进行遍历即可解算出主轴线与图形的交点，此交点即为主轴的端点，进而获得主轴线在图像中的坐标位置、求解出图形的长短主轴、对应图形顶点以及轴向。再与高程数据图融合即可求得主轴端点的坐标以及高程，计算得到目标高度差、纵横比、方向走势等信息参数。通过规定目标区域长短轴比值和高差比的阈值，对同震崩塌滑坡和长条形的地物进行区分，从而解译出同震崩塌滑坡的位置，其中设定的阈值针对不同研究区域应单独考虑确定阈值。
1.2.3 NDVI指数特征
通常，当山体在地震作用下形成同震崩塌滑坡时，由于山体瞬时发生整体失稳，导致失稳区域地表的植被会被大面积破坏，使得植被覆盖度显著降低，与周边环境形成明显差异。同震崩塌滑坡区域与周边环境在植被稀疏覆盖度的显著差异，可以从一定程度上反映出同震崩塌滑坡的分布情况[31-32]。由于植被对近红外波段具有较强的反射性，而在红光波段则具有较强的吸收性，因此，可以利用该特性来区分植被区与非植被区。
目前，归一化植被指数NDVI（Normalized Difference Vegetation Index）是衡量植被稀疏覆盖度常用的度量方式，还可以反映出植物冠层的背景情况[33-35]，如土壤潮湿地面、雪、岩石等。其值表征可见光红光波段与近红外波段的反射率差值与二者之和的比值：
$$ {{NDVI}} = \frac{{{\rho _{{\rm{nir}}}} - {\rho _{{\rm{red}}}}}}{{{\rho _{{\rm{nir}}}} + {\rho _{{\rm{red}}}}}} $$ (11) 式中，
${\rho _{{\rm{nir}}}}$ 表示近红外波段反射率，${\rho _{{\rm{red}}}}$ 表示红光波段反射率。一般来说，NDVI取值在−0.1～0.1的区间内。关于确定NDVI区分植被、土壤、岩石或雪的阈值，目前通常根据经验确定，张华等[36]为了反映NDVI的动态变化，按照等间距法，将NDVI进行分级，分别为：植被覆盖低取值在(0, 0.2]、中低(0.2, 0.4]、中(0.4, 0.6]、中高(0.6, 0.8]和高(0.8, 1]。统计表明，当NDVI处于[−1, 0)一般指示地面覆盖类型为云、水、雪等，这些地物通常对可见光表现出高反射；当NDVI接近0时表示 和 两者几乎相等，指示没有植被覆盖，地表被裸露的土壤或者岩石覆盖；当NDVI处于(0, 1]通常指示该地有植被，其植被覆盖度与该值呈现正相关[37]。本研究取NDVI数值在[−0.1, 0.1]区间内的区域为同震崩塌滑坡区域，从而降低植被覆盖高的区域被误检为同震崩塌滑坡区域的可能性。
2. 实例验证
2.1 研究区概况
以2014年8月云南鲁甸地震龙头山镇区域为研究案例，文中开展同震崩塌滑坡的解译识别方法的验证。该研究区位于云南省昭通市鲁甸县。该区域曾于北京时间2014年8月3日16时30分发生了6.5级地震，震源深度达12 km[38-39]。地震的震区属于高山峡谷地貌，由于正处在雨季，地震引发了严重的滑坡、泥石流、滚石以及堰塞湖等次生灾害，同震崩塌滑坡较为发育。
2.2 数据来源
区内地形环境复杂、地物类型多样、山体破坏显著，存在着多处局部型同震崩塌滑坡物源。原始数据包括：震后鲁甸地区高分一号卫星影像数据及鲁甸地区数字高程模型（http://glovis.usgs.gov/)，影像拍摄日期为2014年8月20日，本研究使用的原始数据参数信息如表1所示。
表 1 原始数据参数信息表Table 1. Information of the input data原始数据 谱段范围/μm 空间分辨率/m 高分一号（GF-1）
卫星影像数据0.45～0.90 2.00 0.45～0.52 8.00 0.52～0.59 0.63～0.69 0.77～0.89 数字高程模型DEM − 30.00 2.3 解译结果分析
本研究以鲁甸地震为实验案例验证局部阈值二值化结合特征策略的同震崩塌滑坡解译方法的可行性。图5为鲁甸县龙头山镇区域地理位图。
![]()
图 5 鲁甸县龙头山镇影像（影像数据来源于高分一号卫星2014年8月20日影像数据）Figure 5. Image of Longtoushan Town, Ludian County (Image from GF-1 dataset, on August 20, 2014)首先，运用局部阈值二值化分割方法对该研究区域全色影像进行50次蒙特卡洛步计算，设定置信概率
${T_P} = 0.80$ ，获得该研究区域的局部二值化分割结果图，如图6所示。经过局部阈值二值化分割之后，高分一号全色影像被分为前景和背景两部分。从分割结果来看，前景中包含较多细小斑块，且存在大量的假阳性崩塌解译结果，与实际情况不符。究其原因，可能是该研究区的全色影像数据中地物整体偏暗，地物间色差并不明显，局部阈值二值化方法不能有效的区分出同震崩塌滑坡目标与周围背景，导致检测结果中存在着大量的杂波噪声。因此，本研究在初步分割结果基础上引入多特征策略融合方法。首先结合研究区DEM数据引入区域坡度信息，获得融合坡度信息的改进二值化局部阈值分割的同震崩塌滑坡识别结果，如图7所示。
从图中可知，引入坡度信息之后，去除了少量地势平坦区域的虚假杂波结果。通过对比分析，前景中仍然包含有长条形地物和植被区导致的假阳性结果，整体识别效果有待进一步提升。
以融合坡度信息的局部阈值二值化分割结果为基础，继续引入主轴轴向特征对同震崩塌滑坡进行识别筛选，获得结果如图8所示。从图中可看出，基于主轴轴向特征的同震崩塌滑坡识别方法进一步去除了细小斑点噪声，一定程度上抑制了长条形地物和细小斑点噪声对识别结果的影响。但本质上依然还存在高植被覆盖区域导致的误报率偏高问题。为了解决这个问题，本研究引入NDVI指数特征对上述结果进一步筛选。以[−0.1, 0.1]为NDVI指数阈值，对上述判定结果进行修正，获得融合NDVI指数特征的同震崩塌滑坡分割效果图，如图9所示。
观察以上的判定结果可得，局部阈值二值化的结果经过融合地形坡度特征，主轴特征以及NDVI指数特征判定后，图像的干扰噪声明显减少。其中融合坡度特征，主要去除坡度低于5%的缓坡区域干扰；针对崩塌体的主轴特征为纵向，而沿山道路、裸露河漫滩等干扰地物主轴特征多为横向的特点，引入主轴特征，实现了排除横向主轴地物干扰的目的；针对高植被覆盖地区发生山体震裂崩塌概率相对较小的特点，引入NDVI指数筛选，主要去除此类高植被覆盖区域对识别结果的影响。
2.4 精度验证
为验证同震崩塌滑坡的识别结果，对光学遥感影像进行了人工目视解译。同震崩塌滑坡在光学遥感影像上具有与非灾害地形不同的影像特征，可以从三个方面进行判断：颜色、色调特征，同震崩塌滑坡会破坏原始地表覆盖，在遥感图像上表现为地表覆盖的颜色连续性发生变化，对于新发生的崩塌山体，山体表面无植被覆盖且物质构成疏松，地表反射率较强，在影像上的色调较浅；图像纹理方面，纹理是图像上色调的变化频率，崩塌体不同部位的形态、色调都有所差异，所以崩塌体在纹理上呈现出无规则状态；形状特征，在遥感图像上，崩塌体通常表现为双沟同源、圈椅、椭圆等特殊的平面形态。
基于以上特征，并结合地形的高程数据，我们对检测区域进行了人工目视解译，同时为保证目视解译结果的准确性，文中以研究区发生震裂崩塌灾害前的影像数据（图10）作为参照进行目视解译。共识别出该区域内16个同震崩塌滑坡，识别结果如图11中的红色区域所示。其中，标记为蓝色的3处道路沿线的小型岩土裸露区域由于其遥感影像特征与同震滑坡的特征相似，无法确定其是同震崩塌滑坡或是施工产生的挖方边坡，且由于时间较久无法现场核实和筛查。因此，为了更合理地说明文中提出方法的有效性，文中将这三处疑似区域作为假阳性地物进行考虑。
![]()
图 10 研究区震前影像（影像数据来源于WorldView-2卫星2011年12月6日影像数据）Figure 10. Pre-earthquake images in the study area (Image from WorldView-2 dataset, on December 6, 2011)![]()
图 11 目视解译结果示意图（红色区域为同震崩塌滑坡，蓝色区域为挖方边坡等疑似地物）Figure 11. Schematic diagram of visual interpretation results (pink polygons are the co-seismic collapses, red polygons are the suspected objects due to excavation)针对上节基于局部二值化结合多特征策略的同震崩塌滑坡识别结果与目视解译结果进行统计量化对比分析，由对比结果可知，基于局部二值化结合多特征策略（地形坡度特征、主轴轴向特征、NDVI特征）的方法准确的识别出了所有同震崩塌滑坡目标，去除了假阳性地物的干扰，提高了检测准确性，有效地提高了同震崩塌滑坡解译识别的精度。
由表2可知，基于全色影像局部阈值二值化和地形坡度特征的同震崩塌滑坡识别真阳率已达100%，而假阳率高达156.3%。主轴轴向特征的引入，进一步将同震崩塌滑坡假阳率降低至112.5%，在地形坡度特征基础上对同震崩塌滑坡误检数量进行了有效控制。NDVI指数特征的引入将同震崩塌滑坡假阳率从112.5%降低至23.1%。可以看出本方法最终的解译结果，无法识别出同期影像中存在的道路挖方边坡等疑似区域，有必要通过现场勘察进行进一步筛查。
表 2 检测解译效果统计Table 2. Detection and interpretation effect statistics可以发现，从所去除假阳性结果的数量上看，利用NDVI指数是更为有效的手段，表明相比地形坡度特征以及解译目标几何形态的主轴特征，本案例中同震崩塌滑坡在植被稀疏覆盖度上的差异更为显著，植被区对同震崩塌滑坡识别的影响较大。但对于其他区域而言，通过NDVI指数特征去除假阳性结果是否最为有效仍有待进一步研究。因此，在实际工作中，应该将地形坡度特征、主轴特征、NDVI特征进行综合考虑，从而尽可能地去除解译结果中假阳性地物干扰，对结果进行修正，提高同震崩塌滑坡的识别精度。
3. 讨论
（1）文中运用多特征融合策略来对同震崩塌滑坡识别结果开展进一步的判定分析，特征策略的获取精度可能对识别结果精度产生一定影响，主要表现为识别结果与真实同震崩塌滑坡边界不统一。因数据资源有限，本文采用的DEM数据分辨率为30 m，其分辨率相对于影像数据（分辨率2 m）较低且不匹配，故在系统计算主轴端点高程时，无法实现1 m级别的高程计算。采用更高精度DEM数据能够较为准确的描述特征位置的高程信息，或许能够进一步提高同震崩塌滑坡解译的精度。因此，在后续的工作研究中将开展尝试运用高精度的DEM数据来参与同震崩塌滑坡识别工作，进一步验证DEM精度对识别效果的影响。
（2）文中选取NDVI的阈值为[−0.1, 0.1]，此阈值是针对鲁甸山体震裂崩塌地区，经过实验比选以及参考相关文献最终确定，对于其他山体震裂崩塌区域并不具有准确适用性，但仍可作为参考，在此阈值基础上加以调整，确定不同区域的NDVI划分阈值。由于时间和资源的限制，本研究并未能够针对不同研究区域进行NDVI特征分析，后续可开展不同研究区域NDVI特征的对比分析。
（3）文中案例研究所采用的云南鲁甸县龙头山镇的同震崩塌滑坡影像，研究区实际覆盖区域面积10.5 km2，大型同震崩塌滑坡共16处。利用文中提出的光学遥感影像多特征融合解译方法能够有效识别该区域内的同震崩塌滑坡。而对于更广区域的遥感影像，考虑到其背景更为复杂，所包含的崩塌滑坡也不仅限于同震所产生的，因此本方法是否有效仍有待进一步验证。一种可行思路是将整幅高分一号影像分割为若干幅10 km2的子图像，并分别对每个子图像利用本文所述方法进行同震崩塌滑坡解译，但文中建议的NDVI指数阈值[−0.1, 0.1]是否适合各子图像，还需要进行研究和论证，这也是下一步研究中所需要重点解决的问题。
（4）研究主要从同震崩塌滑坡的发育地形、几何形态以及植被分布差异三个考量出发，选择地形坡度特征、主轴特征、NDVI指数特征来进行同震崩塌滑坡的遥感影像解译。除此之外，高山峡谷地区地震诱发崩塌滑坡的一些典型规律理论上也可以作为附加特征进行融合解译。例如，许强等[40-41]指出汶川地震诱发滑坡具有“背坡面效应”，即在与发震断裂带近于垂直的沟谷斜坡中,在地震波传播的背坡面一侧的滑坡发育密度明显大于迎坡面一侧。这些前期研究所得到的同震崩塌滑坡分布的规律，也可为遥感影像解译结果的进一步修正提供支撑。
4. 结论
文中针对目前同震崩塌滑坡的光学遥感影像识别方法中存在假阳率偏高的问题，开展了基于多源遥感数据融合的同震崩塌滑坡识别方法研究，相比于传统同震崩塌滑坡识别方法，有效提升了识别精度。具体结论如下：
（1）针对传统的全局阈值二值化导致的识别精度低、前景和背景区分不明显等问题，本研究介绍了基于全色影像局部阈值二值化的同震崩塌滑坡识别方法，利用了该方法对大尺度影像进行分割操作，避免了利用单一的全局阈值无法对影像细节做出较好的分割的问题。
（2）针对局部阈值二值化方法识别假阳率仍然较高等问题，本研究引入了多特征策略组（地形坡度特征、主轴轴向特征、NDVI指数特征）来进行同震崩塌滑坡的识别。通过多特征策略，有效地排除了灰度特性与同震崩塌滑坡相似的平缓区域、长条形地物以及植被覆盖区对识别结果的影响，进一步降低假阳率，提高识别准确性。
（3）文中研究依托2014年云南鲁甸地震龙头山镇区域，对同震崩塌滑坡识别方法进行实例例证。结果表明，基于蒙特卡洛迭代计算的局部阈值二值化结合多特征策略方法，能够有效避免传统二值化方法识别精度低、无法区分同震崩塌滑坡区与非同震崩塌滑坡区的问题，可以有效提升同震崩塌滑坡的识别精度。
-
![]()
图 2 全色影像局部阈值二值化初步检测流程图
Figure 2. Flow chart of collapse identification for panchromatic image based on local threshold binarization
![]()
图 5 鲁甸县龙头山镇影像（影像数据来源于高分一号卫星2014年8月20日影像数据）
Figure 5. Image of Longtoushan Town, Ludian County (Image from GF-1 dataset, on August 20, 2014)
![]()
图 10 研究区震前影像（影像数据来源于WorldView-2卫星2011年12月6日影像数据）
Figure 10. Pre-earthquake images in the study area (Image from WorldView-2 dataset, on December 6, 2011)
![]()
图 11 目视解译结果示意图（红色区域为同震崩塌滑坡，蓝色区域为挖方边坡等疑似地物）
Figure 11. Schematic diagram of visual interpretation results (pink polygons are the co-seismic collapses, red polygons are the suspected objects due to excavation)
表 1 原始数据参数信息表
Table 1 Information of the input data
原始数据 谱段范围/μm 空间分辨率/m 高分一号（GF-1）
卫星影像数据0.45～0.90 2.00 0.45～0.52 8.00 0.52～0.59 0.63～0.69 0.77～0.89 数字高程模型DEM − 30.00 
下载: 导出CSV
-
[1] 范一大,吴玮,王薇,等. 中国灾害遥感研究进展[J]. 遥感学报,2016,20(5):1170 − 1184. [FAN Yida,WU Wei,WANG Wei,et al. Research progress of disaster remote sensing in China[J]. Journal of Remote Sensing,2016,20(5):1170 − 1184. (in Chinese with English abstract) [2] 常昊,张吕. 云南鲁甸Ms6.5级地震震区滑坡易发性分析[J]. 中国地质灾害与防治学报,2017,28(2):38 − 48. [CHANG Hao,ZHANG Lyu. Analysis of Susceptibility causes of landslides triggered by earthquake in Ludian Ms6.5 earthquake region[J]. The Chinese Journal of Geological Hazard and Control,2017,28(2):38 − 48. (in Chinese with English abstract) [3] 曹颖,黄江培,钱佳威,等. 利用时移层析成像方法揭示与2014年云南鲁甸MS6.5地震有关的P波速度变化[J]. 地球物理学报,2021,64(5):1569 − 1584. [CAO Ying,HUANG Jiangpei,QIAN Jiawei,et al. Application of time-lapse seismic tomography based on double-difference tomography to reveal P wave velocity changes related to the 2014 Ludian MS6.5 earthquake[J]. Chinese Journal of Geophysics,2021,64(5):1569 − 1584. (in Chinese with English abstract) [4] 韩继冲,张朝,曹娟. 基于逻辑回归的地震滑坡易发性评价—以汶川地震、鲁甸地震为例[J]. 灾害学,2021,36(2):193 − 199. [HAN Jichong,ZHANG Zhao,CAO Juan. Assessing earthquake-induced landslide susceptibility based on logistic regression in 2008 Wenchuan earthquake and 2014 Ludian earthquake[J]. Journal of Catastrophology,2021,36(2):193 − 199. (in Chinese with English abstract) DOI: 10.3969/j.issn.1000-811X.2021.02.034 [5] 胡华,吴轩,张越. 基于模拟试验的强降雨条件下花岗岩残积土斜坡滑塌破坏机理分析[J]. 中国地质灾害与防治学报,2021,32(5):92 − 97. [HU Hua,WU Xuan,ZHANG Yue. Experimental study on slope collapse characteristics of granite residual soil slope under heavy rainfall[J]. The Chinese Journal of Geological Hazard and Control,2021,32(5):92 − 97. (in Chinese with English abstract) [6] 魏正发,张俊才,曹小岩,等. 青海西宁南北山滑坡、崩塌成因及影响分析[J]. 中国地质灾害与防治学报,2021,32(4):47 − 55. [WEI Zhengfa,ZHANG Juncai,CAO Xiaoyan,et al. Causes and influential factor analysis of landslides and rockfalls in north & south mountain areas of Xining City,Qinghai Province[J]. The Chinese Journal of Geological Hazard and Control,2021,32(4):47 − 55. (in Chinese with English abstract) [7] 龙玉洁,李为乐,黄润秋,等. 汶川地震震后10 a绵远河流域滑坡遥感自动提取与演化趋势分析[J]. 武汉大学学报·信息科学版,2020,45(11):1792 − 1800. [LONG Yujie,LI Weile,HUANG Runqiu,et al. Automatic extraction and evolution trend analysis of landslides in Mianyuan River basin in the 10 years after Wenchuan earthquake[J]. Geomatics and Information Science of Wuhan University,2020,45(11):1792 − 1800. (in Chinese with English abstract) [8] 孙国庆,陈方,于博,等. 2001—2017年尼泊尔中部地区滑坡变化及其影响因素[J]. 中国科学院大学学报,2020,37(3):308 − 316. [SUN Guoqing,CHEN Fang,YU Bo,et al. Landslide change and its influence factors in central Nepal from 2001 to 2017[J]. Journal of University of Chinese Academy of Sciences,2020,37(3):308 − 316. (in Chinese with English abstract) DOI: 10.7523/j.issn.2095-6134.2020.03.003 [9] 许冲,戴福初,陈剑,等. 汶川Ms8.0地震重灾区次生地质灾害遥感精细解译[J]. 遥感学报,2009,13(4):754 − 762. [XU Chong,DAI Fuchu,CHEN Jian,et al. Remote sensing fine interpretation of secondary geological disasters in the hardest hit areas of Wenchuan Ms8.0 earthquake[J]. Journal of Remote Sensing,2009,13(4):754 − 762. (in Chinese) DOI: 10.11834/jrs.20090416 [10] LU P,QIN Y Y,LI Z B,et al. Landslide mapping from multi-sensor data through improved change detection-based Markov random field[J]. Remote Sensing of Environment,2019,231:111235. DOI: 10.1016/j.rse.2019.111235
[11] LI Z B,SHI W Z,LU P,et al. Landslide mapping from aerial photographs using change detection-based Markov random field[J]. Remote Sensing of Environment,2016,187:76 − 90. DOI: 10.1016/j.rse.2016.10.008
[12] STUMPF A,KERLE N. Object-oriented mapping of landslides using Random Forests[J]. Remote Sensing of Environment,2011,115(10):2564 − 2577. DOI: 10.1016/j.rse.2011.05.013
[13] VAN DEN EECKHAUT M,KERLE N,POESEN J,et al. Object-oriented identification of forested landslides with derivatives of single pulse LiDAR data[J]. Geomorphology,2012,173/174:30 − 42. DOI: 10.1016/j.geomorph.2012.05.024
[14] MARTHA T R,KERLE N,JETTEN V,et al. Characterising spectral,spatial and morphometric properties of landslides for semi-automatic detection using object-oriented methods[J]. Geomorphology,2010,116(1/2):24 − 36.
[15] SUN W Y,TIAN Y S,MU X M,et al. Loess landslide inventory map based on GF-1 satellite imagery[J]. Remote Sensing,2017,9(4):314. DOI: 10.3390/rs9040314
[16] LESHCHINSKY B A,OLSEN M J,TANYU B F. Contour Connection Method for automated identification and classification of landslide deposits[J]. Computers & Geosciences,2015,74:27 − 38.
[17] LI Y G,CHEN G Q,HAN Z,et al. A hybrid automatic thresholding approach using panchromatic imagery for rapid mapping of landslides[J]. GIScience & Remote Sensing,2014,51(6):710 − 730.
[18] HAN Z,SU B,LI Y G,et al. An enhanced image binarization method incorporating with Monte-Carlo simulation[J]. Journal of Central South University,2019,26(6):1661 − 1671. DOI: 10.1007/s11771-019-4120-9
[19] HAN Z,LI Y G,DU Y F,et al. Noncontact detection of earthquake-induced landslides by an enhanced image binarization method incorporating with Monte-Carlo simulation[J]. Geomatics,Natural Hazards and Risk,2019,10(1):219 − 241. DOI: 10.1080/19475705.2018.1520745
[20] CASTELLANOS F J,GALLEGO A J,CALVO-ZARAGOZA J. Unsupervised neural domain adaptation for document image binarization[J]. Pattern Recognition,2021,119:108099. DOI: 10.1016/j.patcog.2021.108099
[21] XIONG W,ZHOU L,YUE L,et al. An enhanced binarization framework for degraded historical document images[J]. EURASIP Journal on Image and Video Processing,2021,2021(1):13. DOI: 10.1186/s13640-021-00556-4
[22] 皮新宇,曾永年,贺城墙. 融合多源遥感数据的高分辨率城市植被覆盖度估算[J]. 遥感学报,2021,25(6):1216 − 1226. [PI Xinyu,ZENG Yongnian,HE Chengqiang. High-resolution urban vegetation coverage estimation based on multi-source remote sensing data fusion[J]. National Remote Sensing Bulletin,2021,25(6):1216 − 1226. (in Chinese with English abstract) [23] ASHOK A,RANI H P,JAYAKUMAR K V. Monitoring of dynamic wetland changes using NDVI and NDWI based landsat imagery[J]. Remote Sensing Applications:Society and Environment,2021,23:100547. DOI: 10.1016/j.rsase.2021.100547
[24] 陈安,李景吉,黎文婷,等. 2001—2018年雅砻江流域植被NDVI时空动态及其对气候变化的响应[J]. 水土保持研究,2022,29(1):169 − 175. [CHEN An,LI Jingji,LI Wenting,et al. Spatiotemporal of NDVI in the Yalong River basin from 2001 to 2018 and its response to climate change[J]. Research of Soil and Water Conservation,2022,29(1):169 − 175. (in Chinese with English abstract) [25] SHEN J X,EVANS F. The potential of landsat NDVI sequences to explain wheat yield variation in fields in western Australia[J]. Remote Sensing,2021,13(11):2202. DOI: 10.3390/rs13112202
[26] 岳思聪,赵荣椿,王庆. 基于象素主轴方向灰度变化特征的特征点检测算法[J]. 西北工业大学学报,2008,26(2):162 − 167. [YUE Sicong,ZHAO Rongchun,WANG Qing. Feature point detection using intensity variations along pixel principal orientation axes[J]. Journal of Northwestern Polytechnical University,2008,26(2):162 − 167. (in Chinese with English abstract) DOI: 10.3969/j.issn.1000-2758.2008.02.006 [27] 杨雨奇,高晓光,冯晓毅,等. 基于主轴分析和团块特征提取的ISAR目标检测方法[J]. 西北工业大学学报,2010,28(5):689 − 694. [YANG Yuqi,GAO Xiaoguang,FENG Xiaoyi,et al. A new method for ISAR target detection based on chief axis analysis and block feature extraction[J]. Journal of Northwestern Polytechnical University,2010,28(5):689 − 694. (in Chinese with English abstract) DOI: 10.3969/j.issn.1000-2758.2010.05.010 [28] 郝军保,邵磊. 惯性主轴方向的最佳判别[J]. 连云港职业大学学报,1993,6(1):92 − 95. [HAO Junbao,SHAO Lei. The best discrimination of the direction of inertia spindle[J]. Journal of Lianyungang Technical College,1993,6(1):92 − 95. (in Chinese) [29] 蔡建乐. 用特征矩阵的伴随矩阵求解惯量主轴方向[J]. 大学物理,1995,14(9):21 − 22. [CAI Jianle. A calculation for the principal axes of inertia by adjoint matrix of eigen matrix[J]. College Physics,1995,14(9):21 − 22. (in Chinese with English abstract) [30] 罗斌. 灰度图象的惯性主轴方向特征[J]. 安徽大学学报(自然科学版),1998,22(4):40 − 42. [LUO Bin. Least inertia moment axis of grey scale image[J]. Journal of Anhui University (Natural Sciences),1998,22(4):40 − 42. (in Chinese with English abstract) [31] 周智勇. 基于Landsat遥感影像的围场县植被覆盖时空格局变化[J]. 水文地质工程地质,2020,47(6):81 − 90. [ZHOU Zhiyong. Change in temporal-spatial pattern of vegetation coverage in Weichang County based on Landsat remote sensing image[J]. Hydrogeology & Engineering Geology,2020,47(6):81 − 90. (in Chinese with English abstract) [32] 贺军亮,韦锐,李丽,等. 基于时间序列植被指数资料的承德市植被覆盖时空演变分析[J]. 水文地质工程地质,2020,47(6):91 − 98. [HE Junliang,WEI Rui,LI Li,et al. Temporal and spatial evolution of vegetation cover in Chengde based ontime series NDVI data[J]. Hydrogeology & Engineering Geology,2020,47(6):91 − 98. (in Chinese with English abstract) [33] 杜春雨,范文义. 叶面积指数与植被指数关系研究[J]. 林业勘查设计,2013(2):77 − 80. [DU Chunyu,FAN Wenyi. Research and analysis of the correlation between leaf area index and vegetation index[J]. Forest Investigation Design,2013(2):77 − 80. (in Chinese with English abstract) DOI: 10.3969/j.issn.1673-4505.2013.02.035 [34] 潘霞,高永,汪季,等. 植被指数遥感演化研究进展[J]. 北方园艺,2018(20):162 − 169. [PAN Xia,GAO Yong,WANG Ji,et al. Review on vegetation index using remote sensing evolution[J]. Northern Horticulture,2018(20):162 − 169. (in Chinese with English abstract) [35] 张慧,李平衡,周国模,等. 植被指数的地形效应研究进展[J]. 应用生态学报,2018,29(2):669 − 677. [ZHANG Hui,LI Pingheng,ZHOU Guomo,et al. Advances in the studies on topographic effects of vegetation indices[J]. Chinese Journal of Applied Ecology,2018,29(2):669 − 677. (in Chinese with English abstract) [36] 张华,李明,宋金岳,等. 基于地理探测器的祁连山国家公园植被NDVI变化驱动因素分析[J]. 生态学杂志,2021,40(8):2530 − 2540. [ZHANG Hua,LI Ming,SONG Jinyue,et al. Analysis of driving factors of vegetation NDVI change in Qilian Mountain National Park based on geographic detector[J]. Chinese Journal of Ecology,2021,40(8):2530 − 2540. (in Chinese with English abstract) [37] 武正丽. 2000~2012年祁连山植被覆盖变化及其对气候的响应研究[D]. 兰州: 西北师范大学, 2014 WU Zhengli. The research of the vegetation change and the sensitivity between NDVI and climatic factors in Qilian Mountains from2000to2012[D]. Lanzhou: Northwest Normal University, 2014. （in Chinese with English abstract）
[38] 周庆,吴果. 鲁甸6.5级地震崩滑地质灾害分布与成因探讨[J]. 地震地质,2015,37(1):269 − 278. [ZHOU Qing,WU Guo. Seismic landslides and seismogenic structure of the 2014 Ludian ms6.5 earthquake[J]. Seismology and Geology,2015,37(1):269 − 278. (in Chinese with English abstract) DOI: 10.3969/j.issn.0253-4967.2015.01.021 [39] 田颖颖,许冲,徐锡伟,等. 2014年鲁甸MS6.5地震震前与同震滑坡空间分布规律对比分析[J]. 地震地质,2015,37(1):291 − 306. [TIAN Yingying,XU Chong,XU Xiwei,et al. Spatial distribution analysis of coseismic and pre-earthquake landslides triggered by the 2014 Ludian ms6.5 earthquake[J]. Seismology and Geology,2015,37(1):291 − 306. (in Chinese with English abstract) DOI: 10.3969/j.issn.0253-4967.2015.01.023 [40] 许强, 李为乐. 汶川地震诱发滑坡方向效应研究[J]. 四川大学学报（工程科学版）, 2010, 42（增刊1）: 7 − 14 XU Qiang, LI Weile. Study on the direction effects of landslides triggered by Wenchuan earthquake[J]. Journal of Sichuan University （Engineering Science Edition）, 2010, 42（Sup 1）: 7 − 14. （in Chinese with English abstract）
[41] 许强,李为乐. 汶川地震诱发大型滑坡分布规律研究[J]. 工程地质学报,2010,18(6):818 − 826. [XU Qiang,LI Weile. Distribution of large-scale landslides induced by the Wenchuan earthquake[J]. Journal of Engineering Geology,2010,18(6):818 − 826. (in Chinese with English abstract) DOI: 10.3969/j.issn.1004-9665.2010.06.002 -
期刊类型引用(5)
1. 字汝芬，刘佳佳，王宇鸿，段平，李佳. 云南鲁甸M_s6.5级地震震后滑坡的时空分异特征. 中国地质灾害与防治学报. 2025(01): 73-83 . 
本站查看
2. 杨利亚，俞淑洋，杨静，殷非凡. 基于改进DDRNet网络的遥感影像山体滑坡识别. 北京测绘. 2024(03): 393-397 . 百度学术
3. 唐然，任穗川，范宣梅，许强. 大型红层缓倾岩层滑坡形成机制——以川北断渠滑坡为例. 成都理工大学学报(自然科学版). 2024(04): 649-662 . 百度学术
4. 裴鹏程，黄帅，袁静，张智康. 走滑断层作用下上覆土层的变形破坏机理. 中国地质灾害与防治学报. 2024(06): 115-127 . 本站查看
5. 范宣梅，王欣，戴岚欣，方成勇，邓宇，邹城彬，汤明高，魏振磊，窦向阳，张静，杨帆，陈兰，魏涛，杨银双，张欣欣，夏明垚，倪涛，唐小川，李为乐，戴可人，董秀军，许强. 2022年M_S6.8级泸定地震诱发地质灾害特征与空间分布规律研究. 工程地质学报. 2022(05): 1504-1516 . 百度学术
其他类型引用(1)
计量
- 文章访问数: 1731
- HTML全文浏览量: 1342
- PDF下载量: 285
- 被引次数: 6
邮件订阅
RSS