应用D-InSAR技术识别分析通渭县李店乡滑坡群灾害

王睿博，甘淑，张荐铭

(1.昆明理工大学国土资源工程学院，云南 昆明 650093；2.云南省高校高原山地空间信息测绘技术应用工程研究中心，云南 昆明 650093)

1 引 言

滑坡是中国常见的地质现象，每年由于滑坡灾害造成的损失不计其数。而滑坡群则是大范围的滑坡，其灾害影响更胜于滑坡。2019年9月14日～9月16日，受连续降雨影响，甘肃定西市通渭县常家河镇李家店乡发生特大型滑坡群灾害(图1)。虽未造成人员伤亡，但是灾害造成的农业、企业、供电、通讯、道路、桥梁等经济损失约 1 994.7万元。

图1 无人机航拍影像

合成孔径雷达差分干涉测量(D-InSAR)是指利用同一区域的两幅干涉图像，其中一幅是形变前获取的干涉图像，另一幅是形变后获取的干涉图像，通过差分处理(除去地球表面，地形起伏)来获取地表形变的测量技术。其作为一种先进的地面监测技术，与常规方法相比，在识别大范围地区的地表微小形变方面有着巨大优势，精度最高可达毫米，实现方法主要有3种：二轨法，三轨法以及四轨法。目前，已有许多相关学者将此技术应用于形变观测，文献[1]利用D-InSAR技术监测京津城际高速铁路沉降情况，发现了沉降漏斗，并为铁路方案的比选提供重要依据。文献[2]以青藏铁路拉萨段为试验区，使用D-InSAR技术监测了大范围冻土形变，发现其形变情况非常吻合冻土冻胀融沉的物理变化规律。文献[3]利用D-InSAR技术监测了九寨沟震后高精度形变情况，快速高效地获取震损物源的分布与数量特征。文献[4]利用D-InSAR技术监测了河北峰峰煤矿地区地表沉降情况，其利用多模式雷达，获取了该矿区的沉降量，精确地反映了矿区的地理信息。文献[5]利用D-InSAR技术监测米林滑坡，获取滑坡形变值，证明了D-InSAR技术用于滑坡早期识别与动态监测的可行性。但D-InSAR仍存在许多局限，如受到影像个数、DEM精度、时空间基线等影响，均会造成形变结果精度不同。本实验则利用二轨D-InSAR技术对此滑坡群形变进行观测，并根据形变结果分析不同精度DEM下的滑坡范围，判定灾害原因。

2 研究区概况及数据来源

研究区(图2)位于通渭县李店乡，其地形山高坡陡，土地贫瘠，自然条件严酷。李店乡是一个干旱乡，年降水量 380 mm左右，年蒸发量则在 1 500 mm以上。根据文献[6]研究记载，通渭县受到了1718年通渭7.8级大地震和1920年海原8.5级大地震的影响，在境内诱发了大量的黄土地震滑坡，现有的地形由大面积的滑坡体组成，结构较为松散。从土质方面来看，该地区与文献[7]研究区域地理位置相近，土质为黄土，一般呈浅黄色，具大孔隙，垂直节理发育，湿陷性强，易崩塌，黄土溶蚀微地貌发育，存在不稳定性[6]。黄土湿陷系数随着天然含水率、饱和度、干密度、液限、塑限、塑性指数的增大而减小，而随着孔隙比、压缩系数的增加呈增大的趋势[7]。

本次研究选取的数据来源于哨兵1A(Sentinel-1)卫星，模式为干涉测量宽幅(IW)的SAR影像数据，并使用了与SAR影像相匹配的精确轨道数据(表1)。

同时，本实验采用了不同分辨率精度的DEM，分别是 90 m，30 m以及 12.5 m(图3)，用以比较不同精度下形变量范围。根据文献[8]结论，90 m与 30 m精度的DEM选用SRTM数据——主要是由美国航空航天局(NASA)和国防部国家测绘局(NIMA)联合测量，由SRTM数据中心(http：//srtm.csi.cgiar.org/srtmdata/)获取。12.5 m精度的DEM则是来自于日本的ALOS卫星，其搭载的PALSAR传感器所制作的DEM，精度约为 12.5 m，可满足此次实验需求，由欧空局网站(https：//search.asf.alaska.edu/#/)获取。

Sentinel-1干涉对参数表 表1

图2 研究区范围

图3 3种精度DEM

3 D-InSAR处理过程

本次D-InSAR技术应用实验主要采用二轨法，其基本思想是利用已知的外部DEM来去除地形引起的干涉相位，其关键步骤是外部DEM与主影像的配准[9，10]。因此，本实验利用SARscape软件，运用不同精度的DEM来进行D-InSAR配准，以提高实验精度，判定滑坡区域。二轨D-InSAR处理流程如图4所示。

图4 D-InSAR处理流程图

3.1 基线估算与干涉图生成处理

(1)基线估算

基线估算主要用于评价干涉相对的质量。此步骤主要计算空间基线、时间基线、多普勒偏移等系统参数，用以判断获取的SAR数据影像是否能够做干涉，这是D-InSAR技术进行的基础。空间基线是两个传感器之间相对位置的向量，可以将其分解为两个基线分量：平行基线和垂直基线。其中，垂直基线是判定影像能否做干涉的最重要的一个参数，若其超过了临界基线最小或最大值时，两幅图像则无相位信息，没有相干性，即无法做干涉。垂直基线的长短也反映了干涉测量系统对高程的敏感程度[9]。

本次实验中，将两景SAR影像导入SARscape软件后，对VV极化影像裁剪出研究区域，再对裁剪后的两景影像进行基线估算，得到空间基线长度为 102.356 m，临界基线的最小/最大值为 -5 452.484 m与 5 452.484 m，2π模糊位移为 0.028 m。由此可见，空间基线远小于临界基线，存在相位信息，有相干性，可以满足D-InSAR技术的处理要求。

(2)干涉图生成处理

裁剪后的两幅SLC影像共轭相乘，并分别与3幅 90 m，30 m和 12.5 m精度的DEM进行配准，此过程设置相同的制图分辨率，配准结束后分别生成了3幅精度不同的干涉条纹图。此时还无法用目视判定不同精度DEM配准的影像差别，但经过去除平地效应之后，根据去平干涉图像(图5)，已经可以大致看出，DEM精度不同，去平图像已经有很大差别，但此时去平后的干涉条纹图的噪声还异常明显，必须要进行去噪处理。

图5 部分去平干涉图

3.2 滤波相干处理与相位解缠

(1)滤波相干处理

SARscape提供了3种滤波器，分别是Adaptive滤波器、Boxcar滤波器和Goldstein滤波器。本实验均采用Goldstein滤波器，实验处理得到滤波后的干涉条纹图(图6)以及干涉图的相干系数图(图7)。相干系数图中，相干系数越大，表示相干性越高，反之，则表示变化较大，相干性较低。

图6 部分滤波干涉图

图7 相干系数图

(2)相位解缠

相位的变化是以2π为周期的，当相位的变化量超过2π时，相位就会重新开始，并不断开始循环。因此，相位解缠是对前几个步骤(去平地效应和滤波)处理后的相位进行解缠处理，即把缠绕相位解缠绕后，恢复成真实的干涉相位，目的主要是解决2π的模糊问题，使地形上的线性变化的信息与相位一一对应。此步骤也是InSAR数据处理过程中最重要和最困难的步骤。

SARscape提供了3种解缠方法：Region Growing(区域增长法)、Minimum Cost Flow(最小费用流法)和Delaunay MCF。其中，Minimum Cost Flow是默认的解缠方法，这种方法是将相位解缠问题转化为最小化问题，通过在全局范围内搜索路径和最短枝切来求得最小化问题的最优解。该方法可以应用于规则网络(网格)，也可以用于不规则网络(三角网)[11，12]。用此方法解缠时，优先采用正方形的格网，考虑所有像元，当像元的相干性小于所设阈值时，可对该像元做掩膜处理。相较而言，当解缠相对困难时，Minimum Cost Flow比Region Growing可以取得更好的结果。

通过实验对比，此处解缠选择默认的Minimum Cost Flow方法效果较好。由此得到相位解缠后条纹图(图8)。

图8 部分相位解缠后条纹图

3.3 相位精确修正与形变图示表达

(1)轨道精炼与重去平优化

轨道精炼和重去平是利用选取的GCP控制点，将干涉后的相位和解缠后的相位精确修正。在选取GCP控制点时，其原则应该远离形变区域、解缠错误的相位跃变区域和残余地形相位。因此，平地区域是比较好的选择控制点的位置。将GCP控制点作为稳定的参考点，通过计算可以得到控制点的相位误差，从而将这些误差去除，并重新修改解缠图像头文件中的轨道参数，使得相位更加精确。

轨道精炼主要有轨道优化和多项式优化两种方法，如果选择轨道优化的方法，则至少需要7个控制点，否则会自动改为多项式优化。多项式优化原理是从解缠相位中估算相位斜坡，但要求控制点个数必须与多项式次数所需要的控制点个数相对应，否则多项式次数会自动降低[13]。多项式优化虽然条件较为苛刻，但其健壮性更好，即使是在基线小的情况下也可以使用。

本实验选取Automatic Refinement(自动优化法)，让SARscape自行判定最佳优化方式。

(2)相位转形变和地理编码

此步骤把合成的相位与经过绝对校准的解缠相位相结合，转换为形变数据，再通过地理编码将形变量显示到制图坐标系统中，默认得到的是LOS方向的形变，最终得出形变结果。如图9所示，图中形变量的值可能是由于滑坡形变方向和数据的LOS方向相差较大，从而导致测得数据值与实际形变值相比总体偏小。

图9 形变区域图

4 滑坡识别应用分析探讨

4.1 不同DEM精度的形变监测滑坡识别

将图9中3种不同精度DEM形变结果分别通过ArcGIS软件将形变区域放大分析，并将3种结果叠加，提取出主要沉降与抬升区域(图10(g))，用以辅助分析滑坡。经分析，分别得到三种DEM精度研究区形变图(图10(a)、(b)、(c))。

将ArcGIS中3种形变结果处理为图层文件，转换KML格式，叠加至Google Earth经纬度相同区域，形成Google Earth叠加形变图(图10(d)、(e)、(f))，此图是将平面二维形变信息转化为立体三维信息。

图10 分析图

现结合图10中3种不同精度DEM形变结果及Google Earth地形情况对滑坡区域进行详细分析：

(1)根据90 m精度DEM形变结果，结合Google Earth叠加地形效果，可以看到滑坡群明显滑坡有5处，形变范围较为明显，形变值较其他两种精度DEM所做结果偏大。A滑坡处坡度较缓，有明显的过渡区域，形变量依次由大到小，山顶滑坡形变最大，约为 0.041 m～0.047 m，到山谷有明显抬升部分，约为 0.008 m～0.014 m。B滑坡处山体较陡，过渡区域较少，滑坡形变较为明显。C、D、E三处滑坡与B滑坡情况较为相似，都属于山体较陡区域，形变量变化速度较快。

(2)根据30 m精度DEM形变结果，测出6处滑坡，A、B、D、E、F区域与 90 m精度DEM形变结果在形变范围方面基本相似，形变量范围介于 90 m精度DEM与 12.5 m精度DEM测得结果之间，故此精度形变值可靠性较高，形变量最大处的值约为 0.035 m～0.041 m，抬升处形变值约为 0.007 m～0.012 m。新测得C处滑坡，属于小范围滑坡，其形变量变化范围不大。

(3)12.5 m精度DEM形变结果总体相似于 30 m精度DEM形变结果，但是其形变区域较为零散，形变值相较于其他两种精度DEM测得结果整体偏小，抬升区域不明显，总体判断滑坡较为困难。

4.2 滑坡识别结果分析探讨

通过分析3种不同精度的DEM测出的形变结果，可以大致判断滑坡群主体范围：李店乡东部及东北部方向，以阳山里村西部区域形变量最大，滑坡区域最密集。以 30 m精度DEM显示效果最好，既能显示出所有的形变区域，又能显示出较为突出的形变量，整体相较于其他两种精度DEM测得效果更好，此结果与文献[14]结论相同。

现利用30 m精度DEM测得形变结果，对6处滑坡进行量化分析，分别在坡顶和坡脚各选取一点并连接成线，生成此线上6幅地形剖面形变图，分布如图11(a)所示。

从图11(b)滑坡剖面图可以看出，6处剖面图均符合滑坡特征，坡顶到坡底总体呈现下沉趋势，中部有滑坡物堆积部分，形变量有明显抬升。由形变量及像元个数可以判定出A滑坡的波浪形堆积物大致是从第5像元到第30像元处，约有 0.001 m～0.004 m的抬升，总抬升量约有 0.01 m；B处滑坡堆积物存在于第7像元到第12像元区域，约有 0.005 m抬升；C滑坡堆积物存在于6～10像元，约抬升了 0.013 m；D滑坡相似于C滑坡，堆积物存在于4～14像元，约有 0.012 m抬升；E滑坡抬升区域较少，大致是由于滑坡体未完全滑下，属于整体性糯滑，分布于第5～9像元区域，抬升了 0.008 m；F滑坡堆积物最多，根据图10(b)(e)及图11(a)可以看出，此处为较大U型山谷，四周滑坡均往中部堆积，故剖面显示出较高的抬升量，约有 0.021 m，存在于约第9～26像元，跨度约有17个像元。

图11 剖面分析图

结合图10中(d)(e)(f)，滑坡最大形变处(深蓝色区域)发生在道路左侧，说明滑坡对道路影响较大，结合无人机航拍，符合现场情况。Google Earth叠加影像可以明显看出山体区域两侧均有整体性的滑动形变，根据滑坡的特点，滑坡变形范围呈现的是连续的面状体，通常形状上陡下缓，近似于圆弧形，由图11中可以看出形变部分也呈现面状特点，且滑坡变形部分(蓝色区域)均分布在坡体较陡区域，其周围区域范围内也出现了零散的抬升区域(黄色区域)，在经过三次不同精度DEM的处理后，每次处理结果均有形变，形变量及形变主体范围基本不变。这说明了此处形变在这12天里空间和时间上是连续变化的，因此实验结果符合滑坡特征。

由于滑坡还具有“大雨大滑、小雨小滑、无雨不滑”的特点，滑坡当日，结合天气情况，通渭县李店乡地区降水量达 72.6 mm，已经远超平均日降水量，随着降水入渗土地中，致使斜坡土体含水量增大，斜坡体自重也随之增大，降雨入渗到达滑坡带，使滑坡带饱和，使得土体抗剪强度显著降低，从而引发斜坡失去稳定性并蠕滑变形，最终发生滑坡群灾害。

5 结 论

本此实验以通渭县常家河镇李店乡滑坡群为观测目标，利用二轨D-InSAR技术对滑坡群进行差分干涉测量，获取了不同精度DEM下的形变值、形变范围。通过比较，最终以 30 m精度DEM测出的形变效果最好，以此判定了6处滑坡群主体部分，均发生在李店乡东部及东北部方向，以阳山里村西部区域形变量最大，滑坡区域最密集。根据形变值与形变范围，结合滑坡特点、剖面情况、无人机航拍、土质状况及天气情况，判定滑坡原因是降水量过大引发的灾害。

此次实验只利用D-InSAR技术观测了研究区12天内滑坡群的沉降值，并未研究其时序变化。因此，若要持续监测此区域，获取时间序列上的变化情况，还可以利用PS-InSAR、SBAS-InSAR等技术进一步研究。