唐山马兰庄铁矿露天开采边坡变形监测的GB-InSAR技术

刘作利 刘景玉 申修强 马 辉 田卫明 左荣虎

(1.唐山首钢马兰庄铁矿有限责任公司；2.北京理工大学信息与电子学院；3.北京理工雷科电子信息技术有限公司)

我国有露天矿山约7万座，尾矿库1. 2万座，排土场10万余座，该类边坡一旦失稳，不仅会影响企业生产，更会威胁作业人员的生命安全[1]。因此，通过对边坡变形进行实时监测，并对监测数据进行一系列处理分析，对于准确掌握边坡灾害发生、发展和变化机理，进而实现边坡灾害预警具有重要意义[2-3]。全站仪、测量机器人和水准仪等传统测量技术在变形监测中已经得到了广泛应用，但精度易受气候、时间、通视等条件的限制[4]。

合成孔径雷达(Synthetic aperture radar, SAR)是一种基于微波传感器的雷达，具有全天候和一定穿透性等优点。差分干涉合成孔径雷达技术(Differential interferometric synthetic aperture radar，D-InSAR)是SAR的一个重要应用，在近十余年内得到了快速发展和广泛应用。星载和地基干涉合成孔径雷达技术是干涉合成孔径雷达技术(Interferometric synthetic aperture radar，InSAR)的2种重要形式。星载干涉合成孔径雷达具有获取变形范围大的优势，可应用于大区域地表沉降监测[5-9]，但由于星载干涉合成孔径雷达具有较长重返周期、固定成像姿态和低空间分辨率等不足，导致其在边坡变形监测中效果不理想。地基干涉合成孔径雷达测量技术(GB-InSAR)是星载干涉合成孔径雷达干涉测量技术的有效补充，不仅具有最优观测姿态和连续观测性能，而且具有灵活多变、分辨率高、平台稳定、观测周期短、造价相对低廉等优点，属于非接触测量方法范畴，适用于对危险边坡进行实时变形监测[10-12]。本研究以唐山马兰庄铁矿为例，针对该矿边坡特点，采用GB-InSAR技术进行边坡变形监测，结合监测数据处理分析结果，进一步评价利用该技术进行边坡变形监测的可行性。

1 GB-InSAR技术原理

1.1 SAR技术原理

SAR是利用一个小的真实天线的运动来等效构成一个长天线，从而可以在运动方向上获得一个等效的大孔径(合成孔径)天线，利用该方式，雷达对目标的分辨率将会有所提高。根据合成孔径成像技术原理，设雷达系统发射的信号带宽为B，电磁波在真空中的传播速度为c，则雷达斜距向分辨率为c/2B[13]。当信号带宽为300 MHz时，系统能达到的最大距离向分辨率为0.5 m。若雷达发射的电磁波波长为λ，雷达滑动的最大距离为L，则雷达方位角分辨率为λ/2L。

本研究采用的GB-InSAR雷达系统信号处于Ku波段(12.5～18 GHz)，可以获得的最高方位向分辨率约为4 mrad。通过将SAR与调频连续波技术(Stepped frequency continuous wave,FMCW)相结合，监测区域可以被分割为较多二维小单元(图1)。

图1 GB-InSAR监测空间分辨单元

GB-InSAR雷达系统在获取目标区域数据过程中，由于天线运动轨道固定，对某一区域目标重复监测成像时天线都处于相同的空间位置，因而空间基线为0。雷达信号往返路径之差即为雷达视线方向目标距离的变化量，具体可以表现为相位差[13]。

1.2 地基雷达差分干涉变形测量原理

利用D-InSAR技术可以对SAR图像进行干涉处理并去除干涉图中的非形变相位，实现对地表进行大范围、高精度、全天候变形监测。由于地基雷达主机部分始终在线性导轨上滑动，滑轨的起点和终点位置始终未发生变动，SAR图像数据的轨道是严格重合的，即地基雷达并不存在空间基线。如此，地基雷达和星载雷达的干涉测量方式有较大区别，地基雷达差分干涉处理时，无需去平地效应和去地形效应，因此从理论上分析，差分干涉SAR图像的相位信息仅包含变形信息。GB-InSAR可通过测量视线方向上回波信号的相位变化实现高精度变形测量(图2)，视线方向的形变量Δd与干涉相位φdef关系可以表示为

(1)

式中，λ为雷达信号的波长。

图2 差分干涉雷达变形监测原理

由干涉相位解缠理论[9]可知，干涉图中相邻2个永久散射体(Permanent scatterer,PS)点位之间的相位差不大于1/4λ，即为边坡雷达变形监测的最大形变分辨能力。因此，地基差分干涉SAR图像的变形测量分辨能力与雷达波长密切相关。如本研究采用的Ku波段雷达波长约为1.8 cm，则相邻2个形变点位的最大形变分辨能力为4.5 mm。根据差分干涉雷达的基本理论，目标点的干涉相位可以表示为

φInt=φdef+Δφatmosphere+Δφtrack+Δφnoise，

(2)

式中，φdef为地表变形相位；Δφatmosphere为大气相位；Δφtrack为由定轨误差引入的相位；Δφnoise为由时间去相关等因素引入的随机相位。

当式(2)中的后3项误差较大时将无法获得高精度的地表变形信息。为此，在本研究边坡雷达系统中采用了永久性散射体技术(PS)。该方法首先在SAR图像中根据幅度阈值、幅度离差阈值、相关性阈值等指标优选出高相关性的目标点(PS点)；然后对Δφatmosphere、Δφtrack、Δφnoise进行建模并估计补偿大气误差、轨道误差和噪声误差相位；最后获取PS点的形变信息。数据解算过程的重要环节为相位解缠[9]，特别是在GB-InSAR矿区边坡变形监测数据处理中，直接影响了边坡变形监测精度。进行相位解缠时首先选择1组永久散射体候选子集，通过构建三角网进行插值处理，而后进行空间维和时间维解缠，最终通过处理分析永久散射体间的位移信息来消除大气、轨道等噪声误差，进而获取高精度的变形监测数据。

2 GB-InSAR系统组成与工作流程

2.1 系统组成

GB-InSAR系统的主体部分由高分辨调频连续波(FMCW)雷达、高精度电控位移台组成(图3)。整个系统包括主控计算机设备、远程数据传输系统设备等。调频连续波雷达通过发送线性调频的电磁波信号可以实现快速成像，最快观测一次耗时小于2 min。高精度电控位移台可保障雷达在运行过程中平稳可靠，为雷达合成孔径计算过程的必备外界条件。实时监测软件包括二维和三维2个模式，可以有效地从雷达观测结果中识别出相应的地形位置信息。雷达系统的主控计算机除了安装实时监测软件以外，还可以构建数据库及发布网站，可将雷达观测数据以网站的方式进行实时发布。当雷达与监控室之间无法进行有线网络通讯时，数据远程传输设备可以通过无线方式实时发送雷达观测数据。

图3 边坡监测雷达系统主体部分外观

2.2 系统工作流程

设备工作时，边坡雷达实时处理软件可以实时显示雷达工作状态，可以实现数据实时处理、处理结果实时存储与显示、历史数据实时回放等功能。雷达观测数据通过网络传输至服务器，再利用服务器中构建的数据库对大量的雷达变形监测数据进行统一管理，最终观测结果可以通过多个客户端进行实时显示。本研究雷达系统的变形监测数据可与三维地形信息进行叠加显示，对变形数据可以进行打点曲线分析和实时预警，此外，用户还可以在线查询特定区域在各时期的变形信息。GB-InSAR系统的工作流程如图4所示。

该系统在形变测量的基础上，设计了振动观测环节，通过对目标回波进行计算，从而可以计算出强散射体的振幅和频率值。

3 工程应用

3.1 边坡概况

马兰庄铁矿资源储量为8 912.15万t，SFe平均品位为29.64%，目前露天采矿规模达到 300万t/a，选矿处理原矿能力达180万t/a，精矿产量为80万t/a。该矿位于河北省迁安市马兰庄镇，西距北京约200 km，西南距唐山约80 km。矿区中心的地理坐标为东经118°36'、北纬40°06'。该矿露天矿区呈椭圆状，整个采场长半轴约1 100 m，短半轴约900 m，采场边坡为典型岩质边坡，无植被覆盖，边坡高度落差大，最大开采深度超过200 m；边坡角为38°～47°，在如此陡峭的边坡上安装监测设备难度较大；此外，采场存在雾霾、粉尘，能见度不高。随着该矿露天开采深度进一步增加，边坡稳定性成为制约采场安全生产的重要因素。

图4 边坡雷达系统工作流程

3.2 边坡变形监测方案

根据该矿露天采场的实际条件，结合被监测边坡危险区域的分布特点，将边坡雷达安置于矿坑东坡和西坡地基较稳定处。为使边坡变形监测工作不受雨雪等不利天气条件影响，满足恶劣气象条件下不间断监测要求，在矿坑东侧边坡监测位置建立了活动房(图5)，将GB-InSAR安置于该活动房中进行边坡变形观测，可以有效避免恶劣天气对仪器的影响(图6)。将边坡监测系统服务器设置于矿区机房，监控人员通过客户端可以完成对雷达工作的实时控制以及对监测结果的分析(图7)。GB-InSAR采用24 h不间断监测方式，测量频率为2～ 10 min/次，采用实时处理方式实现对整个边坡表面的连续变形进行实时监测。

图5 雷达活动房

3.3 监测数据分析

矿坑西坡监测时间为2016年9月12日—2016年12月20日。矿坑东边坡监测时间为2016年12月10日—2016年12月20日，监测期分为2段：第1段监测期为2016年12月10日—2016年12月16日；第2段监测期为2016年12月17日—2016年12月20日。矿坑东边坡2次监测的区域完全一致，雷达观测距离为300～900 m，观测范围恰好覆盖整个露天采场的东边坡。GB-InSAR通过主动发射调频连续波信号并接收成像，而后对成像结果进行PS点选择、相位解缠及误差改正，单轨变形结果计算的最短耗时小于2 min，从而实时获取边坡变形信息。矿坑西坡的监测期超过3个月，相对于东坡，尽管监测期较长，但未监测出明显的变形信息，表明在监测期内该区域的边坡基本处于稳定状态。矿坑东坡的监测期仅有11 d，虽然监测期较短，但2个监测时段均监测出矿坑东坡在同一区域内出现形变，该区域为90 m×40 m(长半轴×短半轴)椭圆形，标高由高至低，变形量逐渐增大(图8)。

图6 安装于活动房内的边坡雷达

图7 监控室远程边坡雷达控制

分析图8可知：第1个监测段内的最大变形量达到-150 mm(符号为负，表现为接近雷达方向的视线向形变)；第2个监测段内雷达重新启动，变形结果重新清零，边坡最大变形量约-120 mm，也表现为接近雷达方向的视线向形变。2次观测结果均显示出边坡变形区下部变形速度快、变形量大的特点，2次监测期积累的雷达视线向最大变形量达到-270 mm。

为进一步分析边坡变形特点，针对第2观测期内边坡变形较大的1个点位(图9)进行了打点曲线分析，如图9所示。分析图9可知：该点位变形速率较稳定，观测期内的最大变形量约-120 mm。为更精确地计算该点位的变形速率，利用边坡雷达实时监测软件，在综合分析边坡变形特征及监测误差的基础上，将矿坑东坡第2监测期内的累计变形量与第2个监测期的总监测时间相除，结果如图10所示。分析图10可知：变形速率所显示的位置与变形区域完全一致，局部区域的最大形变速率超过-1 mm/h，同时可以看到变形区域底部的变形量大于上部。

图8 开采东坡形变时间序列

图9 点位变形曲线和变形速率曲线

图10 观测期内平均形变速率

依据本研究DB-InSAR边坡变形监测结果，对边坡变形区进行了现场考察，发现该变形区顶部的边坡岩体产生了明显裂缝，区域上部出现的裂缝显示局部位置变形量超过180 mm(图11)，并且下部变形量大于上部。总体上，雷达监测结果与现场实际情况基本一致。矿山工程技术人员根据本研究边坡变形监测结果，及时对变形区域进行了边坡处理，有效避免了边坡灾害发生。

图11 矿坑东坡现场裂缝

4 结 语

针对唐山马兰庄铁矿露天采场边坡特征，为确保边坡稳定，采用GB-InSAR技术进行了变形监测，在采坑东部边坡发现了一椭圆形变形区域，区域最大累计变形量达到-270 mm，并且该区域底部变形量大于上部，与现场实际情况基本吻合。矿山工程技术人员根据本研究监测结果，及时采取了边坡处理措施，有效避免了边坡灾害发生。研究表明，采用GB-InSAR技术可以对露天采场边坡进行高精度实时变形监测，在一定程度上可以实现边坡灾害预警。

[1] 佴 磊，徐 燕，代树林，等.边坡工程[M].北京：科学出版社，2010.

[2] 郑颖人，陈祖煜，王恭先，等.边坡与滑坡工程治理[M].2版.

北京：人民交通出版社，2010.

[3] 王宏伟，王 攀，王永强，等.马钢凹山露天采场南帮边坡监测及稳定性分析[J].金属矿山，2016(3)：35-37.

[4] 刘 超，高井祥，王 坚，等.GPS/伪卫星技术在露天矿边坡监测中的应用[J].煤炭学报，2010，35(5)：755-759．

[5] 张 洁，胡光道，罗宁波.InSAR技术在滑坡监测中的应用研究[J].工程地球物理学报，2004(2)：147-152.

[6] BERARDINO P,COSTANTINI M,FRANCESCHETTI G,et al.Use of differential SAR interferometry in monitoring and modelling large slope instability at Maratea (Basilicata, Italy) [J].Engineering Geology, 2003,68(1-2):31-51.

[7] 王桂杰，谢谟文，邱 骋，等.D-InSAR技术在大范围滑坡监测中的应用[J].岩土力学，2010，31(4)：1337-1344.

[8] 魏海霞，高照忠，叶长斌.融合D-InSAR 技术与ArcGIS 软件的矿区开采沉陷监测[J].金属矿山，2016(11)：128-131.

[9] 汪 磊，邓喀中，薛继群，等.融合概率积分模型与D-InSAR 的开采沉陷预计[J].金属矿山，2016(2)：160-163.

[10] 杨红磊，彭军还，崔洪曜.GB-InSAR监测大型露天矿边坡形变[J].地球物理学进展，2012，27(4)：1804-1811.

[11] 李如仁，杨 震，余 博.GB-InSAR集成GIS的露天煤矿边坡变形监测[J].测绘通报，2017(5)：26-30.

[12] 王志勇，孟庆颖. InSAR相位解缠算法在矿区沉降监测中的对比与分析[J].遥感信息，2014，29(6)：47-51.

[13] 赵东寅，申其鸿，马海涛，等.国产地基合成孔径雷达监测预警系统在紫金山金铜矿露天采场边坡位移监测的应用[J].中国安全生产科学技术，2015(4)：54-58.