同震崩塌滑坡的光学遥感影像多特征融合解译方法

2022-12-30韩征方振雄傅邦杰吴斌晖李艳鸽李常丽陈光齐

中国地质灾害与防治学报 2022年6期

韩征，方振雄，傅邦杰，吴斌晖，李艳鸽,3，李常丽，陈光齐

（1. 中南大学土木工程学院，湖南长沙 410075；2. 轨道交通工程结构防灾减灾湖南省重点实验室，湖南长沙 410075；3. 重载铁路工程结构教育部重点实验室，湖南长沙 410075；4. 九州大学工学府，福冈 819-0372，日本）

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 方法流程示意图Fig.1 Schematic diagram of the proposed method

1.1 局部阈值二值化

同震崩塌滑坡识别的本质是对光学影像中潜在的物源区域进行分割操作。针对同一幅影像上不同部分具有不同灰度的现象，利用单一的全局阈值无法对影像细节做出较好地分割。因此，本研究采用局部阈值二值化方法对影像进行初步分割，来进行同震崩塌滑坡识别。

局部阈值二值化方法首先利用蒙特卡洛迭代来对影像进行随机分块，对分块后的局部影像采用经典大津（Otsu）二值化方法来区分识别目标的前景与背景。根据蒙特卡洛迭代结果的置信度来判断像元是否属于崩塌体。基于全色影像与地形数据融合的崩塌体识别方法具体操作思路如下：

（1）对H×W大小的影像进行蒙特卡洛随机分块，共设置N步蒙特卡洛计算步，H和W分别指代影像随机分块的高度和宽度。

（2）通过每次计算步获得的蒙特卡洛随机分块图像(大小为H×W)。

（3）对上述图像块进行经典大津（Otsu）二值化阈值处理，获得每个像元在第k次蒙特卡洛计算步的二值化分割结果fk(x,y)，其中k

（4）对每次蒙特卡洛计算步获得的像元分割结果进行二值化概率分布统计。

（5）设置置信度TP为概率阈值，大于TP为前景目标，小于TP则为背景杂波。我们在研究[18]中发现，当TP=0.80具有较好的效果，即假设100 次蒙特卡洛迭代计算中，某像元有80 次以上被判识为目标前景，则该像元在最终结果中被判为同震崩塌滑坡。

基于局部阈值二值化的全色影像同震崩塌滑坡识别方法的具体思路如图2 所示。

图2 全色影像局部阈值二值化初步检测流程图Fig.2 Flow chart of collapse identification for panchromatic image based on local threshold binarization

1.2 多特征策略融合

同震崩塌滑坡的影像灰度特征与裸露的土壤极为相似，但裸露的土壤不一定是同震崩塌滑坡，如裸露的河漫滩、开垦的斜坡地、泥土道路等地物。这些地物很难有效地与同震崩塌滑坡进行区分，致使我们在研究[17]中提出的方法在同震崩塌滑坡解译结果上假阳率偏高，导致检测的准确性受到影响。

为了排除这些误判为同震崩塌滑坡的假阳性地物干扰，研究首先引入地形坡度策略对初步识别结果进行筛选。但是，受数字高程模型DEM 精度的限制，地形坡度策略虽然能够去除少部分地形的干扰，但在河漫滩、坡脚等处仍然可能存在误判。因此，在地形坡度策略的基础上，研究引入了主轴特征对识别方法加以补充。最后，考虑到植被区对最终结果的影响，研究尝试结合以下特征策略（地形坡度、主轴轴向、NDVI指数）对初步的检测结果进行二次筛选，尝试进一步降低同震崩塌滑坡识别假阳率。

1.2.1 地形坡度特征

地形坡度表征每个像元到与其相邻的像元方向上的最大高程变化率。通常来说，坡度值越小，地势越平坦；坡度值越大，地势越陡峭。考虑到同震崩塌滑坡的发育特征，即通常发育在坡度较陡的山体，我们沿用已有研究[18]所建议的5°坡度作为阈值分界，对初步检测结果进行过滤，去掉坡度较缓的区域。

1.2.2 主轴轴向特征

同震崩塌滑坡的滑移方向多为由高到低的滑流，因此，其主轴轴向应当具有明显的方向性。同震崩塌滑坡受重力因素影响，其发生移动一般表现为从山坡到山谷的运动，多为从上到下的直线式滑行，通常轴向特征表现为近竖直向。沿山道路的主轴轴向则多以横向为主，为此通过目标体的主轴轴向可以有效区分同震崩塌滑坡与道路。此外，河岸边裸露的河漫滩的主轴方向也多为横向。因此，可以根据地物的主轴轴向特征降低裸露泥土道路、河漫滩对同震崩塌滑坡检测识别检测的影响。

目标区域主轴轴向的计算过程如图3 所示。

图3 计算主轴特征示意图Fig.3 Schematic diagram of calculating spindle characteristics

设影像中的目标区域为R，一共有N个像素，(x,y)处的像素值为f(x,y)。则目标的质心定义为目标像素坐标的均值(,)：

其中，i为x轴方向的阶数，j为y轴方向的阶数。

目标的惯性矩定义为：

其中，惯性主轴定义惯性矩I(θ)为最小的方向，对I(θ)求导：

令I′(θ)=0可得：

可以发现，式（2）—（7）有两个解

为了区别最大值和最小值，求I(θ)的二阶导数，得

将θ1和θ2分别代入式（9）中进行计算，当I′′(θ)>0时，θi即为所求的惯性主轴方向，主轴轴向示意图如图4 所示。

图4 主轴轴向特征示意图Fig.4 Schematic diagram of spindle axial characteristics

利用主轴特征对解译结果进行修正的步骤中，关键在于如何得到主轴的位置以及确定主轴端点的坐标和高程。文中方法具体如下：首先计算得到目标区域的主轴轴向[28−30]，再对主轴进行遍历即可解算出主轴线与图形的交点，此交点即为主轴的端点，进而获得主轴线在图像中的坐标位置、求解出图形的长短主轴、对应图形顶点以及轴向。再与高程数据图融合即可求得主轴端点的坐标以及高程，计算得到目标高度差、纵横比、方向走势等信息参数。通过规定目标区域长短轴比值和高差比的阈值，对同震崩塌滑坡和长条形的地物进行区分，从而解译出同震崩塌滑坡的位置，其中设定的阈值针对不同研究区域应单独考虑确定阈值。

1.2.3NDVI指数特征

通常，当山体在地震作用下形成同震崩塌滑坡时，由于山体瞬时发生整体失稳，导致失稳区域地表的植被会被大面积破坏，使得植被覆盖度显著降低，与周边环境形成明显差异。同震崩塌滑坡区域与周边环境在植被稀疏覆盖度的显著差异，可以从一定程度上反映出同震崩塌滑坡的分布情况[31−32]。由于植被对近红外波段具有较强的反射性，而在红光波段则具有较强的吸收性，因此，可以利用该特性来区分植被区与非植被区。

目前，归一化植被指数NDVI（Normalized Difference Vegetation Index）是衡量植被稀疏覆盖度常用的度量方式，还可以反映出植物冠层的背景情况[33−35]，如土壤潮湿地面、雪、岩石等。其值表征可见光红光波段与近红外波段的反射率差值与二者之和的比值：

式中，ρnir表示近红外波段反射率，ρ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

2.3 解译结果分析

本研究以鲁甸地震为实验案例验证局部阈值二值化结合特征策略的同震崩塌滑坡解译方法的可行性。图5 为鲁甸县龙头山镇区域地理位图。

图5 鲁甸县龙头山镇影像（影像数据来源于高分一号卫星2014 年8 月20 日影像数据）Fig.5 Image of Longtoushan Town,Ludian County (Image from GF-1 dataset,on August 20,2014)

首先，运用局部阈值二值化分割方法对该研究区域全色影像进行50 次蒙特卡洛步计算，设定置信概率TP=0.80，获得该研究区域的局部二值化分割结果图，如图6 所示。经过局部阈值二值化分割之后，高分一号全色影像被分为前景和背景两部分。从分割结果来看，前景中包含较多细小斑块，且存在大量的假阳性崩塌解译结果，与实际情况不符。究其原因，可能是该研究区的全色影像数据中地物整体偏暗，地物间色差并不明显，局部阈值二值化方法不能有效的区分出同震崩塌滑坡目标与周围背景，导致检测结果中存在着大量的杂波噪声。

图6 局部阈值二值化分割结果Fig.6 Local threshold binarization segmentation results

因此，本研究在初步分割结果基础上引入多特征策略融合方法。首先结合研究区DEM 数据引入区域坡度信息，获得融合坡度信息的改进二值化局部阈值分割的同震崩塌滑坡识别结果，如图7 所示。

图7 融合坡度信息的判定结果Fig.7 Judgment result of fusion slope information

从图中可知，引入坡度信息之后，去除了少量地势平坦区域的虚假杂波结果。通过对比分析，前景中仍然包含有长条形地物和植被区导致的假阳性结果，整体识别效果有待进一步提升。

以融合坡度信息的局部阈值二值化分割结果为基础，继续引入主轴轴向特征对同震崩塌滑坡进行识别筛选，获得结果如图8 所示。从图中可看出，基于主轴轴向特征的同震崩塌滑坡识别方法进一步去除了细小斑点噪声，一定程度上抑制了长条形地物和细小斑点噪声对识别结果的影响。但本质上依然还存在高植被覆盖区域导致的误报率偏高问题。为了解决这个问题，本研究引入NDVI指数特征对上述结果进一步筛选。以[−0.1,0.1]为NDVI指数阈值，对上述判定结果进行修正，获得融合NDVI指数特征的同震崩塌滑坡分割效果图，如图9 所示。

图8 融合主轴轴向特征的判定结果Fig.8 Judgment result of fusion spindle axial feature

图9 融合NDVI 指数特征的判定结果Fig.9 Judgment result of fusion of NDVI index features

观察以上的判定结果可得，局部阈值二值化的结果经过融合地形坡度特征，主轴特征以及NDVI指数特征判定后，图像的干扰噪声明显减少。其中融合坡度特征，主要去除坡度低于5%的缓坡区域干扰；针对崩塌体的主轴特征为纵向，而沿山道路、裸露河漫滩等干扰地物主轴特征多为横向的特点，引入主轴特征，实现了排除横向主轴地物干扰的目的；针对高植被覆盖地区发生山体震裂崩塌概率相对较小的特点，引入NDVI指数筛选，主要去除此类高植被覆盖区域对识别结果的影响。

2.4 精度验证

为验证同震崩塌滑坡的识别结果，对光学遥感影像进行了人工目视解译。同震崩塌滑坡在光学遥感影像上具有与非灾害地形不同的影像特征，可以从三个方面进行判断：颜色、色调特征，同震崩塌滑坡会破坏原始地表覆盖，在遥感图像上表现为地表覆盖的颜色连续性发生变化，对于新发生的崩塌山体，山体表面无植被覆盖且物质构成疏松，地表反射率较强，在影像上的色调较浅；图像纹理方面，纹理是图像上色调的变化频率，崩塌体不同部位的形态、色调都有所差异，所以崩塌体在纹理上呈现出无规则状态；形状特征，在遥感图像上，崩塌体通常表现为双沟同源、圈椅、椭圆等特殊的平面形态。

基于以上特征，并结合地形的高程数据，我们对检测区域进行了人工目视解译，同时为保证目视解译结果的准确性，文中以研究区发生震裂崩塌灾害前的影像数据（图10）作为参照进行目视解译。共识别出该区域内16 个同震崩塌滑坡，识别结果如图11 中的红色区域所示。其中，标记为蓝色的3 处道路沿线的小型岩土裸露区域由于其遥感影像特征与同震滑坡的特征相似，无法确定其是同震崩塌滑坡或是施工产生的挖方边坡，且由于时间较久无法现场核实和筛查。因此，为了更合理地说明文中提出方法的有效性，文中将这三处疑似区域作为假阳性地物进行考虑。

图10 研究区震前影像（影像数据来源于WorldView-2 卫星2011 年12 月6 日影像数据）Fig.10 Pre-earthquake images in the study area (Image from WorldView-2 dataset,on December 6,2011)

图11 目视解译结果示意图（红色区域为同震崩塌滑坡，蓝色区域为挖方边坡等疑似地物）Fig.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]指出汶川地震诱发滑坡具有“背坡面效应”，即在与发震断裂带近于垂直的沟谷斜坡中,在地震波传播的背坡面一侧的滑坡发育密度明显大于迎坡面一侧。这些前期研究所得到的同震崩塌滑坡分布的规律，也可为遥感影像解译结果的进一步修正提供支撑。