2021, Vol. 41
2. 西南石油大学油气藏地质及开发工程国家重点实验室，成都市新都大道8号，610500;
3. 西南石油大学油气空间信息工程研究所，成都市新都大道8号，610500;
4. 武汉大学测绘遥感信息工程国家重点实验室, 武汉市珞喻路129 号, 430079;
5. 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室, 武汉市徐东大街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]。相干点目标的二次差分相位可表示为：
| $ \Delta {\varphi _{{\rm{diff}},i - j}} = \frac{{4{\rm{ \mathsf{ π} }}}}{\lambda }\Delta \upsilon {B_T} + \frac{{4{\rm{ \mathsf{ π} }}{B_ \bot }}}{{\lambda R\sin \theta }}\Delta h + {\varphi _{{\rm{res}},i - j}} $ | (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 |
大光包滑坡为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 |
图 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 |
|
图 8 特征点累积形变量与降雨量关系 Fig. 8 Relationship between cumulative deformation of feature points and precipitation |
本文引入大光包地区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 特征点形变量统计 Tab. 1 Statistics of feature points deformation |
为进一步分析形变漏斗特征，以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，空间基线设为100 m，生成212个干涉对，剔除相干性较差的干涉对后余下183个干涉对，共提取57 225个相干点，所提取的形变速率如图 10所示。对比分析IPTA方法和StaMPS-SBAS方法所得结果可知，两者趋势一致，飞跃前冲式碎屑流区(Ⅱ区)和主滑坡体堆积区(Ⅳ区)均呈现出形变漏斗，且形变数值一致，表明两者结果可靠。
|
图 10 StaMPS-SBAS方法得到的雷达视线向形变速率 Fig. 10 Line-of-sight deformation rate obtained by StaMPS-SBAS method |
本文利用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方法不仅可以减少时空失相关和大气延迟的影响，还可以通过矢量化数据节省计算空间，提高数据处理速度和效率，迭代求解得到的结果也更加可靠，在大型滑坡形变监测中具有较好的可靠性和广阔的应用前景。
| [1] |
Werner C, Wegmuller U, Strozzi T, et al. Interferometric Point Target Analysis for Deformation Mapping[C]. IEEE International Geoscience and Remote Sensing Symposium, Toulouse, 2003
(  0) |
| [2] |
Werner C, Wegmuller U, Wiesmann A, et al. Interferometric Point Target Analysis with JERS-1 L-Band SAR Data[C]. IEEE International Geoscience and Remote Sensing Symposium, Toulouse, 2003
( 0) |
| [3] |
Liu Y Y, Zhao C Y, Zhang Q, et al. Land Subsidence in Taiyuan, China, Monitored by InSAR Technique with Multisensor SAR Datasets from 1992 to 2015[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2018, 11(5): 1509-1519 DOI:10.1109/JSTARS.2018.2802702
( 0) |
| [4] |
李广宇, 张瑞, 刘国祥, 等. Sentinel-1A TS-DInSAR京津冀地区沉降监测与分析[J]. 遥感学报, 2018, 22(4): 633-646 (Li Guangyu, Zhang Rui, Liu Guoxiang, et al. Land Subsidence Detection and Analysis over Beijing-Tianjin-Hebei Area Based on Sentinel-1A TS-DInSAR[J]. Journal of Remote Sensing, 2018, 22(4): 633-646)
( 0) |
| [5] |
刘雨鑫. 利用时序InSAR技术提取格尔木市及周缘地区的形变[D]. 武汉: 武汉大学, 2019 (Liu Yuxin. Surface Deformation Monitoring of Golmud and Surrounding Areas Using the Temporal InSAR Technique[D]. Wuhan: Wuhan University, 2019)
( 0) |
| [6] |
Ferretti A, Prati C, Rocca F. Permanent Scatterers in SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8-20 DOI:10.1109/36.898661
( 0) |
| [7] |
俞晓莹, 姜成岭, 张建, 等. IPTA监测圣佩德罗湾港口地表时序沉降[J]. 测绘科学, 2012, 37(6): 21-25 (Yu Xiaoying, Jiang Chengling, Zhang Jian, et al. IPTA Monitoring Long-Term Series Surface Deformation of San Pedro[J]. Science of Surveying and Mapping, 2012, 37(6): 21-25)
( 0) |
| [8] |
张海波, 李宗春, 许兵, 等. IPTA方法在地面沉降监测中的应用[J]. 测绘科学技术学报, 2016, 33(2): 145-149 (Zhang Haibo, Li Zongchun, Xu Bing, et al. Ground Subsidence Monitoring Using Interferometric Point Target Analysis[J]. Journal of Geomatics Science and Technology, 2016, 33(2): 145-149)
( 0) |
| [9] |
崔圣华, 裴向军, 黄润秋, 等. 大光包滑坡不连续地质特征及其工程地质意义[J]. 西南交通大学学报, 2019, 54(1): 61-72 (Cui Shenghua, Pei Xiangjun, Huang Runqiu, et al. Discontinuities and Engineering Geological Significances of Strong Earthquake-Induced Daguangbao Landslide[J]. Journal of Southwest Jiaotong University, 2019, 54(1): 61-72)
( 0) |
| [10] |
李天涛, 裴向军, 黄润秋. 强震触发大光包巨型滑坡运动特征研究[J]. 水文地质工程地质, 2014, 41(2): 116-121 (Li Tiantao, Pei Xiangjun, Huang Runqiu. A Study of Motion Features of the Daguangbao Large-Scale Landslide Induced by the Wenchuan Earthquake[J]. Hydrogeologyand Engineering Geology, 2014, 41(2): 116-121)
( 0) |
| [11] |
许向宁, 李胜伟, 王小群, 等. 安县大光包滑坡形成机制与运动学特征讨论[J]. 工程地质学报, 2013, 21(2): 269-281 (Xu Xiangning, Li Shengwei, Wang Xiaoqun, et al. Characteristics of Formation Mechanism and Kinematics of Daguangbao Landslide Caused by Wenchuan Earthquake, Sichuan, China[J]. Journal of Engineering Geology, 2013, 21(2): 269-281 DOI:10.3969/j.issn.1004-9665.2013.02.013)
( 0) |
| [12] |
Huang R Q, Fan X M. The Landslide Story[J]. Nature Geoscience, 2013, 6(5): 325-326 DOI:10.1038/ngeo1806
( 0) |
| [13] |
Dai K R, Li Z H, Tomás R, et al. Monitoring Activity at the Daguangbao Mega-Landslide(ChinaUsing Sentinel-1 TOPS Time Series Interferometry[J]. Remote Sensing of Environment, 2016, 186: 501-513 DOI:10.1016/j.rse.2016.09.009
( 0) |
| [14] |
陈磊, 赵学胜, 汤益先, 等. 基于时序InSAR技术的大光包滑坡变形监测[J]. 矿业科学学报, 2016, 1(2): 113-119 (Chen Lei, Zhao Xuesheng, Tang Yixian, et al. Deformation Monitoring of the Daguangbao Landslide with Time Series InSAR Technique[J]. Journal of Mining Science and Technology, 2016, 1(2): 113-119)
( 0) |
2. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, 8 Xindu Road, Chengdu 610500, China;
3. Institute of Petroleum and Natural Gas SpatialInformation Engineering, Southwest Petroleum University, 8 Xindu Road, Chengdu 610500, China;
4. State Key Laboratory ofInformation Engineeringin Surveying, Mapping and Remote Sensing, Wuhan University, 129 Luoyu Road, Wuhan 430079, China;
5. State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academyfor Precision Measurement Science and Technology, CAS, 340 Xudong Street, Wuhan 430077, China;
6. China Railway Eryuan Engineering Group Co Ltd, 3 Tongjin Road, Chengdu 610031, China