甘肃通渭黄土滑坡二维形变时序监测

蒲虹宇，张立峰，何毅，陈宝山，陈毅，何旭

（1. 兰州交通大学测绘与地理信息学院，甘肃 兰州 730070；2. 地理国情监测技术应用国家地方联合工程研究中心，甘肃 兰州 730070；3. 甘肃省地理国情监测工程实验室，甘肃 兰州 730070）

0 引言

黄土滑坡是一种突发性极强的地质灾害，近年来，随着黄土高原地区经济发展、人类活动以及在地震或暴雨等诱发因素作用下常呈大面积高发状态[1]。传统的滑坡监测手段以野外实地调查为主，譬如GPS、水准监测，该方法可识别勘测滑坡性质、范围以及运动规律等，但需要耗费大量的时间和人力，无法快速、准确地获取区域性滑坡信息，阻碍后续工作的进展。随着遥感技术的发展，滑坡体信息的识别从长期以来的点到线，线到面的模式得到改变，提高了精度与效率[2]，其中光学影像以其视点高、视域广、时效好等特点成为滑坡体识别优先手段。然而，光学影像受云雾以及天气影响大，且只能根据滑坡体在影像上的光谱特征与纹理特性进行简单的定性分析与解译。合成孔径雷达干涉测量（interferometric synthetic aperture radar，InSAR）技术具有覆盖范围广、穿透云层、全天候运作、精度高的特点[3−4]，能够定量对滑坡体进行描述，在对大型滑坡追溯、监测与普查方面已经发挥重要作用。在滑坡灾害追溯方面，Dong等[5]采用Sentinel-1 与ALOS-2 PALSAR 数据对发生在2017 年造成80 多人遇难的新磨村滑坡灾害进行了追溯，监测结果表明此次灾害前已有明显加速阶段；戴可人等[6]采用InSAR 方法对2018 年南峪乡滑坡进行形变追溯，发现了滑前明显的累计与加速形变信息，但未获取到临滑前明显加速阶段；侯燕军等[7]利用L 波段升轨ALOS 数据，对兰州市普兰太有限公司进行了有效识别。滑坡灾前形变追溯能够分析引起形变加速的原因，进而了解了该地区滑坡成因机理，为后续监测预警以及预测工作打好基础。

文中将采 用SBAS-InSAR（small baseline sunset In-SAR）技术，融合升/降轨Sentinel-1A 数据对甘肃省通渭县常河镇滑坡灾前进行二维形变追溯，研究其灾前演变过程。结合降雨量信息分析滑坡诱因，并通过D-InSAR技术进行形变追溯佐证，有利于了解滑坡诱导因素、滑坡内部机理以及更好的进行灾后稳定性评估，为未来该区域滑坡监测与预警提供理论支撑。

1 研究区概况

2019 年9 月14 日，受多日降雨影响，甘肃省定西市通渭县常家河镇小庄村发生山体滑坡，滑坡体长约560 m，宽约800 m，平均厚度约20 m，体积约7.74×106m3，属于大型滑坡。经相关文献核查统计[8]，此滑坡灾害损失范围涉及通渭县1 个乡镇17 个村民小组 2 975 人（未有人员伤亡），农业、水利、道路、桥梁等直接经济损失 2347.2 万元。现场调查显示，在县道X087 的23 km+400 m—24 km+400 m 处地面出现大面积蠕滑造成的拉张裂缝，部分农田表面发育裂缝，常家河前进砖瓦厂工人宿舍后墙坍塌、阳坡大桥垮塌等现象[9]。

常家河镇位于甘肃省定西市通渭县西南向，见图1（a）、图1（b）则为Sentinel-1A 数据覆盖范围，其灾后影像如图1（c）所示。研究区在1718 年和1920 年发生两次大地震，通渭县东南部一带形成大面积滑坡隐患，导致这些规模巨大的滑坡体对当今甚至未来地质灾害活动仍有强烈影响[10]。该区域属陇中黄土高原丘陵沟壑区，断裂发育，分布有常家河与苦水河两大河流，因此新构造运动与剥蚀作用强烈。加上区内以粉砾质泥岩和砾质泥岩为主，河流两侧成为滑坡最发育的区域。由于连续降雨是引起黄土滑坡的重要诱因之一，而该区域降水主要集中在6—9 月，且多以连续降雨的形式出现[9]，因此每年7—9 月是该区域滑坡多发期。

图1 研究区位置Fig.1 Location of the study area

2 数据与方法

2.1 数据

在研究区和研究期间内采用欧洲航天局提供的升/降轨渐进扫描(TOPS)模式地形观测Sentinel-1A 场景图像共63 景，其中升轨数据有31 景，其获取日期为2018 年9 月1 日—2019 年9 月8 日；降轨数据32 景，日期为2018 年9 月6 日—2019 年9 月13 日。升/降轨轨道卫星主要参数见表1。在SBAS-InSAR 技术中，文中研究中还采用美国地质调查局公布的航天飞机雷达地形任务SRTM（Shuttle Radar Topography Mission）数据提供的90 m 分辨率DEM(数字高程模型)对地形相关的相位进行了去除和编码[11]。另外，采用中国资源卫星中心陆地观测卫星服务平台提供的降水数据进行了后续分析（http://36.112.130.153:7777/DSSPlatform/index.html）。

表1 Sentinel-1A 卫星SAR 数据主要参数Table 1 Main parameters of Sentinel-1A satellite SAR data

2.2 D-InSAR

D-InSAR 是具有时间与空间上形变的雷达影像做差分处理消除平地相位与地形相位以便获取微小形变信息的技术，根据所使用的雷达影像数量可分为二轨法、三轨法和四轨法[12−13]。实验采用二轨法，主要利用外部DEM 模拟地形相位，然后同两幅影像干涉获取到的具有形变的地形相位进行差分处理，获取最终的变形相位。D-InSAR 方法差分干涉图的每个像素的相位值φ可表示为：

式中：φdef——地表形变信息；

φorb——卫星轨道误差；

φatm——大气噪声；

φdem——高程误差；

φn——噪声误差。

对除地表形变信息以外的其他相位误差进行去除，最终得到的视线向地表形变则可表示为：

式中：∆D——视线向地表形变；

λ——波长/cm。

2.3 SBAS-InSAR

SBAS-InSAR 是一种利用多幅主影像获得短空间基线差分干涉图数据集的时间序列InSAR（TS-InSAR）分析方法[14−15]。主要研究思路是基于该区域合适的时空基线阈值，生成阈值内根据不同主影像自由组合而成的SAR 影像干涉对，随后对高相干点进行相位解缠处理反演出时间序列沉降信息[16]。该方法适用于非城镇区域，在传统InSAR 技术上有效降低了时间与空间失相干的影响。

SBAS-InSAR 技术的主要原理如下。首先将N+1幅SLC（Single Look Complex）影像按时间序列获取，分别为：t0,t1,t2,···,ti,···,tN。然后将符合预先设定的时空间基线阈值的影像进行配对和差分干涉处理，得到M对干涉对。此时，干涉图中每个像素的相位值的组成部分同式（1）。利用轨道信息与外部DEM 数据可逐影像对去除平地与地形效应，然后对高相干性像元的相位采用最小费用流方法进行相位解缠，再根据最小二乘法和奇异值分解（Singular Value Decomposition，SVD）法对每个像元进行形变速率与DEM 误差的第一次求解。由于大气相位的干扰，会影响结果精度，因此在第一次基线估算后，通过时域上的高通滤波（365 d）和空间域上的低通滤波（1 200 m）进行大气相位估算与删除，接着使用相同方法进行第二次形变速率与DEM 误差估算，得到最终的变形结果。SBAS-InSAR 技术在SARscape5.5 中有四个主要步骤，其具体处理流程如图2所示。

图2 SBAS-InSAR 基本流程Fig.2 Basic process of SBAS-InSAR

在连接图生成步骤中，本研究通过设置降轨数据以日期为20190104 影像作为主影像，升轨数据以日期为20181203 影像作为主影像，120 d 时间基线阈值，20 m垂直基线阈值这些参数，将其他所有影像都根据主影像进行精配准与干涉图生成，文中所使用升/降轨数据的具体获取日期与时空基线关系如图3 所示；第二步干涉处理升/降轨数据则分别生成了205 组与142 组干涉对，并采用Goldstein 滤波方法提高了信噪比，在此步骤中，编辑去除相干性不好与解缠结果不好的干涉对；第三步SBAS 技术反演中，通过选择形变大小接近0值且位于山脊与远离形变区的GCP 控制点共38 个进行了轨道精炼与重去平，并利用SVD 方法获得地表形变速率；第四步地理编码则是对SBAS 反演结果进行地理坐标投影，选择虚拟值去除可得到更精确的结果。

图3 升降轨SAR 影像时空基线图Fig.3 Spatiotemporal baseline image of SAR imagery of ascending/descend orbit

2.4 基于升降轨的滑坡二维形变获取

随着利用InSAR 技术应用于滑坡监测研究的不断深入，科研者们发现同一区域在升/降轨InSAR 监测结果存在较大的差异，甚至完全相反[17]，这是因为InSAR技术一般只能得到LOS 方向上的一维形变，很难对滑坡体真实形变特征进行描述。再加之SAR 影像存在的叠掩、透视收缩、阴影等几何畸变对LOS 向形变结果产生影响[18]，因此融合多个InSAR 监测结果，利用信息冗余，将主影像异常像元值用对应地理位置副影像正常像元值进行替代，从而对丢失信息进行补偿[19]，能在不同维度对滑坡体形变进行研究，获得更为准确的形变特征，满足滑坡监测的工程需求。假设朝向雷达卫星方向的形变为正，远离雷达卫星运动方向的形变为负[6]。根据雷达卫星成像几何可知，垂直方向dh、南北方向dn以及东西方向de可构成雷达LOS 向的形变[18]，如式（3）所示。

式中：θ——雷达视线方向的入射角/（°）；

α——雷达卫星航向角/（°）。

由于常规InSAR 对垂直方向形变监测最为敏感，东西方向形变次之，而南北方向对于LOS 向形变量的贡献最小，因此视线向上的南北方向形变分量可以忽略不计，所以卫星LOS 向形变量与InSAR 监测的形变量的关系可写为:

3 结果与讨论

3.1 SBAS-InSAR 形变追溯

文中基于上述SBAS-InSAR 技术，获取了甘肃通渭县常河镇从2018 年9 月—2019 年9 月升/降轨视线方向上的形变，并基于苦水河建立600 m 缓冲区对形变结果进行裁剪（图4）。图中，形变为负值表示地物位移远离卫星视线方向，用蓝色代表；形变为正值则表示地物位移靠近卫星视线方向，用红色代表；黑色椭圆则圈定出滑坡体所在区域。升轨向形变场结果显示，升轨数据获取的滑坡体最大年平均速率为−24.47 mm/a；降轨形变场结果显示，降轨数据获取的滑坡体最大年平均形变速率为18.63 mm/a。升/降轨数据由于不同的观测几何，其获取到的雷达影像的后向散射值不同，进而导致其探测到的形变范围与形变量级有一定的差异。对于该滑坡体，升/降轨数据获取的形变范围与形变量级差异较小，且都显示出西南走向的形变几何，表明两个轨道结果相干性好且吻合度高，可进行升/降轨数据联合分析，对单轨道探测不完全的区域进行互补。根据该滑坡体的DEM 与升/降轨LOS 形变结果，可以推测其滑动方向大致沿东西向分布，因此利用东西向形变结果与垂直向形变结果能够更好的反应此次滑坡形变过程。基于上述二维分解方法，联合升/降轨LOS 向形变结果获取了常家河镇区域的二维形变。二维形变图东西向形变场显示，该滑坡体最大年平均形变速率为29.94 mm/a，整体向东移动；垂直向形变场显示该滑坡体最大年平均速率为−18.25 mm/a，垂直向下形变。两形变场方向对比分析显示，该区域形变量级大的地方多分布于苦水河附近，且东西向的形变量级大于垂直向。

图4 SBAS-InSAR 不同方向的形变结果Fig.4 Different direction of SBAS-InSAR deformation result

为了更好的反应该滑坡体运动过程，用滑坡体大致边界裁剪出东西向与垂直向InSAR 结果，做出形变时序图（图5），其中越呈现红色的区域代表沿着该方向抬升的形变越大，越凸显蓝色的区域则表示沿着该方向沉降的形变越大。二维时序图最终结果显示，在2018 年9 月到2019 年9 月这一年里，该滑坡体最大累计东向形变量为33.53 mm，最大累计西向形变量为18.36 mm，最大累计上升形变量为10 mm，最大累计下沉形变量为15.99 mm。其中，西向最大形变分布在滑坡体最北端，苦水河右侧区域，东向最大形变区域分布在县道X087 与苦水河中间区域，即苦水河左侧与县道右侧区域，可以说东西向最大形变量的分布规律符合苦水河侵蚀河岸所造成的坡体下滑现象。垂直上升形变区分布在滑坡体周边区域，以滑坡体后缘为主，而该区域附近存在落水洞、地下暗河以及少数裂缝，连续降雨沿裂缝、落水洞等集中通道快速下渗致使滑带界面处出现泥化、软化现象，并逐步促进落水洞等与滑带及地下水的连通。随着时间的推移，降雨入渗产生的地下水不断集聚，出现滞水现象，致使地下水位上升[7,21]。文献[20]指出，地下水位与地面沉降呈正相关。因此，该滑坡体后缘区域因为地下水位的上升，导致在垂直向上的形变出现灾害发生前的抬升现象。垂直下沉形变区域分布在苦水河附近，且距苦水河越近蓝色越明显，可以说垂直向的形变分布规律也满足苦水河侵蚀作用造成的坡体下移现象。王浩杰等[9]将此次滑坡定性为推移-牵引式滑坡，根据灾害发生前垂直向形变结果来看，抬升最大累计量区域（10 mm）位于坡体后缘，滑坡下沉最大累计量区域（15.99 mm）则位于苦水河西侧坡体前缘，满足推移-牵引式滑坡灾害发生规律。对二维向形变进行交叉分析，该滑坡体最先发生形变的区域处于县道X087 与苦水河之间的耕地上，这可能与降水以及灌溉引起的水位变化相关[22]。随后，在东西向时序图上显示，分布在苦水河附近的区域开始形变且形变量级逐渐与耕地上的形变量级融合成为形变中心；垂直向时序图上显示，分布在苦水河附近的区域形变量级则逐渐取代耕地上的形变量级成为形变中心。在滑坡发生前，垂直向与东西向形变较大的区域基本介于县道X087 与苦水河之间，且两方向上的最大累计形变都沿苦水河分布，分别可达到−15.99 mm 与33.53 mm。为验证苦水河是否对该区域形变有影响，选取县道X087 与苦水河距离相近的三条水平线对垂直方向上20190908 期SBAS结果进行剖面分析（图6）。由图6 可以明显看出，剖面图具有以苦水河为中心的沉降漏斗，距其48 m 缓冲区内累计沉降量都小于0 mm；在县道X087附近区域形变都较小，距离苦水河120 m 时，区域内形变值大幅度降低。因此，可以推测该滑坡发生的诱发因素之一是苦水河对河岸的侵蚀作用。

图5 二维形变时序形变图Fig.5 Sequence diagram of two-dimensional deformation

图6 垂直向剖面图Fig.6 Vertical profile

据报道，此滑坡发生前有连续降雨现象。徐辉等人[23]研究表明滑坡灾害发生前10 d 特别是前5 d 的降雨是导致滑坡发生的关键因素，这与王浩杰等[7]研究的此次通渭滑坡结果相吻合，表明降雨可能是此滑坡的诱发因素。然而，降雨是否对此次滑坡有长时间序列的影响未知，因此文中选取二维时序图上远离苦水河的四个点的时序统计结果与月降水量进行分析，探究月降雨对此滑坡的影响。按照近苦水河的垂直距离的大小依次对四个点命名，a、b、c、d 四点距苦水河的垂直距离分别为395 m、146 m、111 m 以及103 m，其统计结果如图7 所示。根据降雨量多少，将研究区累积形变时间段以2019 年5 月为分界点分为两个阶段，2018 年9 月—2019 年5 月降雨量较少将其标记为Ⅰ阶段，Ⅱ阶段（2019 年5—9 月）则为降雨量较多月份。从东西向形变结果折线图与时序点可以看出，Ⅰ阶段形变较为稳定是形变累积阶段；在2019 年5 月后的Ⅱ阶段形变加速，四个点的形变斜率都普遍增大，以距离苦水河最近的d 点形变加速最为明显。从垂直向形变结果与时序点可以看出，Ⅰ阶段中，b、c 两点在2018 年9—10 月形变斜率大于另外两点，在2018 年11 月—2019 年2 月四个点形变都较为平稳，2019 年3—4 月，b、c、d 三点形变斜率增大，而在2019 年5 月之后，随着降雨量的增加，四个点形变呈现出直线的加速阶段。据王浩杰等[7]的研究，c 点处于的台缘区域横向裂缝和陡坎十分发育，土壤裸露，遇水更容易软化，水土保护能力低于位于耕地处的a、b、d 三点。我们知道裂缝、落水洞处区域，连续降雨会引起地下水位抬升，进而引起向上的形变，c 点虽然有着地下水位抬升引起的抬升，但其位于的台缘区，有临空条件，裸露土壤软化导致的下滑趋势更大，因此c 点在垂直向中Ⅱ阶段与其他三点呈相反趋势加速。时序统计结果表明，此滑坡体的灾前形变过程经历了两个阶段：形变累积（缓慢形变）阶段和加速形变阶段。将加速形变阶段同该时段降雨统计结果交互分析，2019 年5—9 月降水量增加的同时，相应的滑坡变形也明显增大，呈现明显相关性。因此，可以说降雨为此滑坡发生的诱发因素之一。李鹏程[24]研究表明，通渭县内河流为降雨型补给河流，汛期水量大且集中在6—9 月，来势凶猛，对河岸的侵蚀作用更加强烈。因此，在垂直向上累积沉降量的形变累积阶段，出现距离苦水河进的b、c 两点形变斜率变大的现象（2018 年9—11 月），这也进一步说明苦水河对河岸的侵蚀作用是此滑坡发生的诱发因素之一。

图7 时间序列形变与降雨量交互结果Fig.7 Interaction between time series deformation and rainfall

3.2 内符合精度验证与D-InSAR 结果佐证

SBAS-InSAR 结果显示，此次滑坡灾害并非偶然，该区域在滑坡发生前一年时间里一直存在形变。为了验证SBAS-InSAR 结果的可靠性，利用内符合精度方法与D-InSAR 技术进行验证。D-InSAR 技术中升轨数据选择时间为2019 年9 月8 日与2019 年9 月20 日，降轨数据获取时间为2019 年9 月1 日与2019 年9 月25 日。

陈有东等[2]指出，相干性系数越高，所得干涉相位的条纹越清晰，干涉相位观测量越可靠，且平均相干系数在0.5 以上就能提取到有用的形变信号。因此，文中对这一年中不同季节的干涉相对进行了平均相干系数计算，结果如图8 所示。在图8 中，不同季节的干涉对平均相干系数均在0.5 及以上，说明了文中实验结果可靠。利用D-InSAR 技术进行佐证，两个轨道的D-InSAR结果如图9 所示。因为滑坡体滑动后，形变量级大，出现未能解算出的区域，而这区域位于县道X087 与苦水河之间，与SBAS 技术追述的灾前形变最大区域相吻合。相关报道指出[6]，此次滑坡灾害致使道路出现多处裂缝与阳坡大桥区域垮塌；王浩杰等[7]研究则表明，县道X087 在灾害发生前无裂缝，阳坡大桥灾害发生前只有一条裂缝，均处于稳定状态。在文中实验，SBAS 追溯的结果显示，一年间道路县道X087 与阳坡大桥区域形变较小；D-InSAR 升轨结果显示县道X087 灾后形变最大为13.37 mm，阳坡桥区域形变为8 mm，降轨结果则表明县道X087 灾后形变最大值是11.01 mm，阳坡大桥区域最大形变值是2 mm，这两种技术的监测结果均与实际情况相吻合。综上所述，此次SBAS-InSAR 结果有很好的可信度。

图8 升/降轨不同季节相干性图以及相干性系数Fig.8 Coherence diagrams and coherence coefficients in different seasons of ascending/ descend orbits

图9 升/降轨D-InSAR 形变结果Fig.9 D-InSAR deformation results of ascending/descend orbit

4 结束语

文中基于SBAS-InSAR 技术和升/降轨Sentinel-1A雷达卫星影像获取了甘肃省通渭县滑坡灾前的形变，建立二维模型分析此滑坡形变特征，探讨影响滑坡发生的因素，并用内符合精度验证方法与D-InSAR 技术对SBAS-InSAR 技术的可靠性进行了验证。结果显示：

（1）通渭滑坡从2018 年9 月起便存在形变，随着时间的推移形变量逐渐增大。至滑坡发生前东西向与垂直向最大累计形变量达33.53 mm 与−15.99 mm。

（2）二维时序累计图显示，此滑坡灾害发生前区域内的形变中心沿着苦水河分布，且形变严重区域基本位于县道X087 与苦水河之间，苦水河对于河岸的侵蚀作用可能是此次滑坡诱发因素之一。

（3）将时序结果与降雨量数据联合分析，滑坡区域内远离苦水河的时间序列形变可分成两个阶段：缓慢形变（形变积累）阶段与加速形变阶段。在加速形变阶段，形变趋向与降水量变化趋向高度一致，因此可以表明降水也是此次滑坡滑动的主要诱发因素。

（4）由于二维形变的日期是2018 年9 月6 日—2019 年9 月8 日，距离滑坡发生尚有6 d 时间，因此不能得到临滑阶段明显的加速形变期，这说明基于Sentinel 数据的InSAR 技术运用在滑坡长期监测预警中可能存在时间分辨率不足等缺点，因此需结合其他雷达影像、技术等完备滑坡监测预警体系。