基于Sentinel-1A数据的地表形变监测

刘智学，薛东剑*

（1.成都理工大学 地球科学学院，四川 成都 610059）

我国西南地区受印度板块向亚欧板块俯冲影响，是地震活动发生极频繁且强烈的区域，地震带来的次生地质灾害不可计数。在所有地质灾害中，滑坡是最常见的类型，占地质灾害总量的70%以上[1]。由2008年汶川地震（里氏震级8.0 Ms、矩震级8.3 Mw）引发的大光包滑坡是世界上规模最大的滑坡点之一，其内潜藏的大量安全隐患足以对当地造成数十年的影响。因此，为避免震后次生灾害威胁，需要对高风险区域进行长期监测。

早期的滑坡形变探测方法包括大地测量法、GPS法、自动伸缩计法和分布式光纤法等[2]；但这类以点探面的监测方法很难全面反映滑坡区域的形变特征，易导致分析结果与实际情况间存在较大误差。合成孔径雷达（SAR）是一种有源成像系统，基于其发展的DInSAR地表形变监测技术可大面积连续监测地表 cm级甚至mm级的形变[3]，在地震形变场提取、火山活动以及矿区沉降监测等方面卓有成效，如侯丽燕[4]等基于多次DInSAR处理研究了地震序列发震断层特征，加深了对周边活动断裂性质的认识；Wempen J M[5]利用DInSAR技术定量分析了矿区开采初期和随着采矿方向推进所造成的地面沉降，为沉降动力的深入研究提供了理论依据。然而，传统的DInSAR技术更注重短期形变监测，在解决失相干现象时略有不足。因此，长时间序列差分干涉测量方法应运而生，主要包括适用于城市、乡镇等相干目标密集区域的永久散射体干涉测量（PS-InSAR）和同时适用于城市区域和自然表面的小基线干涉测量（SBAS-InSAR）方法，后者采用了更多的主影像来构造干涉像对，有效避免了长时间基线带来的失相干，在地表沉降与滑坡形变监测中应用广泛。程海琴[6]利用时序InSAR技术分析了龙门山地区震后的滑坡灾害点时空分布特征，定量揭示了滑坡物质迁移特征和震后滑坡体的地表形变特征；陆超然[7]等分别利用差分干涉和SBAS方法进行了大范围地灾隐患点的探测，并总结了二者的优劣；莫玉娟[8]等基于SBAS技术提取了小金县的地表沉降 趋势。

本文分别利用DInSAR和SBAS-InSAR技术对包含研究区域的22景哨兵一号（Sentinel-1）C波段数据进行处理，提取分析了该滑坡点的视线向累积形变量、时序形变数据和年平均形变速率，并对比了两种数据提取方法的形变结果，为大光包滑坡稳定性探测工作提供了数据参考。

1 研究区概况

大光包滑坡（31e38′N、104e07′E）位于四川省绵阳市安州区，四川盆地西北与山地区域交界处，西北方向有大面积山林覆盖，地势落差大，平均海拔相对较高；地质结构复杂多样且降水量较大，地貌构造包括平坝、丘陵以及低中型山脉3种类型。该区域地处龙门山断裂带，地质灾害较多，尤以地震为最，中小型地震频繁发生（图1）。在本文选取的雷达影像数据的成像时间内，四川地区共发生36次里氏4.0级以上的浅源地震，其中包括11次不低于里氏5.0级的中强地震。滑坡分别距离2019年12月9日安州区发生的里氏4.6级地震和2020年2月3日成都市青白江区发生的里氏5.1级地震只有18 km和85 km。

图1 研究区位置示意图 （审图号：GS(2019)3333）

区域内亚热带湿润季风气候显著，全年降雨量足，冬春与夏秋季变换时，旱、雨季交替出现，中小型河流交织密集，水资源充沛。为分析研究区旱、雨季降雨量变化与地表形变速率变化之间的线性关系，本文制作了滑坡周边100 km范围内的月均降雨量数据柱状图，如图2所示。

图2 研究区降雨量数据

2 数据来源与研究方法

2.1 数据来源

Sentinel-1卫星所提供的C频段雷达遥感影像能提取大范围、高精度的地面运动数据，是欧洲委员会和欧空局针对欧洲航天哥白尼计划的第一颗对地观测卫星[9]。该卫星为双星座卫星系统，成像重访周期最短可达6 d，且所有数据均免费提供给公众。本文采用Sentinel-1A卫星干涉宽幅模式数据进行研究。

2.2 DInSAR技术

DInSAR技术基于重复获取的同一区域后向散射信号，通过共轭相乘获取干涉相位[10]，获取相干性良好的雷达复数影像差分干涉两次成像之间的相位差信息，再去除干涉相位中包含的地形起伏相位信息，减少大气效应、噪声等的影响，得到最终形变结果。根据DInSAR处理时对不同时相SAR影像的数目要求，常用的 DInSAR 技术可分为两轨法、三轨法和四轨法[11]3种，本文采用两轨法进行处理。该方法将研究区域前后不同时间形成的两景雷达影像作为干涉像对来进行差分干涉处理。由于其获得的干涉相位中既包含形变相位，又包含地形相位信息（忽略大气扰动的影响），因此需以实验数据成像区域的先验DEM数据为外部参考，对相同空间几何条件下的地形相位进行模拟，并将其从合成干涉数据中去除，这也意味着DEM数据的准确性将直接影响DInSAR的结果。其数学表达 式为：

式中，ΔΦD为差分相位，是模拟得到的地形相位。

2.3 SBAS-InSAR技术

Berardino P等为提高差分干涉测量的精度，利用SBAS-InSAR技术来解决失相干等问题。该方法基于不同时间段的单视数据提取相干目标[12]，从而获取大范围时序形变规律。其优势在于，对数据量要求不高，通过链接多个长时间序列影像，构成小基线集合，进而增加数据利用率，最终获得mm级形变速率结果。SBAS-InSAR的数据处理流程如图3所示。

图3 SBAS-InSAR的数据处理流程

在SBAS-InSAR测量方法中，将研究区某一时间段内的N景数据按预设条件进行配对组合，生成最多M组像对，即

为获取合适的干涉像对进行形变结果提取，需设置一个合适的时空基线阈值来对数据进行选择。由于Sentinel传感器平台重访周期较短，为提高形变探测精度、避免完全失相关，本文将最大临界基线阈值设置为45%。像对组合经过干涉工作流处理后，将生成解缠后的干涉相位图，则第i组像对生成的干涉图中任一像素的干涉相位为：

式中，x、r分别为像素点的方位向和距离向坐标；A、B分别为第i组像对的主辅影像。

为更准确地获取真实相位值，解算结果中(-π,π)的相位值，解缠后相位表达式为：

式中，tn、tn-1分别为该组像对两幅影像的成像时间；vn为对应的平均相位速率。

加入控制点文件对解缠后的结果进行轨道精炼，用以校正干涉相位以及解缠后所得的相位信息，如恒定相位的去除。本文输入雷达坐标系下的控制点进行优化处理，将上述结果描述为矩阵形式，即

得到一个M×N矩阵，为避免矩阵秩亏，采用奇异值分解[13]最小二乘处理后，再对各成像间隔内速度积分，即可得到研究区缓慢形变的最终形变时间序列。SBAS时序累积形变如图4所示。

图4 SBAS时序累积形变

3 研究结果与分析

为了方便对DInSAR形变数据进行综合分析，将每组像对的累积形变量制作成图表（图5），从而直观反映形变结果。将每组影像像对间的累积形变值与地震、降水量数据进行对比，可以看出，它们之间存在一定的相关关系，如滑坡区在11月—次年3月步入旱季，降雨量较少，暂不考虑雨水引发的滑坡体滑动，主要查看地震数据，这期间四川省内共发生浅源可感地震11次，其中2020年2月的5.1级地震距离滑坡点仅不到30 km，较之以往冬季的形变量更大，但也没有产生明显滑动，说明中小型浅源地震的发生对滑坡体影响较小，推测地震发生频次更高、且距离滑坡体更近时，滑坡形变数据才会产生更明显的变化；滑坡点周边在2018年3月—10月未发生5级以上的中型地震，可将地震对滑坡滑动的影响忽略不计，主要考虑降水量变化，随着雨季降水量的明显增加，滑坡也开始表现出更大的形变数据波动，受滑坡体附近水流冲刷作用影响，滑坡部分地区有一定程度的 抬升。

图5 DInSAR形变结果的累积形变区间

SBAS-InSAR测量得到的研究区近两年累积形变量（图6）表明，滑坡点及其周边位置虽一直处于动态变化中，但其在固定时间段内的总体抬升与沉降量相互平衡，因此其一维视线方向上的形变值变化并不明显，重大地质灾害发生的可能性较小。

图6 SBAS累积形变量

为了更深入地了解滑坡点附近的地表形变，本文选取滑坡上的A、B、C、D、E五个位置，将其平均沉降速率绘制成曲线图进行分析，如图7所示，可以看出，D、E点有少量抬升，B、C点有沉降现象， A点附近则近于平稳。根据滑坡点周围地质条件推测，坡体高处的风化剥蚀堆积物导致了D、E点的抬升，而B、C点受河流冲刷影响，缓慢下降。

图7 部分点位平均形变速率/(mm/a)

4 结 语

针对大光包滑坡实验区地形地貌复杂，实地勘测地理位置信息难度较大、危险性较高的特点，本文通过非接触传感器获取的遥感影像来研究地震和降水量对滑坡体产生的影响。利用Sentinel-1A卫星影像获取实验区相位数据，下载并处理了2018年3月— 2020年2月的22景干涉宽幅模式的单视复数影像；利用DInSAR和SBAS-InSAR技术计算并提取了滑坡点在该时期内的阶段形变量、时序累积形变量和年平均形变速率；结合旱、雨季降雨量变化以及周边地区浅源地震发生频次对滑坡形变速率变化规律进行了总结。研究结果表明，滑坡体在实验数据成像时间段内的形变量不大，且波动较小；滑坡区域的最大抬升约为10 cm，累积沉降量最大不超过6 cm，年平均形变速率趋于稳定，最高不超过100 mm/a，说明大光包滑坡已处于较稳定的状态；在春、冬交替之际，随着降雨量的减少，河流冲刷作用减弱，形变趋势也稍有减缓；周边地区中小型浅源地震的发生没有对滑坡点造成显著危害。