矿区地表沉陷D-InSAR监测分析与应用

孙建，王道顺，姚文举,3*，李进峰，甄洪帅,3，张方龙

(1.山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队)，山东 兖州 272100；2.高分辨率对地观测系统济宁数据与应用中心，山东 兖州 272100；3.自然资源部采煤沉陷区综合治理工程技术创新中心，山东 兖州 272100；4.梁山县自然资源和规划局，山东 梁山 272100)

0 引言

地下煤炭资源采出后，会导致开采区周围岩体原始应力平衡发生破坏，造成应力重新分布和覆岩移动，引发地表沉陷，从而进一步破坏耕地、建(构)筑物和其他基础设施，严重制约矿区和城市的可持续发展，而对开采沉陷区进行实时动态监测则是研究矿区开采沉陷规律、预防灾害发生和解决矿区环境问题的重要手段，因此开展矿区地表沉陷监测具有重要意义[1-2]。

传统的矿区地表沉陷监测方法主要通过布设岩移观测站，利用GNSS、全站仪、水准仪等常规测量手段定期进行沉陷观测，此类方法监测精度虽然可观，但监测点位少、布设的监测点容易遭到破环，难以整体反映矿区开采的地表沉陷特征，此外如果在观测条件恶劣的地区还容易对测量人员构成安全威胁[3]。合成孔径雷达差分干涉测量技术(Differential Interferometric Synthetic Aperture Radar, D-InSAR)作为一种新兴的空间对地观测技术，通过非接触测量的方式克服了传统监测方法的不足，可全天候连续获取大范围区域的微小形变信息。2004年，吴立新[4-5]在国内首次分析了D-InSAR技术在煤矿开采沉陷变形监测中的优势，并将该技术成功应用于唐山市及开滦矿区。2008年，刘广等人[6]利用ENVISAT和JERS-1雷达数据对河北峰峰矿区进行了地表沉陷监测，分析了不同波段条件下干涉相位特性。2014年，王志勇等人[7]以济宁矿区为研究对象，采用ALOS PALSAR雷达数据对地下采煤引起的地表沉降进行了精细化分析与解译，证明D-InSAR在矿区地表开采沉陷监测中的精度已达到cm级。2017年，郭山川等人[8]利用两轨法D-InSAR技术对黄土高原的大同市南郊区进行了采煤沉陷监测，并与开采沉陷预计结果进行了验证，分析了结果差异原因，此后众多学者不断将D-InSAR技术应用于矿区地表开采沉陷监测研究当中[9-14]。

本文在分析双轨D-InSAR获取地表形变原理的基础上，对两景高分辨率RADARSAT-2超宽精细(Extra Fine,XF)模式雷达影像进行差分干涉处理，获取了霄云矿区2016年12月18日—2017年2月25日开采沉陷区域和沉陷值，对监测结果和精度进行了分析，证明了D-InSAR技术在矿区地表开采沉陷监测中的有效性。

1 D-InSAR技术原理

D-InSAR技术是通过对同一地区前后跨越形变期获取的SAR单视复数影像干涉得来的相位差，进而获取地面形变信息的技术。根据在差分干涉处理过程中去除地形相位方法的不同，可以将D-InSAR测量形变分为二轨法、三轨法、四轨法等[15]。二轨法所需数据量少，干涉处理简单，依靠外部DEM数据模拟去除地形相位，由此获取地表形变信息。随着高精度覆盖全球陆地地表数字高程模型(Digital Elevation Model, DEM)数据的免费使用，二轨法已经成为D-InSAR获取地表形变的常用方法。

图1为D-InSAR技术的基本原理图。在雷达两次成像前后，地物点P发生形变后的位置为P'，形变位移矢量为△d。在不考虑其他相位误差的条件下，前后两景影像的干涉相位可表示为：

图1 D-InSAR基本原理图

由于干涉条件的限制，D-InSAR的干涉相位φint中通常包含不同的相位成分，其中φref为参考相位，φtopo为地形相位，将这两项从干涉相位去除后即可得到形变相位φdef，因此可以得到：

(2)

由式(2)可以看出，D-InSAR技术监测得到的形变量是在雷达视距方向的投影，数据处理流程见图2。

图2 D-InSAR数据处理流程图

2 研究区与数据源

2.1 研究区概况

霄云煤矿位于山东省济宁市金乡县霄云镇境内，隶属于济宁矿业集团有限公司。矿井于2012年8月建成投产，可采储量4791.5万t，主采3煤赋存稳定，可采范围内平均厚度3.49m，倾角一般为10°～21°，采用走向长壁采煤法，综合机械化一次采全高采煤工艺。截至2017年3月，研究区共包括7个工作面，工作面开采情况见表1，区域概况见图3。

表1 研究区工作面信息

图3 研究区域概况图

2.2 数据源

矿山开采沉陷具有大量级、非线性等特点，当地表形变梯度较大，超过InSAR技术所能测量的梯度范围时，容易造成干涉图干涉条纹杂乱和相位混叠，导致解缠错误，从而难以获取准确的开采沉陷信息[16]。InSAR监测的最大形变梯度可表示为[17]：

dmax=λ/μ

(3)

式中：λ—雷达波长，μ—影像像元分辨率。可以看出，InSAR可监测的最大形变梯度与所用雷达数据参数相关，雷达波长越长、像元分辨率越高，能监测到的形变梯度越大，越适用于矿区开采沉陷监测，因此本次实验选择加拿大RADARSAT-2 XF的C波段5m高分辨率SAR影像作为监测数据源，该影像在矿区开采沉陷监测中具有显著优势[18]。外部DEM数据选择使用日本太空发展署(JAXA)免费分发的30m分辨率ALOS World 3D数据用以模拟并去除地形相位。所用影像对参数见表2。

表2 干涉SAR影像参数列表

3 数据处理与分析

3.1 数据处理过程

双轨D-InSAR方法可以快速获取监测区域的形变信息，但需要保证监测区域具有良好的时间和空间相干性，若所用数据的空间基线超过临界基线值，会导致完全失相干，无法获取有效监测信息，因此在数据处理前一般需要对所用数据进行基线估计，以便判断所用数据能否达到理想的监测结果。由表2可知，两景影像时间基线为72d，空间基线仅为-23.96m，远远小于其空间临界基线值，满足影像干涉要求。

实验选择2016年12月18日影像作为主影像，2017年2月25日获取的影像为辅影像，根据图2的数据处理流程进行差分干涉处理。首先对两景影像进行配准处理，并保证配准精度小于1/8个像元，生成原始干涉图(图4a)；在距离向和方位向按照1∶2的视数比进行多视处理，去除平地相位和地形相位后生成差分干涉图(图4b)；由于数据成像和处理过程中不可避免会受到噪声影响，实验选择自适应的Goldstein滤波方法对差分干涉图进行滤波处理(图4c)，减弱噪声相位，并生成相干系数图(图4d)；为了恢复真实相位值，利用最小费用流方法进行对干涉相位进行解缠处理，设置相干解缠阈值0.2，逐像元进行相位解缠。

图4 差分干涉处理结果图

在远离形变区域、无相位解缠错误且相干性良好的区域选择19个控制点，利用三次多项式模型去除残余的恒定相位和相位坡道，轨道精炼后的相位中误差为-0.0057rad，标准差为0.1614rad；去除残余相位后将相位换算为形变，经过地理编码后生成监测区域的沉陷结果(图5)。

图5 研究区沉陷量图

3.2 监测结果分析

从图5可以看出，由于煤炭开采活动影响，研究区域在2016年12月18日—2017年2月25日内发生明显沉陷现象，共分布有3处沉陷盆地，沉陷盆地空间分布特征与采空区分布实际情况吻合，证明了D-InSAR技术在矿区地表沉陷监测中的有效性。

为了进一步对监测结果进行分析，实验设置沉陷等值线初始值-0.005m，间距0.01m，绘制沉陷等值线见图6。

图6 研究区沉降等值线图

可以看出，西侧沉陷盆地沉陷范围最大，中间位置的沉陷盆地沉陷量最大，东侧沉陷盆地沉陷量最小。结合图3采空区情况分析，a区域监测时间段内1311工作面自西向东开采，北侧1309工作面为已停采区，1311工作面的采动引起了1309已采区上覆岩层“活化”而产生二次破裂，造成沉陷盆地中心位于两工作面中间位置且东部沉陷量明显大于西部，最大沉陷量为0.118m，而1309工作面自西向东开采停采5个月，处于沉陷衰退阶段，东侧区域沉陷量小于0.015m；b区域两个沉陷盆地已接连成片，1304工作面的开采同样引起了1302-1、1302-2和1306已采区的二次变形并使沉陷量和沉陷范围扩大，最大沉陷量达到0.174m，1312工作面自东向西开采，仅停采2个月，最大沉陷值-0.064m，因此仍能呈现出沉陷盆地形态。

值得注意的是，根据研究区的煤层开采厚度和其它地质采矿信息，本次研究监测到的工作面开采下沉量显著小于理论值。分析其原因认为，本次实验监测的时间区间为2016年12月18日至2017年2月25日(共69天)，未完整覆盖1311和1304工作面的全部开采时段，监测时段之后地表仍处于移动变形活跃期，因此监测实验得到的整体下沉值偏小。由于缺少实测数据，本次实验采用内符合评估的方式对监测结果精度进行评价，通过统计形变量及其中误差的大小和分布来验证监测结果的可靠性。

从表3、图7可以看出，监测结果的沉陷量峰值为-0.001m，沉陷分布接近于0m，符合大范围稳定的实际情况；沉陷精度峰值-0.026m，平均值0.011m，符合D-InSAR技术的cm级监测精度要求。

表3 精度统计表

图7 监测沉降量统计分布图

4 结论

本文利用D-InSAR技术对覆盖研究区域的两景RADARSAT-2雷达数据进行差分干涉处理，获取了该区域2016年12月18日至2017年2月25日时间段内的沉降信息，并对监测结果和精度进行了分析。

(1)D-InSAR监测到研究区域共分布有3处沉陷盆地，监测到的最大沉陷量达到-0.174m。

(2)研究区域的监测结果沉陷分布特征与工作面开采实际位置高度吻合，沉陷量大小和沉陷范围与工作面开采情况相符。

(3)监测区域沉陷量的统计分布峰值趋近于0m，内符合精度误差平均值0.011m，整体沉陷趋势符合实际情况，证明D-InSAR技术在矿区地表沉陷监测中具有良好的应用效果。