基于时序InSAR的矿区滑坡前地表运动特征分析

2022-09-21贺黎明裴攀科吴立新张香凝

东北大学学报（自然科学版） 2022年9期

贺黎明, 裴攀科, 吴立新, 张香凝

(1. 东北大学资源与土木工程学院，辽宁沈阳 110819; 2. 中南大学地球科学与信息物理学院，湖南长沙 410083)

近年来,随露天矿产资源的大量开发利用，矿区发生的滑坡灾害越来越多，造成巨大的破坏和经济损失[1].分析滑坡演化过程在风险预防中起关键作用[2]，其中对滑坡前地表运动特征进行分析显得尤为重要，该特征可以帮助预测边坡变形的发展趋势[3]，对确保露天矿山安全生产具有重要意义.

传统的滑坡形变监测方法一般需知滑坡位置等先验知识，主要方法有实地测量法、北斗/GNSS监测法、分布式光纤法等.这些监测方法都是将局部特征点作为数据采集目标，从局部映射到整体进行分析，无法全面且精准反映整个研究区域内滑坡形变的情况，不能应用于人员难以进入且处于高风险的易滑坡区，容易造成“测区未滑，滑区未测”的现象.随遥感技术的发展，大量新的滑坡变形监测手段涌现并得到广泛应用，如雷达差分干涉测量[4]、无人机航空摄影测量[5]、三维激光扫描[6]等，但每一种方法都有其自身的局限性和适用范围[7]，如何精准预测滑坡并降低滑坡灾害所造成的经济财产损失及人员伤亡仍是一个亟待解决的问题.利用时序InSAR获取矿区地表变形场信息[8-9]，对矿区滑坡前地表运动特征进行分析，结合先验知识对滑坡进行早期预测、预警是一种重要且有效的方法[10].

实验以鞍山市鞍千哑巴岭露天采场边坡为研究对象，基于44景Sentinel-1雷达影像，对鞍千哑巴岭露天采场2019年11月25日发生滑坡的西南侧边坡进行了遥感监测.利用多轨SAR数据提取边坡形变信息，对不同轨道SAR数据集反演结果进行可靠性分析，提出评价不同轨道SAR数据集监测结果可靠性的判别依据.在此基础上，发现鞍千哑巴岭露天采场西南侧边坡滑坡前地表运动异常加速的变化特征，并对诱发滑坡的因素进行分析.研究结果表明，多轨时序InSAR技术可为矿区地质灾害预测预警提供可靠科学依据[11]，对于保障矿山安全生产具有重要意义.

1 研究区域

研究区域为隶属于鞍山市鞍千矿业有限责任公司的哑巴岭露天采场，采场位于辽宁省鞍山市东部，距离鞍山市中心的直线距离为12 km，生产规模为1×107t/a.该采场的范围为：北纬41°05′50″～41°06′30″，东经123°07′50″～123°08′40″，矿体长约1 900 m，平均水平厚度约120 m，海拔标高约为100～200 m，矿体倾角一般大于75°，倾向为NE或SW，如图1所示.由于长期持续开采，哑巴岭采场对比2008年投产以来发生了很大变化，其大部分区域已经变成深凹露天开采区，形成很多易滑坡危险区.

图1 鞍千哑巴岭露天采场地理位置图Fig.1 Location map of the Anqian Yabaling open pit

2 数据集与方法

2.1 SAR数据集

研究使用哨兵1A/B(Sentinel-1 A/B)C波段(波长约为5.6 cm)SAR影像来监测获取鞍千哑巴岭采场边坡的时间序列形变图.针对覆盖整个鞍千哑巴岭采场的研究区，在欧空局网站收集了2019年6月7日至2019年11月22日共29景干涉宽幅模式(IW)的升轨(T25，T98)单视复数(single look complex，SLC)影像和2019年6月6日至11月21日共15景干涉宽幅模式(IW)的降轨(T105)单视复数影像，影像时间与空间分布分别如图2和图3所示.

图2 SAR影像时间分布图Fig.2 SAR images time distribution

图3 实验数据覆盖范围Fig.3 Coverage range of the experimental data

2.2 数据处理方法

由于SAR卫星影像数据的限制，我国将InSAR技术应用于滑坡预测领域中的时间较晚，近些年随Sentinel卫星影像数据免费公开使用及国产雷达卫星投入使用，InSAR在该领域的研究进展迅速[12].将干涉技术与合成孔径雷达(SAR)技术结合形成的时序InSAR技术为地质灾害所引起的大规模位移监测提供了新手段，可通过不同时间从SAR传感器收集的影像相位信息来监测细微的地表形变，通过多个不同轨道提供的多视角结果可更全面地评价监测结果的可靠性.

数据处理分为三大步骤，分别为基础资料收集、SBAS-InSAR数据处理、滑坡识别与滑坡前地表运动特征分析，技术路线如图4所示.

图4 技术路线流程图Fig.4 Technical flow chart

SBAS-InSAR技术即小基线集技术，是基于早期小基线差分干涉图集反演地表形变的方法所演变而来的一项技术.SBAS-InSAR技术的基本原理是基于一定数量同一研究区域的SAR影像，通过处理大量短时空基线干涉像对，应用奇异值分解(singular value decomposition, SVD)得到相干点的位移速率和位移时间序列.

SBAS-InSAR算法基于N+1幅单视复数SAR影像，通过限定最大时间和空间基线的阈值组合成大量短基线像对以减弱时空去相关现象.一方面，由于像对空间基线较小从而避免了部分地形因素对残余相位的影响；另一方面，在时间基线域中将干涉图分离成几个独立的子集合，只要干涉图子集合在时间域上重叠，就可以正确分离出对应相位值[13]，从而使用最小二乘法得到每一个子集合的形变时间序列.在M个干涉图中选择一组共同的相干像对，通过相位展开从干涉相位中分离出绝对相位值，从而将所有像素都映射到一个具有高相干性和先验信息的参考基准点.基于所选择的高相干点目标建立观测方程，再将所有的小基线子集合通过奇异值分解法进行解算，从而得到所需要的时间形变序列结果.SBAS-InSAR技术将每一组D-InSAR所得到的形变结果作为SBAS-InSAR处理的观测值，这大大增加了时间采样率，同时通过约束基线的方法以达到时空去相关的目的，有效提高了处理效率.

3 基于时序InSAR的矿区滑坡前形变监测

通过处理实验数据得到了三组不同轨道影像数据集所获取的哑巴岭露天采场视线(line of sight，LOS)的时间序列形变结果，如图5～7所示.其中值为负值(红色)表示远离卫星的地面运动，而正值(蓝色)则表示朝向卫星的地面运动.三组不同轨道数据集所获取的InSAR监测结果表明哑巴岭采场西南侧边坡存在异常变形，变形范围不仅局限于边邦陡坡，下面将依次对每一幅时间序列图进行详细分析.

图5 哑巴岭露天采场升轨1(T25)数据集LOS向形变时间序列图Fig.5 Time series graphs of deformation in LOS direction of the Yabaling open-pit slope ascending 1(T25)dataset

哑巴岭露天采场升轨1(T25)数据集所获取的LOS向形变时间序列图如图5所示，整个监测期间采场西南侧边坡中心区域出现了较大量级和范围的沉降.根据意大利IFFI工程所改进的滑坡活动性分级体系[14]，充分证明该区域为发生滑坡的潜在危险区，可初步根据监测结果快速识别滑坡危险区，为进一步分析矿区滑坡前地表运动特征提供研究范围.

哑巴岭露天采场升轨2(T98)数据集LOS向形变时间序列如图6所示，可知在整个监测期间，哑巴岭采场西南侧边坡中心区域随时间逐渐出现了较大量级和范围的沉降，最终在2019年11月3日至11月15日的监测结果中达到最大形变值.值得注意的是，最后一幅影像的拍摄时间距离滑坡事件发生还有10天，按照这种沉降模式的发展趋势，后续的10天中其最大沉降值将会继续增加.

图6 哑巴岭露天采场升轨2(T98)数据集LOS向形变时间序列图Fig.6 Time series graphs of deformation in LOS direction of Yabaling open-pit slope ascending 2(T98)dataset

哑巴岭露天采场降轨(T105)数据集所获取的LOS向形变时间序列如图7所示，在整个监测期间，哑巴岭采场西南侧边坡中心区域的LOS向形变值变化范围及量级均较小.采场西南侧边坡中心区域整体变化趋势与升轨数据集的监测结果相反，形变值多为正值(蓝色)，表现为朝向卫星的地面运动，无法通过该数据集的监测结果直接快速判定易滑坡危险区的范围，需要对该数据集的监测结果再进行详细分析.

图7 哑巴岭露天采场降轨(T105)数据集LOS向形变时间序列图Fig.7 Time series graph of deformation in LOS direction of Yabaling open-pit slope descending(T105)dataset

综上所述，升轨1(T25)与升轨2(T98)数据集所获取的鞍千哑巴岭露天采场的形变场基本一致，在采场的西南侧边坡中心区域出现了较大范围和量级的地表形变位移.同样，在两组升轨(T25，T98)数据集的监测结果与降轨(T05)数据集所获取的监测结果中也发现了部分区域在形变范围和量级上存在较大差异.为了定量分析该采场滑坡前的地表形变时空演化规律，需要对这三个不同轨道数据集时序InSAR处理后的结果再进行分析.

提取跨越整个异常形变区域的三条剖线AA′，BB′及CC′如图8所示，分别在每一条剖线上选取5个特征点(共计15个特征点)对其进行时间序列形变分析.

图8 剖面线及特征点地理位置分布图Fig.8 Geographical location of the profile lines and feature points

1) 纵剖面AA′：升轨1(T25)、升轨2(T98)及降轨(T105)纵剖面AA′特征点位移时间序列图如图9所示.

升轨1(T25)和升轨2(T98)数据集的结果显示所有特征点均表现出了持续下沉的趋势，且该断面的中心区域P04，P07两个特征点的累计形变值最大，如图9a,9b所示.从图9c的降轨(T105)

图9 升轨1(T25)、升轨2(T98)及降轨(T105) 纵剖面AA′特征点位移时间序列图Fig.9 Ascending1(T25), ascending 2(T98)and descending(T105)displacement time series diagram of the characteristic pointsalong the longitudinal section AA′(a)—升轨1(T25); (b)—升轨2(T98); (c)—降轨(T105).

数据集所获取的形变结果可知，该断面所有特征点的形变范围较小，整体呈线性变化，且变化趋势很小.采矿所引发的地表形变与时间常常呈非线性关系.从图9b升轨2(T98)中的LOS向形变时间序列中不难发现该断面所有的特征点在2019年6月12日至2019年10月10日(红色虚线处)均处于一个缓慢线性下沉的趋势，但在2019年10月10日(红色虚线处)至2019年11月15日忽然出现了加速沉降的现象，这一形变演化过程对应于矿区单点沉降中的初始沉降、加速沉降两个过程.说明在发生滑坡前该断面所有的点均存在一个加速下沉的特征，这符合采矿引发地表形变的规律，可以作为分析矿区滑坡前地表运动特征的参考依据.

2) 纵剖面BB′：升轨1(T25)、升轨2(T98)及降轨(T105)纵剖面BB′特征点位移时间序列如图10所示.

图10 升轨1(T25)、升轨2(T98)及降轨(T105) 纵剖面BB′特征点位移时间序列图Fig.10 Ascending1(T25), ascending2(T98)and descending(T105)displacement time series diagram of characteristic points along the longitudinal section BB′(a)—升轨1(T25); (b)—升轨2(T98); (c)—降轨(T105).

升轨1(T25)和升轨2(T98)数据集的结果显示所有特征点均表现出持续下沉的趋势，但位于纵剖面两端的两点P02和P04的累计形变值较小，中心区域P05，P08，P11的累计形变值较大，说明形变的中心区域位于特征点P05，P08，P11之间.同样的，在图10b升轨2(T98)中不难发现P05，P08，P11在2019年6月12日至2019年10月10日(红色虚线处)期间，所有特征点均处于一个缓慢线性沉降的趋势，但在2019年10月10日(红色虚线处)至2019年11月15日，P05，P08，P11这三个特征点忽然出现了加速沉降的现象，说明在发生滑坡前该断面中心区域的点均存在一个加速下沉的特征.从图10c降轨(T105)数据集所获取的形变结果可知，在整个监测期间，该断面所有特征点的形变范围相比纵剖面AA′的范围有所扩大，整体仍呈线性变化，且变化趋势很小.

3) 纵剖面CC′：升轨1(T25)、升轨2(T98)及降轨(T105)纵剖面CC′特征点位移时间序列如图11所示.

升轨1(T25)数据集的结果显示，所有特征点的累计沉降量均小于纵剖面AA′与BB′，如图11a所示.升轨2(T98)数据集的结果显示，位于纵剖面CC′上所有的特征点在整个监测期间均以线性形式进行沉降，但沉降速率及沉降范围相比纵剖面AA′与CC′都要小，未能发现存在加速下沉的特征点，如图11b所示.降轨(T105)数据集的结果显示，在整个监测期间，该断面所有特征点仍呈线性变化，且变化趋势很小，如图11c所示.以上结果表明，纵剖面CC′已经是形变范围的边界区域，属于线性稳定沉降区域.

图11 升轨1(T25)、升轨2(T98)及降轨(T105) 纵剖面CC′特征点位移时间序列图Fig.11 Ascending1(T25), ascending2(T98) and descending(T105)displacement time series diagram of characteristic points along the longitudinal section CC′(a)—升轨1(T25); (b)—升轨2(T98); (c)—降轨(T105).

通过对纵剖面AA′、BB′与CC′的分析以及从升轨2(T98)数据集所获取的所有特征点的形变时间序列图来看，不难发现P01，P04，P05，P07，P08，P10，P11，P13这8个特征点在整个监测期间存在滑坡前的加速时期(2019年10月10日至2019年11月15日)，如图12所示.图中，红色点表示这些出现异常加速现象的点，白色点表示在整个监测期间表现为线性匀速沉降的点.

图12 不同沉降模式特征点分布图Fig.12 Distribution map of characteristic points of different settlement modes

以线性行为为特征区域，是极限平衡状态下的不稳定性指标.没有迹象表明这种行为可以持续多长时间，或可能会演变为加速行为或稳定行为.在这种情况下，如果这些稳定区域围绕加速区域，就极有可能受到加速区域移动的影响[15].从2019年10月10日至2019年11月25日(滑坡事件发生时)在哑巴岭矿区滑坡位置顶部后缘位置存在一个扇形加速沉降区域，如图13所示.这成为直接诱导滑坡事件发生的主要因素，根据这一现象可以圈定出一个可能发生滑坡的范围，即易滑坡区.

图13 易滑坡区示意图Fig.13 Sketch map of the landslide area

4 InSAR多轨道观测值之间的可靠性检验

为了解释来自三组不同观测轨道InSAR测量结果之间的可靠性，比较分析了从三组不同的哨兵1A / B(Sentinel-1 A / B)数据集获取的鞍千哑巴岭采场西南侧边坡形变区域中所有像素的InSAR平均形变测量值，整个过程分为以下两个步骤：

1) 将时序InSAR反演所获得的LOS向形变场重采样为100 m×200 m的网格，从而达到最小化地理定位误差(解决测量之间的空间不匹配问题)的目的；2)异常形变区域(蓝色)光学影像如图14所示，选择了三个轨道共有的局部形变区域中所有形变监测数据不为空值的像素共计324个.异常形变区域所有像素平均形变时间序列图如图15所示，统计了三条不同轨道数据集在该形变区域内所有像素LOS向平均形变值的时间序列，并分别进行线性拟合.根据线性拟合的结果可以发现：升轨1(T25)与升轨2(T98)的判定系数R2分别为0.89与0.98，基本表现出相同的沉降趋势，而降轨(T105)与升轨1(T25)和升轨2(T98)趋势线呈相反关系.这些现象表明从不同数据集获得的不同沉降模式彼此相关.升轨的两个数据集监测所得到的结果在2019年6月至2019年11月之间表现出相似的线性运动，使这两个数据集所得到的监测结果具有可信度.

图14 异常形变区域光学影像Fig.14 The optical image of the abnormal deformation area

图15 异常形变区域所有像素平均形变时间序列图Fig.15 Time series graph of the average deformation of all pixels in the abnormal deformation area

同时不难发现，只有升轨2(T98)数据集的监测结果显示了在滑坡前存在一个突然加速的区域使得整个滑坡事件得以发生.对整个形变区域平均形变时间序列线性拟合结果进行分析，升轨2(T98)数据集的判定系数R2为0.98，这说明整个大的形变区域是以线性模式进行均匀沉降，这符合滑坡事件发生的基本沉降规律，所以针对鞍千哑巴岭露天采场，以整个大范围的平均形变时间序列的线性拟合效果的差异来判断不同轨道数据集之间监测结果的可靠性是合理的.

比较相同或不同轨道之间时序InSAR反演结果是定量检查测量一致性的常见做法[16].本次实验所观察到的三组数据集之间的差异可能源于采样不足、地理参考不确定、季节性偏差、长时间间隙、不同的影像获取时间及不同的轨道参数等因素.在本次实验数据集中观察到的升轨和降轨所获取的结果之间产生相反的变化趋势.原因主要包括：三组不同轨道数据集拥有不同的入射角及航向角，其对研究区域真实形变的敏感度不同，由于采矿所引发的地表沉降导致其LOS向形变中可能会包含水平形变的影响，这些因素共同造成了差异性.

5 矿区滑坡因素分析

影响岩质边坡稳定性的因素可分为两个方面，即内部与外部因素.内部因素为哑巴岭采场西南侧岩质边坡的一些固有属性，而外部因素通常会引起边坡失稳从而发生滑坡，如工程爆破、降雨、地震等.尽管目前没有特别有力的证据表明鞍千哑巴岭露天采场在2019年11月25日的滑坡是由于采场西南侧边坡的不稳定因素而引起的，但考虑到该矿区的特殊性质，强烈和频繁的开采活动及矿石爆破工作有可能构成了边坡失稳的主要外部因素.

孔隙压力是滑坡的主要驱动因素[17]，季节性滑动速率的变化与降水驱动的瞬时孔隙压力的增加有关.鞍山地区位于辽河平原的东南部边缘，暖温带大陆性季风普遍影响当地气候.这使得该地区四季分明、干冷同期、雨热同季.鞍山地区2019年平均降水量约为686 mm，其中85%的降雨集中在5月至9月，但该年11月份的降水量却异常增多，占据非雨季降水量的34%.2019年11月25日位于该地区的哑巴岭露天采场边坡发生一次大型滑坡.辽宁鞍山53339号区站2019年11月的日降水量数据和鞍山地区11月份的平均温气与天气数据如图16所示.

图16 鞍山地区2019年11月部分气象数据Fig.16 Part of the meteorological data of Anshan area in November 2019

由图可知：11月12日开始到滑坡前陆陆续续一直在下雨或下雪，11月19日气温升高，在之前积雪融化的基础上又开始降雨.由于前期降水持续时间较长，降水量也很大，尤其是在发生滑坡的前两天，直接导致哑巴岭西南侧边坡裂隙中进入大量的雨水，长时间大范围的降水导致软弱层(滑面)形成，诱导本来处于线性稳定沉降的坡面发生滑坡事件.