DS-InSAR技术在复杂艰险山区滑坡形变监测中的应用前景研究

Study on the Application Prospect of DS-InSAR Technology in Landslide Deformation Monitoring in Complex and Dangerous Mountainous Areas

摘要: 分布式目标雷达干涉(DS-InSAR)技术既继承了传统InSAR技术观测范围广、自动化程度高等技术优势，又能提高复杂艰险山区时序InSAR相干点密度，可应用于滑坡形变监测。本文利用覆盖康定附近的Sentinel-1A升轨影像，基于KS统计检验方法和自适应滤波，提取出了位于岩屑区、裸露土坡的DS点，并进行滑坡形变分析，实验结果为：该滑坡区域在2017年3月~2018年10月的监测时间内，年平均沉降速率约为10.9 mm/y，最大沉降量达到21.0 mm；且发生较大沉降时间段主要集中于每年5~8月，与该区域集中降雨时间相符。实验表明DS-InSAR技术对复杂艰险山区进行滑坡形变监测具有重要的参考价值。

Abstract: Distributed scatter InSAR (DS-InSAR) not only has the advantage of large scope of observation and high degree of automation in traditional InSAR technology, but also can increase the coherence point density of time series InSAR in complex and dangerous mountainous areas. Therefore, it could be used into landslide deformation monitoring. This paper used Sentinel-1A ascending images covering some areas of Kangding city as the experiment data. DS points located in debris and bare soil slope area were extracted by applying KS statistical test method and adaptive filtering. And then landslide deformation analysis is carried out. The result shows that the average annual subsidence velocity of the landslide area is about 10.9 mm/y and the maximum settlement is 21.0 mm. The large subsidence period is mainly concentrated from May to August each year, which coincides with the time of concentrated rainfall in this area. Experiment shows that DS-InSAR technology has important reference value for landslide deformation monitoring of complex and dangerous mountain area.

1. 引言

复杂艰险山区山高谷深，地层岩性复杂多变，地质构造复杂、深大断裂发育，新构造运动活跃、地震频繁强烈，气候恶劣多变，致使内、外动力地质作用十分强烈，具有“显著的地形高差”、“强烈的板块活动”、“频发的山地灾害”、“脆弱的生态环境”等四大特点。在复杂艰险山区工程施工要面对极端艰险复杂的地形与地质环境，高原地区频发的地震、滑坡、泥石流等山地灾害给重点工程建设及后续运营造成巨大安全风险，本文就DS-InSAR技术在复杂艰险山区滑坡形变监测的应用前景开展相关研究工作。

卫星合成孔径雷达干涉(InSAR)技术具有观测范围广、自动化程度高等技术优点，已被证明是滑坡监测的有效手段 [1] [2]。在其基础上发展起来的永久散射体雷达干涉测量(Persistent Scatterer InSAR, PSI)技术受时空失相干和噪声影响较小，能够为形变场提供可靠的观测结果 [3]。但PSI技术只针对永久散射体点目标(如房屋、铁塔等人工建筑物)进行分析，在复杂艰险山区进行滑坡形变监测时，人工建筑物稀少，PS点几乎没有 [4]。除永久散射体以外，还存在地物性质较为一致且在一定时间内保持稳定的散射体(即分布式散射体，DS点)，能为复杂艰险山区InSAR形变监测提供更丰富的信息。本文引入DS-InSAR技术，开展复杂艰险山区滑坡形变监测，对岩屑区和裸露土坡进行分析，识别出分布式散射体(Distributed Scatter，DS点)，增加用于形变监测的点目标数量，提高滑坡形变监测结果的可信度。

2. 技术路线

在PS-InSAR技术的基础上，基于点目标的时间序列分析(IPTA)、相干点时域分析(TCP InSAR)等技术的提出，进一步推动时序InSAR技术的发展 [5] [6]。然而，对于PS点目标稀少的区域，这些算法均无法起到形变监测的作用。尤其是川藏铁路沿线山高谷深，引起侧视成像雷达阴影和叠掩，可识别的PS点几乎没有，致使PSI技术滑坡形变监测结果空间分辨率低，相位解缠效果差 [7]。Ferretti等学者于2011年提出了SqueeSAR技术，以分析识别DS点，来增加自然地表的点目标数量。SqueeSAR技术利用KS算法探测出分布式散射体，再进行建模解算和位移提取，并通过对山区SAR数据进行SqueeSAR实验证明了其有效性。为提高川藏地区滑坡识别和位移监测的能力，本文引入分布式散射体(如岩屑区、裸露土坡等)的识别分析，具体的DS点提取流程如图1所示。

Figure 1. The flow chart of extracting DS points

图1. DS点提取流程

探测DS点时，釆用KS检验方法为每个像素点确定统计同质点(Statistical Homogeneity Point, SHP) [8]。如果一个像素点的统计同质点大于所设阈值，则将其定义为DS备选点(DSC)。对DSC进行自适应滤波，逐点恢复相位时间序列，提高DSC点的相干性。通过计算相干系数从DSC点集中筛选出真正的DS点，得到较高空间密度的形变场。

KS检验基于通用分布的拟合优度检验方法，采用经验分布函数确定两像素之间是否具有一致性 [9]。设总体分布F的独立同分布样本为 $X_1,X_2,\cdots ,X_k$，其经验分布函数为：

$F_N\left(z\right)=\{\begin{array}{l}0,x<X_{\left(1\right)}\\ \frac{k}{N},X_{\left(k\right)}\le x\le X_{\left(k+1\right)}\\ 1,X_{\left(N\right)}\le x\end{array},k=1,\cdots ,N-1$ (1)

式中， $X_{\left(1\right)}\le X_{(2)}\le \cdots \le X_{(k)}$ 为从小到大重新排序后得到的序列。

对于两个不同的像素点而言，判断其振幅一致性所对应的零假设和备选假设为：

$H_0:F_N^i\left(z\right)\equiv F_N^j\left(z\right)\leftrightarrow H_1:F_N^i\left(z\right)\ne F_N^j\left(z\right)$ (2)

式中， $F_N^i\left(z\right)$ 和 $F_N^j\left(z\right)$ 分别代表两个不同像素点i和j的经验分布函数。

Ferretti等通过引入柯尔莫哥洛夫距离D_N进行上述假设检验：

$D_N=\sqrt{\frac{N}{2}}{}_i{}^{\sup }\left|F_N^i\left(z\right)-F_N^j\left(z\right)\right|$ (3)

式中D_N的累积概率密度函数为：

$P\left(D_N\le t\right)=H\left(t\right)=1-2{\displaystyle {\sum }_{N=1}^{\infty }{\left(-1\right)}^{N-1}\cdot {\text{e}}^{-2\cdot N^2\cdot t^2}}$ (4)

所对应显著性水平α为：

$\alpha =1-H\left(t\right)$ (5)

一般情况下，取拒绝域为0.05，当统计检验量大于0.05时，接受零假设，两个不同像素点i和j分布相同，为一对SHP；当统计检验量小于0.05时，接受备选假设，认为两点不是一对SHP。

确定了点P旳SHP之后，即可判定此点是否为DSC。判定过程中需设定SHP的数目阈值，SHP数目大于给定阈值才是DSC。对影像进行逐点分析之后，即可得到全部的DSC。然后对影像进行自适应滤波，提高DSC点的相干性。滤波的本质，是针对备选点进行逐点计算，提高备选点的相干性，同时提供非常好的可视化效果，去除了影像中大部分斑点噪声。依次计算每个DS备选点的相干系数值，并设定经验相干性阈值，判断相干性较高的DSC为DS点。

得到所需DS点集后，采用分步三维解缠算法进行相位解缠，即可获取DS点相位时间序列 [10]。相位解缠后，把除形变相位外的其他相位成分都归为形变提取的残余误差项。为了最终得到形变相位，需要先扣除相位值中的残差项。残差项可以分为空间相关与空间非相关部分。其中，残差项空间非相关部分可以被当作噪声，通过空间域低通滤波的方法扣除；对于残差项空间相关部分可以通过最小二乘迭代的方法和先时间域高通滤波、再空间域的低通滤波的方法扣除。最后，将DS点形变相位转化为DS点形变时间序列 [11] [12]。

3. 实验分析

实验利用覆盖康定附近复杂艰险山区的Sentinel-1A升轨影像数据，影像获取时间为2017年3月13日至2018年11月15日，使用合适的参数生成干涉影像集合，采用KS统计检验方法提取DS点集，使用分步三维相位解缠方法进行时空相位解缠，进而提取形变时间序列。经实验分析，发现朋布西乡存在滑坡隐患灾害点。

朋布西乡是四川省甘孜藏族自治州康定市辖乡，位于康定西南部，管辖面积4236平方公里，人口0.3万，S215省道经此辖域，是道孚前往塔公草原的重要路段。实验选用覆盖朋布西乡的49景Sentinel-1A升轨影像(雷达波长5.6 cm，多视处理后影像空间分辨率30 m)，影像数据基本信息如表1所示。

Table 1. Image information

表1. 影像基本信息

通过进行SLC影像的生成、轨道校正、配准及裁剪等预处理工作，裁剪后的影像覆盖范围如图2所示，实验影像覆盖范围31.2 × 14.3 km²。该区域山高谷深，有植被覆盖，使用DS-InSAR技术增加滑坡监测相干点数量，提高大面积区域滑坡隐患点探测的成功率。

Figure 2. Image coverage area

图2. 实验影像覆盖区域

以2018年1月7日影像为主影像，将其余48景影像一一与主影像干涉，生成48组干涉对，时空基线图如图3所示。选取的主影像在时间序列上靠中间位置，使所有影像时间基线合适，保证影像的干涉质量。

Figure 3. Baseline and time interval

图3. 时空基线图

实验提取DS点的方法为KS统计检验方法，主要包括以下关键步骤：1) 通过选取合适的窗口大小，运用KS统计检验方法逐个选取SHP点，搜索窗口大小为15，设定平均阈值100对点进行筛选，得到DSC点集；2) 对得到的DSC点集进行滤波，减弱斑点噪声效应，达到增强相干性的目的。设置像元相干系数阈值0.2，从DSC点集中进一步筛选得到DS点。通过KS统计检验方法共选取得到DS点89672个，相位解缠后得到相位解缠结果如图4所示。

Figure 4. Phase unwrapping diagram

图4. 相位解缠图

通过分析干涉相位的空间相关部分和空间非相关部分，逐步分离雷达侧视角误差、大气延迟和轨道误差，最终提取得到雷达视线方向(LOS向)年平均形变速率，如图5所示。从图5可以看出，提取得到的DS点几乎全是雷达扫射到的裸露土坡，为常见的分布式散射体。B区域山坡发生了沉降，年平均形变速率约为−10.9 mm/y，而其余坡体未发生类似B区域的沉降变化，标记B点为探测到的滑坡隐患点。

Figure 5. Annual average deformation rate in LOS direction

图5. LOS向年平均形变速率

Figure 6. Time-series of deformation

图6. 形变时间序列

把第一景影像作为零形变起始点，B点在2017年3月至2018年11月间的形变时间序列如图6所示。从B点形变时间序列可以看出，B区域发生较大幅度的沉降主要发生在2017年4月30日至2017年8月28日和2018年4月25日至2018年8月23日两个时间段，最大沉降量为21.0 mm，从2018年8月23日后，B区域沉降逐渐趋于平稳。已有研究资料表明，康定市多年平均降雨量804.5 mm，主要集中在5～9月，夏季降水量达621 mm [13] [14]。强降雨使得大量的地表水渗入地下，使得软弱层不断受到侵蚀，导致岩土体的抗剪强度降低；土体的重度逐渐增大导致斜坡土体聚合力的减弱；土体裂隙中产生静水压力，增加坡体的滑动推力，从而降低坡体稳定性，引发滑坡沉降。从沉降量时间序列中可以发现，发生较大幅度沉降的时间段主要分布于每年5~8月。随着雨季来临，降雨量增大，开始出现滑坡且沉降现象加剧，在经过雨季之后，随着降雨量的减少，以及滑坡区域进行自然修复，沉降现象得到缓解。

4. 结论

本文以复杂艰险山区为实验区域，利用覆盖康定附近复杂艰险山区的49景Sentinel-1A SAR影像，采用KS统计检验方法提取出DS点集，并基于提取到的DS点相位空间相关性，运用时序InSAR技术进行滑坡形变监测。实验结果为：朋布西乡存在滑坡隐患灾害点，2017年3月至2018年11月，其年平均形变速率约为−10.9 mm/y，最大沉降量达到21.0 mm；且发生较大沉降时间段主要集中于每年5~8月，与该区域集中降雨时间相符。实验表明，DS-InSAR技术对复杂艰险山区进行滑坡形变监测具有重要的参考价值。