我国冶金矿山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可表示为
(2) |
其中:αA是dA, h与北向的夹角;α0是雷达的航向角。设传感器入射角为θ,A点的LOS向形变dA, los由dA, v和dA, h-los在LOS向投影的矢量和组成,可表示为
(3) |
将式(1)、式(2)代入式(3)可得:
(4) |
因此,在降轨数据条件下,A点的LOS向形变可转换为坡向形变:
(5) |
类似地,根据升轨数据的几何成像(图 1b)可知,dA, h在LOS所在铅垂面的投影dA, h-los可表示为
(6) |
A点的LOS向形变dA, los由dA, v和dA, h-los在LOS向投影的矢量和组成,表示为
(7) |
将式(1)、式(6)代入式(7)后可以得到dA和dA, los的关系,与式(4)相同。因此,在升轨数据条件下,A点的LOS向形变转为坡向形变的公式与式(5)相同[10]。在本文实验中,三组数据的航向角(以北向起始顺时针为正)α0和入射角θ如表 1所示。由于大孤山铁矿矿坑西北帮岩体破碎,工程地质条件复杂,为了保证安全生产,开采坡度角β约为37°,坡度方向与北向的水平夹角与点目标所在位置有关.
大孤山铁矿位于我国辽宁省鞍山市, 是鞍山钢铁集团鞍山矿业公司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组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幅差分干涉图.
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.时间序列结果表明, 持续大量降水会导致边坡形变量和形变速率的增大, 这也是在雨季提高测量机器人的监测频率、缩短监测周期的主要原因.为了保障安全生产, 应在每年雨季对西北帮高陡边坡实施稳定性监测.
工程地质资料(图 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) 通过研究同期的降水数据及大孤山工程地质条件, 对变形原因进行分析, 表明露天矿的稳定性主要受岩体结构、岩性与降水的影响.
[1] |
廖明生, 王腾. 时间序列InSAR技术与应用[M]. 北京: 科学出版社, 2014. (Liao Ming-sheng, Wang Teng. Time series InSAR technology and application[M]. Beijing: Science Press, 2014.) |
[2] |
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 |
[3] |
Ferretti A. Nonlinear subsidence rate estimation using permanent scatters in differential SAR interferometry[J]. IEEE Transactions on Geoscience & Remote Sensing, 2000, 38(5): 2202-2212. |
[4] |
Berardino P, Fornaro G, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]. IEEE Transactions on Geoscience & Remote Sensing, 2002, 40(11): 2375-2383. |
[5] |
Hooper A. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches[J]. Geophysical Research Letters, 2008, 35: L16302. DOI:10.1029/2008GL034654 |
[6] |
Amelung F, Galloway D L, Bell J W, et al. Sensing the ups and downs of Las Vegas:InSAR reveals structural control of land subsidence and aquifer-system deformation[J]. Geology, 1999, 27(6): 483-486. DOI:10.1130/0091-7613(1999)027<0483:STUADO>2.3.CO;2 |
[7] |
Lanari R, Mora O, Manunta M, et al. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42: 1377-1386. DOI:10.1109/TGRS.2004.828196 |
[8] |
Galloway D L, Hudnut K W, Ingebritsen S E, et al. Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope Valley, Mojave Desert, California[J]. Water Resources Research, 1998, 34(10): 2573-2585. DOI:10.1029/98WR01285 |
[9] |
Zhao C Y, Zhang Q, Ding X L, et al. Monitoring of land subsidence and ground fissures in Xi'an, China 2005—2006:mapped by SAR interferometry[J]. Environmental Geology, 2009, 58: 1533-1540. DOI:10.1007/s00254-008-1654-9 |
[10] |
Wei L H, Zhang Y, Zhao Z G, et al. Analysis of mining waste dump site stability based on multiple remote sensing technologies[J]. Remote Sensing, 2018, 10(12): 20-25. |
[11] |
Fariba M, Bahram S, Masoud M, et al. Monitoring surface changes in discontinuous permafrost terrain using small baseline SAR interferometry, object-based classification, and geological features:a case study from Mayo, Yukon Territory, Canada[J]. GIScience & Remote Sensing, 2018, 56(4): 485-510. |
[12] |
刘治保. 鞍山矿山区域生态修复与景观再生研究——以鞍山大孤山铁矿矿山区域为例[J]. 美与时代(城市版), 2017(3): 63-64. (Liu Zhi-bao. Regional ecological restoration and landscape regeneration of mining areas in Anshan—a case study on the Dagushan Iron Mine[J]. Beauty and Time(City Edition), 2017(3): 63-64.) |
[13] |
Huber M, Gruber A, Wendleder A, et al.The global TanDEM-X DEM: production status and first validation results[C]//International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences.Melbourne, Australia, 2012: 45-50.
|
[14] |
The TanDEM-X 90m digital elevation model[EB/OL].[2018-10-31].https://geoservice.dlr.de/web/dataguide/tdm90/.
|
[15] |
王忠红, 刘华艳, 洪利, 等. 深凹露天矿边坡破坏成因分析及治理建议[J]. 矿业工程, 2008, 6(5): 23-24. (Wang Zhong-hong, Liu Hua-yan, Hong Li, et al. Analysis of reasons causing slope damage and its prevention in deep mining open pit mine[J]. Mining Engineering, 2008, 6(5): 23-24.) |
[16] |
胡颖鹏.大孤山铁矿露天地下协同开采方案研究[D].沈阳: 东北大学, 2015. (Hu Ying-peng.Study on open-pit and underground collaborative mining scheme of Dagushan Iron Mine[D].Shenyang: Northeastern University, 2015. ) |
[17] |
张文钦, 金锡古. 大孤山高陡边坡稳定性关键技术研究[J]. 有色矿冶, 2012, 28(3): 12-15. (Zhang Wen-qin, Jin Xi-gu. Research of the key technologies for the high steep slope stability in Dagushan Iron Mine[J]. Nonferrous Mining and Metallurgy, 2012, 28(3): 12-15.) |