滑坡具有隐蔽性强、破坏程度大、发生频率高和分布广泛等特点。贵州省位于中国西南地区，地貌类型主要为高原、山地、丘陵和盆地^[1]，气候温暖湿润、多云多雨、植被茂密，这会增大滑坡的隐蔽性和多发性，也会增加光学遥感成像和解译滑坡的难度。《2020年度贵州省地质灾害防治方案》提出“空天地”一体化的地灾隐患调查，使用InSAR开展广域地表形变监测，识别、提取出更多地质灾害隐患点。

卫星InSAR利用微波频段的相位信息探测雷达视线向cm级甚至mm级微小地表形变，具有广覆盖、全天候、全天时、高重访和主动遥感的优势。由于传统D-InSAR受时空去相干和大气扰动等限制，以PS InSAR^[2]和SBAS InSAR^[3]为代表的时序InSAR应运而生，可以在很大程度上解决D-InSAR的应用限制。在山区环境中，分布式散射体(distributed scatterers，DS)数量远大于永久散射体数量，因此SBAS InSAR比PS InSAR更适用于山区场景的滑坡形变监测，能够显著减小SAR影像的数量需求。

2014年以来，中高分辨率Sentinel-1A数据免费向用户开放，极大推进了InSAR在我国的普及和应用。针对重大地质灾害隐患的早期识别，我国学者相继提出空-天-地“三查”技术体系和形态-形变-形式“三形观测”的技术思路，强调InSAR在地质灾害监测中的普查作用，使得InSAR在地质灾害调查中逐渐得到认可^[4-5]。本文以贵州省最近一次滑塌体量大且致灾严重的鸡场镇滑坡为研究背景，获取2018-07~2019-07覆盖贵州省水城县鸡场镇滑坡发生前的31期升轨和30期降轨Sentinel-1A数据，结合30 m分辨率的SRTM DEM数据，使用SBAS InSAR开展滑坡发生前的地表形变监测和隐患早期识别研究。该研究可分析和评估SBAS InSAR在中国西南山区地表形变监测中的应用效果，为贵州省以及中国西南山区滑坡隐患调查与早期识别提供技术参考。

1 研究区域及数据

贵州省位于中国西南地区，地貌类型以高原、山地、丘陵和盆地为主，其中山地和丘陵占92.5%^[1]。贵州省年平均降雨量超过1 000 mm，主要集中在夏季。凹凸不平的地貌类型和时域集中的强降雨特征易引发滑坡、泥石流等地质灾害。2019-07-23水城县鸡场镇发生大型滑坡(图 1白色轮廓区)，该滑坡体属于构造侵蚀、剥蚀型中山地貌，滑坡在高海拔地带发生，演变为沿两条沟道分布的长滑移距离型泥石流，整个滑坡区域长约1 300 m，前后缘高差约460 m，体积约为191.2×10⁴ m³，属于典型的高位远程滑坡，滑坡体和碎屑流的主滑方向约为26°^[6]。

研究区范围为104°35′~104°49′E、26°12′~ 26°27′N(图 1红框)。实验采用欧空局Sentinel-1A渐进式地形观测模式干涉宽幅成像模式VV极化SAR数据集，右视侧视角约为39°。覆盖鸡场镇滑坡的Sentinel-1A数据有升轨第128轨、降轨第164轨，升轨成像时段为2018-07-25~2019-07-20，共31期；降轨成像时段为2018-07-27~2019-07-22，共30期，均为灾前影像。SAR数据集空间覆盖范围如图 1蓝框所示。辅助数据有POD精密定轨星历数据和分辨率为30 m的SRTM DEM数据，其中POD数据用于精化轨道坐标，SRTM DEM数据用于SAR数据精配准、去地形相位和地理编码等。

2 研究方法

本研究使用FaSHPS算法^[7]统计、提取同质像元作为DS候选点，设定干涉对集的平均相干阈值筛选同质像元作为DS点。由于研究区植被较多，去相干情况较为严重，因此使用像对干涉相干性最优且连通的原则生成像对干涉连接图^[8]。使用StaMPS软件进行SBAS InSAR处理，反演地表形变场。整体技术流程如图 2所示。

贵州省水城县地形起伏较大、植被茂密，跨季影像相互配准难度较大。本文使用POD精密定轨星历数据精化轨道坐标，计算辅影像相对于主影像的初始偏移量，使精度达到像元级。在此基础上，使用强度互相关算法从主影像、辅影像中分区块提取若干同名像点，对每个区块同名像点方位向和距离向的偏移量进行最小二乘拟合，求得配准多项式系数，配准精度达到0.1像元。配准多项式系数确定后，使用三次卷积内插法将辅影像重采样到主影像空间上。

假设共有成像日期为[t₀, t₁, …，t_N]的N+1景SAR影像，并且所有SAR影像均已精配准。通过设定时间、空间基线阈值或使用干涉相干性最优且连通的原则，得到M个干涉像对，则干涉像对数M与影像数N+1满足关系式：

$ \frac{N+1}{2} \leqslant M \leqslant \frac{N(N+1)}{2} $

(1)

本研究采用干涉相干性最优且像对连通的原则组建干涉像对连接图：首先从SAR影像中选择具有代表性的小区域(如高相干地物和低相干地物各占1/2)生成全部像对连接图以及相干图(数量为$\frac{N(N+1)}{2}$)，然后对生成的所有相干图按均值大小由高到低排序并顺次连接干涉像对，直到所有干涉像对连通为止，最终形成M个干涉相干性最优且像对连通的连接图。该方法能够有效提高干涉像对集的整体干涉相干性。

影像t_A和影像t_B生成的第i幅干涉图(t_A < t_B)在雷达坐标系下坐标点(x, y)的干涉相位差为δφ_i(x, y)=φ_tA(x, y)-φ_tB(x, y)。使用时空滤波器估计并去除大气相位和高程相位差后，M个干涉相位差组成的矩阵为δ_φ= Bφ，其中B为M×N矩阵。用两个时刻之间的平均速率代替其中的未知相位，得到Bv=ΔΦ。使用SVD算法对B进行逆运算，求取v的最小二乘解，再对各时段内的形变速率进行积分即可得到时间序列形变量，详细公式见文献[9]。

为消除Sentinel-1A数据在山地、丘陵等地区的阴影、叠掩现象，本文对研究区成像时段相对一致的升轨和降轨数据分别进行独立处理，获取视线向形变速率(正值表示地表向靠近卫星方向发生形变，负值表示向远离卫星方向发生形变)。根据InSAR成像几何可知，视线向形变量在水平方向只能分解出东西方向分量，无南北方向分量，这也是卫星InSAR成像几何的固有缺陷。水平形变量V_E-W和垂直形变量V_V表达式为：

$ V_{E-W}=V_{\mathrm{LOS}} \cdot \sin \alpha, V_V=V_{\mathrm{LOS}} \cdot \cos \alpha $

(2)

式中，V_LOS为雷达视线向形变速率，α为降轨和升轨入射角。

为了更好地分析坡体形变特征，将V_LOS投影至斜坡方向：

$ V_{\text {slope }}=\frac{V_{\text {LOS }}}{\cos \gamma} $

(3)

式中，γ为雷达卫星视向线方向和斜坡坡度方向间夹角，详细公式见文献[9]。

3 数据处理与结果分析 3.1 SBAS InSAR处理

使用POD精密定轨星历数据精化Sentinel-1A轨道坐标，生成单视复数(single look complex, SLC)影像，并按照研究区范围裁剪SLC数据。使用干涉相干性最优且像对连通的原则分别对31期(ID: 0~30)升轨影像和30期(ID: 0~29)降轨影像进行干涉像对连接，分别生成64个干涉对和68个干涉对(图 3)。

对研究区SLC数据进行SBAS InSAR处理：1)选择升轨、降轨第1期(ID: 0)影像作为超级主影像，多视视数设为方位向×距离向=1×4，使用SRTM DEM进行主辅影像精配准、差分干涉处理，计算相干图、幅度图，分辨率分别为13.9 m×14.5 m(升轨)和14.0 m×14.5 m(降轨)；2)使用Goldstein算法对差分干涉图进行滤波处理；3)设定相干系数阈值为0.2进行掩膜，使用MCF算法对差分干涉图进行相位解缠；4)设定平均相干阈值为0.41，升轨共提取26 038个DS点，降轨共提取50 049个DS点；5)在DS点上建立和解算观测方程，使用时空滤波器去除大气相位和高程残差相位，使用SVD算法解算DS点的时间序列形变；6)使用SRTM DEM对反演的地表形变和幅度图进行地理编码。

3.2 结果分析与讨论

升轨、降轨Sentinel-1A影像SBAS InSAR反演的2018-07~2019-07研究区地表雷达视线向平均形变速率结果如图 4所示。DS点空间分布存在3种情况：1)西向阳坡DS点多数存在于升轨Sentinel-1A影像上，少量存在于降轨影像上；2)东向阳坡DS点多数存在于降轨Sentinel-1A影像上，少量存在于升轨影像上；3)地势较为平坦的DS点同时存在于升轨和降轨Sentinel-1A影像上。上述现象是由卫星雷达右视侧视成像造成，可通过融合升降轨时序InSAR形变场提升地势起伏较大山区、丘陵的形变监测能力。

鸡场镇滑坡体(图 4紫色区域)在发生滑移前植被覆盖度高，DS点较稀疏，上沿及周边DS点形变速率约为0，说明鸡场镇滑坡发生前滑坡体上沿及周边区域并未出现大型地表形变，表现为稳定状态。查询鸡场镇滑坡发生前两周天气状况可知(表 1)，鸡场镇滑塌体在滑坡发生当天有强降雨，在滑坡发生前两周内有11 d出现不同程度降雨。因此初步判断鸡场镇滑坡是连续降雨、强降雨诱发的突发性滑坡，已超出Sentinel-1A的12 d重访周期监测能力，这与Wang等^[1]的研究结论一致。

在研究区内划定出A、B、C、D、E五个形变区，主要分布于北盘江两侧斜坡体上(图 4黑框)。将各形变区升降轨视线向形变沿竖直方向、东西水平方向分解，沿坡向合成，结果如图 5所示。竖直形变量向上(抬升)为正值，向下(沉降)为负值；水平形变量向西为正值，向东为负值。各形变区沿坡向形变量及方向示意图中箭头线段长度与坡向形变量成正比。由图 5等高线分布情况可知：1)山坡坡度相对较陡处坡向形变量相对较大；2)地势平坦地区DS点密度明显大于地势陡峭处DS点密度；3)地势平坦地区主要表现为竖直形变，地势陡峭处同时存在竖直形变和水平形变。5个形变区的潜在威胁目标如图 5所示：1)形变区A。形变分布在斜坡体、采矿区或加工厂，形变可能与地下采矿和矿物加工的抽排水有关，主要威胁对象有工厂、居民住房、学校和道路；2)形变区B。形变分布在斜坡体，形变可能由斜坡体失稳造成，主要威胁对象有居民住房和道路；3)形变区C。形变分布在斜坡体、加工厂，部分斜坡体具有明显的老滑坡痕迹，形变可能由矿物加工的抽排水或斜坡失稳引起，主要威胁对象有居民住房、工厂、道路和河流；4)形变区D。形变分布在斜坡体，形变可能由斜坡体失稳造成，主要威胁对象有居民住房和道路；5)形变区E为露天采矿区，附近有煤炭仓储中心、煤炭加工厂和发电厂等，形变可能与采矿和矿物加工的抽排水有关，主要威胁对象有居民住房、道路和附近工厂。

4 结语

使用SBAS InSAR对升轨、降轨Sentinel-1A数据进行处理，提取贵州省水城县研究区5个明显形变区。结合Google Earth推断引起地表形变的主要原因有斜坡体失稳、地下/露天采矿和矿物加工的抽排水等，潜在威胁对象包括居民住房、学校、工厂和道路等。鸡场镇滑坡发生前升轨、降轨Sentinel-1A的SBAS InSAR形变场并未出现明显形变信息，且在滑坡发生当天及前两周内有11 d出现不同程度降雨，因此初步判断鸡场镇滑坡是由外部因素(如降雨)触发的突发性滑坡。

SBAS InSAR能有效应用于山区滑坡隐患早期识别和形变监测，但对于外部因素(如降雨)触发的突发性滑坡，卫星雷达存在重访周期较长的局限性，如Sentinel-1A为12 d，Sentinel-1A/B双星为6 d。因此，随着卫星数量增加，可通过多星组网构建卫星星座，缩短卫星雷达的重访周期，提升突发性滑坡隐患的早期识别和形变监测能力。本研究在贵州水城县滑坡隐患调查中所用技术、方法及实验结论可为贵州省以及中国西南山区滑坡隐患调查与早期识别提供参考。

致谢: 感谢英国利兹大学Andy Hooper教授提供StaMPS软件，感谢中山大学蒋弥教授提供FaSHPS软件，感谢欧空局提供Sentinel-1A数据，感谢美国航空航天局提供SRTM DEM数据。