时序InSAR技术通过提取观测时间段内具有稳定散射特征的高相干点集合，并依据不同相位信号分量的时空分布特性来实现误差分离，从而获得经过地形改正以及轨道、大气误差去除的形变速率场，代表性技术主要有永久散射体(persistent scatterer interferometry，PSI)技术^[13]和小基线集(small baseline subsets，SBAS)技术^[14]，具体细节可参见文献[13-14，16-17]。

本文采用的时序InSAR处理流程^[18]综合了PSI和SBAS两种方法，保留了PSI技术处理的全分辨率和高精度的特点，以及SBAS通过构成短基线集来获得高相干干涉对的处理方法，同时不预设形变模型来进行参数估计，这样能获得更多相干点，从而更可靠地提取地表形变量。该处理流程如图 1所示，从缠绕的差分干涉图出发，无需多视处理，利用距离向频谱滤波找出在短时间干涉像对中失相干很小的点，即滤波后相位微失相干(slowly-decorrelating filtered phase，SDFP)点。利用低通滤波器从干涉相位观测值中分离出空间相关的相位分量，包括地表形变相位、大气湍流混合相位、部分轨道误差以及空间相关的地形误差；由视角引起的空间不相关的相位，如空间不相关的地形误差和相位中心误差，可以用残余相位与垂直基线的线性关系拟合估计去除。去除掉这两部分相位之后得到像素的去相关噪声，根据噪声水平的大小设定阈值筛选出SDFP点。对SDFP上的缠绕相位进行时空三维解缠，恢复相位绝对值。解缠后，利用高通滤波器滤除大气误差，再减去空间相关和空间不相关的误差相位，就可以解算出最终的形变量。

图  1  滑坡隐患识别监测时序InSAR处理流程

Figure 1.  Workflow of Identification and Monitoring Potential Landslide Using Time-Series InSAR

InSAR测量得到的形变量是地表三维形变矢量(包含南北方向分量$d_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}}$、东西方向分量$d_{\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{t}}$和垂直方向分量$d_{\mathrm{u}\mathrm{p}}$)在雷达LOS上的投影^[19]。以右视成像为例，视线方向形变$d_{\mathrm{L}\mathrm{O}\mathrm{S}}$可以表示为：

式中，${\alpha }_h$和$\theta$分别代表飞行指向角和名义入射角。

以卫星视角30°、飞行方向角345°的情况为例，如图 2(a)所示，假设InSAR测得的LOS形变为30 mm(缓慢变形滑坡体常见的年平均形变量)，使用3 mm作为形变测量的标准差，理论上存在着无数种不同方向的三维形变矢量LOS投影量等于30 mm。3 mm的LOS测量精度对应于垂直方向的精度约为3.5 mm，对应于零多普勒面内的水平形变测量精度为6 mm。零多普勒面方向的水平形变可以看做是其他方向形变在该方向上的投影。零多普勒方向的测量精度最高，为6 mm，换算到其他各个方向的精度如图 2(b)所示。同时可以看出在此成像几何情况下，3 mm的LOS形变测量精度转换到南北方向精度降低约8倍，可以达到23.2 mm。传统InSAR技术不具有沿卫星飞行方向形变的测量能力，所以在卫星飞行方向$\sigma =\mathrm{\infty }$。值得注意的是，3 mm的LOS形变测量精度在城市地区通常可以达到，但是在复杂山区滑坡监测中，受相干性、大气延迟等因素影响，很难达到。滑坡研究中一般比较合理的精度要求在10 mm左右^[20]，此时对应的南北方向水平形变测量精度接近80 mm。

图  2  InSAR LOS观测对滑坡三维运动灵敏度分析示意图

Figure 2.  Diagram of the Sensitivity of InSAR Line-of-Sight Observation to the Real Motion of Landslides

在山区滑坡监测中，不同坡体的形变方向差异很大。假设滑坡体沿着最大地形梯度方向(坡度$\phi$、方向角${\alpha }_{\mathrm{a}\mathrm{s}\mathrm{p}}$)发生位移，考虑到InSAR测量结果是滑坡体真实三维形变LOS上的投影，结合某一轨道的卫星成像几何参数(入射角$\theta$和飞行方向角$\alpha$)和区域地形，可以得到该轨道下InSAR形变测量结果$d_{\mathrm{L}\mathrm{O}\mathrm{S}}$与该区域真实形变$d_{\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}$之间的定量化关系：

这样就可以估算得到不同坡体上InSAR测量的灵敏度，直观地判断不同坡体形变的真实量级和危险程度，更可靠地捕捉区域内的形变信号，实现对灾害隐患区域的准确探测。

本文选择四川茂县岷江河谷典型区域作为研究区域，从沟口乡开始，沿岷江向上延伸至石大关乡，区域范围以及地形如图 3所示。茂县位于四川省西北部的阿坝藏族羌族自治州，地处青藏高原向四川盆地过渡地区，地貌以高山峡谷为主，地势从西北向东南倾斜。发源于岷山南麓的岷江由北向南贯穿四川省全境，途经松潘、茂县等地，从都江堰流入成都平原，最终在宜宾汇入长江，是长江的重要支流之一。岷江河谷两岸地势陡峭，大多呈V形。研究区域气候类型属于高原季风气候，气象条件受海拔影响显著，暴雨、泥石流等灾害集中多发于春夏季节。

图  3  研究区域地形图

Figure 3.  Topographic Map of the Study Area

茂县地区地质构造复杂，分布有岷江断裂带、松坪沟断层及石大关断层，复杂的地质构造使得区域内岩石高度破碎，地层出露以砂岩、板岩、千枚岩等软质岩体为主。研究区处于南北地震带中段的龙门山地震带和松潘地震带上，地震活动频繁。本文使用了两组Sentinel-1A/1B时间序列数据集，分别为升轨和降轨观测获取，其基本参数如表 1所示，各自的时空基线分布如图 4所示。图 4中每个圆圈代表一景Sentinel-1A/1B影像，每条黑线连接的两景SAR影像组成一个干涉像对进行差分干涉处理。

图  4  Sentinel-1数据干涉对组合图

Figure 4.  Interferometric Combinations for Sentinel-1 Data

采用时序InSAR中的StaMPS方法对覆盖茂县岷江河谷区段的升轨和降轨Sentinel-1数据进行了处理，提取了2014—2017年该地区LOS的地表年平均形变速率场，叠加在谷歌地球

(Google Earth)影像上(见图 5(a)、图 6(a))。图 5、图 6中红色表示沿着LOS远离卫星运动，蓝色表示地物朝向卫星运动。探测到的相干点主要集中在坡体中下部，靠近坡顶的区域受陡峭地形和植被覆盖的影响相干点较为稀疏。从形变的空间分布来看，升轨数据集探测到的形变区域较多，大多集中在岷江河谷的西岸山坡上，最大LOS形变速率超过50 mm/a。降轨数据集中只有擦耳岩村的区域形变较为明显。对比同一区域不同轨道观测获取的形变场，不难看出不同观测几何的数据在同一山区存在明显的灵敏度差异。采用前述的滑坡形变监测灵敏度估算方法，结合参考地形数据，计算得到该区域InSAR LOS形变测量灵敏度分布(见图 5(b)、图 6(b))。升轨数据处理结果中箩兜寨到两河口一线密集分布的形变区在降轨数据集中几乎探测不到，这说明降轨观测数据对于该区域的形变测量灵敏度几乎为零。结合探测灵敏度分布图可以更准确地解译InSAR结果，判断高风险形变区域。

图  5  Sentinel-1升轨数据集时序InSAR形变测量结果

Figure 5.  Time-Series InSAR Results of Ascending Sentinel-1 Data

图  6  降轨数据时序InSAR形变探测结果

Figure 6.  Time-Series InSAR Results of Descending Data

在区域形变探测结果和灵敏度分析的基础上，将该区域划分为沟口乡-飞虹乡、水草坪村-两河口、两河口3大块区域分别开展分析。

茂县沟口乡位于两条峡谷的交汇处，图 7给出了从Sentinel-1升轨数据集中提取的沟口乡-飞虹乡区域年平均形变速率图及其对应区域的谷歌地球光学影像。飞虹乡擦耳岩村所在的坡体、鸡公寨所在山体以及沟口乡政府对面坡体这3块区域可以探测到明显的形变信号。其中鸡公寨滑坡形变量级最大，LOS最大形变速率接近50 mm/a，而沟口乡对面坡体的最大LOS形变速率仅为25 mm/a左右。结合图 5的形变测量灵敏度图可以发现，鸡公寨滑坡在升轨数据集中的形变测量灵敏度在0.8以上，升轨数据集的观测几何非常有利。而茂县沟口乡政府对面的坡体朝向偏南北向，在图 5(b)中该区域的灵敏度仅为0.4左右，实际沿坡面方向的形变速率远比InSAR结果的LOS形变速率大。

图  7  沟口乡-飞虹乡区域Sentinel-1降轨数据形变情况

Figure 7.  Deformation Map from Ascending Sentinel-1 Data in the Region from Goukou Township to Feihong Township

飞虹乡政府对面擦耳岩村所在的斜坡也是呈近南北朝向，图 7(a)的升轨数据集处理结果中LOS形变速率为10 mm/a左右，但是其对应的形变测量灵敏度仅为0.2左右。图 8给出了擦耳岩滑坡区域降轨数据集的形变结果，在降轨数据集中的灵敏度约为0.4~0.5。降轨数据集LOS最大形变速率超过50 mm/a，对应实际沿坡面运动速率超过100 mm/a。图 8中同时给出了该斜坡上一个点对应的降轨时间序列形变量，可以看出该区域时间序列形变呈现周期性的变化，2015年和2017年都出现过形变放缓的趋势。

图  8  飞虹乡擦耳岩滑坡Sentinel-1降轨数据形变图

Figure 8.  Deformation Map from Descending Sentinel-1Data of Ca'eryan Landslide in Feihong Township

图 9(a)是高分辨率光学卫星影像中擦耳岩滑坡的大致形态，图 9(b)是从茂县飞虹乡政府所在位置拍摄的擦耳岩滑坡体的现场实景照片。从图 9中可以看出整面山坡呈现圈椅状的典型滑坡地貌特征，斜坡上部坡度较陡，下部较缓，坡体上存在大量小型崩塌、落石留下的痕迹。滑坡体后缘可见滑动后形成的陡壁，前缘受岷江河水长期冲刷，坡度较陡，形成临空面，存在小型滑塌堆积挤压河道的现象。

图  9  飞虹乡擦耳岩村斜坡

Figure 9.  The Slope in Ca'eryan Village, Feihong Township

茂县水草坪村到两河口之间的河谷西侧山体存在连续几个大型的不稳定斜坡。这一段岷江河道呈西北-东南走向，切割山体形成的坡体朝向与升轨Sentinel-1数据的雷达视线方向近似平行，成像几何有利于观测，探测灵敏度均在0.8以上，形变速率图如图 10所示。其中位于水草坪村和箩兜寨的两个滑坡在雷达LOS形变量最大，均超过50 mm/a，并且滑坡边界在形变速率图中较为清晰，该区域岷江东岸的斜坡无明显的整体变形信号。

图  10  水草坪村-两河口区域哨兵1号升轨数据集形变速率图

Figure 10.  Deformation Velocity from Descending Sentinel-1 Data in the Region from Shuicaoping to Lianghekou

图 11是岷江西岸坡体形变探测结果的放大图，其中A、B、C、D 4个区域分别取自4个不同的滑坡体且4块不稳定的区域是同一片山坡被沟壑分割而成的不同坡体。从2014—2017年Sentinel-1升轨数据集观测到的形变结果可以看出，B所在的滑坡体面积和LOS形变速率最大，滑坡后缘位于距离山顶约1/3处，前缘临江，长期受岷江河水冲刷，最大形变速率超过50 mm/a，出现在靠近坡脚的位置。A区域的滑坡体在面积和形变速率方面都比B区域滑坡小，前缘不临江。水草坪村和箩兜寨分别位于A和B所在的滑坡体上，受到滑坡灾害的直接威胁。

图  11  水草坪村和箩兜寨滑坡形变探测结果放大图

Figure 11.  Enlarged View of the Deformation Map inShuicaoping Village and Luodou Village

C和D所在的两个滑坡体位于箩兜寨滑坡的上游，靠近两河口位置，其中C区域滑坡形变速率较大，超过50 mm/a，并且滑坡前缘临江。D区域的滑坡体相比于C区域形变范围和形变速率都更小。C、D两个滑坡体区域附近均没有居民点，唯一的威胁对象是沿着岷江河谷伸展的213国道。需要注意的是，C、D区域在升轨数据集中的形变探测灵敏度比A、B区域低，在升轨LOS相同大小的形变速率情况下，C、D区域的真实沿坡面形变速率更大。

图 12给出了A、B、C、D对应位置测量点集合的平均时间序列形变量。这4个区域2014-10—2017-07的累积形变量均超过100 mm，其中B区域累积形变超过150 mm，在2016-01前后形变放缓，随后加速，在2017-05之后又再次出现加速；A、C区域大体上呈现线性形变特征；D区域则是在2015-03之后出现短暂的形变减缓。

图  12  水草坪村和箩兜寨A、B、C、D 4个区域时间序列形变

Figure 12.  Time-Series Deformation of Regions A, B, C and D in Shuicaoping and Luodou Villages

2019年11月初，对B区域所在的箩兜寨滑坡进行了实地考察，图 13中两张照片分别是在图 11白色三角形标注的P₁和P₂点处拍摄的。P₁点靠近箩兜寨滑坡的顶部，可以清晰地观察到滑坡后缘因坡体长期下滑形成的陡坎。在P₂点处，山坡上通向箩兜寨山顶的水泥道路有多处受坡体局部下滑推挤影响出现明显的破损痕迹。

图  13  箩兜寨滑坡现场考察照片

Figure 13.  In-situ Photos of the Landslide in Luodou Village

两河口区域位于多条峡谷交汇处，通往黑水县的302省道在黑水河汇入岷江处与213国道汇合。该区域的不稳定斜坡主要从升轨数据形变探测结果中识别得到，图 14给出了两河口区域的升轨数据LOS形变速率图。

图  14  两河口区域Sentinel-1升轨数据集形变速率图

Figure 14.  Deformation Velocity from Ascending Sentinel-1 Data in the Region of Lianghekou

除了黑水河北岸白布村所在的南北向坡体外，岷江两岸山体的坡脚处分布着一些小型滑坡体。形变的范围和形变速率都相对较小，Sentinel-1升轨数据LOS形变速率大约在30 mm/a左右。白布村所在的山坡呈南北向，同时坡体表面沟壑多，局部地形十分复杂，形变探测的灵敏度差异也较大，所以在图 14的形变测量结果中表现为一部分区域朝向卫星运动，一部分区域远离卫星运动。此外，在时序InSAR结果中，白布村滑坡体的大片中心区域无法提取出有效的测量点。根据现场考察的情况看，该坡体上植被覆盖很少，推测这部分地区形变梯度很大，导致了严重的失相干现象。

四川茂县岷江河谷地区受地质条件、地形地貌、地震活动、强降雨等多种因素的影响，发育分布有众多的滑坡灾害隐患，对当地居民的生命财产安全和国道省道等基础设施的安全运营构成了严重威胁。本文收集了46景升轨和42景降轨的Sentinel-1数据，利用时间序列InSAR技术在岷江河谷沟口乡到石大关乡一线区域开展了滑坡隐患识别，共识别出岷江及其支流黑水河沿岸大小滑坡隐患20余处。同时对不同区段的形变测量结果进行了详细分析，其中针对飞虹乡到箩兜寨一线的多个滑坡分别提取分析了时间序列形变，并实地开展了现场调查核实，验证了InSAR探测识别结果的有效性。

本文研究结果表明，利用时序InSAR技术在坡度和坡向各异的复杂高山峡谷地区开展滑坡隐患识别和形变监测时，需要综合考虑SAR数据的成像几何条件与区域地形特征的共同影响分析不同轨道的InSAR测量灵敏度，避免造成对于滑坡体真实形变量级的误判。