基于DInSAR技术的矿区累计时序形变监测研究*

王 瑞，薛玉玲，万士桢，郑玉琳 ，张 婷

(1.赣南科技学院，江西 赣州 341000；2.中国矿业大学 环境与测绘工程学院，江苏 徐州 221116;3.江西理工大学 应用科学学院，江西 赣州 341000)

0 引言

我国矿产资源比较丰富，矿产资源开采在带来可观的经济效益的同时，也给生态环境以及人类的生产生活带来了一定危害。采空区上部岩层在重力作用下会发生一系列变形，造成矿区地表下沉，甚至引发一系列环境地质问题。

随着卫星遥感技术的发展，合成孔径雷达干涉测量(InSAR)技术逐渐成为我国形变监测的主要手段。该技术是一种全新的空间大地测量技术，其具有全天观测、实时观测且观测范围广的优点，能够获取高精度、高分辨率、大范围的观测信息，且无需在地面布置观测点，因此该技术被广泛应用于地表及各类地物形变监测。在监测地震同震变形场时，差分合成孔径雷达干涉测量(DInSAR)技术能够获得雷达视线方向上的形变量[1-2]。

地下煤炭开采常会引起地表沉陷，从小规模小尺度的沉陷逐渐发展到区域性大尺度沉陷，沉陷量及范围主要取决于地质条件、岩石力学性质以及开采技术方法[3]。地表沉陷并不是在短时间内发生的，而是随着开采面的增大而逐渐显现的，因此持续监测是必要的。地下煤炭开采引起的大尺度形变是一个复杂的时空过程，通过地表监测对矿区地质灾害的防治具有重要意义。针对采煤引起的地表沉陷，推广程度最高、应用范围最广的地表沉陷信息获取技术为常规观测站法，即在拟开采工作面上方沿工作面走向和倾向每隔 5～25 m布置一系列观测点，在工作面开采的不同阶段采用高精度水准仪、全站仪或GPS进行实地监测，经数据处理获得沿煤层走向和倾向的下沉量及变形情况。该方法是逐点监测，空间分辨率低、效率不高；在长期监测过程中，观测站易受矿区地表复杂条件影响。

随着测量技术的发展，InSAR技术为解决区域大范围沉陷监测问题提供了一种有效方法，其可以提取存档数据分析沉陷前后的规律[4-7]。本文以西部某矿区的工作面为研究对象，利用Radarsat-2图像作为实验数据，通过图像解译获取沉陷盆地边界。

1 InSAR技术原理及DEM修正

1.1 InSAR技术

InSAR技术是一种图像处理技术，可以从SAR传感器的两个天线干涉通道生成数字高程模型(DEM)和相对相干性。随着InSAR技术的发展，其应用领域不断扩大，如通过同一区域两颗卫星通道测量卫星和地球表面之间的相位差，此相位差的产生原因主要有：①两次或多次卫星轨迹之间的位置差异，但是差异远小于卫星到地球表面的距离，因此可以应用于地形地貌形变监测；②两次或多次卫星采集观测到的区域位移，此位移可能由地下资源的开采、地震及火山运动等因素造成。因此，当SAR系统对同一区域进行两次或多次观测时，如该区域的几何位置相对于传感器发生了变化，则视为地表产生了形变。

SAR系统在采集数据过程中，反射信号会受噪声、地形、大气及地表运动等因素的影响。DInSAR原理及地表形变几何关系示意图见图 1。

图1 DInSAR原理及地表形变几何关系图

地表形变干涉相位的计算步骤为：干涉相位φ表示同一个目标P由卫星S1和S2在两次不同采集的时间内，P点发生了形变而移至P′点。干涉相位φ受多个因素的影响，其计算式如下[8]：

(1)

(2)

φInt=φTopo+φFlat+φDefo+φAtm+φNoise，

(3)

(4)

h=H-R1cosθ，

(5)

h′=H-R2cosθ，

(6)

式中：λ为雷达卫星的波长，R1、R2分别为卫星两次通过目标P、P′时的距离，φS1、φS2分别为主卫星和副卫星位置的干涉相位，φTopo为受地形影响的相位，φFlat为参考面因素引起的相位，φDefo为沿视线向(LOS)地表形变因素引起的相位，φAtm为大气延迟因素引起的相位，φNoise为噪声因素引起的相位。式(4)、式(5)和式(6)揭示了干涉相位φInt与高程h、h′之间的函数关系。已知天线位置的高H、基线长度B、基线水平角α、往返双程相位差N，可得到雷达成像视角θ，再计算形变前后点P的高程h和h′，最后获得形变量。

采用GAMMA的传统连续干涉测量方法解算有效的影像。ZEBKER等[9]的研究成果表明，为能够准确有效解算出雷达视线的形变量，单位像素所能解算出的形变量为波长的一半。RadarSat-2的波长为0.056 m，故干涉图所能解算出的最大单位像素理论变形值为0.028 m。当任意两幅影像之间相邻单位像素的形变大于波长的二分之一时，则可能无法在时域内正确地解缠。

1.2 DEM修正

何秀风等[10]的研究结果表明，在采用DInSAR技术进行沉陷解算时，因垂直基线不会为0，所以DEM精度对两种通法的解算具有关键性作用。因此本文在12.5m DEM的基础上，融合无人机(Unmanned Aerial Vehicle,UAV)获取的开采前地表DEM数据，分别对由原始DEM和UAV获取的DEM数据进行标准差计算，再据此计算两种融合DEM的权重[11-12]。融合模型如式(7)和式(8)所示，融合结果如图2所示。

图2 融合DEM

(7)

hr=Jhy1+(1-J)hy2，

(8)

式中：J为数据融合权重，取值范围在0～1；s1和s2为融合DEM标准差；l为融合影像误差的相关系数；hr为融合后的DEM数据；hy1和hy2为融合前两种原始DEM数据。

2 研究区域概况及数据

2.1 研究区域概况

本文以西部某矿区为研究对象，其位于内蒙古自治区鄂尔多斯市。矿区为不规则多边形，南北长14.40 km，东西最宽6.77 km，标高1 150～1 420 m，地表起伏较大。矿区煤田广阔、资源丰富、煤质优良；煤层埋深较小，约为150～456 m。

2.2 数据采集

利用SAR影像作为实验解算数据，雷达影像采用加拿大太空总署发射的RadarSat-2(搭载C波段传感器，重访周期24 d)采集差分干涉测量的影像数据。干涉影像的采集日期及相关参数如表1所示。

表1 SAR影像参数

3 解算结果及精度分析

3.1 解算结果分析

该矿区某工作面自2018年7月12日开始回采，至2018年10月25日回采完毕，共耗时105 d，工作面总长度约1 263.7 m。工作面地表环境复杂，地势起伏较大，地表土壤属于沙地，表面有许多被雨水冲刷而成的沟壑，传统测量方法难以监测地表变形。因此，本文采用DInSAR技术构建完整沉陷盆地模型，用以分析沉陷规律。

图3为相邻两期图像之间的时间序列形变图。由图3可知，受2S201工作面东侧的2S202工作面开采的残余形变影响，2S202工作面自2018年6月9日至2019年1月11日地表出现形变，形变量随着时间的推移不断减小，此后逐渐消失；2S201工作面于2018年7月12日开始回采，随着工作面的推进，地表开始出现形变，并且形变沿着工作面回采的方向移动；2018年7月12日至2018年10月25日为工作面强回采期，地表形变SAR影像失相干较为严重，仅能识别出较小形变；随着回采的结束，地表形变逐渐减小，DInSAR技术可识别出沉陷盆地边缘，最大识别沉陷量出现在2018年11月24日至2019年1月11日期间，为0.144 m。

图4为2018年7月12日至2019年4月16日期间的时间序列累计形变图。由图4可知，在此时期内视线方向所能识别的累计最大沉陷量为0.342 m，随着工作面的推进和时间的推移，受回采的影响沉陷盆地逐渐增大。图4a-图4c为回采时期的形变图，该时期内地表沉陷量较大，SAR影像失相干严重，识别沉陷量较小。

图3 相邻两期图像之间的时间序列形变图图4 DInSAR累计沉陷监测

3.2 精度分析

在2S201工作面上部布置监测点进行精度验证。采用三角高程中间观测方法进行数据采集，根据坐标位置提取DInSAR时序累计沉陷量。从监测点中选取具有代表性的10个监测点进行误差统计，结果如表2所示。

表2 DInSAR误差统计 单位：m

DInSAR技术无法真实解算大量级形变，不适合大沉陷监测，沉陷越小的区域误差越小。由表2可知：对于小尺度形变，DInSAR技术的精度可在厘米级；当沉陷量增大时，误差为分米级和米级。2S201工作面自南向北回采，地表由南向北沉陷量递减，因此在累计沉陷量监测中北部大于南部。DInSAR技术对于沉陷盆地周边小区域的沉陷识别较为敏感。

4 结论

本文采用DInSAR技术提取因井下煤炭开采引起的地表变形数据，在数据解算过程中引入多源DEM数据，计算每种DEM数据标准差并建立不同权值融合模型。提取西部某矿区某工作面2018年7月12日至2019年4月16日的完整地表形变场进行分析，结果表明，在地表条件特殊且地面观测站较少的情况下，DInSAR技术可以有效分析地表形变规律，其虽受相干影响而无法识别大尺度真实形变，但可以识别沉陷边界。