基于时序InSAR的滑坡早期识别——以贵州省水城县为例

廖 军，吴彩燕，王立娟，谭秋焰，朱新婷

（1.西南科技大学环境与资源学院，四川绵阳 621010；2.四川省安全科学技术研究院，四川成都 610045）

引言

滑坡是指岩体、碎屑、土壤沿斜坡自上而下滑动的现象，也指斜坡在重力作用下沿着坡体内部一个或者多个面做剪切运动的现象［1～3］，在我国的山区和峡谷地带广泛分布，尤其是西南高山区和三峡库区，发生频率最高，造成的损失最严重［4～5］。常规的滑坡识别手段主要有野外实地勘查和遥感影像识别，野外实地勘查能准确圈定滑坡要素，但其工作既费时又低效，并且勘测区域受到地形的极大限制［6］；遥感影像识别滑坡大多是利用专家经验知识通过光学遥感影像的光谱、空间、纹理等信息进行的［7～8］，尽管已经取得了相当不错的成就，但光学影像容易受到天气因素的限制如大雾、降雨等，并且在监测和识别滑动不明显的潜在滑坡能力不足［9～10］。而新出现的合成孔径雷达干涉测量（InSAR）技术不仅能精确获取地表微小形变信息，且具有全天候、全天时、时效性高等优点，成为了国内外学者研究的热点［11～14］。

近年来，InSAR技术逐渐应用于滑坡灾害识别，并取得了一些成功的案例。廖明生等以三峡库区滑坡为对象，使用高分辨SAR数据对其进行监测，成功获得了滑坡发生的时间、位置及形变情况，证明了高分辨率雷达数据有效应用在滑坡监测中［15］；Tofani V等［16］采用永久散射体干涉测量（Persistent Scatterer，PS-InSAR）技术对意大利北部的Santo Stefano d′Aveto滑坡进行分类与监测，结果表明PS技术结合地面数据能够准确判别滑坡的边界和运动状态；Bianchini等［17］以意大利的Calabria地区为研究区，利用PS-InSAR技术对该地区的滑坡进行了识别及滑坡编目的更新；赵超英等［18］利用多源SAR数据来识别黑方台地区的潜在滑坡，结合实地调查和光学影像验证了识别结果的可靠性，并分析典型滑坡的失稳模式；Dong等［19］使用相干散射体（Coherent Scatterer，CS-InSAR）技术更新了丹巴县甲居滑坡的活动边界，绘制了更详细的滑坡位移图。现有的研究虽然在一定程度上可以对潜在滑坡进行识别，但是这些方法大多是从雷达视线方向（LOS）角度出发进行识别研究，虽然可以识别滑坡发生的时间和地理位置，但LOS方向与地表和山坡之间存在一定的夹角，不能直接反映滑坡的变形程度，也不能得到滑坡变形的真实时间序列，难以识别滑坡前后的具体变形趋势。

本研究以滑坡灾害多发的水城县为研究区，通过改进的SBAS-InSAR技术，获得水城县坡度方向形变信息，结合水城县雷达可视分区及谷歌影像剔除无效观测值，利用Anselin Local Moran′s I对获取的监测值进行异常值分析和聚类处理，获取低值聚集区域，最后结合谷歌影像的地表特征，实现水城县潜在滑坡的识别，用于指导该区域滑坡灾害的预测预警。

1 研究区概况及数据源

1.1 研究区概况

水城县位于贵州省六盘水市，邻接云南省，在水系和构造运动作用下，西北地势较东南高。该区属于亚热带湿润季风气候地区，降雨集中在5～10月份，其中6～8月降雨量占全年53.03%［20］，县内河网密闭，水量集中，是造成境内地质灾害的主要诱因之一。水城县地质构造复杂，存在多个构造变形区，地层岩性以二叠系、三叠系以及第四系为主。

水城县境内滑坡广泛分布在飞仙关组、峨眉山组及龙潭组等地层中，该地层中节理发育，岩层破碎，极易在人类活动或降水等因素下发生滑坡等地质灾害［20～21］。

图1 研究区地理位置Fig.1 Geographical location of the study area

1.2 数据源介绍

文中采用的数据源包括欧空局新一代C波段的中高分辨率卫星Sentinel－1A SAR数据，选取了2018年1月～2018年12月共20景降轨数据，入射角为39.5°，访问周期为12 d，拍摄模式为IW，空间分辨率为5×20m（Rg×Az），极化方式为VV，并使用POD精密定轨星历数据（Precise Orbit Ephemerides，POD）进行轨道矫正，定位精度优于5 cm；数字高程模型（DEM）为30 m分辨率的SRTM DEM（Shuttle Radar Topography Mission，SRTM）；从Google Earth中获取的研究区光学遥感影像数据；野外调查资料。

2 研究方法

本研究选取2018/1～2018/12共20景Sentinel－1A SAR降轨数据，利用永久散射体干涉测量（PS-InSAR）技术提取获得相干性高且相对稳定的散射体用做地面控制点（Ground Control Point，GCP），引入到小基线集干涉测量（SBAS-In-SAR）技术处理流程中，获取长时间序列的地表形变信息（LOS向），并转换为坡度方向形变信息；结合水城县地形可视性分区及谷歌影像，剔除不可视区（阴影、透视收缩、叠掩）、高密度林地及灌木林、河流等不可靠的形变信息，保留其他区域可靠形变信息；对获取的可靠形变信息采用Anselin Local Moran′s I指数进行空间异常值分析和聚类处理，获取低值聚集区域，并结合谷歌影像的地表特征，从而对潜在滑坡进行识别。滑坡识别路线如图2。

图2 潜在滑坡识别流程Fig.2 Landslide identification process

2.1 结合PS-InSAR技术的GCP点选择

在SBAS-InSAR处理流程中，选取稳定的地面控制点（GCP）进行轨道精炼和相位偏移的计算，可以有效地估算和去除残余相位及可能存在的斜坡相位，进而消除平地效应，对相位转形变的可靠性起着重要作用［16－18］。

由于人工选点会对形变速率引入误差，文中利用PS-InSAR技术的优势［22］，依次设置振幅离差阈值（≤0.25）、相干性阈值（≥0.8）以及空间域的低通滤波和时域的高通滤波，减小干涉相位中的大气延迟相位误差，得到可靠的稳定散射点的形变信息，并对该结果设置速率阈值选取在－0.1～0.1 mm的稳定点作为GCP点，通过地理坐标到SAR斜距坐标系的转换，即可应用于在SBAS处理流程中。

振幅离差指数（D）是依据SAR图像上同一位置像元的强度值在时间序列上的离散特性进行目标的识别，其计算公式如下［23］：

式中：μ是时序上的某个点的振幅平均值；σ是时序上的振幅的标准差。

2.2 SBAS-InSAR技术原理及方法

小基线集干涉测量（SBAS-InSAR）是在差分干涉测量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）和永久散射体干涉测量（PS-InSAR）发展起来的一种时序监测方法，该方法能有效消弱常规D-InSAR中大气效应，时空失相干等因素对监测结果的影响，同时弥补了PS-InSAR技术不适用于相干点稀疏地区的不足［10，16～17］。在2002年，Berardino等［24］首次提出了小基线集时序监测方法，SBAS-InSAR能够获取非线性的形变结果，在滑坡形变监测中具有良好的适应性。

小基线集的原理是通过设置时间、空间基线的阈值，将多期SAR影像进行自由配对组合生成多个子集，每个子集由超主图像连接，从而形成一个大集合，保证在每个小集合中SAR影像间的基线小，集合间的SAR影像基线大，从而提高数据的采样率，克服空间基线过长导致的空间失相干影响。配合外部高精度的DEM消除地形相位，生成一系列差分干涉图，对差分干涉图进行筛选得到高相干性像元，然后利用地面控制点进行轨道精炼，剔除低相干性和解缠错误的干涉像对，采用奇异值分解法和最小二乘法对多个小基线集联合求解，从而得到研究区的时间形变序列及地表平均形变速率［25－26］。

大多数滑坡都是沿坡面滑动的，而小基线集干涉测量法只能监测雷达视线的形变，不足以反映滑坡的真实变形。因此，为了准确地描述滑坡真实的形变情况，将LOS向形变速率转换为沿坡度方向形变速率［27－28］：

式中：V s l o p e为沿坡度方向的形变速率；V los为LOS向的形变速率；β为雷达视线方向与斜坡坡面的夹角，即视坡夹角；φ为斜坡的坡度/（°）；α为斜坡的坡向/（°）；θ为雷达入射角/（°）；αs为卫星轨道方向与正北方向的夹角，即卫星飞行的方向，升轨数据为负，降轨数据为正。

2.3 聚类和异常值分析原理

采用Anselin Local Moran′s I指数对SBAS-InSAR技术获取的形变干涉点进行空间域的异常值分析和聚类分析［29～30］，通过计算要素与其他所有要素作为相邻要素的数值关系，从而获得数据集中有统计显著性的高值（热点）和低值（冷点），自动对形变点中的异常数据进行聚合，确定适当的空间分析范围，并校正数据的空间相关性，得到形变聚类点置信水平p值和得分，其相关的计算公式如下：

统计量z I i得分的计算方法如下：

式中：I i为形变点i的Anselin Local Moran′s I指数值；x i、x j分别为点i，j的属性；n为选定的窗口内要素的总数；是对应属性的平均值；w i,j是要素i和j之间的空间权重；S2和S分别为为选定的窗口内像元的方差、标准差；μ0为总体均值；Z I i为形变聚类点i的得分；p i为形变聚类点i的显著性置信水平。

3 数据处理

3.1 雷达可视性分析

雷达遥感通过向地面发送脉冲信号，接收地面反射的回波信号，实现对地物的观测。尽管具有全天时、全天候、大面积等无法比拟的优势，由于SAR采用侧视成像方式，其入射角（θ）与实际地形存在夹角，使得SAR影像会出现阴影、透视收缩、叠掩3类常见的几何畸变［27］。根据Sentinel－1A卫星降轨条件下的几何参数及水城县DEM提取的地形参数（坡度、坡向）可获得水城县可视性分区方式，并获得水城县可视性区域（见表2）。

表2 SAR影像可见性分区Table 2 SAR image visibility zoning

Sentinel－1A降轨影像在水城县的可视性区域分布如图3所示，可见区域约占35.76%，低敏感区域约占24.08%，不可见区域约占40.16%。其中，可见区域以主要分布于东南、西南区域；不可见区域主要分布在中心区域。

图3 研究区雷达可视性分区Fig.3 Radar visibility zoning of in study area

3.2 基于改进SBAS技术的数据处理

文中使用精密轨道数据（POD）和30m分辨率的SRTM DEM进行影像轨道误差和地形相位校正。为了提高时间、空间相干性，保证干涉图质量，将时间基线设置为120 d，空间基线阈值为极限基线的5%（5 162.12×0.05=258.106 m），使用3D解缠。超级主影像为2018年7月3日，共生成了138个干涉像对，如图4所示。通过查看生成的相干系数图、去平地干涉图、相位解缠图，剔除相干性较差的干涉对，最终使用了115个干涉相对进行后续处理。

图4 数据集时空基线分布Fig.4 Spatiotemporal baseline distribution of the dataset

图5（a）是去平地后的相位干涉图，水城县由于植被覆盖密集，地形起伏较大，导致地形相位去除效果一般，为了减小植被等因素的影响，选取了Delaunay MCF方法进行相位解缠；为了获得较高的地面分辨率，将方位向与距离向的视数比设置为1：4；为了获得更为清晰干涉条纹，选取了Goldstein滤波进行噪声去除，滤波后的干涉图结果如图5（b），可看出滤波后噪声误差减小，干涉条纹的清晰度提高。

图5 研究区部分干涉处理结果图Fig.5 Results of partial interference processing in the study area

3.3 坡度方向可靠形变信息获取

基于改进的SBAS-InSAR技术，获取了2018年水城县LOS方向的地表变形速率，对于西南山区，LOS向的形变速率不足以真实反映斜坡的真正形变情况，因此将LOS向形变速率转换到坡度方向形变速率。当β接近90°时，cosβ近乎为0，导致V slop e趋于无穷大。为了对产生的异常值进行剔除，Herrera等以cosβ=±0.3为阈值，即V sl ope≤3.33V los。设定0＜C＜0.3时，C=0.3；-0.3＜C＜0时，C=-0.3［27－28，31－32］。V slope表示沿坡度方向的位移速率，块体移动方向只能沿斜坡向下，故剔除V slope为非负的形变点。

基于水城县可视性分区结果，先剔除了由透视收缩、叠掩和阴影不可信区域的形变信息，保留了好可视区域和低敏感区域中的形变信息。由于水城县南部地区常年覆盖着茂密的植被，尽管处于可见区域内，但C波段地表穿透能力有限，即便SAR数据时间间隔短，也会导致失相干，得到有效监测值较少且监测精度有限。此外，滑坡在河流及平坦的地区发生的可能性较低。因此，在可视性分区结果的基础上，结合谷歌影像，将密度较高的林地、灌木林、河流以及坡度小于5°不可靠区域的形变信息进行剔除，保留其它区域可靠形变信息，得到水城县坡度方向年均形变速率（图6）。

图6 研究区2018年坡度方向年均形变速率Fig.6 Average annual deformation rate of the slope direction in study area in 2018

4 结果与分析

4.1 水城县变形监测结果分析

从图6中可以看出整个水城县北部的形变速率较为突出，各地形变速率不均，存在多个形变中心，集中分布在陡坡上。区域1形变较为严重，年平均速率在－75～0 mm/a，矿区开采处的变速率最明显，形变速率达到了－142.3 mm/a；区域2整体地形起伏较大，形变程度也较为突出，其年平均形变速率在－35～0 mm/a之间；区域3受到雷达监测制约以及植被覆盖影响，获取的监测结果有限，形变速率大多在－10～0mm/a，少部分区域形变速率达到了－20 mm/a左右，整体较为稳定。

由上述分析可知，研究区内的变形区具有分布不均匀的特点，主要是人类活动因素（矿区开采、修筑道路）及差异化的地理环境因素（坡度、坡向、植被覆盖率等）造成的。根据已有资料分析，区域1分布着较多的矿山，受矿区开采的影响，导致该区域整体形变量较大，易导致滑坡的发生。区域2地形起伏大、山体陡峭、植被覆盖率低，受地形地貌等地理环境因素的影响，其山体边坡的稳定性会明显低于区域3，发生滑坡的可能性相对较高。

对获取的形变信息进行统计分析，形变速率在－15～0 mm/a的形变点占据总数的40.5%，形变速率在－35～－15 mm/a的形变点占据总数的31.5%，形变速率在－55～－35 mm/a的形变点占据总数的13.5%，形变速率大于－55 mm/a的形变点占总数14.5%。从统计分析结果来看，形变速率大部分集中在－35～0 mm/a，其中形变速率小于－30 mm/a所占比例最大，在一定程度上能够表明水城县地表整体较为稳定，形变主要集中在水城县的北部区域。。

根据上述对水城县坡度方向形变速率分析，结合国内外众多学者针对滑坡事故做出的研究成果［13－18，22－23］，文中选取V slop e=15 mm/a为速率阈值，即V slope＞15 mm/a的区域易发生滑坡事故。

4.2 水城县滑坡识别及验证

4.2.1 滑坡识别

通过Anselin Local Moran′s I对坡度方向形变点进行空间域异常值分析和聚类处理，会得到HH、HL、LH、LL四种聚类类型，其中HH表示高值聚集，HL表示低值包含高值异常，LH表示高值包含低值异常，LL表示低值聚集。通过获得低值聚集区域，结合滑坡在谷歌影像上的形状、分布、微地貌特征，确定滑坡位置及边界。通过本文的方法最终共识别出了14处潜在滑坡，滑坡详细信息见表3，滑坡识别结果分布见图7，局部滑坡识别结果见图8。

图8 部分滑坡识别结果图Fig.8 Some landslide recognition results

表3 研究区滑坡识别结果Table 3 Landslide identification in the study area

根据表3、图7分析，识别出的滑坡有一半以上是由于剧烈的人类工程活动诱发导致，85%的滑坡类型为土质滑坡，在强降雨的情况下更容易出现失稳，结合H1、H5、H9的时序形变分析，在6～9月份形变趋势处于急剧增加的态势，是因为该期间迎来降雨高峰期，雨水富集在松散土层和相对防水的基岩接触面之间，增加了土壤的自重，同时软化了滑带土，使其内聚力和内摩擦角减小，使得坡体的稳定性降低，为土层滑坡的形成奠定了基础。

图7 研究区滑坡识别分布图Fig.7 Distribution map of landslide identification in study area

4.2.2 验证分析

根据表3中潜在滑坡的高程值和坡度值进行分析，95%以上的滑坡分布在高程大于1 000 m，坡度大于20°的区域，借助前人对水城县滑坡分布规律特征的研究成果进行对比验证［20］，表明识别的潜在滑坡与水城县滑坡分布规律特征相吻合，在一定程度上说明了识别结果的可靠性。

此外，通过贵州省地质灾害排查项目及野外实地调查获得了水城县的22处滑坡，用于对本文识别的潜在滑坡进行验证。验证结果：识别获得的14处潜在滑坡，有11处潜在滑坡位置和范围与野外调查的滑坡较为吻合（占识别结果的78.6%），有3处潜在滑坡范围与野外调查的结果存在一定偏差（占识别结果的21.4%），部分识别结果见图8，另有8处滑坡未能识别。

根据图7，图8可知，对于成功识别的滑坡点位置，大多是位于坡度和形变程度较大的区域。未能成功识别的8处滑坡点，但也均监测到了部分形变区域，其未能识别的主要原因是：这些区域大多处于雷达不可见和不敏感区域，获得有效形变点较少；部分地区植被覆盖严重，为避免监测结果受大面积植被因素的影响，剔除了该区域的形变点，减少了该区域的有效时序形变点位提取量，导致无法准确识别。

5 结论与讨论

5.1 结论

本文以滑坡多发的水城县为研究区，探讨了一种基于InSAR技术对潜在滑坡识别中的应用。使用欧空局提供的Sentinel－1 A SAR数据，利用改进的SBAS-InSAR方法获取地表坡度方向形变信息，通过Anselin Local Moran′s I进行聚类分析，获得高异常值区域，并结合滑坡的形状、分布、微地貌特征成功识别出了14处潜在滑坡。将识别的潜在滑坡与野外勘探资料相核实，验证了该时序InSAR用于滑坡隐患识别的有效性及优势。

本研究受限于单一的降轨数据，对地形可视性差的区域无法进行有效识别，因此无法完全提取研究区内潜在滑坡体。文中是对整个县城的进行时序监测，在处理数据上所花费时长较大，为了更快速的获取识别结果，可利用升降轨数据直接对滑坡高易发区进行监测识别，从而有针对性地调查和监测高概率潜在的滑坡范围，达到防止滑坡灾害发生的目的。

5.2 讨论

由于受到SAR系统监测方式的制约，在对西南高山区域的滑坡识别应用中，导致可探测的形变区域面积较少，难以全方位对潜在的地质灾害隐患进行监测识别，限制了通过InSAR技术对潜在地质灾害识别的应用。本研究获取数据类型单一、有效观测区覆盖范围较小，导致部分滑坡点未能进行监测，因此可以选用联合升轨和降轨观测减少观测盲区，提升滑坡的有效探测率。

现有的民用雷达卫星通常工作于X波段（波长3.1cm）、C波段（波长5.6cm）或L波段（波长2.3.6cm），但是能大量获取的免费SAR数据为C波段的Sentinel－1 A数据。对于InSAR形变测量来说，雷达波长决定测量灵敏度，即波长越短，监测精度越高。但是在高山峡谷区域，X波段和C波段的波长较短，难以穿透植被冠层，容易导致时间失相干；使用波段更长、穿透力更强的L波段进行观测时，数据往往具有更好的相干性，在西南山区滑坡识别应用中更具有优势。