基于InSAR 的煤矿采空区地表形变监测与预警

王凤云，陶秋香，陈 洋，韩 宇，郭在洁

（山东科技大学 测绘与空间信息学院，山东 青岛 266590）

近年来，能源需求不断提升，煤炭开采规模逐渐扩大，经过大幅度煤层开采后，地下会形成采空区，容易引发采空区的地表形变，即地表下沉、裂缝和地表抬升等现象，损坏农用耕地、建筑群等基础性设施；矿震是一种非天然地震，是由于地下持续开矿挖井，形成大面积空洞，导致采空区上方性质不同的层状岩层有规律性地弯曲、离层、断裂下沉，从而引发岩层内部储存的弹性能释放的现象，容易引起大面积地表塌陷，危害人民的生命财产安全[1-3]。2011 年11 月3 日，河南义马煤业千秋煤矿因矿震引发重大冲击压事故；2016 年4 月19 日，北京房山发生2.7级矿震；2019 年6 月9 日，吉林长春龙家堡矿业发生2.3 级矿震，导致井下当班作业人员被困；2020 年12 月15 日，陕西榆林市榆阳区发生2.6 级矿震，导致煤矿停产。因此，在提高矿区生产效率的同时，保证煤矿安全开采，迫切需要对矿区地表形变进行有效监测和稳定性分析，建立矿震预警机制，预防潜在矿震等灾害的影响[4-5]。

合成孔径雷达干涉测量InSAR（Interfero -metry Synthetic Aperture Radar）技术因其拥有全天候、全天时、高精度、高空间分辨率及覆盖面广等优点逐步用于矿区的地表形变监测[6]。随着该技术的不断完善，繁衍出DInSAR（Differential InSAR）、SBAS In－SAR（Small BAseline Subset InSAR）等多种地表形变监测手段。许多研究人员将InSAR 技术应用于采空区地表形变监测，姚佳明等[7]选用升、降轨L 波段PALSAR-2 影像数据，利用InSAR 技术对煤矿采空区开展了短期动态地表沉降监测，结合研究区开采信息对煤矿采空范围及开采时间进行反演，验证了InSAR技术对煤矿采区反演的合理性与可靠性；何荣兴等[8]介绍了采空区灾害类型及不同类型的相应案例，分析了采空区灾害发生的一般规律和特征，提出了尽量选择不产生采空区的采矿方法；栗明明等[9]利用39景Sentinel-1 降轨数据，采用SBAS InSAR技术获取隧道周边采空区地表形变发展过程，证实SBAS In－SAR 技术能够获得采空区的毫米级沉降结果。

为深入探究InSAR 技术在采空区地表形变的监测能力及对矿震的预警能力，验证该技术的监测精度，以近期发生矿震的山东某煤矿为研究区，依据矿震发生时间，选取2020-09-01—2020-12-31 覆盖该煤矿的11 景C 波段Sentinel-1A SAR 影像，分别采用DInSAR、SBAS InSAR 技术获取该时段内的矿区形变监测结果，并结合实际进一步分析时间演变生成的累计形变。

1 研究区与数据源

1.1 研究区概况

山东济宁市境内的鲁西煤炭生产基地是全国煤炭生产基地之一，矿产资源丰富，经过长期持续高强度开采，导致大面积土地沉陷，矿震多发频发，生态环境遭到严重破坏。因此，为保护耕地及矿区生态环境，减少由于煤矿采空区岩层的移动变形导致的地表塌陷以及地表建筑损伤倒塌现象，避免或减轻矿震灾害，迫切需要有效监测矿区采空区地表形变情况，查清采空区沉陷现状，加大对采煤塌陷地区的监测监管力度[10]。

2020 年12 月23 日，山东省济宁市曲阜（35.54°N，116.92°E）发生M2.4 级矿震，震中位于已停采的济宁某煤矿采空区。此次矿震，无开采工作面，井下安全，地表无塌陷，无地表建筑物和人员财产损失，各项情况正常。该煤矿位于兖州市以东约15 km，曲阜市西南约10 km，陵城镇附近。

1.2 SAR 卫星数据

Sentinel-1 卫星，载有C 波段合成孔径雷达，具备多种工作模式，是欧洲航天局哥白尼计划（Global Monitoring for Environment and Security，GMES）中的地球观测卫星。该卫星实现单、双极化等若干种不同的极化方式，能提供连续、全天候的雷达影像，拥有高重访频率、高覆盖能力以及极好的时效性和可靠性[11]。选取了覆盖矿区的11 景C 波段、VV 极化的Sentinel-1A 升轨数据，时间跨度为2020-09-01—2020-12-30。SAR 影像的具体参数见表1。

表1 SAR 影像的具体参数Table 1 Specific parameters of SAR images

为满足本次研究，去除地形相位，选用了由美国太空总署（National Aeronautics and Space Administration, NASA）和国防部国家测绘局（National Imagery and Mapping Agency, NIMA）联合测量的地面分辨率为90 m、平均精度16 m 的SRTM3-DEM[12-13]。

2 InSAR 技术基本原理

2.1 DInSAR 原理

DInSAR 技术是通过对同一地区不同时间的2 幅SAR 影像进行差分干涉处理获取干涉相位，利用多图像重复干涉或者引入外部DEM 模拟地形信息去除地形相位，从而获取地表微量形变的测量技术[14-15]。DInSAR 获取的相位可以表示为[13]：

对相位φdef进行相位解缠处理，得到真实形变相位φreal，然后提取雷达视线向上的地表形变量△Rtow，表示为[16]：

2.2 SBAS InSAR 原理

假设雷达传感器在同一研究区不同时刻获取n+1 幅SAR 影像，通过给定时间和空间阈值，生成m 幅差分干涉图。用ta、tb（ta>tb）时刻获取的SAR 影像生成去除地形相位的第k 幅差分干涉图，则第k幅影像在方位-距离像元坐标系（x，r）中的差分干涉相位δφk（x，r）为：

式中：φ（tb，x，r）、φ（ta，x，r）为tb、ta时刻相位；d（ta，x，r）、d（tb，x，r）为视线向累计形变量。

去除地形残余相位、大气延迟相位以及各种噪声相位后，地表形变的平均地表形变速率νT为：

则相位为：

式中：Ej、Sj为主、从影像获取时间；νk为k 时刻像元形变速率。

由此，定义A（j，k）=tk-tk-1，且A 是1 个m×n 的秩亏矩阵，得到矩阵方程：

通过奇异值分解法和最小二乘法可求出平均形变速率相位值，得到地表累计线性形变量[17-18]。

2.3 数据处理主要步骤

数据处理流程图如图1。

图1 数据处理流程图Fig.1 The primary data processing procedure

1）DInSAR 。DInSAR 基于干涉相位获取地表形变信息，关键技术主要包括：主、辅影像预处理、影像配准及重采样、干涉图生成、基线估计、地形相位去除、差分干涉图滤波、相位解缠、地理编码等，通过相位-形变转换，得到研究区的地表形变信息。

2）SBAS InSAR。选择2020 年9 月13 日的影像为主影像，其余为辅影像，建立连接图，设置合适的时间基线和空间基线阈值，将满足时空基线阈值条件的2 幅影像进行差分干涉处理，得到时序差分干涉图，影像滤波后选取高相干像元，并进行相位解缠，选择无残余地形条纹且远离形变区的地面控制点进行相位修正去除相位偏移，估算形变速率和残余地形，利用二次解缠优化干涉图，进行大气滤波估算，去除大气相位，得到时间序列上最终位移结果。

3）对比DInSAR 和SBAS InSAR 得到的各成像时刻的形变量，研究2 种InSAR 技术对采空区地表形变的监测和矿震预警能力。

3 实验结果分析

3.1 DInSAR 形变监测结果

结合Sentinel-1A 影像的获取时间，以12 d 为1个监测时段，利用该技术获取各时段矿区地表形变信息及累计形变信息。

3.1.1 特征点DInSAR影像叠加图分析

根据DInSAR 监测到的地表形变分布，选择形变特征明显并且能够监测到形变数据的9 个特征点进行数据提取与研究分析。9 个特征点分布与DInSAR各时段地表形变影像叠加图如图2，9 个特征点分布与DInSAR 累计地表形变影像叠加图如图3。

由图2、图3 可以看出：

图2 9 个特征点分布与DInSAR 各时段地表形变影像叠加图Fig.2 Overlaying charts of 9 characteristic points and DInSAR subsidence images at each period

图3 9 个特征点分布与DInSAR 累计地表形变影像叠加图Fig.3 Overlaying charts of 9 characteristic points and cumulative subsidence images by DInSAR at each period

1）DInSAR 监测到的地表形变中心（35.533°N，116.918°E）与官方发布的矿震中心“D”（35.54°N，116.92°E）相距约417 m。官方给出的矿震中心位置仅保留了小数点后2 位，是1 个大致位置，且“D”处于DInSAR 监测的地表形变范围之内，2020-09—2020-12 期间，“D”处地表经历了“抬升-下沉-抬升”。因此，DInSAR 监测到的地表形变中心在官方发布的矿震中心位置的误差范围内。

2）在第1（2020－09－01—2020－09－13）、第4（2020－10－07—2020－10－19）、第6（2020－10－31—2020－11－12）3 个时段内矿区上覆地表存在大面积抬升，在第2（2020－09－13—2020－09－25）、第3（2020－09－25—2020－10－07）、第5（2020－10－19—2020－10－31）3 个时段缓慢下沉。在第7 监测时段（2020－11－12—2020-11-24）内，地表出现-8～-10 mm 之间的不规则的大面积下沉，监测区域内最大形变量达到-11.5 mm。在第8 监测时段（2020-11-24—2020-12-06），矿区中心下沉，周边部分区域发生抬升，第9 监测时段（2020-12-06— 2020-12-18）矿区中心仍持续下沉，周边区域大面积抬升，地表高度差逐渐增大，地表结构变形，导致在第10 个监测时段（2020-12-18—2020-12-30）内发生矿震。

3.1.2 特征点DInSAR 形变数据分析

9 个DInSAR 特征点形变数据见表2、表3，9 个DInSAR 特征点形变数据图如图4。

表2 9 个DInSAR 特征点各时段形变量Table 2 Settlement data at each period of 9 characteristic points by DInSAR

表3 9 个DInSAR 特征点累计形变量Table 3 Cumulative settlement data of 9 characteristic points by DInSAR

由表2、表3、图4 分析可得：

1）该采空区工作面已停止开采，且附近无正在开采的工作面，因此监测初期（2020-09-01—2020 -11-12）研究区域的上覆地表形变趋于稳定。至矿震发生前夕（2020 年11 月下旬到2020 年12 月中下旬），监测点中最大累计形变量达-27.1 mm。此次矿震发生时，周边环境感受到异常，但未造成地表塌陷，因此监测到形变值较小。

2）图4 中，9 个监测点在前6 个监测时段（2020-09-01—2020-11-12）内，地表缓慢下沉和缓慢抬升交替性出现，且抬升范围在0.5～7.5 mm 之间，下沉数值小于-7 mm；在第7 个监测时段（2020-11-12—2020-11-24）时，9 个特征点均监测到形变的急剧变化，无抬升，最小下沉数据为-7.8 mm，最大下沉数据达-11.2 mm，9 个特征点在该时段的下沉数值均高于前6 个时段的下沉数值；从第7 个监测时间段到第10 个监测时间段（2020-11-12—2020-12-30），采空区持续下沉。综合来看，2020-11-12—2020-12-30，采空区经过形变量急剧变化后依然持续下沉，经过长时间的下沉累计，采空区内部承受不住地表重力，于2020 年12 月23 日发生矿震。

图4 9 个DInSAR 特征点形变数据图Fig.4 Settlement data diagrams of 9 characteristic points by DInSAR

3.2 SBAS InSAR 形变监测结果

9 个特征点分布与SBAS InSAR 累计形变影像叠加图如图5。9 个SBAS InSAR 特征点累计形变量见表4，9 个SBAS InSAR 特征点形变速率见表5，9个SBAS InSAR 特征点形变速率图如图6。

图5 9 个特征点分布与SBAS InSAR 累计形变影像叠加图Fig.5 Overlaying charts of 9 characteristic points and the cumulative subsidence images by SBAS InSAR at each period

图6 9 个SBAS InSAR 特征点形变速率图Fig.6 Settlement rate diagrams of 9 characteristic points by SBAS InSAR

表4 9 个SBAS InSAR 特征点累计形变量Table 4 Cumulative settlement data of 9 characteristic points by SBAS InSAR

表5 9 个SBAS InSAR 特征点形变速率Table 5 Settlement rate of 9 characteristic points by SBAS InSAR

由图5、表4、表5、图6 分析可得，

1）根据WGS84 坐标系的参数进行计算，SBAS InSAR 监测到发生地表形变的中心位置（35.532°N，116.918°E）与官方发布的矿震中心“D”的距离在460 m 左右。同样，“D”处于SBA InSAR 监测的地表形变范围之内，在2020 年9 月至2020 年12 期间，“D”附近地表经历了“抬升-下沉-抬升”。因此，SBAS InSAR 监测到的地表形变中心在官方发布的矿震中心位置的误差范围内。

2）SBAS InSAR 监测结果显示采空区工作面早已停止开采，上覆地表仍在持续下沉，地表形变活动一直在发生。矿震发生前夕（2020 年11 月下旬到2020 年12 月中下旬），SBAS InSAR 监测到矿区最大累计形变达到-41.3 mm。因此，为防止由于持续形变影响地表结构而产生矿震等强地表活动，防止危害附近人民安全和破坏生态环境，对采空区的监测是必要的。

3）图6 中，9 个监测点在前7 个监测时段（2020-09-01—2020-11-24）内，形变速率在0～-0.4 mm/d之间平缓浮动；到第8 个监测时段（2020-11-24—2020-12-06），监测点的地表形变速率增大，其中3、5、6 3 个监测点的形变速率有较明显的增大，形变速率在数值上分别增大了0.14、0.13、0.21 mm/d；从第8 个监测时间段到第10 个监测时间段（2020-11-24—2020-12-30），监测点的形变速率持续增大，形变量也持续增加。可见，该矿区采空区经过长时间的形变累计，内部承受不住地表重力而发生矿震。

为更进一步分析该采空区形变变化，对该区域进行剖面分析，考察其在时间序列上的形变变化。提取时序累计形变量绘制的采空区地表形变剖面图如图7。

图7 采空区地表形变剖面图Fig.7 Settlement profile of surface

图7 清楚地反映出该采空区在时间域的形变量变化情况。由图7 可见，随着时间的推移，该采空区形变量在逐渐增加，2020 年11 月24 日，最大累计形变量达到-30.4 mm，自此之后，形变量开始增大。由此可以推断，该采空区由于长期持续下沉，导致2020 年11 月24 日后，形变速率加快，地表结构加速破坏，导致矿震的发生。

4 结 语

1）受矿震影响前，DInSAR、SBAS InSAR 技术均监测到矿区上覆地表缓慢变化，监测到的地表形变分别在-0.2 ～-6.9、-1.0 ～-5.2 mm 之间；矿震发生前夕（2020 年11 月下旬到2020 年12 月中下旬），2种InSAR 技术都监测到采空区上覆地表形变发生了较明显的变化，形变持续加快，形变量持续增大，地表持续下沉，DInSAR 监测到最大累计形变量达到-27.0 mm，SBAS InSAR 监测到矿区最大累计形变达到-41.3 mm；矿震发生期间，矿区地表仍持续下沉。

2）DInSAR 技术采用相邻成像时刻的2 幅SAR影像两两差分干涉处理获取相邻时刻之间的地表形变信息，处理耗时，低相干点的监测精度不高；SBAS InSAR 技术采用的是时间序列的差分干涉处理得到各成像时刻的累计地表形变信息，处理过程复杂，无法得到低相干点的地表形变信息；二者各有优缺点。

3）2 种InSAR 技术的监测结果均体现出煤矿采空区经过持续形变，会引发地表形变加剧，地表结构遭到破坏，内部承受不住地表重力，导致矿震发生。研究结果对DInSAR、SBAS InSAR 技术应用于矿震预测提供一定的参考，但尚需结合其它更多的矿震实例做进一步的深入研究。