董继红， 马志刚， 梁京涛， 刘 彬， 赵 聪， 曾 帅， 鄢圣武， 马晓波

(1.四川省地质调查院，稀有稀土战略资源评价与利用四川省重点实验室，成都 610081； 2.四川省智慧地质大数据有限公司，成都 610081； 3.四川省国土空间生态修复与地质灾害防治研究院，成都 610081)

0 引言

滑坡作为频繁发生、破坏强烈的全球自然灾害之一，对于人民生命财产安全造成了巨大的威胁，同时也制约着当地经济的发展[1-2]。四川省位于中国西南山区，地质环境条件复杂，构造发育，是地质灾害高易发区，同时也是我国地质灾害防治工作的重点区域[3]。川西片区不仅地质构造复杂，活动断裂和强震也极为发育，被多条断裂带穿越，同时该区域海拔变化剧烈，降雨充沛，使得该区域地质灾害尤为频发。2017年6月24日，四川省茂县叠溪镇新磨村发生高位隐伏滑坡，造成13人死亡，73人失踪[4-6]； 2020年8月21日雅安市汉源县富泉镇发生山体滑坡，造成7人死亡2人失踪。

面对频发的地质灾害，合成孔径雷达干涉测量技术(interferometric synthetic aperture Radar，InSAR)作为一种大范围地表形变监测的新技术，因其非接触、大范围、空间覆盖范围广、监测精度高等优势，已经广泛用于地震、地面沉降、冰川、地裂缝、滑坡等地质灾害调查与监测研究中，并取得了很好的效果[7-10]。赵超英等[11]基于InSAR技术对甘肃省黑方台区域的黄土滑坡利用多源合成孔径雷达(synthetic aperture Radar,SAR)数据进行滑坡编目，长时间序列监测，并以新塬2号滑坡体为例进行失稳模式识别研究，将识别与监测采用的方式方法进行区分与介绍; 刘星洪等[12]利用光学遥感解译技术与InSAR技术对巴塘—芒康段的地质灾害进行了调查，并通过野外地质调研、地理信息系统(geographic information system,GIS)空间分析和工程类比等工作，证明综合遥感技术方法在青藏高原高山峡谷区的应用，文中通过对光学与InSAR结果对比进行分析，讨论了InSAR技术对隐患的识别效果; 韩冬建等[13]利用小基线集(small baseline subset,SBAS)技术和干涉图堆叠技术基于不同波段的SAR数据对樟木口岸在尼泊尔地震引发的山体滑坡进行识别和形变特征研究，成功识别出了13处滑坡隐患，文中只是对不同波段的SAR数据识别能力进行了对比，其中SBAS技术只是用来监测典型滑坡隐患，并没有探讨SBAS技术对滑坡隐患的识别能力差异性。

Stacking技术利用对解缠图进行相位堆叠，可以快速获取地表形变速率，广泛用于滑坡隐患的早期识别，SBAS技术通过时间域和空间域的解算，可以有效消除对地观测过程中各种误差对结果的影响，在获取地表形变速率的同时可以获取地表形变点长时间序列监测信息。因此本文选择在植被密集区域，开展不同时序InSAR技术在滑坡隐患识别数目和准确率方面对比研究，探讨这2种技术用于滑坡隐患识别效果的优劣，并选择典型滑坡隐患变形特征进行对比分析。

1 研究区概况及数据源

1.1 研究区概况

雅安位于四川盆地与青藏高原的过渡地带，其西侧边缘地带的山地海拔高达5 000余m，东侧低山丘陵地带海拔仅1 000 余m，区内地形高差极大。龙门山断裂带、鲜水河断裂带、安宁河断裂带三大活动活动断裂带均穿过雅安地区，导致区内构造活动强烈，地震频发。同时，雅安是我国内陆降雨最为充沛的地区之一，素有“雨城”之称。极为复杂的地质条件以及充沛的降雨构成了雅安地区脆弱的地质环境，加之日益增多的人类工程活动，使得雅安成为全国地质灾害最易发的区域之一[14-15]。

1.2 数据源

本文收集了覆盖研究区的Sentinel-1A卫星SAR影像数据，其覆盖范围如图1所示。由于研究区域面积较大，升轨数据需要4个轨道数据(图1绿色框)，降轨数据需要2个轨道数据(图1蓝色框)可以完全覆盖研究区域，时间范围为2018年10月—2020年10月，单个轨道是91景数据。表1列出了本文研究中所使用的SAR数据集的基本参数。

图1 研究区域示意图Fig.1 Schematic diagram of the study area

表1 Sentinel-1A卫星参数与数据选用时间范围Tab.1 Sentinel-1A satellite parametersand data selection time range

使用AW3D30数字表面模型来消除InSAR干涉处理过程中的地形相位，同时该数据也被用来辅助SAR影像进行地理编码及计算SAR数据的叠掩与阴影区域。Sentinel-1卫星的精密轨道星历(precise orbit ephemerides,POD)数据被用来辅助Sentinel-1数据的预处理和基线误差改正[9]。

同时在数据处理中对于与高程不相干的随机大气扰动误差，本文利用李振洪教授团队研制的通用型卫星雷达在线大气改正系统(generic atmospheric correction online service for InSAR, GACOS)来去除[16]。

2 时序InSAR技术数据处理流程及结果分析

φint=φfla+φtop+φdef+φnois+φatm+φorb，

(1)

式中：φfla为平地相位；φtop为DEM误差引入的残余地形相位；φdef为地表形变相位；φnois为噪声相位；φatm为大气延迟相位；φorb为轨道误差引起的相位。因此最终想从φint中获取到形变相位φdef，需要经过相位解缠及一系列的误差去除处理。

基于国内外研究[17-21]并结合本文的研究内容，本文采用SBAS技术与相位堆叠技术(Stacking)进行滑坡隐患识别对比分析研究。SBAS技术是由Berardino等[22]于2002年提出的，是针对覆盖同一区域的多幅SLC(single look complex)数据，按照规范设置一定的垂直基线阈值和时间基线阈值进行组合，通过这种形式可以有效降低因垂直基线和时间基线过大引起的失相干现象，同时利用空间滤波提高相干性，通过牺牲分辨率提高信噪比的方式进行多视处理，进一步减弱失相干对地表监测的影响； 利用差分干涉计算获得差分干涉图和相干图，对差分干涉图进行相位解缠，并对解缠相位进行时间域和空间域的形变量估算，减去趋势误差； 通过时空域滤波来去除大气误差的影响，从而更高精度地获取相对于第一幅SAR影像的累积形变时间序列及平均形变速率。

Stacking技术是由Sandwell等[23]在1998年提出的，该技术原理是对多幅干涉图进行加权平均解算，以达到削弱空间上不相关的噪声的影响，在本文的数据处理中对基线精化之后的解缠图进行相位堆叠解算，通过地理编码获取其在地理坐标系下的平均形变速率图[24]。图2为本文数据处理流程图。通过不同时序InSAR技术获取的地表形变速率图，依据形变阈值，结合DEM、解缠图、差分干涉图及光学影像等信息进行滑坡隐患的早期识别解译。

图2 时序InSAR技术数据处理流程Fig.2 Time-series InSAR technology data processing flow

2.1 基于Stacking InSAR技术地表形变监测结果分析

在雅安市数据处理中，在距离向和方位向采用10×2的多视比，垂直基线不大于150 m，时间基线不大于48 d进行干涉对组合，然后对所有干涉对开展差分干涉、自适应滤波、相位解缠处理及大气改正。解缠方法采用最小费用流法，解缠参考点选择在相干性高、且稳定的区域，利用GACOS在线大气改正系统对解缠图进行残余大气延迟误差改正[25]，最后对解缠图进行筛选，筛选标准依据相干性较好、解缠误差较少，对筛选之后的解缠图进行Stacking解算，获取了覆盖研究区域沿雷达视线向的形变速率图，如图3所示。从图3中可以看出升轨数据结果在石棉县、荥经县出现较多空白区域，在整个实验区的左半部分出现较多的解缠误差，主要是为了让植被更加密集的区域(荥经县、汉源县、石棉县等)可以获取更多的观测信息，通过对覆盖该区域的平均相干图的相干性进行查看，选择大于相干性阈值面积占比80%以上的值作为相干性阈值，对低于相干性的区域进行掩模处理。然而在这些区域相干性过低，使得存在较多的解缠误差； 因为研究区域不能被单景Sentinel-1数据所覆盖，这也就使得两景之间影像边界处存在颜色的差异性，后期在图幅拼接的时候未做校正处理，但是对滑坡隐患的识别、解译并不影响。将石棉县区域结果放大，可以很明显发现因升降轨数据卫星飞行方向不同，石棉县等地区域的结果出现明显的差异性，正是由于升降轨数据收集处理分析，避免了单个轨道数据对隐患点的漏判。

2.2 基于SBAS InSAR技术地表形变监测结果分析

在SBAS数据处理过程中，参数设置与Stacking计算过程设置基本一致。通过时间域和空间域的解算获得了沿卫星雷达视线向方向的地表形变速率结果(图4)，其中负值表示远离卫星飞行方向，正值表示接近卫星飞行方向。

通过解译图4发现升轨数据SBAS结果中存在较大的空白区域，主要集中在宝兴县、天全县和荥经县，结合该区域植被覆盖情况，分析认为这与地表植被覆盖情况相关，造成这些区域在时间域上相干性不连续造成空缺。同时一些区域存在部分解缠误差，分析原因是在相位解缠过程中为了顾及更多的信息，基于相干性设置的掩模阈值过低所导致。

在对高速公路的桥梁进行养护的过程中，其实不仅仅只针对于桥梁和高速公路本身，还要兼顾高速公路周围的生态环境以及生活服务，因为它们是一个整体，不可割裂开来，只有做到全面养护才能保证高速公路桥梁的持久利用。

通过对比Stacking结果与SBAS结果发现，在SBAS结果中山脊处存在部分空白区域，以宝兴县降轨结果表现最为明显，主要原因是在SBAS处理过程中剔除了叠掩、阴影区域。同时升轨数据SBAS结果较升轨数据Stacking结果有较大范围空白，主要是由于空白区域在时间域上相干性不连续所造成的。SBAS结果明显好于Stacking数据结果，分析原因主要是Stacking技术只是简单的进行相位加权平均，对一些残余大气误差、DEM误差等未进行去除，而SBAS技术通过时间域和空间域的解算，能有效地消除或削弱解缠粗差、大气误差以及DEM误差等因素的影响，因此后者结果优于前者。

通过以上对比分析，从地表形变监测结果和处理效果来分析，在雅安市SBAS技术获取结果明显优于Stacking技术获取结果。

3 时序InSAR结果对比分析

3.1 汉源县区域地表形变监测结果分析

选取雅安市汉源县区域单独展示对比分析。使用升轨数据获取了其年平均形变速率图(图5)。发现Stacking结果与SBAS结果具有较大的差异性，从Stacking结果中并不能看到明显形变信息，结合SBAS解译的结果在Stacking结果中识别出2处形变区域，并且主要集中在城镇区域，且形变信息不明显，受大气影响严重。在山区出现许多嘈杂点和颜色拉伸现象，如图5(a)中黄色方框所示，分析原因一方面是这些区域相干性差，解缠阈值过低，存在解缠误差； 另一方面是叠掩现象的存在，经地理编码之后，叠掩区域出现拉伸现象。图5(b)为利用SBAS技术获取的结果，可以看到整个结果干净，但是前期处理中设置的相干性掩模阈值过低，因为剔除了一些误差，所以结果中没有出现左图对应位置的嘈杂点，但是在山区还是表现出一些效果差的区域； 同时对叠掩阴影区域的剔除，使得SBAS结果中在一些地形起伏区域出现空白，但是结果更加精确。基于SBAS结果初步可以解译出5处明显的疑似隐患形变点(红色圆圈所示)，中间圆圈所圈的3处隐患点，一处形变信息与另外两处相反，主要是该结果为雷达视线向形变信息。

选择图5(b)左边第一个隐患点，即红色箭头所示，利用SBAS InSAR技术进行时序分析，获取了该隐患点在时间域上的累积形变曲线图，如图6所示。

图6 疑似隐患点累积形变Fig.6 Cumulative deformation of suspectedhidden danger points

由图6可知，该点从2018年10月—2020年10月一直处于不断变形状态，并且目前处于加速变形期，累积形变量超过了0.07 m，佐证了此点为隐患点。结合图3—6对比分析发现在雅安区域SBAS技术更加适合于滑坡隐患的早期识别及监测工作。

3.2 典型滑坡隐患识别对比分析

选取一经野外验证的滑坡隐患点进行不同时序InSAR技术识别对比分析，该点位于宝兴县陇东镇先锋村。初步通过利用升轨数据基于Stacking技术及SBAS技术处理获取了形变速率图，结合光学影像进行解译，获取了该滑坡的边界信息，如图7所示。从SBAS结果中可以发现明显形变信息，形变区与滑坡范围基本一致，最大形变速率超过了-80 mm/a，最大形变区域位于滑坡体左侧及前部区域，同时将解译的边界信息叠加至Stacking结果上，发现在Stacking结果中变形信息并不明显，与SBAS结果变形速率最大区域所对应位置一致。对比可知该隐患点的SBAS技术识别效果明显优于Stacking技术的识别效果，SBAS结果中形变信息空缺区域是因为在时间域上相干性不连续所造成的。因为该隐患点在降轨数据中处于叠掩区域，故没有展示降轨结果。

图7 典型滑坡隐患解译信息

从光学影像上可知解译区植被分布茂密，可见多处明显小规模溜滑迹象，土地裸露。滑坡呈舌状，两侧以冲沟为界，前缘临沟，外凸较明显，无新近堆积，后缘可见早期滑坡壁，坡体上分布有较多农户。

通过野外调查发现该点为已知点滑坡隐患： 唐包滑坡，地理位置为E102°42′53.1″，N30°26′52.7″。滑坡变形区域主要集中在滑坡前缘右侧和中部，变形强烈区域与InSAR识别结果吻合性较好(图8)。图8(a)为滑坡上半部分照片轮廓，从图中可以看到植被发育茂密，植被以灌木为主，乔木次之，可以看到多处滑塌区域，土地裸露。P1位于滑坡前缘右侧，附近居民点院坝开裂长约10 m，宽3～4 cm，房屋也见多条裂缝存在，宽约5～6 cm不等； 公路错开25 cm，裂缝宽约30 cm，长约5 m，局部下陷25 cm。P2位于滑坡前缘中部，地面可见网状拉裂，宽3～8 cm，公路外侧见明显裂缝，路面破坏长度约30 m。

图8 典型滑坡野外调查照片

3.3 雅安市滑坡隐患识别对比分析

基于不同时序InSAR技术获取的结果进行滑坡隐患识别解译，并通过野外进行验证结果的准确性，统计了雅安市不同InSAR技术识别滑坡隐患分布情况(图9)。如图9所示，图中蓝色“共同识别”为基于Stacking技术和SBAS技术均识别出来，并经过野外调查验证为正确的隐患点； “SBAS识别正确”为单独由SBAS技术识别出来并经野外调查验证为正确的隐患点，“SBAS识别错误”为仅由SBAS技术识别并经野外调查验证为错误的隐患点； “Stacking识别正确”为单独由Stacking技术识别出来并经野外调查验证为正确的隐患点，“Stacking识别错误”为仅由Stacking技术识别并经野外调查验证为错误的隐患点。

图9 雅安市滑坡隐患识别分布图Fig.9 Distribution map of landslide hiddendanger identification in Ya’an City

对识别的隐患点进行统计分析(图10)。 可知基于Stacking技术共识别出隐患点83处，经野外验证、核查，共计有48处隐患点被识别出来。Stacking技术识别准确率为57.8%。对识别错误的隐患点进行分析发现，造成隐患识别错误的原因一方面是植被茂密，相干性较差，解缠过程中为了获取更多的信息设置较低的相干性掩模阈值，使得在结果中有较多的解缠误差存在； 另一方面是Stacking技术只是单纯对解缠图进行加权平均，结果中存在较多的残余误差。

图10 雅安市InSAR技术识别地质灾害统计图Fig.10 InSAR technology in Ya’an City toidentify geological hazards statistical map

基于SBAS技术共识别出隐患点54处，经野外验证、核查，共计有38处隐患点被正确识别出来。SBAS技术识别准确率为70.3%。对识别错误的隐患点进行分析发现，造成隐患识别错误的原因一方面植被茂密，相干性较差，造成结果中有一些异常值； 另一方面由于干涉图在时间域相干性不连续，对一些隐患点不能提取时序分析判识，造成误判，以及由其他因素引发的变形，经野外核查不能归为滑坡隐患。

基于Stacking技术和SBAS技术共同识别的隐患点为44处，经野外验证、核查，共计有38处隐患点被两种技术均被识别出来，识别准确率为86.4%。

经过上述讨论可知，在植被密集区域，地形条件复杂，基于InSAR技术进行滑坡隐患识别，影响和限制因素较多，在顾及识别隐患点数目、准确率下，建议SBAS技术与Stacking技术相结合的形式用于该区域地质灾害隐患识别中。

4 结论

本文选择在植被密集区开展利用2018年10月—2020年10月升降轨Sentinel-1数据进行基于Stacking技术和SBAS技术滑坡隐患早期识别对比分析实验，以雅安市为研究区。并选取典型滑坡体结合野外验证资料进行时间序列分析。主要有以下结论：

1)利用Sentinel-1数据基于时序InSAR技术用于雅安市滑坡隐患早期识别具有较强的可行性，通过获取的结果经对比分析发现升降轨数据共同监测可以有效避免滑坡隐患的漏判，增加观测信息。

2)通过统计发现基于Stacking技术识别的滑坡隐患数目最多，结合野外验证资料发现基于SBAS技术识别的滑坡隐患准确率最高。

3)通过对时序InSAR技术识别的滑坡统计分析发现，在顾及识别隐患点数目、准确率等情况下，建议在植被密集区域采用SBAS技术与Stacking技术相结合的形式开展滑坡隐患识别。