基于多轨道SBAS方法的露天矿高陡边坡形变监测

魏恋欢, 封秋月, 毛亚纯, 刘善军

(东北大学 资源与土木工程学院, 辽宁 沈阳 110819)

我国冶金矿山80%的矿石来源于露天开采,随着开采深度的不断增加,我国大多数露天矿已进入深凹开采模式,导致矿坑边坡高度及边坡角不断增大,边坡稳定性随之降低.对于深凹露天矿而言,高陡边坡的稳定性是影响矿山安全生产的重要因素.露天矿边坡滑移灾害的发生表面上是突发的,但其过程是逐步酝酿和发展的,在灾害发生前会存在长期的缓慢形变积累.因此,采取有效手段实时监测高陡边坡的稳定性，研究其形变特点及成因、及时治理潜在的险情，是矿山安全生产的重要保障.

目前,露天矿高陡边坡形变监测一般采用全站仪或GPS等现场测量手段,存在效率低、人力成本高、仅能获得少量点位形变、无法提取大范围分布特征等不足.近年来,随着合成孔径雷达(synthetic aperture radar,SAR)卫星与数据处理技术的不断发展,时序InSAR(time series interferometry SAR)作为一种新型的形变测量手段,具有低成本、高精度等优点,已成功应用于沉降、滑坡、地震、火山活动等地表形变监测中[1].时序InSAR从一组时间序列SAR影像中提取不受时空去相干影响的目标像元,每一个像元的干涉相位φ可分解为如下几个部分:φflat，φtopo，φdef，φnoise.

平地效应相位φflat和地形相位φtopo是由于SAR卫星多次成像时所处位置不同、地球椭球和地表高程变化造成的;φdef为形变相位;φnoise主要是由轨道误差、大气效应等造成的噪声相位[1].针对这些像元开展相位分析,可达到分离各种分量并解算精细形变的目的,实现毫米级的形变监测精度.目前已有的时序InSAR算法主要分为基于点目标的永久散射体法(permanent scatterers InSAR,PS-InSAR)和基于分布式目标的小基线子集法(small baseline subset,SBAS)两大类[2-7].相比于PS-InSAR,SBAS可以有效利用所有时空基线较短的高质量差分干涉图,提高了数据利用率和点密度,更加适合进行长时间的形变监测.

受限于SAR卫星固有的侧视成像特点,时序InSAR得出的是目标像元沿卫星视线(line of sight,LOS)方向的一维形变信息,是真实形变在LOS方向的投影,因此需要由LOS向形变求解真实形变.在大范围沉降监测中,通常认为地物只存在垂直形变,不存在水平形变,简单地用LOS方向的形变结果除以影像入射角的余弦值即可得到垂直形变结果[8-9].然而,在露天矿高陡边坡形变监测中,由于矿坑边坡滑移通常都是沿坡面梯度方向的,而不同空间位置的坡度和坡向差别很大,无法直接根据入射角求解[10].

本文针对上述问题,介绍了SBAS法的基本原理,面向深凹露天矿高陡边坡形变监测需求,提出了一种基于地形的坡向形变求解方法,利用SBAS法分析了覆盖鞍钢集团大孤山铁矿的3组Sentinel-1数据,提取了3组LOS方向形变参数,将3组LOS方向形变分别转换到坡向,与同时期的测量机器人监测数据及降水量进行了对比分析,揭示了大孤山铁矿矿坑边坡的稳定性及影响因素.

1 小基线子集法

作为时序InSAR方法中的一个典型代表,SBAS是在PS-InSAR的基础上发展而来.PS-InSAR从一组SAR影像集中选取一个公共主影像,将所有从影像与主影像干涉得到时序干涉图集.公共主影像的使用使一些时空基线均较大的低质量干涉图参与了PS点目标提取与形变解算,导致人工地物稀少的区域点目标密度偏低,难以精确估算形变分布信息[2-3].SBAS针对这一问题,采用了多主影像干涉模式,根据时空基线分布特征将SAR数据集划分为多个子集,每个子集内时空基线均较短,各子集间通过少量长时空基线的干涉对相连,可以保证数据的高时间采样率和高点密度,降低时空去相干的影响.SBAS的多主影像干涉模式增加了干涉图数量,避免了因干涉图数量少导致的方程秩亏问题,短基线干涉组合可以降低数字高程模型(digital elevation model,DEM)误差对形变量测精度的影响,有利于直接解算研究区域的形变[6].因此,本文采用SBAS法对大孤山铁矿的高陡边坡进行时序InSAR分析.

在SBAS法中,首先根据数据集的时空基线分布特征选择超级主影像,将所有从影像都以超级主影像为参考进行粗配准和精配准.其次,根据多主影像策略生成小基线差分干涉图集.再次,因为干涉图中获得的相位是真实相位对2π的余数值,其取值范围是在-π到π之间的，称为“缠绕相位”，不能直接转换为形变值。通过最小费用流算法对序列差分干涉图进行解缠，将“缠绕相位”加上2π的整数倍，即可恢复为“真实相位”，称为“解缠相位”。用奇异值分解(singular value decomposition,SVD)法将解缠后的差分干涉图结合起来生成一个包含大气、地形和形变信号的时间序列.然后,通过二次解缠和滤波消除地形、大气信号,最终得到LOS方向形变速率和形变量[4,11].

2 边坡坡向形变求解方法

以图1a中的A点为例，假设其沿坡度方向的位移为dA，可按如下公式分解为水平分量和垂直分量:

(1)

式中：β为坡度；dA,v与dA,h分别为dA的垂直与水平分量。SBAS提取的LOS方向形变值dA,los实际上是dA,v和dA,h在LOS向的投影之和。其中dA,v在LOS向的投影可直接根据传感器的入射角θ计算。为了求解dA,h在LOS向的投影，需要先将dA,h投影到LOS向所在的铅垂面[10].由于升轨数据和降轨数据的LOS向不同,需要分别讨论.这里“升轨”是指卫星由南向北飞行,“降轨”是指卫星由北向南飞行.受到地球自转的影响,升轨时卫星实际前进方向为北偏西,相应的LOS方向在水平面的投影为东偏北；降轨时卫星实际前进方向为南偏西,相应的LOS方向在水平面的投影为西偏北.

根据降轨数据的几何关系(图1a)可知，dA,h在LOS所在铅垂面的投影dA,h-los可表示为

dA,h-los=dA,h·sin(αA-α0).

(2)

其中：αA是dA,h与北向的夹角；α0是雷达的航向角。设传感器入射角为θ，A点的LOS向形变dA,los由dA,v和dA,h-los在LOS向投影的矢量和组成，可表示为

dA,los=dA,h-los·sin(θ)-dA,v·cos(θ).

(3)

将式(1)、式(2)代入式(3)可得：

dA,los=dA·cos(β)·sin(αA-α0)·sin(θ)-dA·sin(β)·cos(θ).

(4)

因此，在降轨数据条件下，A点的LOS向形变可转换为坡向形变：

(5)

类似地，根据升轨数据的几何成像(图1b)可知，dA,h在LOS所在铅垂面的投影dA,h-los可表示为

dA,h-los=dA,h·sin(αA-α0-π)=-dA,h·sin(αA-α0).

(6)

A点的LOS向形变dA,los由dA,v和dA,h-los在LOS向投影的矢量和组成，表示为

dA,los=-dA,h-los·sin(θ)-dA,v·cos(θ).

(7)

将式(1)、式(6)代入式(7)后可以得到dA和dA,los的关系，与式(4)相同。因此，在升轨数据条件下，A点的LOS向形变转为坡向形变的公式与式(5)相同[10]。在本文实验中，三组数据的航向角(以北向起始顺时针为正)α0和入射角θ如表1所示。由于大孤山铁矿矿坑西北帮岩体破碎，工程地质条件复杂，为了保证安全生产，开采坡度角β约为37°，坡度方向与北向的水平夹角与点目标所在位置有关.

表1 时序SAR数据信息

3 研究区及数据概况

3.1 研究区概况

大孤山铁矿位于我国辽宁省鞍山市,是鞍山钢铁集团鞍山矿业公司5大露天铁矿之一,是我国典型的短深型深凹露天矿,也是亚洲最深的露天铁矿之一[12].大孤山铁矿于1970年转入深凹露天纵向开采,形成了东西长1 700 m、南北宽1 520 m、垂直深度达348—398 m左右的椭圆形露天采坑[12].矿区年降雨量平均为720 mm,一年中以6—8月降雨量最多,合计438 mm,9—11月降雨量次之,为151 mm.

如图2所示,大孤山铁矿的采坑(黄框区域)位于矿区西部,开采技术为高陡边坡开采,露采工作呈阶梯形向地下延伸,其东侧及东北侧的“倒L”形区域为排土场.该排土场大部分区域为岩石松散裸露堆积.排土场的东南侧是大球场尾矿库.大球场尾矿库属于山谷型尾矿库,沟口向北,其余方向均为山脊,尾矿库坝的筑坝方式是排土场渣土分层碾压筑坝.对大孤山铁矿的高陡边坡的稳定性进行持续监测至关重要,是确保矿山安全生产的关键.

3.2 数据概况

本文收集了覆盖大孤山铁矿的3组Sentinel-1数据,空间覆盖范围如图3所示,详细参数如表1所示.本文采用的DEM数据为德国TanDEM-X计划获取的DEM,空间分辨率为3 rad·s,绝对高程精度约1 m[13-14].由于SBAS法使用的差分干涉图时空基线都相对较短,DEM误差对差分干涉质量的影响相对较小.

图4中展示了3组Sentinel-1数据的获取时间以及SBAS法生成的差分干涉组合.其中,升轨Sentinel-1A数据生成了145幅差分干涉图；轨道编号为P105的降轨Sentinel-1B数据生成了250幅差分干涉图；轨道编号为P3的降轨Sentinel-1B数据生成了229幅差分干涉图.

4 结果分析

4.1 监测结果分析

SBAS法提取的大孤山铁矿采坑区域的LOS向平均形变速率如图5所示,其中红色色标表示地物目标向远离传感器方向移动,蓝色色标表示地物目标向靠近传感器方向移动；黄色箭头指示的Az方向表示卫星飞行方向,LOS方向表示卫星视线方向.交叉验证表明,三组Sentinel-1结果具有较高一致性,监测结果较可靠.大孤山铁矿高陡边坡的主要形变区域位于采坑西北帮,已在图5中用绿色多边形标出,该区域的大部分地物目标LOS向形变速率均超过了-60 mm/年,最高达到了-221 mm/年.

采坑西北帮的形变可大致分为A，B，C三个小区域.其中,A,B区域在升轨结果中为远离传感器方向形变,在两组降轨结果中均为朝向传感器方向形变.由于边坡的真实形变一般是沿坡向朝向坑底的,SBAS结果仅为真实形变在LOS向的投影.考虑到升降轨数据的LOS方向不同,升降轨结果存在差异是合理的.C区域在三组结果中都显示为远离传感器方向形变,主要是由于该区域位于矿坑边缘部位、坡度相对较小且坡向接近南向,因此升降轨结果中的差异不明显.

图6中的Google Earth历史影像显示,2018年9月10日前A区域发生了大规模的岩体移动,岩体的移动方向为沿坡面指向坑底方向,符合露天开挖导致的滑移特征.根据本文提出的方法,求解了位于A区域滑坡体上缘P1点的坡向形变.图7展示了P1点的坡向形变量曲线,为了方便对比,三组数据都以2017年6月5日为参考,截止2019年6月初累计坡向滑移量均超过了-600 mm.

大多数露天矿的滑坡均起因于降水引起的沿软弱结构面的滑动,每年的雨季都是滑坡灾害发生最频繁的时期.降水对边坡稳定性的影响较大,水的浸泡、渗透作用可以降低软弱面的抗剪强度,增加下滑力,从而导致边坡破坏[15].图7显示,2017年和2018年夏季持续大量降水之后,P1点在三组结果中均表现出了加速形变,入冬后逐渐减速.不过,总体而言P1点的加速形变不大,未出现严重的边坡失稳问题.

在现场布设了2个监测点(图6中P2，P3),位于B区域,于2017年8月13日至2018年11月1日期间使用测量机器人进行了连续监测.图8所示的是P2和P3点的坡向形变时间序列,为了方便对比分析,3组SBAS结果都以测量机器人第一次观测时刻(2017年8月13日)为参考.监测结果表明,3组SBAS结果与测量机器人监测结果存在高度一致性,P2点和P3点形变规律相似,2017年8月13日至2019年6月初累计坡向形变均超过了-500 mm,均在每年夏季大量降水后产生了一定的加速形变,入冬时均有减速,未发生明显失稳问题.

P2点和P3点的SBAS监测结果与测量机器人监测结果之间的均方根误差如表2所示,其中S1B_P105数据的均方根误差最小,分别为17.4 mm和19.2 mm.时间序列结果表明,持续大量降水会导致边坡形变量和形变速率的增大,这也是在雨季提高测量机器人的监测频率、缩短监测周期的主要原因.为了保障安全生产,应在每年雨季对西北帮高陡边坡实施稳定性监测.

表2 P2，P3点均方根误差

4.2 工程地质条件分析

工程地质资料(图9)显示,大孤山铁矿露天采坑的东端帮与南帮以花岗岩为主,岩性及构造条件良好,总体较为稳定.西帮(包括西北帮和西南帮)可分为典型的4个区域[16].西南帮为混合岩区域,也是典型的块裂体工程地质模型,但节理发育情况好于西北帮,且断层较少,稳定性较好.西北帮的岩体较为破碎,工程地质条件复杂,其稳定性直接影响矿山的安全生产.

西北帮可大致分为A,B,C三个区域.A区域为品位较低的楔形矿体,其西南侧以F15断层为界,东北侧以F14断层与B区域接触,F14断层周边的破碎带稳定条件较差,威胁着采场北帮的铁路运输；B区域的岩体以混合岩为主,其稳定性受几条较大的断层控制,属于典型的块裂体工程地质模型[17]；C区域的岩体以千枚岩为主,为典型的碎裂和散体工程地质模型,节理裂隙非常发育,岩石强度较低,在雨水的作用下容易发生滑坡现象[17].虽然B,C区域在Google Earth历史影像中并未出现明显的边坡滑动,但从岩土结构上看,易发生形变破坏的地段多为千枚岩地段和混合岩地段(B,C区域),这主要是由于该类岩体整体完整度低,坚固程度差所造成的.考虑到C区域位于矿坑边缘且所占范围逐渐减小,尽管千枚岩的抗压、抗剪强度均较低,易产生形变破坏,其稳定性对矿山安全生产的威胁不大[17].因此,B区域是威胁大孤山铁矿安全生产的潜在滑坡区,应持续进行形变监测.

5 结 论

1) 利用SBAS法处理了多轨道时序SAR数据,对大孤山露天矿的高陡边坡进行了稳定性监测.

2) 提出了一种基于地形特征的LOS向形变到坡向形变的解算方法,得到了采坑西北帮的坡向形变结果,与实测数据存在高度一致性.

3) 通过研究同期的降水数据及大孤山工程地质条件,对变形原因进行分析,表明露天矿的稳定性主要受岩体结构、岩性与降水的影响.