基于SBAS-InSAR技术的煤矿开采沉陷监测与分析

doi:10.11873/j.issn.1004-0323.2022.4.1021

基于SBAS-InSAR技术的煤矿开采沉陷监测与分析

张香凝^,¹, 贺黎明^,¹, 刘翠芝¹, 王兴杰¹, 唐永亮², 何蓉花¹

1.东北大学资源与土木工程学院测绘工程系，辽宁沈阳 110819

2.辽宁省地质环境监测总站，辽宁沈阳 110000

Monitoring and Analysis of Coal Mining Subsidence based on SBAS-InSAR Method

Zhang Xiangning^,¹, He Liming^,¹, Liu Cuizhi¹, Wang Xingjie¹, Tang Yongliang², He Ronghua¹

1.School of Resources and Civil Engineering，Northeastern University，Shenyang 110819，China

2.Liaoning Geo -environmental Monitoring Station，Shenyang 110000，China

通讯作者: 贺黎明（1984-），男，江西萍乡人，副教授，主要从事全球卫星导航定位及雷达遥感相关理论和应用研究。E⁃mail: heliming@mail.neu.edu.cn

收稿日期: 2021-04-15 修回日期: 2022-07-14

基金资助:

国家自然科学基金项目. 41974028
中央高校基本科研业务专项资金. N2201013

Received: 2021-04-15 Revised: 2022-07-14

作者简介 About authors

张香凝(1996-)，女，辽宁沈阳人，硕士研究生，主要从事合成孔径雷达干涉测量方面的研究E⁃mail:332814679@qq.com , E-mail：332814679@qq.com

摘要

地下煤炭资源大量开采导致的地表形变，引发严重的安全和环境隐患，雷达干涉测量技术是高精度、大范围地表形变监测的重要手段之一。以辽宁省沈阳市沈北新区蒲河煤矿为例，采用SBAS-InSAR技术探测2018—2019年矿区地表形变结果，获取了采煤引起地表形变的时空分布特征，结合采场所在区域的地质条件和变形诱发因素，利用数值模拟技术对观测形变结果进行模拟分析，进而讨论了蒲河煤矿地面沉降在时间和空间上的变形规律和机制。InSAR形变监测结果显示，开采区域内存在两处沉降漏斗，且数值模拟结果与InSAR形变观测值分布规律一致，反演结果接近实际情况，可为相关部门制定地面沉降防治措施提供科学依据。

关键词： 煤矿 ; 地面沉降 ; SBAS-InSAR ; Okada模型

Abstract

The surface deformation caused by large-scale mining of underground coal resources results in serious safety and environmental hazards. Radar interferometry technology is one of the important means of high-precision and large-scale surface deformation monitoring. This paper takes Puhe coal mine in Shenbei New District， Shenyang City， Liaoning Province as an example， SBAS-InSAR method is used to detect the surface deformation results in the mining area from 2018 to 2019， and the temporal and spatial distribution characteristics of surface deformation caused by coal mining are obtained. Combined with the geological conditions and deformation inducing factors of the mining area， numerical simulation technology is used to simulate and analyze the observed deformation results， and then the deformation law and mechanism of land subsidence in time and space of the Puhe coal mine is discussed. The InSAR deformation results show that there are two subsidence centers in the mining area， and the numerical simulation results are consistent with the observed InSAR deformation results， and the inversion results are close to the actual situation， which could provide scientific information for disaster prevention of land subsidence in the mining region.

Keywords： Coal mine ; Ground subsidence ; SBAS-InSAR ; Okada model

PDF (6151KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

张香凝, 贺黎明, 刘翠芝, 王兴杰, 唐永亮, 何蓉花. 基于SBAS-InSAR技术的煤矿开采沉陷监测与分析. 遥感技术与应用[J], 2022, 37(4): 1021-1028 doi:10.11873/j.issn.1004-0323.2022.4.1021

Zhang Xiangning, He Liming, Liu Cuizhi, Wang Xingjie, Tang Yongliang, He Ronghua. Monitoring and Analysis of Coal Mining Subsidence based on SBAS-InSAR Method. Remote Sensing Technology and Application[J], 2022, 37(4): 1021-1028 doi:10.11873/j.issn.1004-0323.2022.4.1021

1 引言

在我国，煤炭在一次性能源中占有主体地位，长期大量的煤炭资源开发，对土地资源、植被生长、水资源和地表构筑物等造成了巨大损失和不利影响。因此，在地下煤炭开采过程中对矿区进行地表形变监测、分析及预测显得尤为重要^［1-2］。

近年来，短基线集干涉测量技术^［3］（Small baseline subsets，SBAS-InSAR）越来越广泛的应用在采煤、采油、地下水开采以及地震或火山运动等引起的地表形变监测中。SBAS-InSAR技术由于时空间基线较短，避免了空间和时间失相干^［4］，与其他时序InSAR技术相比，具有干涉图相干性较高、运算效率较高以及可获取高精度形变监测结果等优点，已在矿山形变监测方面取得了众多成果。2011年尹宏杰等^［5］应用SBAS-InSAR技术对湖南冷水江锡矿山进行形变监测，揭示了矿区沉降漏斗的发展和演化过程；2015年赵伟颖等^［6］利用SBAS-InSAR技术获取了徐州矿区采动区的累积沉降量，通过对比水准测量结果，两者相关性较高；2019年肖亮等^［7］选取21景SAR数据，应用SBAS-InSAR技术获取了湖南常德某典型钻井水溶岩盐矿的时序形变数据，证实了SBAS-InSAR技术应用于钻井水溶岩盐矿山沉降监测的可行性；2020年栾元重等^［8］应用19景Sentinel-1 SAR数据利用SBAS-InSAR技术获取了菏泽市某煤矿2016—2017年间一处较大沉降漏斗的沉降速率，通过对比水准测量结果，证明SBAS-InSAR监测结果精确可靠，且该研究应用沉降速率和时序沉降量准确分析了各个采煤工作面的沉降情况。

1985年，日本学者提出的Okada模型在地震、火山、滑坡以及矿山地面沉降等地质灾害的形变模拟中得到广泛应用^［9］。2012年吴勤伟等^［10］简化弹性半空间位错模型并结合蒙特卡罗法反演了煤矿塌陷机理；2021年任文静等^［11］利用SBAS-InSAR方法结合位错模型成功对店坪煤矿形变进行了监测和反演。

以上研究应用InSAR技术进行矿区地表形变高精度监测，取得了较好的研究结果，但将InSAR技术与矿区煤层参数反演相结合，分析地表形变的时空分布规律及其变形机制等方面有待进一步研究，实验以辽宁省沈阳市沈北新区蒲河煤矿作为研究区域，应用SBAS-InSAR技术对覆盖该区域的大量SAR影像进行处理，获取时间序列的地表形变信息，并将形变监测结果作为观测数据，利用非线性贝叶斯反演算法^［13-14］进行煤层参数反演，利用简化的Okada模型模拟煤层的张性分量在垂直方向上引起的地表形变，并对研究区域的地表沉降分布、动态变化和沉降成因进行分析，以期更加清楚地掌握蒲河煤矿沉降漏斗的时空演化规律及其变形机制，为矿区合理防治地质灾害提供理论依据和技术支撑。

2 研究区域与数据

辽宁省沈阳市沈北煤田煤层顶底板岩性均属软岩层，顶板破碎、易冒顶，工程地质条件及岩体稳定性较差。煤层开采后，上覆岩体产生裂隙或断移，与此同时，地下水、地表水通过裂隙向采空区渗漏，加速上覆岩体破坏，在井下采空区产生顶板冒落，形成冒落带，进而引起岩层和地表形态产生移动，易引发地面塌陷、地面沉陷。图1为研究区域卫星遥感影像及SAR影像覆盖情况。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 研究区域分布图

Fig.1 Distribution of study area

选用C波段的Sentinel-1 SAR数据，时间范围为2018年1月至2019年8月（共42景），SAR卫星参数如表1所示，外部DEM数据为德国TanDEM-X任务框架中获得的分辨率为90 m的DEM。

表1 SAR卫星数据参数

Table 1 SAR satellite data parameters

序号	参数	Sentinel-1
1	轨道方向	升轨
2	入射角	37.1°
3	航向角	-166.2°
4	时间范围	20180111—20190810
5	影像数量	42

新窗口打开| 下载CSV

3 SBAS-InSAR基本原理

2002年，Lanari等^{［3，12，18］}提出小基线集InSAR （Small Baseline Subset InSAR，简称SBAS-InSAR）方法，其基本思想是按照集合内基线距小、集合间基线距大的原则，将所有影像按照一定的时空基线阈值进行配对获得若干个集合，利用奇异值分解（SVD）方法将多个小基线集进行联合求解。

将N+1幅覆盖相同区域的SAR影像图，根据空间基线和时间基线的条件组合干涉对，得到M幅差分干涉图，为简化模型，忽略去相干因素、高程误差、大气误差等因素，则第j幅差分干涉图中像元的相位值为：

\begin{array}{l} δ φ_{1} (r, x) = φ (t_{A} . r . x) - φ (t_{B} . r . x) \\ \approx \frac{4 π}{λ} [d (t_{A} . r . x) - d (t_{B} . r . x)] \end{array}

其中： $λ$ 为雷达的波长； $d (t_{A}, r, x)$ 和 $d (t_{B}, r, x)$ 分别为像元在时间 $t_{A}$ 和 $t_{B}$ 沿雷达视线方向的形变。对差分干涉图相位进行相位解缠， $δ φ_{i} (i = 1, . . ., M)$ 为相对于解缠参考点的相位值。主、辅影像对应的时间序列分别为：

\begin{array}{l} I M = [I M_{1}, \dots, I M_{m}], \\ I S = [I S_{1}, \dots, I S_{m}] \end{array}

若主、辅影像按时间序列排列即为 $I M_{j}$ > $I S_{j} (j = 1, \dots, M)$ ，方程组含N个未知数的M个方程，差分干涉图中相位表示如下：

δ φ = A φ

若基线集中含多个子集，A^TA为奇异矩阵。如果有N个不同的子基线集，那么矩阵A的秩为N-L+1。对矩阵A进行奇异值分解得到：

A = U \sum V^{T}

其中：U为M*N正交矩阵；∑的对角素为奇异值 $σ_{i} (i = 1, . . ., N)$ 。式（4）的最小二乘范数解为：

\hat{φ} = V [\begin{matrix} Σ^{- 1} & 0 \\ 0 & 0 \end{matrix}] U^{T} δ φ

将相位求解转化为相位变化速度求解，则待求参数向量表示为：

V^{T} = [\begin{matrix} V_{1} = \frac{φ_{1} - φ_{0}}{t_{1} - t_{0}}, \dots, \\ V_{N} = \frac{φ_{N} - φ_{N - 1}}{t_{N} - t_{N - 1}} \end{matrix}]

根据奇异值分解，去除DEM误差和大气相位，可求解各时间段内的平均速度 $V$ ，在时间域上积分便可得到时序形变量。

4 SAR数据处理结果与分析

基于42幅Sentinel-1B降轨影像、TanDEM数据，利用SBAS-InSAR技术对沈阳市蒲河煤矿研究区域进行了数据处理，设置时间基线小于72 d，空间基线小于2%，进行时空间基线组合，最大空间基线为109.6 m。

采用SBAS-InSAR技术进行处理，距离向与方位向的分辨率为20 m×20 m，生成干涉相对，根据像对的连接关系，引入外部90 m TanDEM数据对影像对进行配准、去平，应用Goldstein滤波降噪以及最小费用流进行解缠生成216幅干涉图，通过GCP控制点来去除和估算残余相位和解缠相位^［17］，之后对干涉图进行去平和精炼，再次解缠和滤波，提高相位解缠的可信度，设置相干系数阈值为0.2，得到相干性较高的干涉解缠相位，并在此基础上计算时间序列地表位移，最后进行地理编码获取WGS-84坐标系下的地表累积形变量及形变速率^［15］。

图2为获取的蒲河煤矿以2018年1月11日为参考时间的形变时间序列图，每幅图左上角为SAR影像的获取时间，其中负值表示地面下沉，正值表示地面抬升，红色区域表示沉降严重区域。从结果可清晰地看到有两两个沉降漏斗出现在崔公堡东部和大望村北部，在接近两年的时间里，矿区两处沉降中心的最大形变量分别达-157.26 mm和-193.1 mm。东南处的沉降中心在大望村东北处；西北处的沉降中心在崔公堡东部。根据时间序列地表形变的发展规律来看，沉降中心持续扩大且沉降量级不断变大，及时、准确的进行地面沉降监测对减小矿区开采对周边环境带来的影响尤为重要。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 LOS向时间序列地表形变图（时间范围20180111—20190530）

Fig.2 Surface deformation map of Los time series （Time range 20180111—20190530）

2018年1月11日—2019年8月10日期间的平均沉降速率如图3（a）所示，矿区整体平均沉降速率集中分布在>-40 mm/a，且矿区两处较为明显的沉降漏斗（蓝色椭圆处）最大沉降速率分别为-88.74 mm/a和-100.68 mm/a，图3（a）中三角点位为3处最大沉降速率点分布位置，为进一步分析沉降漏斗的发育状况，选取两个沉降漏斗A、B区域，并沿蒲河矿区分布做一条剖线分析其在时间序列上的沉降变化情况，提取剖线上的时序累计沉降量以及沉降速率结果，绘制折线图。如图3（b）所示，下沉曲线保持良好的连续性，两个沉降漏斗在时间范围内的沉降量在逐渐增加，沉降漏斗A最大累计沉降量达-157.26 mm；沉降漏斗B累计沉降量有两个峰值-193.1 mm、 -183.98 mm，且剖线上的沉降结果显示沉降漏斗B包含的两个小漏斗有随着煤层开采而发育的趋势。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 形变速率场、沿剖线各期沉降量以及沉降区内开裂的建筑物图

Fig.3 The figure of deformation rate field， settlement amount in each period along the profile line and cracking building in the settlement area

由于沈北蒲河煤田煤层顶底板岩性属软岩层^［16］，岩体稳定性较差，煤层开采后，在地质环境条件、采矿工程条件等影响因素的作用下，易引起岩层和地表形态的移动，产生地面塌陷或地裂缝等地质灾害。结合野外调查数据，沉降漏斗附近区域的崔公堡、大望村的部分道路和房屋均出现严重开裂如图3（c）所示。由此可见，矿区塌陷不仅带来了巨大的经济损失，而且严重威胁居民的人身安全，因此需要针对矿区塌陷采取合理措施，及时、准确地进行地面沉降动态监测，以便最大程度地减少矿区开采对周边环境带来的影响。

5 InSAR形变场模拟

为进一步分析沉降漏斗的发育状态，实验以SBAS-InSAR形变结果作为观测量，利用非线性贝叶斯反演算法进行煤层参数的反演，应用Okada组合模型进行形变模拟，对该地区的地表沉降分布以及动态变化进行分析。

1985年，日本学者Okada^［9］通过对已有的弹性半空间断层引起地表形变研究成果的分析，提出了点源和有限矩形面源的断层位错引起地表形变的通用表达式，后来被人们称为Okada断层位错模型。1992年，Okada对该断层位错模型理论进一步完善。

由于煤矿开采引起的形变主要体现在垂直方向上，假定煤层倾角为0，仅分析张性分量引起的地表形变，简化后Okada模型的主要参数有6个，包括错动面沿走向的长度L、沿倾向的宽度W、错动面几何中心的深度d、错动面与正北方向之间的夹角（走向φ）、错动面的几何中心M点在地表的投影坐标（x₀，y₀）。实验应用以上6个参数对沉陷漏斗进行反演和分析。

由于蒲河煤矿有多个煤炭开采层，考虑Okada模型的参数自由度和特征以及蒲河煤矿的两个沉陷漏斗分布情况，假设煤田地表形变是由地下3个有限矩形源共同作用下引起的，一个有限矩形源位于矿区沉降漏斗A中心处。从图3（b）中的剖线来看，由于沉降漏斗B有两处较为明显的沉降中心，故另两个有限矩形源分别分布在沉降漏斗B的两个沉降中心处（如图3（a）虚线矩形框为矩形源分布处），使用Okada组合模型对煤层参数进行反演，最佳拟合参数如表2所示，表中X表示笛卡尔坐标系下煤层中心的横坐标；Y表示笛卡尔坐标系下煤层中心的纵坐标。

表2 蒲河煤矿Okada组合模型最佳拟合参数

Table 2 The best fitting parameters of Okada combination model in Puhe coal mine

序号	参数	最佳拟合参数
序号	参数	组一	组二	组三
1	X/m	-15.99	644.91	303.95
1	Y/m	-24.01	-1080.1	-939.95
2	长度/m	615.99	646.42	640.73
3	宽度/m	352.91	229.36	320.6
4	深度/m	482.09	568.67	481.29
5	走向/°	24.52	11.01	147.68
6	张量/m	-0.59	-0.04	-0.67

新窗口打开| 下载CSV

第一组Okada组合模型反演中，对SBAS-InSAR形变监测结果建立独立坐标系，坐标原点为（123.548° E， 41.995° N），投影到笛卡尔坐标系下煤层中心的平面坐标为（-15.99 m E，24.01 m N），煤层深度482.09 m，有限矩形长为615.99 m，有限矩形宽为352.91 m，走向角（以北方向顺时针旋转到断层面的角度）为24.52°，位错面开度为-0.59 m；第二组Okada组合模型反演中，对SBAS-InSAR形变监测结果建立独立坐标系，坐标原点为（123.548° E， 41.995° N），投影到笛卡尔坐标系下煤层中心的平面坐标为（644.91 m E，-1 080.1 m N），煤层深度568.67 m，有限矩形长为646.42 m，有限矩形宽为229.36 m，走向角（以北方向顺时针旋转到断层面的角度）为11.01，位错面开度为-0.04 m；第三组Okada组合模型反演中，对SBAS-InSAR形变监测结果建立独立坐标系，坐标原点为（123.548° E， 41.995° N），投影到笛卡尔坐标系下煤层中心的平面坐标为（303.95 m E，-939.95 m N），煤层深度481.29 m，有限矩形长为640.73 m，有限矩形宽为320.6 m，走向角（以北方向顺时针旋转到断层面的角度）为147.68°，位错面开度为-0.67 m。即该区域地表形变是由分别位于东-15.99 m，北24.01 m，482.09 m深的有限矩形源，东644.914 m，北1 080.10 m，481.29 m深的有限矩形源以及东303.95 m，北939.95 m，481.29 m深的有限矩形源共同变化引起的。蒲河煤矿埋深约530 m^［20］，地面海拔高度约100 m，煤层深度约为430 m左右，本文采用的3个有限矩形源，对蒲河煤矿的煤层深度模拟结果较为准确，与实际煤层深度较为接近。从形变模拟得到的3个有限矩形源的分布来看，崔公堡采区有限矩形源形状接近正方形，沉降漏斗的形变走向在横向和纵向均有发育，而大望村采区受两个沉降漏斗共同作用，且根据煤炭开采的实际情况，两个有限矩形源走向、深度均不同，两者为中心位置水平分布且边缘处有些许重叠的长方形，故该采区的形变偏横向发育。

图4（a）为InSAR观测得到的解缠相位图，图4（b）为Okada组合模型模拟InSAR解缠相位图；两个结果较为接近，较为完整地模拟出两个沉降漏斗中心区域，且对于沉降漏斗中心区域外也得到了较好的模拟效果。图4（c）为解缠相位的观测形变与模拟形变之间的残差。从图中可以看出，整个研究区域组合模型模拟地表形变的残差分布较小，精度较高。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 InSAR结果与模拟结果对比图

Fig.4 Comparison of InSAR results and simulation results

对比SBAS InSAR形变结果和Okada组合模型模拟的形变结果，即提取图4中剖线A1A2上的形变量如图5所示，红色实线为SBAS-InSAR观测形变值，蓝色虚线为模拟形变值，二者拟合程度较好，且现象一致，并在同一位置显示了漏斗式沉降结果。对形变残差进行统计和拟合，如图6所示，整体残差以0为中心近似呈正态分布，95%的残差分布在［-50 mm， 50 mm］，和InSAR形变观测值分布规律一致，且残差曲线形态较高瘦，均值1.47 mm，标准差为16.76 mm，说明本次模拟残差控制较好，反演结果接近实际情况。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 剖线A1A2模拟和观测形变对比图

Fig.5 The comparison of simulated and observed deformation of profile A1A2

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 残差分布直方图

Fig.6 The histogram of residual distribution

6 结论

基于2018年1月—2019年8月间的42景Sentinel-1 SAR数据，使用SBAS-InSAR技术以沈阳市蒲河煤矿为研究区域进行地表形变监测，以该形变结果作为反演的观测量，使用Okada组合模型对煤层参数进行反演，并将观测结果与模拟结果进行对比分析，主要研究结论如下：

（1）研究区域主要分布两个沉降漏斗，形成较为明显的沉降中心，最大累计形变量分别达 -157.26 mm和-193.1 mm，对应的平均沉降速率达 -88.74 mm/a和-100.68 mm/a。通过野外调查可知，沉降漏斗对应于煤矿开采位置，该区域由于长期的煤田开采导致地表沉降，在沉降漏斗附近形成了沉陷水域，致使区内民宅及公共建筑物严重受损，道路及地下管线等市政设施遭到破坏，农田退化，土地资源受损严重，造成了巨大的经济损失。

（2）对比InSAR形变结果和Okada组合模型模拟的形变结果，二者呈现较好的一致性，Okada组合模型对研究区域内煤层参数的反演结果较为可靠，残差分布直方图显示整体残差以0为中心近似呈正态分布，且残差较小，充分证明应用3个有限矩形源的Okada组合模型对蒲河煤矿的煤层参数反演结果较好。

本文通过形变模拟反演验证了SBAS-InSAR技术监测地下煤层开采导致的地表形变精度较高，然而基于单轨数据的SBAS-InSAR结果在监测多煤层开采引起的复杂地面沉降存在一定误差，在后续的研究中可应用多种InSAR技术结合或多轨道SAR数据融合的方法得到更高精度的形变监测结果^［21-22］。Okada组合模型模拟中通过假设多个有限矩形源对煤层参数进行反演，结果较为可靠，可为后期蒲河煤矿地面沉陷防治计划制定及沉降漏斗发育特征估计提供科学依据，此外，基于数值模型的煤田开采参数反演也可为开采变形机制研究提供参考。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Liu

Xiaoling

， Wei

Aolin

， Wang

， et al.

Analyses of development characteristics and formation of geological disasters in coal mine area of northern Shaanxi

［J］.The Chinese Journal of Geological Hazard and Control，2016，27（4）：70-73.