中国是世界上地震活动最频繁和受损最严重的国家之一，如2008年“5·12”汶川特大地震、2010年青海玉树地震、2013年四川雅安地震、2017年四川九寨沟地震等都造成了十分严重的生命财产损失^[1-2]。而在地形地貌复杂的山区，地震滑坡是地震引发的多种次生地质灾害中最为常见且破坏力最强的灾害^[3]，其造成的生命财产损失可占一次地震的一半以上^[4]。因此准确快速地获取地震滑坡灾害的范围、位置等信息，对开展地震应急救援和灾后重建等工作具有重要意义。

目前，地震滑坡信息主要通过人工现场勘查及遥感影像目视解译等方法来获取^[5]，然而现场勘查常因震后地面交通阻塞受到限制，需要耗费大量的人力物力，效率低耗时长，而目视解译方法则过于依赖专家经验，工作量过于庞大，无法满足地震应急救援的工作要求^[6]。近年来，随着高分辨率遥感影像的不断发展，越来越多的自动/半自动提取算法如变化检测方法、面向对象提取方法等被广泛应用在滑坡信息的提取识别中^[7-11]。Zhao等^[12]针对尼泊尔地震，利用基于像元的变化检测方法对Landsat-8数据进行滑坡体提取，并取得精度较高的提取结果。Li等^[13]基于航空遥感影像结合变化检测和水平集算法分别提出基于边缘的ELSE和基于区域的RLSE半自动化滑坡识别算法，RLSE提取结果较好，准确率达90%，漏检率约20%。Martha等^[14]综合利用滑坡体的光谱、空间和形状特征，使用面向对象的算法实现对滑坡体的半自动化提取。闫琦^[15]结合超像素分割和显著性检测提出一种基于面向对象的滑坡信息提取算法，结果表明该算法具有较高的召回率和准确率，但对于面积大、信息量多的影像，检测性能相对不足。然而，基于像元的变化检测虽然能够快速提取变化信息，但是其没有考虑到几何、纹理等特征信息，从而引起很多非滑坡体地物的误判；相反，面向对象的滑坡提取算法虽然充分利用遥感影像丰富的地物信息，提高了滑坡识别的精度，但是由于地物信息繁乱复杂，选取合适的特征以及分割参数的选取都需要人工干预，自动化程度较低，且时效性较差。

因此，本文基于变化检测和面向对象提取算法，深入分析地震滑坡的遥感影像特征，研究一种结合变化检测和面向对象算法的地震滑坡信息快速提取方法。该方法在基于像元的变化检测基础上，结合滑坡的纹理、几何和空间结构特征，筛选误判像元，并利用改进的区域生长算法优化滑坡体提取范围，从而实现地震滑坡信息的高效率、高自动化提取，为震后应急救援提供信息技术支撑。

1 研究区与数据 1.1 研究区域概况

本文选取2017年8月8日九寨沟7.0级地震(震中北纬33.2°，东经103.82°)的核心地区九寨沟县漳扎镇(北纬32°54′~33°25′，东经103°38′~104°41′)为研究区(图 1)。漳扎镇位于四川省北部，隶属于四川省阿坝藏族羌族自治州九寨沟县，地处青藏高原和四川盆地过渡地带，全镇面积1306km²，地势南高北低，地形以山地为主，平均海拔约为3400m。研究区位于塔藏断裂和虎牙断裂附近^[16]，是本次地震受灾最严重的地区之一，同时漳扎镇辖区内还包含有著名的九寨沟风景区，也受到地震的严重破坏^[17-18]。本次地震在研究区内引发了大量的地震滑坡，从其中选取2处地震滑坡分布最为密集的地区作为本次实验样区，如图 1所示。其中，样区1为九寨天堂洲际大饭店北部大型滑坡体；样区2为九寨沟景区内熊猫海附近滑坡群。

1.2 研究数据

本文选取研究区震前2017年7月29日和震后2017年9月7日的Sentinel-2A卫星遥感影像作为实验数据。Sentinel-2A卫星是欧空局(https://sentinel.esa.int/)对地观测任务计划的一部分，携带的多光谱成像仪，覆盖包括可见光、近红外等13个光谱波段，同时相较于同类型的Landsat-8卫星，Sentinel-2A卫星有着更高的空间分辨率(10m)，以及更短的重返周期(5d)，这使得Sentinel-2A卫星能够更好地完成地震滑坡信息快速提取的任务。

2 研究方法

本文针对地震滑坡信息提取的主要方法和流程如图 2所示。首先对震前震后遥感影像数据进行几何校正、辐射校正等数据预处理，并通过自动云掩膜提取算法得到去云的遥感影像；接着利用投票法对3种基于像元的变化检测阈值结果进行整合，得到滑坡体的初始掩膜；然后基于多特征参数，对初始掩膜进行逐层修正，剔除非滑坡体，再利用改进的区域生长算法对滑坡边界进行优化，最后通过形态学运算得到更加精确的滑坡范围。同时研究以变化向量分析法(change vector analysis，简称CVA)作为对比试验方法，将两种方法得到的结果与目视解译结果进行对比分析，从而对地震滑坡信息快速提取方法的提取精度做出评估。

2.1 数据预处理

本文下载的Sentinel-2A遥感影像数据均为L1C级产品，已经经过正射校正和几何精校正，几何配准误差小于0.5个像元，满足变化检测算法所需的配准精度要求^[19]，所以只利用ENVI5.4对遥感影像进行辐射定标、FLAASH大气校正、图像镶嵌、图像裁剪等预处理后得到10m可见光和近红外影像。云掩膜自动提取算法是依据波长越短的电磁波越容易受到大气中的水汽和气溶胶的散射作用影响，所以云在波长较短的波段中反射率大于波长较长波段这一现象^[20-21]，选取Sentinel-2A中波长最短的波段Band1-Coastal aerosol(60m)进行实验分析。在该波段中，云层的灰度值很高，与非云区域的差异度比其他波段更为突出，故本研究利用OTSU自动阈值分割法, 得到基础云掩膜M₁。但是考虑到Band1空间分辨率仅为60m，生成的掩膜对于后续研究中使用的10m可见光和近红外波段的误差较大，所以引入只利用可见光波段计算的HOT指数^[22]，其计算公式为

$ \mathrm{HOT}=B-0.5 R, $

(1)

式中：B为蓝波段(10m)辐射值，R为红波段(10m)辐射值。但是HOT指数缺陷明显：蓝色建筑物会被识别出来，同时会对一些高亮度无云地区产生高估^[23]。所以本文利用改进的区域生长算法(RGA，在2.4节中详细介绍)，将云掩膜M₁与HOT指数结合起来得到更加精确的云掩膜。

2.2 变化检测提取

基于像元的变化检测方法是直接对遥感影像灰度值进行分析处理，该方法容易实现，可以更好地保留遥感影像的真实变化信息，快速有效地提取地震滑坡信息。因此本文选取最为实用的3种变化检测方法，即差值法、比值法和变化向量分析法对震前震后遥感影像进行变化信息提取。差值法和比值法在实现原理上十分相似，通过将配准后的两时相影像上各波段的对应像元值相减或相除，构造差异影像，并通过OTSU自动阈值算法计算出变化阈值，来提取变化区域，使用两种方法分别在红绿蓝三波段进行变化提取，共得到6个滑坡体掩膜。变化向量分析法是用一个空间向量表示遥感影像的多波段变化，在每个像元位置得到一个同时具有变化方向和变化大小的变化向量，根据比较变化向量的模和设定的阈值判断该像元点的变化信息，该方法能综合红绿蓝三波段信息，只得到1个滑坡体掩膜。而后在得到的7个滑坡体掩膜的基础上，引入基于混合思想的变化检测方法^[24]，将目标区域分别使用多种变化检测算法，利用投票法综合分析7种方法的阈值分割结果，将所有结果的缺点最小化，从而得到准确率较高的滑坡体掩膜。

2.3 逐层特征信息筛选

目前国内外已有很多学者在面向对象的滑坡体提取研究中提出滑坡在遥感影像上的特征^[25-29]，总结归纳为：在光谱特征上，滑坡体比周围的植被和陈旧裸土的亮度相对较高，这是变化检测算法能够提取滑坡信息的最主要原因；在形状特征上，震后滑坡多呈现为舌状活片状，也有椭圆形、梨形、马蹄形等，普遍没有统一的形态，而大部分人工地物如建筑物、耕地、道路等经常呈现为较为规则的形状；在纹理特征上，震后滑坡整体纹理较为细腻，而与之易混淆的耕地则较为粗糙，所以可以辅助形状特征，有效地排除耕地的干扰。

依据上述总结的地震滑坡的特征信息，对投票法得到的变化检测结果逐层分析，剔除非滑坡体，从而得到最终的滑坡提取掩膜。算法具体包括以下几个步骤：

1) 将投票法综合分析得到的滑坡体掩膜作为输入T₁，震前震后多光谱影像作为辅助数据。

2) 选取多个光谱特征参数，如归一化植被指数NDVI、归一化水体指数NDWI等，利用光谱指数对震区进行初步筛选，得到结果T₂；

3) 对T₂做4邻域的连通区域划分，并计算所有连通区域的形状参数，同时基于震后影像和连通区域计算纹理参数，利用Otsu自动阈值确定算法，得到各参数的最佳阈值；

4) 利用形状参数阈值，排除T₂中形状规则的建筑物、道路，得到结果T₃；

5) 利用纹理参数阈值，排除T₃中耕地，最终得到更加准确的滑坡体掩膜。

由于逐层特征信息筛选是建立在已得到的变化检测的结果上设计的，所以在此只选取针对性强、效果显著的特征参数进行干扰排除。光谱特征参数采用归一化植被指数NDVI和归一化水体指数NDWI，公式如下：

$ {\rm{NDVI}} = \frac{{{\rm{NIR}} - {\rm{R}}}}{{{\rm{NIR}} + {\rm{R}}}}, $

(2)

$ {\rm{NDWI}} = \frac{{{\rm{G}} - {\rm{NIR}}}}{{{\rm{G}} + {\rm{NIR}}}}, $

(3)

式中：NIR为近红外波段反射率，G为绿波段反射率，R为红波段反射率。这2个指数被广泛应用在众多植被和水体提取研究中，能够很好地增强植被和水体的信息，并弱化其他干扰噪声，从而将滑坡体掩膜中的植被和水体误判剔除。形状特征参数采用长宽比r和紧致度c，公式如下：

$ r=\frac{l}{w}, $

(4)

$ c=\frac{l \times w}{a}, $

(5)

式中：l为区域外接矩形的长度，w为区域外接矩形的宽度，a为区域面积；r为区域外接矩形的长宽比，能够很好地表征区域的线性程度，从而排除道路信息的误判；c为区域最小外接矩形面积与区域面积的比值，比值越接近1，说明该区域的形状越接近矩形。由于人工地物往往较为规则，因此，本文使用矩形度这一特征参数去掉人工地物等错误信息。纹理特征参数采用区域方差来描述，公式如下

$ \sigma^{2}=(N-1)^{-1} \cdot \sum\limits_{i=1}^{N}\left(P_{i}-\bar{P}\right)^{2}, $

(6)

式中：N为区域中的像元个数，P_i是第i个像元的反射率，$ \bar P$是区域内像元的平均反射率。

由于经过上述方法筛选剔除，只剩余部分耕地作为误判地物，而耕地区域与滑坡区域的纹理差异较大，用区域方差完全可以区分，故没有选取计算更为复杂、耗时更长的灰度共生矩阵作为纹理特征参数。

2.4 改进的区域生长算法

为了进一步完善滑坡范围，并排除一些误判的薄云干扰，本文使用改进的区域生长算法(regional growth algorithm, 简称RGA)对特征信息筛选后的结果进行优化^[10]。RGA是一种应用广泛的基于区域的图像分割方法，而在本文中，通过对RGA的改进，使其能够结合背景影像和滑坡体掩膜优化滑坡提取范围，具体流程如图 3(a)所示。首先，规定初始种子点是震后滑坡掩膜中各个连通区域的第一个像元。其次，将掩膜图像中孤立的像元判定为噪声直接剔除。然后，将变化向量法提取的变化检测结果与滑坡掩膜叠加在一起，并使用变化向量检测结果中的差异值建立每一个初始生长区域对应的生长标准，其置信区间为$\left[m_{i}-\sigma_{i}, m_{i}+\sigma_{i}\right] $，其中，$ {{m_i}}$与$ {{\sigma _i}}$分别是第i个区域的均值与标准偏差。区域生长算法应用时像元的生长过程如图 3(b)所示，相同的颜色表示该区域内像素的灰度值和属性值有着相似的特征。如果一个临近种子点的像元值处在置信区间之内，那么这个像素就被标记为滑坡，同时如果该像素本身已被判定为云掩膜区域，则该区域会被标记为云从而被剔除。以此迭代，当没有像元被发现并满足这种条件时迭代终止。

2.5 算法精度评估

为评估快速提取算法得到的滑坡掩膜，本文将基于Google Earth、GF-1和GF2等遥感影像的人工目视解译提取的滑坡范围作为验证数据，并利用混淆矩阵分别计算准确率A(accuracy)、漏检率M(miss rate)、总体精度OA(overall accuracy)以及Kappa系数等精度评价指标，公式如下：

$ A = \frac{{{\rm{TP}}}}{{{\rm{TP}} + {\rm{FN}}}}, $

(7)

$ M = \frac{{{\rm{FP}}}}{{{\rm{TP}} + {\rm{FP}}}}, $

(8)

$ \mathrm{OA}=\frac{\mathrm{TP}+\mathrm{TN}}{N}, N=\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}, $

(9)

$ \begin{array}{l} {\rm{Kappa }} = \frac{{{\rm{OA}} - {P_{\rm{e}}}}}{{1 - {P_{\rm{e}}}}}, {P_{\rm{e}}} = \\ \frac{{({\rm{TP}} + {\rm{FN}})({\rm{TP}} + {\rm{FT}}) + ({\rm{TN}} + {\rm{FT}})({\rm{TN}} + {\rm{FN}})}}{{N \times N}}, \end{array} $

(10)

式中:TP(ture positive)为把实际滑坡正确提取为滑坡的像元数；FN(false negative)为把实际滑坡错误提取为非滑坡的像元数；TN(ture negative)为把实际非滑坡正确提取为非滑坡的像元数；FP(false positive)为把实际非滑坡错误提取为滑坡的像元数。A代表提取到的地震滑坡结果准确的比率；M指提取结果中未提取出来的滑坡像元占所有实际发生滑坡像元的比率；OA是指所有提取正确像元之和与总像元数的比值；Kappa系数能够准确反映最终提取的精度。

3 实验结果 3.1 滑坡体提取结果

依据本文上述方法流程，对整个研究区的Sentinel-2A遥感影像进行地震滑坡信息的快速提取，并将本文方法提取结果与CVA的提取结果进行对比，研究区整体和样区1、2的滑坡提取结果分别如图 4~图 6所示。通过对比图 4(c)和图 4(d)发现，本文算法在研究区整体提取效果良好，大部分滑坡体都被准确识别出来。与CVA提取结果相比，本文算法提取滑坡范围与人工目视解译结果更为一致，尤其在大面积滑坡边缘处更为明显；而在厚云的边缘处，本文算法的误检像元更少，其主要原因是改进的区域生长算法不仅可以完善滑坡范围，还将与云掩膜有重叠的区域判定为薄云，从而从提取结果中剔除。而观察图 5(c)、5(d)和图 6(c)、6(d)样区1、2中的两种方法的提取结果对比，可以更加清晰地看出，相较CVA方法本算法对滑坡边界进行了优化，同时将大部分邻近的零碎斑块整合成一个完整的滑坡体，尽管也将少部分夹杂在滑坡区域间的非滑坡像元误判为滑坡，但整体来看提高了提取精度。

3.2 提取结果精度评估

为了更直观地对两种算法进行评估，本文算法与CVA算法的滑坡提取精度评估结果如表 1所示。对比表 1数据可以看出，在两个样区中，本文算法的提取准确率分别为86.65%和87.26%，均低于CVA的准确率93.83%和96.57%，其原因主要是上述的在进行区域生长和形态学运算时，将一些夹杂在滑坡区域间的非滑坡像元误判为滑坡，从而使得准确率有所降低，但幅度并不十分明显。在研究区整体中本文算法的准确率达到68.33%，高于CVA算法的43.14%，但仍然低于在样区1、2中的准确率，分析图 4(c)、4(d)发现原因为本文为了尽量减少提取结果的人为干预而采用的自动云提取算法，虽然将大部分云区域掩膜掉，但是仍有部分面积小的薄云很难自动识别出来，而云在遥感影像中的特征与滑坡体极为相似，经常被误识别为滑坡，所以导致整体的准确率偏低。再观查漏检率，本文算法基本保持在20%~25%之间，而CVA算法则在50%左右，证明本文算法在对滑坡范围的优化上有着明显的改进，但仍然有一些面积较小或亮度低的滑坡未被准确提取出来，其中面积小的滑坡在区域生长算法和形态学运算后很容易被当做噪声剔除，而亮度较低的滑坡在初始的变化检测提取时没有被作为变化像元检测出来。而观察整体精度发现，不论本文算法还是CVA在样区均基本在90%左右，而研究区整体则达到98%，这主要是由于未发生改变的像元和经过算法处理后剔除的非滑坡体变化像元较多，使得整体精度较高。最后观察表 1中两算法的Kappa系数可以看出，本文算法在样区1、2和研究区整体分别达到0.78、0.76、0.70，明显高于CVA的0.59、0.57、0.42，说明本文算法的提取精度处于较高水平，满足地震应急救援时的提取精度需求，可以准确地提取出地震滑坡信息。

3.3 算法效率评价

由于本研究是为满足地震应急需求所提出的滑坡提取算法，所以在保证较高提取精度的前提下，需要尽可能地合理优化算法流程，选取合适的特征参数，从而减少算法的运行时间，以精准快速地提取地震滑坡信息。本研究的各算法均在Python平台上实现，电脑配置为Intel Core i7-7700HQ@2.80 GHz，16 GB内存，Windows 10操作系统。实验结果表明，在处理4081像素×5637像素的10m分辨率遥感影像时，完成几何校正、辐射校正、图像配准、镶嵌裁剪等预处理过后，去云/阴影阶段耗时39s，基于投票法优化的变化检测提取方法耗时155s，基于区域生长算法和逐层特征分析的优化方法耗时158s，本文算法流程整体耗时不超过6min，仅比传统的CVA算法多不到3min，所以可以认为本文提出的地震滑坡信息提取算法基本满足地震应急的时效性需求。

4 结论

本文提出一种结合变化检测和面向对象算法的地震滑坡信息快速提取方法，通过分析基于像元的变化检测算法在滑坡提取中的特点，并总结归纳地震滑坡的光谱、形状、纹理等遥感影像特征，从而剔除非滑坡体变化的提取范围，最后利用改进的区域生长算法对滑坡边界进行优化，从而更加精确地提取滑坡信息。

本文以“8·8九寨沟地震”引发的地震滑坡为例，基于Sentinel-2A卫星遥感影像数据，利用本文提出的方法提取震区滑坡信息，精度评价结果表明算法提取整体精度均在90%以上，Kappa系数均超过0.7，研究区内大部分滑坡均被准确识别，同时提取的滑坡范围也与人工目视解译结果较为一致。实验结果证明本文提出的地震滑坡快速提取方法满足地震应急救援时精度需求，可以高效准确地提取出地震滑坡信息，从而为开展地震应急救援和灾后重建等工作提供空间信息支撑。

本文提出的地震滑坡信息提取方法从影像去云阶段开始，就已经尽量减少人为干预，提高算法自动化程度，但是在特征筛选阶段阈值的选取仍依靠个人经验，因此，在后续的研究中会进一步提高算法的自动化程度，从而更加简单高效地提取滑坡信息，以更好地满足地震应急救援的准确性和时效性的需求。