基于干涉点目标分析方法的大光包滑坡形变监测

于 冰 王 杨 马德英 蒋荣乾 张 过 周志伟

1 西南石油大学土木工程与测绘学院，成都市新都大道8号，6105002 西南石油大学油气藏地质及开发工程国家重点实验室，成都市新都大道8号，6105003 西南石油大学油气空间信息工程研究所，成都市新都大道8号，6105004 武汉大学测绘遥感信息工程国家重点实验室，武汉市珞喻路129号，4300795 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室，武汉市徐东大街340号，430077 6 中铁二院工程集团有限责任公司，成都市通锦路3号，610031

干涉点目标分析(interferometric point target analysis, IPTA)方法[1-2]是对PS-InSAR技术的发展，其通过建立相位模型，计算并分离PS点的各相位分量，再根据结果不断精化模型，迭代求解研究区的形变信息。IPTA方法在地表沉降、基础设施形变等方面应用较为广泛[3-5]，但在滑坡形变监测方面应用较少。

本文利用C波段Sentinel-1A影像数据，采用IPTA方法对大光包滑坡进行形变监测，提取该区域2016～2018年的形变信息，并结合全球降水测量(global precipitation measurement, GPM)计划日均降水量数据，对特征点进行时序形变分析，将该方法获得的结果与StaMPS-SBAS方法的处理结果进行对比，分析IPTA方法应用于大型滑坡形变监测的可行性与潜力，为大光包滑坡的监测与预警提供技术和信息支撑。

1 IPTA方法基本原理

IPTA方法的核心思想是对相干点目标进行处理，通过构建和更新二维线性相位回归模型，迭代求解相邻相干点的线性形变速率增量和高程改正值，直至逼近线性形变速率和高程误差的最优解[1-2]。在IPTA处理中，根据影像幅度和相位信息，计算每幅影像的光谱相干系数和时序振幅稳定性，综合选取高质量的相干点目标[6]。光谱相干系数可表征光谱多样性的大小，聚焦良好的相干点目标具有与分布式目标不同的光谱特性，相似的相干点目标具有相似的光谱特性，表现为低光谱多样性，光谱相干系数较高。振幅稳定性可表征观测目标对雷达波后向散射的稳定程度，相干点目标具有较高的散射稳定性。本研究中时序振幅稳定性阈值设为1.4，光谱相干系数阈值设为0.5，分别获得2个候选相干点目标点集，然后对二者进行合并，剔除同名点，共获取32 758个初始相干点。

相干点的差分干涉相位主要包含形变、残余地形、轨道误差、大气延迟和噪声等相位。IPTA以二维线性相位回归模型为基础，即在空间域对相干点与局部参考点的差分干涉相位进行二次差分，以削弱大气延迟和非线性形变的影响；在时间域对差分干涉相位进行回归处理来分析相位变化过程[7]。相干点目标的二次差分相位可表示为：

(1)

式中，φdiff,i-j为相干点目标间的差分相位差，Δv、Δh分别为建立连接关系的点目标的线性形变速率增量和高程改正值增量，BT、B⊥分别为差分干涉图的时间基线和空间垂直基线，λ为雷达中心波长，R为雷达至地面目标的斜距，θ为雷达入射角，φres,i-j为残余相位，由轨道误差相位φorb_res,i-j、大气延迟相位φatm,i-j和噪声相位φnoi,i-j组成。

在IPTA中引入残余相位标准差来判断回归分析的最优解，并采用逐步回归、迭代改进的方式不断计算修正值，更新线性形变相位、高程相位和残余相位，计算所有相干点目标间的线性形变速率增量和高程误差改正值增量，选定一个相对稳定点作为参考点，从而获取每个相干点相对于参考点的线性形变速率和高程改正值。根据IPTA官方说明书，残余相位标准差小于1.2的相干点可满足精度要求。由于本文研究区地形起伏较大，且无先验信息，因此将第一次回归分析的高程改正值阈值设为±2 000 m，形变速率阈值设为±150 mm/a。较大的搜索范围会解算出较大的形变速率和地形误差，后续逐步改进阈值，缩小范围。在最终的回归分析中，高程改正值阈值设为±500 m，形变速率阈值设为±80 mm/a。最终输出的相干点目标相位标准差均低于1.2，符合实验精度要求。

从差分干涉相位中扣除线性形变和高程误差相位可得到残余相位。其中，残余轨道误差在InSAR中表现为基线误差，可通过基线精化改正[8]；噪声在空间维为高频信号，在空间域对残余相位进行低通滤波可去除；大气延迟为时间维的高频信息，非线性形变为时间维的低通成分，在时间域进行低通滤波可得到非线性形变。最后将线性形变与非线性形变进行组合，计算得到每个相干点的形变时间序列。IPTA数据处理流程如图1所示。

图1 数据处理流程Fig.1 Flow chart of data processing

2 研究区概况与实验数据

2.1 研究区概况

大光包滑坡为2008年汶川地震触发的最大规模滑坡，汶川地震后，该地区海拔由3 047 m降至2 964 m，其长达1.8 km的剪滑光面和高达690 m的滑坡堰塞坝引起国内外学者的广泛关注[9-11]。大光包滑坡位于四川省绵阳市安州区高川乡，中心坐标为104°07′0.2″E、31°38′19.8″N，处于汶川地震发震断层上盘，距地震中心85 km，与地震地表破裂带相距4.5 km[9]，其地理位置如图2所示。

图2 大光包滑坡地理位置Fig.2 Geographical location of Daguangbao landslide

大光包滑坡地质结构复杂，根据滑坡要素、滑坡堆积体形态、长度和延伸方向、岩性分布，可将滑坡体分为12个区域(图3)。Ⅰ为侧向前冲式碎屑流区，表层物质易随白果林沟冲刷流动形成堰塞湖；Ⅱ为飞跃前冲式碎屑流区，表层物质向川林沟流动形成堰塞湖；Ⅲ为左侧滑坡扰动体，岩体破碎严重，呈碎块状；Ⅳ为主滑坡体堆积区，由大光包顶抛射降落于该处，主要为泥盆系沙窝组白云岩；Ⅴ为沟道碎屑流区，多为黑色灰岩碎块，沿黄洞子沟向下流动堆积；Ⅵ为滑坡后缘二次滑塌体，在滑动时与主滑坡体分开后停留于此；Ⅶ为坡面碎屑流区，岩性主要为灰黑色灰岩；Ⅷ为雨水冲刷滑坡表层物质形成的沟口冲积扇；Ⅸ为白云岩和板岩构成的滑坡左侧断壁；Ⅹ为滑坡后缘断壁倒石锥，岩体松散破碎；Ⅺ为滑坡右侧主滑动面，滑床为震旦系灯影组；Ⅻ为滑坡后缘断壁，出露大光包粉砂岩至白云岩的地层剖面[10-11]。

图3 大光包滑坡分区Fig.3 Zoning of Daguangbao landslide

研究区地势西高东低，地形起伏大，地貌以高山为主，具有逆冲-推覆-滑脱-走滑的构造特点[9]。由于研究区常年降雨，地质环境脆弱，该区在接下来数十年内将会处于不稳定状态，再次发生山体滑坡的风险高于地震前[12-13]，对其进行形变监测具有重要意义。

2.2 实验数据

本文采用覆盖实验区的宽幅模式(IW)升轨Sentinel-1A影像，幅宽为250 km，空间分辨率为5 m×20 m，极化方式为垂直同向(VV)极化，入射角为33.83°。研究时间跨度为2016-01～2018-12，影像获取时间间隔为12 d，共计65幅。干涉组合采用单一主影像模式，空间基线阈值设定为100 m，图4为65幅影像的成像时间与空间基线分布。本文选用的DEM数据为2000年美国宇航局主导测量的SRTM-3 DEM数据，分辨率为90 m，用于地理编码和差分处理以去除地形相位。

图4 时空基线分布Fig.4 Distribution of spatiotemporal baselines

3 监测结果与分析

图5为采用IPTA方法得到的大光包滑坡地区雷达视线向年形变速率。从图中可以看出，监测期间大光包滑坡处于形变状态，滑坡体左部较为稳定，形变速率为-15～2 mm/a。飞跃前冲式碎屑流区(Ⅱ区)和主滑坡体堆积区(Ⅳ区)呈现明显漏斗状形变，表面物质松散不稳定，受川林沟、黄洞子沟和雨水冲刷影响，形变量大于其他区域，雷达视线向最大形变速率达-50.590 mm/a，最大累积形变量达-170.726 mm。

图5 雷达视线向形变速率Fig.5 Line-of-sight deformation rate

以72 d为时间间隔输出2016-01-13～2018-12-28滑坡的形变时间序列，结果如图6所示。由图可知，监测时段内大光包滑坡右侧形变漏斗逐渐发育，监测时间到684 d时，形变漏斗已初具规模，雷达视线向累积形变量为-111.543 mm。随着时间推移，形变漏斗的范围逐渐扩大，到2018年年底，形变漏斗长度达767 m，宽度达687 m，雷达视线向累积形变量为-170.726 mm。

图6 大光包滑坡形变时间序列Fig.6 Deformation time series of Daguangbao landslide

在大光包滑坡形变漏斗中心和边缘区域选择6个特征点进行时间序列形变分析，特征点位置如图7所示，形变结果如图8所示。可以看出，6个特征点均呈线性滑动趋势，但形变速率不同。A、B、C点位于形变量较大的漏斗中心，所在位置的坡度分别为49.73°、53.31°和51.49°，整体形变量和形变速率逐年增大，形变速率约为-48.387～-47.361 mm/a，累积形变量约为-141～-138 mm；D、E、F点位于形变量较小的漏斗边缘区域，所在位置的坡度分别为33.67°、25.80°和31.06°，形变发育较为稳定，形变速率约为-17.788～-9.219 mm/a，累积形变量约为-57.432～-31.117 mm。

图7 特征点位置Fig.7 Location of feature points

本文引入大光包地区2016-01-01～2018-12-28的日均降水量数据，结合特征点时序形变进行比较，并统计雷达重访周期内暴雨发生前后特征点的形变量，结果如表1所示。从图8可以看出，大光包地区的降雨量具有明显的季节性特征，每年5～10月雨水充足，8月常有强降雨发生，滑坡上松散物质受雨水冲刷影响，碎石块在滑坡表面滑动[14]，引起累积形变量增大。2016-08-22降雨量为118.815 mm，2016-08-16～09-09特征点的最大形变增量为-10.829 mm；2017-08-22降雨量为80.400 mm，2017-08-11～23特征点的最大形变增量为-5.877 mm；2018-06-25降雨量为71.849 mm，2018-06-19～07-01特征点的最大形变增量为4.858 mm；2018-07-10降雨量为 60.782 mm，2018-07-01～25特征点的最大形变增量为2.162 mm。2018年强降雨发生后，滑坡上松散物质受雨水影响，碎石块在滑坡表面形成新堆积体并伴有滑坡体吸水膨胀抬升[14]，引起特征点形变量减小。综上表明，降水是影响滑坡稳定性的主要因素之一。

表1 特征点形变量统计

为进一步分析形变漏斗特征，以D点为起点，途经A点，F点为终点，作剖面线进行研究。从上到下提取该剖面线共74个点的累积形变量，绘制以144 d为间隔的剖面线形变量时序图(图9)。

图9 剖面线形变量时序Fig.9 Time series deformation of the profile

研究区地形起伏较大，最大高差达245 m，在监测期间内该坡面线累积形变量逐年增大，滑坡顶部形变量较小，形变集中在滑坡中部碎屑流区。该区域坡向为东南向，坡度为45°～50°，松散堆积体易受到降雨及黄洞子沟冲刷影响而产生较大形变，监测时段内最大形变量达-167 mm。

采用StaMPS-SBAS方法对Sentinel-1A影像数据进行处理，时间基线设为50 d，空间基线设为100m，生成212个干涉对，剔除相干性较差的

干涉对后余下183个干涉对，共提取57 225个相干点，所提取的形变速率如图10所示。对比分析IPTA方法和StaMPS-SBAS方法所得结果可知，两者趋势一致，飞跃前冲式碎屑流区(Ⅱ区)和主滑坡体堆积区(Ⅳ区)均呈现出形变漏斗，且形变数值一致，表明两者结果可靠。

图10 StaMPS-SBAS方法得到的雷达视线向形变速率Fig.10 Line-of-sight deformation rate obtained by StaMPS-SBAS method

4 结 语

本文利用IPTA方法对升轨Sentinel-1A影像数据进行处理，获取大光包滑坡2016～2018年的形变时间序列，并结合GPM日均降水量数据进行特征点累积形变分析，得到以下结论：

1) 大光包地区在2016～2018年发生明显形变，形变主要集中在飞跃前冲式碎屑流区和主滑坡体堆积区，呈线性滑动趋势，雷达视线向形变速率在-50.590～19.175 mm/a之间，最大累积形变量达-170.726 mm。

2) 结合GPM日均降水量分析发现，大光包滑坡形变与降雨量变化具有一定相关性。研究区每年5～10月雨水充足，该时间段应结合降水数据，加强滑坡监测，作好相应防范措施，防止滑坡泥石流等灾害发生。

3)IPTA方法和StaMPS-SBAS方法均探测出飞跃前冲式碎屑流区和主滑坡体堆积区出现形变漏斗，形变量级一致，且形变整体趋势相同，证明了2种方法的有效性。

4) 作为PS-InSAR方法的发展，IPTA方法不仅可以减少时空失相关和大气延迟的影响，还可以通过矢量化数据节省计算空间，提高数据处理速度和效率，迭代求解得到的结果也更加可靠，在大型滑坡形变监测中具有较好的可靠性和广阔的应用前景。