高分辨率时序InSAR在矿区开采沉陷监测中的应用——以霄云矿区为例

甄洪帅，张红月，姬宗皓，殷幼松

(1.自然资源部采煤沉陷区综合治理工程技术创新中心，山东 济宁 272000；2.山东省鲁南地质工程勘察院，山东 兖州 272100；3.济宁市采煤沉陷地治理中心，山东 济宁 272000；4.环球星云遥感科技有限公司，浙江 湖州 313200)

0 引言

霄云矿区位于济宁市金乡县霄云镇，于2008年4月开工建设，区内地势平坦，地面标高为35.8～44.9m，总体趋势为西高东低。土地肥沃，地面村庄较多。区内地层为走向近EW、倾向北的单斜构造，井田北至刘吴庄断层，西为曹马集断层，东距鱼台县城约20km，南起各煤层露头。该矿区自2011年底到2015年中由水准测量观测到的最大沉降量超过300mm。为了获取霄云矿区地表沉降的现状资料以服务于后期矿区治理方案的规划设计，研究该矿区采煤沉陷的空间分布及时序演化特征具有重要意义。然而，传统观测手段虽然具有精度高、监测量级大等特点，但在矿区面域形变监测应用中存在着经济代价高、观测周期长、采样点密度难以反映矿区整体形变特征等突出问题[1-7]。Lidar技术可以通过高程变化从而提取矿区形变特征，但费用昂贵，在大规模工程化应用过程中存在一定局限性[8-10]。利用合成孔径雷达干涉测量(InSAR)技术方法进行采煤沉陷区的沉降监测在德国、澳大利亚和波兰等国研究较早[11-13]，近年我国相关学者应用该技术也对矿区沉陷开采做了相应的研究及应用工作[14-17]。而InSAR技术监测矿区沉陷的应用效果取决于诸多因素，主要包括：SAR影像空间分辨率、影像数量及时间采样频率、雷达波长、矿区工作面地物特征及其开采特点等。就空间分辨率这一要素而言，中等分辨率SAR影像在解译过程中容易造成相位混叠、相干点目标密度较小等问题，对整个矿区形变场空间形态及量级程度的刻画存在不足，继而影响矿区沉陷范围线的提取精度[18-22]。鉴此，该文以加拿大RADARSAT-2超精细宽(ExtraFine，XF)模式5m高分辨率时序影像为数据源，利用时序InSAR技术对济宁市霄云矿区开展沉陷调查监测，为后期矿区综合治理及生态修复工作的顺利开展提供切实可靠的现状资料。

1 SAR影像选择

影像源的选择上主要根据当前主流SAR卫星存档数据可利用情况，并结合采煤沉陷区地面沉降监测对SAR数据源幅宽大、分辨率高、重访周期稳定及数据获取可靠等方面的参数需求，最终确定使用RADARSAT-2的超精细宽高分辨率SAR存档数据(ExtraFine，XF)，实现工作区全覆盖、高精度、分年度的监测，确保监测区因时间差异而形成的塌陷形变可以在InSAR结果中有很好的体现。同时，为了平衡影像成本和监测效果，选用每年度时间上分布较为均匀的5期影像，2016--2018年共30景XF时间序列数据进行观测，SAR数据集基本参数见表1。

表1 SAR数据集基本参数

2 时序InSAR分析技术

2.1 差分相位模型

时序InSAR分析是D-InSAR监测技术的应用拓展，通过对不同时刻获取的同一区域SAR影像上受失相干现象影响较小的像元(高相干点目标)开展相位稳定性分析，以实现长时间跨度地表沉降信息的高精度获取。在经过时间序列SAR影像同名点位的配准、从雷达坐标系向地理坐标系的编码转换、不同时刻影像多基线组合选取及利用外部DEM对干涉相位去除地形影响后，时序差分相位模型可表示为：

φx,i=φdef,x,i+Δφε,x,i+φα,x,i+Δφorb,x,i+φn,x,i

式中：φx,i代表第i幅差分干涉对上第x个像元的差分干涉相位；φdef,x,i为雷达视向上的形变相位；Δφε,x,i为外部参考DEM不准确引入的残余地形相位误差；φα,x,i代表大气延迟误差；Δφorb,x,i为卫星轨道数据不准确引入的相位误差；φn,x,i为其他噪声相位，如：热噪声、体散射变化误差及数据处理中引入的误差等。差分相位模型的构建是各种时序InSAR分析方法的解算基础。

2.2 监测方法

当前，根据研究对象的不同，从差分干涉相位中获取形变成分主要有2种处理思路：一种是假定形变区符合某种时序演化规律，先解算出主体形变成分，然后从残余相位中分离出剩余部分，两者之和便是形变信息；另一种则是通过从差分干涉相位中逐步剥离其他干扰相位成分进而获取总的形变信息。与城市缓慢地表沉降可采用前一种时间线性形变模型描述不同，由煤炭开采引起的形变发育过程较为复杂、沉陷量级大，且各工作面开采时段并不一致，故应该采用第二种处理思路，即：不假定时间形变模型的空间相关性InSAR分析方法[23-24]。其处理步骤包括：时序SAR影像配准、干涉图生成、地形相位去除、初选相干点目标、空间相关和非相关视角误差估算、终选相干点目标、相位解缠和大气相位获取等，形变分析过程如下。

假设φdef,φα.Δφorb在一定空间尺度上具有相关性，Δφε.φn在此距离上非空间相关，且在x以为中心、L为半径的圆内其均值为零。以第x像素半径L内所有高相干点目标的相位平均值可表示为：

式中：

由DEM误差造成的相位差与垂直基线成正比例：

φε，x,i=B⊥,x,i·Kε,x

式中：Kε,x为比例常量，改写为：

Kε,x是唯一与基线相关的因子，可通过最小二乘法估算。相位稳定性评定因子表述如下：

式中：φ＇ε,x,i为由Kε,x估计不准确而引起的残余误差，由于相邻相干点目标间空间相关DEM误差相差并不大，故该项对解缠误差的影响较小。对上式进行三维迭代相位解缠处理，并利用大气、轨道和形变信号三者之间时空频谱特性的差异，实现形变信号和噪声信号的有效分离。

3 矿区开采沉陷提取

3.1 干涉组合选取

对于时序InSAR分析而言，可通过两种组合模式构建干涉对，一种是单一主影像的放射状网络；另一种则为短时间、空间基线的小基线集组合。前一种模式一般在影像数量大于25期以上时采用；而后一种模式能够以更高的空间密度提取到相干点目标，对矿区空间形变特征的反应更加明显、直观。该文受限于数据量少且影像间隔时间较长，采用小基线集干涉组合方式，将空间垂直基线300m、时间基线200d作为阈值，所选组合情况如表2所示。

表2 研究区图幅干涉组合时空基线

3.2 差分干涉处理

RADARSAT-2XF模式影像单视情况下方位向分辨率为2.5m左右，而距离向则为5m。为抑制相干斑噪声且同时保持高分辨影像的精细化监测优势，采用距离向和方位向1∶2的窗口对配准后影像进行多视处理，在提高相位信噪比的同时，使得方位向和距离向分辨率接近一致，改善影像判读效果。此外，由于煤炭开采区多位于农田等植被覆盖度高、相干性低的区域，需要对去除地形贡献后的相位进行滤波。

3.3 研究区沉陷监测结果

在对影像解译处理过程中需要注意两个方面：一方面是大气条件和解算传播误差的影响。大气一般在空间1～3km呈现出较强的相关性，且随着距离的增大，相关性逐渐减弱；而解算传播误差同样和距离有关，空间范围越大，误差积累越明显。鉴此，为有效控制大区域时序InSAR分析过程中大气条件及解算传播误差的影响，以霄云矿区边界外扩1～3km为界，对获取的时序差分干涉相位进行解算分析。另外一方面则是需要排除地下水开采这类由居民生产活动产生的地面沉降现象。与矿采形变不同，地下水开采不会在短期内发生剧烈的地表形变，因此其相位条纹不会过于明显；此外该类现象多发生在人口聚集地，可通过形变场与村庄空间分布的叠加分析予以判别。综合上述考虑后，得到2016—2018年研究区范围内的沉陷演化过程，如图1和2所示。

图1 霄云矿区时序形变演化过程

图2 霄云矿区累计沉降量及其等值线

4 结果分析

时序InSAR技术获取的监测结果既可以从空间上分析矿区开采引起的沉陷对于周边环境的影响，尤其在“三下压煤”开采情况下地面沉陷对于建、构筑物的安全稳定构成的威胁，同时也可以与开采面的工作时段相结合，从而综合矿区开采工艺、地质条件等因素对沉陷发育规律进行研究，进一步加深矿区沉陷过程的认知，为开采沉陷控制方案的制定提供辅助资料。

4.1 研究区地质采矿条件

霄云煤矿采用立井开拓，井底水平为-790m，采煤方法为综采放顶煤采煤法，开采标高为-430m～-1500m[25]。其中，首采区地面影响范围内有共有11个村庄，分布较为零散；二采区范围内有9个村庄。该区为全隐蔽式华北型石炭二叠纪含煤地层，地层自上而下分别为第四系、古近系、二叠系、石炭系、奥陶系。区内总体为向北倾伏的单斜构造，地层倾角10°～21°，浅部陡深部缓；断裂构造较发育，共发现断层19条，其中落差大于100m的2条；50～100m的2条；20～50m的11条，断层走向多为NE向和EW向，且主要断层控制程度均已达到可靠或较可靠。矿区煤系地层平均总厚度260m，含煤19层。可采、局部可采煤层3层(3、12下、16上)，平均总厚度5.37m。3煤层厚2.25～7.69m，平均3.98m，稳定全区可采；12下煤层厚0.34～1.28m，平均0.71m，局部可采；16上煤层厚0.60～0.74m，平均0.68m，较稳定局部可采。

4.2 空间影响分析

霄云矿区井田存在较厚的冲积层，在进行工作面开采时，地表下沉量相比于基岩情况下要大，由于厚冲积层的移动传递，会使得浅层土体固结压缩时形变较为均匀，进而反映到地表时会相对平缓，不会引起较大量级的集中形变；此外，开采煤层埋藏较深，矿区内煤层顶板以砂岩为主且较为稳定，石炭-二叠系硬岩超过150m，条带开采条件下上覆岩层与硬岩岩梁将整体缓慢下沉，使得开采引起的地表下沉、变形值会比较均衡[26]，这些从监测结果中得到了很好的体现。霄云煤矿主采3煤层，采深约-790m，采用条带式或隔面开采，因此，开采沉陷并不存在非常集中的大量级变形。从图3中可以看出，在监测时段内，地面沉陷主要集中在2015—2018年度的开采工作面周边，而2012—2014年的开采面已经基本上达到了稳沉状态，不再有明显的下沉趋势。

4.3 开采沉陷规律分析

霄云矿区共划分为6个采区，2012—2018年的开采工作面位置分布如图3中蓝色框体。对矿区内不同开采时段工作面选择点位绘制沉陷形变曲线以探究地面沉陷发育规律，所选取的6处具有代表性点位空间位置如图3中A—F所示，形变序列演化曲线见图4。

图3 矿区开采面空间沉降影响范围

经分析，图4中的6处形变曲线主要体现出：①InSAR形变演化过程与开采时段引起的沉陷吻合度高。图4a位置沉陷由2015年工作面的开采所导致，从曲线发育趋势可以看出，存在明显的收敛趋势，而结合空间沉陷范围可以看出，大部分在2012—2014年间的开采面周边已经达到稳沉状态，不存在突出形变迹象；而在图4b中可以看到，在2017年工作面尚未处于开采状态时，地面点位周边较为稳定，而随着开采工作的进行，随之发生快速下沉过程，这些点位的沉降迹象都与矿区沉陷发育规律较为一致。②距离开采工作面的远近不同，地表沉陷曲线演化形态非常相似，这一点如图4b和图4c所示。两处点位虽然在时序沉降数值以及曲线的缓和程度上存在一定差异，但不难看出，在不受邻近开采工作面影响情况下，两者曲线形态保持一致；同时，距离开采工作面越近，曲线的变化程度越为剧烈，沉陷量值越大。③多个邻近开采工作面共同作用对于沉陷曲线的形态存在显著影响。为了便于对比分析不同工作面、不同开采时间对开采沉陷的影响，将D，E和F点绘制于一张图上，如图4d所示。从图4d中可以看出，D点2016年开采引起的形变仍然处于中后期发育过程中。虽然E和F点位都主要受2018年开采工作面影响，但可以明显看出这两处点位变形过程并非遵循先稳定而后快速沉陷的演化规律，这主要是因为这两处点位在2018年工作面进行开采之前，周边还存在2016年的开采面，由于2018年开采时该处尚未达到稳沉状态，因此出现了2016年中期到2017年初的缓慢沉陷过程，而到了2017年基本处于沉陷发育的末端。在2018年随着新的工作面进行开采，导致地面沉陷现象又重复发生，便出现了上述的形变发育混叠过程。

5 结论

该文利用15期RADARSAT-2XF模式高分辨率数据，以济宁市霄云矿区为例，采用基于空间相关性分析的时间序列InSAR方法对其采空区沉陷现状进行了解译获取，结果表明：

(1)高分辨率SAR影像与时序InSAR技术的结合可以实现高精度的矿区开采形变信息的准确获取，具有监测效率高，经济代价小、可延续性好等突出优势，在大范围采煤沉陷区沉陷监测方面极具应用潜力。

(2)高分辨率时序InSAR分析方法可以很好捕捉到开采沉陷区的时空形变特征，通过对时序相位成分的分析，在获取沉陷区累计沉降量的同时，还可以得到每景影像获取时刻的变形信息，为矿区地面沉陷演化过程的分析提供了资料。

(3)受当前雷达卫星运转条件以及传感器自身成像要求的限制，InSAR技术在矿区沉陷监测中的应用效果受到了一定程度的制约，但可以预见，随着SAR卫星及传感器的迭代更新及InSAR技术方法的不断发展，其在矿区沉陷监测方面将会发挥更加突出的作用。