2017年8月8日，四川省阿坝州九寨沟县发生7.0级地震，震源深度为20 km，地震造成了严重的人员伤亡与财产损失。截止至8月13日，地震造成25人死亡，17多万人受灾，7万多间房屋不同程度受损(中国地震局)。此次地震诱发的滑坡数量多，大部分滑坡分布在九寨沟景区附近，为景区的恢复重建带来了巨大的困难。因此，在震后需要开展滑坡灾害点普查，编制滑坡目录，调查滑坡的空间分布特征，为震区景区恢复与防治、滑坡危险性分析、滑坡预报提供信息基础。

常规的滑坡普查方法主要有人工实地调查与光学遥感影像解译两类方法。地震造成的滑坡分布范围广、数量众多，山区地形复杂且地域广阔，依靠常规的人工野外实地调查耗时耗力，且地震诱发的滑坡大部分分布在坡度相对陡峭的山坡，危险性极大。遥感技术具有宏观、快速、覆盖范围广的特征优势，能快速地获取震后灾区的遥感影像。近年来，国内外学者提出了一批基于光学影像的半自动滑坡快速识别算法，包括震后单景影像的滑坡信息提取方法(李松 等，2015；Martha 等，2010；Lu 等，2011；Stumpf和Kerle，2011；Mondini 等，2013)与震前震后影像变化检测识别方法(Singh，1989；朱林，2016；赵祥 等，2009；李松 等，2010)，实现了滑坡信息的有效识别。由于地震后震区天气环境恶劣，光学影像受云雾天气影响难以有效成像，因此常导致滑坡信息识别不全面。

合成孔径雷达SAR (Synthetic Aperture Radar)卫星技术不受云雾等恶劣天气的影响，可全天时、全天候地成像，突破了光学影像的局限性，逐渐成为地震应急与震害评估的主流遥感技术手段。相关学者在基于SAR图像的滑坡识别方法方面开展了大量探索。Chorowicz等(1998)基于滑坡体地貌特征对SAR图像中的滑坡进行解译。Bovolo和Bruzzone (2005)利用基于小波变换的多尺度变化检测方法提取了SAR图像中的滑坡体信息。Rott和Nagler(2006)利用雷达差分干涉技术以及边坡形变分析技术评估了SAR图像中的滑坡灾害。陈莹等(2011)利用灾害前CBERS-02B数据和Radarsat-1 SAR数据相结合的方法，提取了汶川地震北川县城的滑坡灾害分布，证明了方法的有效性。随着SAR传感器的发展，SAR从原先的单一工作模式逐渐发展成为多波段、多极化、多角度的工作模式，为地震滑坡识别应用提供灵活性保障(郭华东和李新武，2011)。然而，当前的识别方法主要利用SAR图像的纹理特征、极化特征等单一特征，并未融合SAR图像中的多源特征，尤其是全极化SAR图像数据中的特征，使得滑坡点普查精度不高，不能满足实际应用的需求。

与单极化SAR数据相比，全极化SAR数据包含目标的信息量更加丰富，获取的目标极化散射特性与目标本身的形状结构、物理属性之间有本质的联系，能够更准确地识别地物(陈劲松 等，2004)。王兴玲等(2014)基于机载全极化SAR图像，采用Bayes决策的方法提取了北川境内的滑坡灾害体分布信息。现有的基于震后单时相极化SAR影像滑坡体识别研究较少，绝大部分研究主要集中于震害建筑物的提取。在特征应用方面，缺乏极化特征的挖掘研究，同时对纹理特征的应用也有限(Zhai 等，2016)，滑坡识别准确度不高。

针对此问题，本文以九寨沟地震为例，采用中国首颗分辨率为1 m的C频段多极化高分三号卫星数据，在深入分析滑坡全极化影像数据特征的基础上，研究BP神经网络自动分类技术，构建集极化特征、纹理特征、地形特征等特征的滑坡自动识别模型，实现地震滑坡点的普查，为九寨沟地震景区恢复重建及科考提供参考依据。

8月8日九寨沟地震震中位于岷江断裂、塔藏断裂和虎牙断裂附近(图1)。李渝生等(2017)基于野外调查的地表变形，推测可能发震断层位置如图1中虚线所示。震后短时间内获取了不同模式的SAR卫星数据，本文采用高分三号全极化数据(获取时间：2017年8月9日)，几何纠正处理后影像空间分辨率为8 m。影像彩色合成图像由HH(红色)、VV(绿色)、HV(蓝色)进行彩色合成(图2)。采用Lee滤波方法对数据进行滤波来降低图像噪声影响。为验证方法的有效性，利用震后无人机获取的航空正射影像(获取时间：2017年8月12日)、高分二号光学影像作为参考对比影像，无人机影像空间分辨率为0.14 m，高分二号光学影像空间分辨率为0.8 m。

图 1 九寨沟及其周围区域构造格架图 Figure 1 Tectonics frame of Jiuzhaigou and its adjacent area

图 2 研究区高分三号彩色合成图像 Figure 2 Color composite of GF-3 image Research area

地震滑坡具有破坏性强、破坏范围大的特点，通常造成地表覆盖植被的破坏。地震滑坡形成的废墟由乱石、泥土、杂物等无序堆积而成，散射特征杂乱无序，随机性较大，主导散射机制为体散射。但是受某些大体积杂物或石块的影响，在极化SAR影像中常常表现出一些不规则的亮点分布，使其功率比相同体积的植被冠层要高。

为分析地震滑坡体、居民区、植被在高分三号影像中的散射特征差异，在研究区的高分三号影像中绘制剖面线，统计剖面线上不同地物类型像元的特征值，用来分析不同类型地物目标在高分三号影像中的特征表现。3类地物的HH、VV和HV共3个极化通道剖面线图中(图3)，剖面线数值(图3中的Y轴)是剖面线经过的每行像素的强度均值，而每行像素的相对位置是图3中的X坐标。由此可见，HH和VV极化通道上的强度是非常相近的，而HV极化通道的强度均大于HH和VV。对于HH、VV、HV极化，滑坡体的散射强度均高于植被与居民区的散射强度。从剖面线的形状可以发现，植被与居民区分布区域，其3个极化通道上的散射强度随像素位置变化均没有太明显的抖动，然而滑坡体的散射强度随像素位置的变化抖动非常厉害，有显著的高峰和低谷，这个现象说明，植被、居民区的散射更加均匀化，更稳定，而滑坡体的散射特征随像元位置的变化浮动较大。滑坡体组成物质复杂，因此其散射也由多部分组成，包括二面角强散射区、背景区及弱散射区等。加之滑坡体的斑点效应也较突出，因此滑坡体分布区域的散射特征在各极化通道中的强度是浮动变化的。

图 3 不同地物类型散射特征剖面线分析 Figure 3 Analysis of scattering characteristics section lines of different ground features

极化SAR数据提供了观测目标的散射机制分析，利用Yamaguchi极化目标分解方法可获取研究区的奇数次(单次)、偶次、螺旋体和随机散射机制强度图。在研究区的高分三号影像中，采用与图3相同位置的剖面线来统计经过像元的散射分量值，进而分析不同地物目标类型在高分三号影像中的散射分量强度(图4)。对比偶次散射、奇次散射和随机多次散射，滑坡体区域特征值均大于植被和居民区分布区域，三者的奇次散射和偶次散射分量均高于随机多次散射分量。滑坡体分布区域的奇次散射和偶次散射强度整体比植被与居民区强很多，而偶次散射成分，三者相当。

图 4 不同地物类型散射机制剖面线分析 Figure 4 Profile analysis of scattering mechanism of different ground objects

基于上述滑坡散射特征分析，本文提出了以高分三号全极化SAR数据为数据源的BP神经网络地震滑坡提取方法。首先对高分三号数据进行地理编码、滤波、极化分解等数据处理，在此基础上提取影像的极化特征和纹理特征。随后，将震后获取的高分二号卫星数据与高分三号卫星数据进行精配准，并在高分二号影像上选取典型的滑坡点样本作为分类训练样本。最后，综合利用极化特征、纹理特征、地形特征信息，基于选取的训练样本，构建BP神经网络模型，实现全区域滑坡体的普查(图5)。

SAR在成像过程中不可避免的产生斑点噪声，斑点噪声的存在严重影响了SAR图像的质量。数据预处理主要针对高分3号全极化影像数据，其目的是最大程度地抑制斑点噪声影响，提高目标识别的准确度。数据预处理主要包括图像滤波、地理编码、几何校正。

图 5 滑坡普查技术流程图 Figure 5 Technical flowchart of landslide general investigation

极化SAR图像获取之后，需要对图像进行解译，揭示隐含在图像中的目标信息，为下一步的决策提供依据。极化目标分解理论为理解目标和提取目标信息提供了一种不同于其他遥感信息提取的手段。极化分解本身不会带来关于地表的新信息，但是可从中提取被成像地区的特征并分析产生后向散射的原因，从而更好地辅助极化SAR图像解译(王超 等，2008)。目标的极化特性与其形状结构有着本质的联系，相对完整地描述了雷达目标的电磁散射特性，可反映目标表面粗糙度、对称性和取向等其他雷达参数不能提供的信息，是完整刻画目标特性所不可或缺的内容。本实验中采用Yamaguchi分解方法获取奇数次(单次)、偶次、螺旋体和随机散射机制特征信息。

SAR图像中含有丰富的纹理特征，图像中滑坡体的纹理特征与居民区、植被区的存在差异，采用灰度共生矩阵来计算空间纹理特征，计算的纹理特征主要包括均值、方差、同质性、异质性、熵、角二阶矩、相关性、对比度。

地形特征主要描述滑坡的微地貌信息，通过数字高程模型DEM(Digital Elevation Model)计算获得，主要包括高程和坡度信息。

训练样本用来训练构建多个特征集合，为后续大区域信息识别提供分类特征规则集。基于精配准的高分2号影像选择滑坡训练样本，选择特征明显、分布范围较大的滑坡体图像块作为训练样本，最小训练样本图像块为4个像元，最大为64个像元，训练样本的选取在图像中尽可能均匀分布，本实验中选择的训练样本个数为11。选取图像块覆盖范围的各特征均值作为训练样本的输入特征。

BP(Back Propagation)神经网络是1986年由Rumelhart和McCelland为首的科学家小组提出，是一种按误差逆传播算法训练的多层前馈网络，是目前应用最广泛的神经网络模型之一(赵一博 等，2013)。BP网络能学习和存储大量的输入—输出模式映射关系，而无需事前揭示描述这种映射关系的数学方程。它的学习规则是使用最速下降法，通过反向传播来不断调整网络的权值和阈值，使网络的误差平方和最小。BP神经网络模型拓扑结构包括(图6)：输入层(input layer)、隐层(hide layer)和输出层(output layer)(徐黎明 等，2013)。

图 6 BP神经网络模型拓扑结构(骆剑承 等，2001) Figure 6 BP neural network model topology (Luo, et al.，2001)

(1)输入层与输出层节点数的确定。输入层的节点通常与参与的特征空间的维数是一致的。本实验参与的特征主要包括极化特征、纹理特征和地形特征，特征维数为14，因此，输入层的神经元数设置为14。输出层的节点数目与分类结果数据一致，本实验分为3类，因此输出层的节点数目设置为3。

(2)隐含层节点数确定。BP神经网络主要特征是具有非线性映射的功能，意味着最基本的三层BP网络也可完成任意的N维到M维的映射，如果选取的隐含层数目过多，会导致网络变得复杂，从而增加训练时间，导致网络陷入瘫痪，影响信息提取结果。本实验中隐含节点数目的确定参考屈晓辉等(2007)的公式来确定，隐含节点数目= $ \dfrac{{A \times B + \dfrac{1}{2}A \times \left({{B^2} + B} \right) - 1}}{{A + B}}$ ，式中，A为分类数，B为特征向量维数。根据此公式，该网络中的隐含节点数目设置为20。

(3)阈值、学习速率和动量因子的确定。神经元输入、输出函数的阈值θ决定了S型函数曲线变化的梯度。θ值越小，曲线越陡，函数随自变量变化变得敏感，易发生振荡；当θ值越大时，曲线平坦，学习速度加大，但学习过程线性化，对网络的识别能力造成影响。通过多次分析与实验，θ的初始值设定为0.25。学习速率影响算法的收敛程度，因此学习速率不能设置过大，但是也不能太小，如果太小，会造成训练时间很长(王旭 等，2000)，本文中学习速率设置为0.04。动量因子可滤除学习过程中的高频振荡，提高收敛速度(朱矩波和马士林，1998)，文中动量因子最大迭代次数设置为1500次。

综合利用极化特征、纹理特征与地形特征，基于获取的建筑物样本，构建BP神经网络模型，模型参数设置详见第4节，实现滑坡体的提取，识别得到的滑坡体分布(图7)。研究区共识别滑坡841处，面积约为7.64×10⁴ m²。滑坡分布距离震中较近，主要分布在距震中约10 km范围内。识别的滑坡多数面积较小，以中小型滑坡为主。叠加道路矢量数据发现，滑坡沿道路较为发育，形成的堆积物易造成道路的堵塞，这在无人机影像可清晰地发现(图8)。分析地形特征发现，滑坡沿河谷较为发育，尤其是在九寨沟景区的熊猫海附近。

图 7 滑坡普查结果分布图 Figure 7 Landslide distribution map

图 8 无人机航拍影像 Figure 8 Images of unmanned aerial vehicle

为定量地评估方法的有效性，将识别结果与高分二号光学影像和无人机航空影像目视解译结果进行对比，因震后获取的高分二号光学影像中含云量较高，部分区域未有效成像，而九寨沟地震后，中国地震局地壳应力研究所携带的无人机拍摄的航空影像获取了被云遮盖的影像部分，两者的结合能完整地还原震后场景。同时为了验证多特征应用信息识别的有效性，在基于相同方法的基础上，仅利用单一纹理特征、极化特征对滑坡体分布信息进行提取(图9)，并对其识别精度进行评价。

图 9 滑坡普查结果分布图 Figure 9 Landslide distribution map

定量评估指标选择用户精度UA(User’s Accuracy)，制图精度PA(Producer’s Accuracy)，错分误差CA(Commison Error)，漏分误差OE(Ommison Error)，总体精度OA(Overall Accuracy)与Kappa系数，指标利用混淆矩阵计算获得。不同特征应用的滑坡体识别精度定量评估表明(表1)。综合多特征的BP神经网络方法的滑坡体识别总体精度为92.8%，相较于利用单一极化特征、纹理特征方法的识别准确度要高。基于单一纹理特征与单一极化特征方法的滑坡体总体识别精度相差不大，证明两类特征均能有效地识别滑坡体信息。仅利用极化特征、纹理特征方法的滑坡体制图精度比多特征综合应用的精度要高，而用户精度较低，这说明了仅利用单一的特征存在较大的错分误差，而加入地形特征之后，有效地降低了滑坡体的错分误差，证明地形特征能较好地剔除非滑坡信息。

表 1(Table 1)

表 1 滑坡识别精度分析 Table 1 Accuracy analysis of landslide identification results 特征 类别 PA/% UA/% CA/% OE/% OA/% Kappa系数 极化特征、纹理特征、地形特征 滑坡体 64.5 82.9 17.0 35.5 92.8 0.715 纹理特征 滑坡体 70.5 67.6 32.4 29.5 80.0 0.63 极化特征 滑坡体 76.7 74.6 25.4 23.2 80.5 0.65

表 1 滑坡识别精度分析 Table 1 Accuracy analysis of landslide identification results

采用无人机航空正射影像来对比分析滑坡提取的准确性，图7中为利用高分三号全极化数据提取的滑坡分布，从航空无人机拍摄的滑坡可见(图10)，基于高分三号全极化SAR影像可较完整地提取滑坡灾害体的分布。高分三号影像中，部分山脊被误分为滑坡，这是由于山脊在图像上散射特征较为明显，与滑坡体形成的亮斑混淆在一起，通过加入地形特征，即滑坡体在坡度上的空间分布特征可排除部分错分的滑坡信息，消除山脊带来的影响。

地震应急过程中震区灾情的调查通常分为灾情普查与灾情详查两部分。灾情普查主要满足震后初期灾情分布范围与破坏程度的需求，为救援力量分配、灾情初步评估提供依据，这个阶段要求影像的覆盖范围大。灾情详查主要用于震后灾情空间分布范围、位置、破坏程度的调查，可为精细灾情评估、应急精准服务、震后科考提供信息，此阶段要求灾情普查的遥感影像具有较高的空间分辨率。

图 10 局部滑坡分布区域方法对比图 Figure 10 Enlarged contrast map of local landslide distribution area

无人机技术具有机动、快速获取灾区灾情信息的能力，获取的无人机影像空间分辨率高，可实现滑坡体的详查。基于无人机影像采用提出的多层次面向对象分析的方法提取了漳扎镇附件区域的滑坡空间分布(图11)，该方法可有效地实现滑坡灾害的提取，且空间分布的范围、受灾程度等信息均可清晰的展现，满足地震滑坡详查的目的。但无人机获取的影像与卫星影像相比，尽管其空间分辨率高，但其覆盖的范围较小，在震后应急与灾情评估的短时间内，亟需了解整个震区灾情的空间分布位置与范围，从而计算灾害空间分布密度，掌握灾区受灾范围，从而为滑坡灾害危险性分析与预测提供基础，进而为地震应急救援辅助决策提供服务。在此实际需求背景下，需要的是获取大范围的灾情分布情况，即实现灾害空间分布的详查。高分三号影像具有覆盖范围广的特征优势，同时在云、雾等恶劣天气条件下也有效地成像，在利用其极化特征与纹理特征的基础上，可实现大范围区域性的灾情调查。同时，覆盖部分区域的高分辨率光学遥感图像也可为滑坡体样本的选取及后续精度的评估验证提供数据支撑。

图 11 无人机影像滑坡提取结果 Figure 11 Landslide extraction results using unmanned aerial vehicle image

针对光学遥感技术易受云雾等恶劣天气条件影响难以有效成像、信息无法识别的问题，提出了融合多维特征的神经网络模型的高分三号极化SAR数据滑坡信息识别方法，通过与仅利用单一纹理特征与极化特征的信息识别方法进行实验对比分析，结果表明，本文提出方法的识别精度得到了显著的提高，证明了方法的有效性。

本文提出的方法融合了极化SAR图像的极化特征、纹理特征以及地形等辅助特征数据，构建多维特征神经网络模型，避免了利用单一特征造成识别信息冗余的难题。基于选取的样本，挖掘样本数据的特征，构建滑坡识别特征规则，避免了特征阈值人工选择的问题，实现信息的自动化识别。

基于国产高分三号极化SAR卫星数据，分析了不同地物类型目标在图像中的特征差异，并实现了高精度地滑坡空间分布识别，检验了高分三号卫星数据在地震应急服务中的能力，也促进国产卫星数据在防震减灾中的应用。

本文建立的方法融合了极化特征、纹理特征和地形特征，但提出的分类方法以像元为单位，且特征的选取是人工进行的，在高分辨率遥感图像信息识别方面仍存在一定的局限性。未来研究中，需引入极化SAR图像面向对象与深度学习方法相结合的方法，进一步提高信息识别准确度。

志 谢 此次实验的数据获取来源于中国资源卫星应用中心，在此表示衷心的感谢！