基于短基线InSAR技术的滑坡监测研究——以四川茂县为例

彭艺，吴云龙，杨玉，袁华清

(1.中国地震局地震研究所，中国地震局地震大地测量重点实验室，湖北 武汉 430071；2.防灾科技学院，河北省地震动力学重点实验室，河北 三河 065201)

1 引 言

我国疆土辽阔、物质资源丰富，全国各地存在着各种各样的地质条件，高原、山地、盆地、林地等各种地形遍布全国。这些各不相同的地质背景，加上各地不同的气候影响，致使各种地质灾害的发生，这大大影响着我国的民生建设和经济发展。而其中，地面沉降作为一个典型的环境地质灾害，不仅反映了沉降区域的地表变化趋势，也表征了相应区域的地质变化规律，同时也可以在此基础上探究其他灾害，如地震、泥石流等的成因和发展趋势。因此，合理利用现有的科学监测方法，选定较为典型的研究区域进行监察和预测，同时对监测结果进行分析并制订相应的应对措施，是一项十分重要的研究。传统滑坡沉降的监测主要依赖于人工勘测例如三角高程测量，即使是使用了GNSS等手段获取数据，仍然存在着小范围、低精度、高成本以及分辨率低和数据处理周期长等各种缺陷[1]。随着卫星遥感领域的进步和发展，基于遥感手段的滑坡沉降监测越来越成为一种重要的方法。

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)[2]技术是一种星载测量技术，相比于传统的测量方式，将Sentinel-1卫星影像数据经相干处理得到沉降区域的地表形变信息，从而直观地表征研究测区的地表沉降情况，具有全天候、全天时、覆盖广泛、分辨率高和精度高的优点，目前在山体滑坡、三角洲地区沉降[2]、地震火山、矿区沉降以及城市地下水开采引发的地面沉降监测等各方面已经得到广泛应用。长安大学的赵延岭利用干涉图堆叠技术和SBAS技术对长江三峡区域的滑坡进行了探测与识别，获取形变信息，分析滑坡成因，做出合理预测，整体效果比较理想[3]。InSAR技术主要分为两种：传统的差分干涉测量方法(D-InSAR)和时间序列差分干涉测量，其中D-InSAR技术的构建原理是对覆盖同一地区的两幅干涉图像(一幅形变前的干涉图像，另一幅形变后的干涉图像)进行差分处理以获取目标区域地表形变[4，5]。在早期的地表形变监测中，D-InSAR技术发挥着非常重要的作用，但是随着科学技术的发展，其对时间和基线去相干过程中产生的误差、干涉图噪音、大气效应和外部DEM误差的消除越来越无法达到精度[6]，为了更好地提高影像的相干性，在此基础上提出SBAS-InSAR小基线技术，选择最短基线对影像对进行配准，提高了数据处理的精度，同时小基线技术对于DEM的数据精度要求相较于其他InSAR手段并不高，是目前应用较为广泛的一种地表形变监测手段[7，8]。

为研究SBAS-InSAR技术在地面沉降监测中的应用，本文以四川省茂县为例，利用SBAS-InSAR技术研究该地山体滑坡的发生以及对该地地形产生的影响，进而结合当地地质气候条件推测我国西南地区的地形变化，再而分析我国全国范围内山体滑坡的趋势和影响，对经济文化发展的影响也可以因地制宜地相应得到预测和解决。

2 研究区域概况

四川省茂县是四川省阿坝藏族羌族自治州的一个下辖县，该县处于四川省的西北部，地质上处于青藏高原和川西平原的交界地，地势变化大，总体呈现西北高东南低的趋势。茂县地貌以高山峡谷为主，海拔高低悬殊比较大，气候受到西风环境和印度洋西南季风影响，局部气候较为复杂，多种气候状况都有。常见各种气候灾害，再加上位于岷江上游，水资源丰富，春夏季节经常会有暴雨、冰雹、泥石流等自然灾害发生，同时也位于地震多发区，极端天气相对较多，出现大型滑坡的情况也相对较多。

除当地的地理气候环境容易导致山体滑坡等的自然灾害的发生，茂县的地质条件特殊，这也决定着当地地质灾害的发生。根据青藏高原的整体地势走向，自西向东将青藏高原东缘分为青藏高原地貌区、龙门山高山地貌区和四川盆地地貌区[9]。而导致这些地貌产生的原因是地下的地质构造，青藏高原东缘地质构造自西向东分别是松潘-甘孜褶皱带、龙门山造山带和四川盆地[10]。本文选择的四川省茂县位于青藏高原地貌区和龙门山高山地貌区的交界处，即松潘-甘孜褶皱带和龙门山造山带的交界[11]，周围发生过多次较大的历史地震。地质活动活跃，地形起伏大，再加上岷江上游水源的冲刷，这里形成陡峭的峡谷地势。此外，茂县的地层岩石分布着砂岩、板岩、泥页岩和碳酸盐等，土质易变质变形，在地下构造变化的过程中极其容易受到受影响，岩石易碎，容易产生地面沉降，如果遇到雨水冲刷，山体滑坡和泥石流无法避免[11，12]。再加上地势崎岖不平、交通不便，在自然灾害出现后经常无法第一时间进行救援，而这里土地、矿产和生物资源极其丰富，历史上也是兵家必争的要塞之地，所以四川茂县的山体滑坡现象是我们必须时刻监测预警并制定合理措施的一项工作。

历史上，四川茂县发生过两次十分严重的山体滑坡灾害，对这里造成了非常大的影响。1933年8月25日，四川省茂县叠溪发生的7.5级强烈大地震诱发山体滑坡，几乎将这个位于四川省西部的羌族聚居区全部掩埋，只剩一座破败不堪的城隍庙，据不完全统计，死亡人数约有 2 500人，造成严重损失。2017年6月24日，茂县叠溪镇新磨村的山体突发高位垮塌，被埋100多人，还造成了 2 km的堵塞河道，1 600余米的道路掩埋，以及光缆受损有 3 km长的损害，基本中断基础通信。经初步分析，这是一起由强降雨诱发的山体滑坡，再加上茂县地处地震断裂带以及汶川地震对山体的影响，在连日降雨的情况下发生了山体滑坡。

3 研究方法

针对四川茂县的山体滑坡监测，本文采用合成孔径雷达(Interferometric Synthetic Aperture Radar，InSAR)的技术进行测量，相对于一些需要人工进行操作的三角高程测量、GNSS测量等，InSAR技术具有全天候、全天时、覆盖广泛、分辨率高和精度高的优点，是目前对山体滑坡进行监测精度最高也最合理的方法。在2002年，Berardino等人在InSAR技术的基础上提出了小基线技术，获取地表形变以研究大尺度上的形变[8]，后经Hopper等人逐步优化，逐步提高了传统InSAR技术时间采样率，减弱空间失相干造成的影响和地形对差分的影响，增强了形变监测的可靠性，在地表形变监测中发挥着非常重要的作用。

本文选取了2017年3月15日～6月19日期间四川茂县以及周边区域的8景Sentinel-1A数据，利用Hooper[13]提出的SBAS-InSAR数据处理模型进行处理，对影像数据进行匹配，然后进行干涉处理、滤波处理和解缠，得到最小基线构成的三角网后，再进行反演得到地表形变的速率，最后提取时间序列曲线，并进行时序分析，找到滑坡和疑似滑坡的区域。

3.1 算法模型

在进行反演之前，要先对影像的相位进行计算。根据卫星影像获取的N+1幅SAR影像，形成M幅干涉图，生成干涉相位，再利用外部数据去除地形相位影响，根据最小基线原则，利用三维相位解缠的方法进行解缠，之后进行滤波处理以消除各种影响。

在相应时间间隔内，使用同一传感器获取覆盖同一地区的N+1幅影像，根据集合内基线距最小、集合间基线距最大的原则，通过给定阈值，可生成干涉图M幅(M和N满足(N+1)/2≤M≤[N(N+1)]/2)，第i幅干涉图任意像元(x，y)处得到的干涉相位可表示为：

фobs，i(x，y)=фtopo，i(x，y)+фdefo，i(x，y)+фatm，i(x，y)+фnoise，i(x，y)

(1)

其中фobs表示缠绕相位，фtopo表示地形相位，фdefo表示DEM误差引起的相位，фatm表示大气噪声引起的相位，фnoise表示一起热噪声等其他相位，这些相位都可以通过滤波来去除。

之后再经过图像匹配、干涉处理、去除地形相位、解缠计算、滤波处理后进行最小二乘计算，得到地表形变速率。其中最小二乘计算的公式如下：

AVlos=φobs

(2)

Vlos=A+φobs

(3)

(4)

提取形变时序后，在各相位点上联立方程组，利用最小二乘进行解算。利用式(2)求矩阵A+，使用的是矩阵的奇异值分解法(SVD)求其广义逆，即A=USVT。再利用形变量的差和时间差得到各个时间点的形变速率矩阵即式(4)，带入式(3)可得到形变速率的矩阵[1，2]，即通过A+=VS+UT求得A+，最终求得Vlos[6，8，13，14，15]。

3.2 计算流程

本文使用SBAS-InSAR模型，对卫星影像图完成干涉处理、相位解缠和滤波处理，最后得到形变速率并分析其形变规律[16]。即先选择的主影像，将副影像与其进行配准，并进行干涉处理，参照DEM高程数据去除地形相位影响，根据小基线原则，生成小基线集干涉图，设定阈值后筛选SDFP候选点，通过估算其稳定性，选取稳定性高的点以去除估算误差，在将此时得到的形变差利用三维相位解缠的方法进行解缠；最后进行滤波处理去除噪声，提取形变时序，利用最小二乘法，同时利用矩阵的奇异值分解法建立方程组运算。具体流程如图1所示：

图1 SBAS处理数据流程图

4 实验结果

本文选取2017年3月15日～6月19日期间茂县8幅Sentinel-1A影像数据，根据茂县地理位置103°39′46″E，32°4′47″N，裁剪至合适大小，基于SBAS-InSAR模型，将各影像两两连接匹配，生成连接图和时间基线，并自动定义超级主影像为20170315，再将各像对之间的匹配之后会对相距较远的两个数据进行解缠，得到最小基线的影像对配准三角网(图2)。

图2 影像对连接图、时间基线分布图、解缠后三角网

将解缠后的影像对进行干涉处理、滤波处理和3D解缠。参照拼接好的茂县区域的DEM高程数据，选取阈值为0.2进行DEM法干涉处理，同时选择高斯滤波去除噪声提高影像对相干性，并使用最小费用法进行3D解缠。本文选择相干性较好的影像对20170315-20170327的相干系数图、干涉图、解缠后影像图进行展示(图3)，并剔除相干性不好的影像对。

图3 20170408-20170514相干系数图、干涉图、解缠后影像图

选取相干性较好的影像对20170408-20170514，选取45个GCP点进行轨道精化和重去平以提高影像对的配准精度，选取GCP点时遵循尽量均匀分布，在未发生形变区域处人工选取的原则。经两次反演，分别进行二次解缠，以达到优化干涉像对、去除DEM残余误差的目的，以此反演出地表形变速率和残余地形，以及使用大气滤波以去除大气信息、估算形变组分。将两次反演的SAR坐标系下地表信息转换至WGS-84坐标系下，通过内插和平滑的方法进行地理编码，得出地表形变平均速率，经密度分割后选取测区特征点绘制时间序列曲线分析其走势，并结合地质气候条件进行分析。

5 结果分析

将测区的地表形变速率进行密度分割(图4)，该区域最大抬升速率大约为 50 mm/a，即图中红色的部分，而且分布在这里的点非常少，说明出现抬升的地方相对来说比较少。同时也可以看出，最大沉降速率可达到 -100 mm/a，如图4所示蓝色部分分布的点较多，说明这段时间内出现了比较严重的滑坡。

图4 形变速率分布图

同时选取4个地势较高地特征点和1个地势较低地特征点绘制时间序列曲线(图5，图6)，分析监测区域的地面沉降状况。5个特征点的地表在3月～4月初整体形变十分微小，基本没有变化，甚至会有微小的地表抬升，但是从4月～5月初，出现了明显的高程值变化，特征点处呈现较陡的曲线走向，说此处发生了一定程度的山体滑坡，不过在5月中旬，5个特征点的高程变化值又集中在0的上下，说明5月中旬这里没有发生大的地势形变事件，即没有影响较大的自然灾害发生，而在5月末～6月中旬这段时间，时间序列曲线却出现了骤变，如图前4个特征点的时间序列曲线急速下降为负值，第5个点的时间序列曲线迅速上升，并且按照曲线的走向，在6月19日之后仍会按照这个趋势下降和上升，说明此段时间内此地发生了较为严重的山体滑坡，并且随着时间的增加，山体滑坡会继续发生，而且甚至可能会有更加严重的自然灾害发生，可能会对附近居民的生活以及道路产生十分严重的影响。

图5 特征点1-4时间序列曲线

图6 特征点5时间序列曲线

四川茂县的山体形变发生上述地势变化，与其所处的地势条件和该处的气候变化有着非常大的关系。这里地处四川盆地西北部，山地较多，四周分布着海拔为 1 000 m～3 000 m的青藏高原、云贵高原、巫山、大巴山和大娄山等几座高山，周围群山环绕，西北边缘是很长的龙门山脉，中间地势低矮，可明显分为边缘山地和中间盆地两部分。同时，四川盆地也是我国著名的红层盆地，其地表岩石主要是由紫红色砂岩和页岩组成，它们极易风化发育成紫色土，土壤中含有丰富的钙、磷、钾等营养元素，土壤极其肥沃，但是土层浅薄，母岩疏松，非常易于崩解，很容易在大雨天气随着雨水被冲走造成泥石流，严重的甚至可以造成山体滑坡。除此之外，这里属于亚热带季风性湿润气候，夏季平均温度24℃～28℃，极端高温达到36℃～42℃，年平均降水量在 1 000 mm～1 300 mm之间，但是年内降水分布不均，整体气候呈现冬干、春旱、夏涝、秋绵雨的分布趋势，70%～75%的雨量都集中在6月～10月的夏季，高温多雨的夏季很容易发生连续暴雨天气，由此造成山体滑坡的大型自然灾害现象。

所以正是由于四川茂县处于这样的地形和气候条件中，才会出现如图所示的时间序列的变化。图中3月份5个特征点的高程变化值基本都集中在0的上下，会有微小的形变，因为四川盆地的地表岩石易于风化，而此时正值春天，相对来说，风比较大，再加上紫色土土质疏松，容易崩解，岩石的风化以及土壤的流失都是这里地势发生变化的原因。根据2017年全年茂县月降水量(图7)分析，4月份开始进入夏天，雨水天气逐渐增多，处在山脉地区再加上土壤疏松，泥土随着雨水就从地势高的地方流失到地势较低的地方，也就产生了图中所示的前4个特征点高程变化值是负值，而第5个特征点的高程变化值则为正的情况，且明显从图中分析得到，第5个特征点的地势相较于前4个都比较低。5月中旬，根据图中茂县平均气温分布图以及降水量分布图，可以发现5月份茂县的平均温度并不算高，而且降水量也相对较少，因为5月17日～18日这段时间，茂县遭受了大风、冰雹自然灾害的袭击，所以地势没有大规模的变化。

图7 2017年茂县月降水量和平均温度图

5月中旬的冰雹等极端天气结束之后，正式进入夏季，在亚热带季风性湿润气候的影响下，气温持续上升，降水量也开始持续增加。此后降水集中，强降雨天气持续，在雨水的冲刷下，疏松的土壤开始发生滑坡，所以5月中旬以后，地势高的4个特征点的高程变化值发生骤变且为负值，而地势较高的第5个特征点则发生了抬升。因为当地地势变化大，历史地震致使茂县岩体结构发生破坏和裂缝，在持续的强降雨的冲刷下，具有单面山斜坡结构的后缘岩层被拉断[12]，发生了顺层滑动，诱发持续的泥土流失。截止到本次试验研究的6月19日，特征点持续发生沉降和抬升，并且按照时间序列曲线，在未来的时间，该地区还会按照这个趋势继续沉降和抬升，而且曲线变化很陡，很有可能发生更为严重的自然灾害，如大型山体滑坡。而事实也确实如此，2017年6月24日6时左右，茂县叠溪镇新磨村山体突发了高位垮塌，道路和村民被掩埋，是一起十分严重的自然灾害事件，对当地的居民造成了十分严重的人身和财产伤害。

6 结 论

经过上述的实验过程，对2017年3月～6月之间四川茂县的卫星影像数据，基于小基线技术，对卫星影像进行了干涉，滤波和解缠等，再进行两次反演得到其形变速率，最后对其进行地理编码，找到特征点的时间序列，通过不同时间段的高程值变化来研究该处的形变情况，并同时结合当地全年月平均降水和气温数据分析其成因，最终得到的结论有：

(1)四川省茂县在2017年3月15日～6月19日期间，地表发生了一定程度的沉降，最大沉降速率和抬升速率分别达到 -100 mm/a和 50 mm/a，且地面沉降随着降水和温度的变化而改变。

(2)3月茂县附近整体地表形变不大，直至5月中旬，出现一定幅度的沉降，到5月中旬出现停止沉降的现象，之后到6月，沉降急速增加，并且有按照这个速度持续沉降的趋势。

(3)结合2017年茂县月平均降水和温度数据，以及当地地质条件和构造，分析得出：处在地质构造带交界处的茂县地质条件特殊，岩体结构已遭受破坏，持续的降水容易诱发山体滑坡的发生。

以上就是基于小基线技术对3月～6月四川茂县卫星影像数据的形变速率相干性分析，并结合当年月平均温度和降水数据以及当地的地质构造背景分析，得出了此处可能会在将来发生大型山体滑坡的结论，整体吻合当地具体滑坡变化。依据SBAS-InSAR技术对地质活动活跃的地区进行定期监测，并结合当地地质条件和降水数据做出分析和预测，方法合理且结果可信，对城市的发展和规避一些自然灾害带来的损失具有重要意义。