刘翠芝¹, 王兴杰¹, 贺黎明¹, 唐永亮²
1. 东北大学 资源与土木工程学院, 辽宁 沈阳 110819;
2. 辽宁省地质环境监测总站, 辽宁 沈阳 110032
收稿日期：2021-07-12
基金项目：国家自然科学基金资助项目(41974028)；中央高校基本科研业务费专项资金资助项目(N180104014)。
作者简介：刘翠芝(1965-)，女，黑龙江绥化人，东北大学副教授。

摘要：近距离多煤层开采引起的采空区地面沉降在当前工矿区沉陷中占有较大的比例, 而天基InSAR技术在工矿区大范围形变监测方面具有显著优势.以2018年1月至2019年8月间的43景C波段Sentinel-1B SAR数据为基础, 通过时序InSAR技术对近距离多煤层开采的康平煤田进行了大范围沉降监测.以Okada模型为基础, 结合近距离多煤层分布特征, 采用叠加理论的原则, 提出了一种基于多源模型的近距离多煤层参数反演方法.通过定量计算弹性半空间下双源模型和三源模型近距离多煤层参数反演的精度, 发现三源模型精度更高, 更符合小康矿近距离多煤层开采的实际特征.本研究可以为基于InSAR技术的近距离多煤层开采沉降监测及参数反演提供参考.
关键词：近距离多煤层参数反演时序InSAR多源Okada模型形变模拟
Parameter Inversion Method of Close Multi-coalseam Based on Time Series InSAR
LIU Cui-zhi¹, WANG Xing-jie¹, HE Li-ming¹, TANG Yong-liang²
1. School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China;
2. Liaoning Geological Environment Monitoring Station, Shenyang 110032, China
Corresponding author: HE Li-ming, E-mail: heliming@mail.neu.edu.cn.

Abstract: Ground subsidence caused by close multi-coalseam mining occupies a large proportion in current mining activities. Due to the significant advantages of space-based InSAR technique in large-scale deformation monitoring, the subsidence caused by the close multi-coalseam mining in Kangping coalfield was monitored based on the 43 scenes C-band Sentinel-1B SAR data from January 2018 to August 2019. Considering the distribution characteristics of the close multi-coalseam, a close multi-coalseam parameter inversion method using superposition theory based on multi-source Okada model was proposed. By quantitatively calculating the accuracy of inversion results of dual-source model and three-source model of close multi-coalseam in elastic half space, it is found that the three-source model has a higher accuracy, which is consistent with the actual characteristics of the close multi-coalseam of the Xiaokang mine. This study provides some references to the subsidence monitoring and parameter inversion of close multi-coalseam mining.
Key words: close multi-coalseamparameter inversiontime series InSARmulti-source Okada modeldeformation modeling
煤矿开采带来的采空沉陷是我国诸多资源开采型城市面临的严重的地质灾害之一^[1], 由于煤采区开采的持续性和区域多煤层地质结构的复杂性, 地表沉降主要表现为非线性和大梯度的特征, 给当地的经济发展和人民生活带来了极大的危害, 严重制约资源型城市的可持续性发展^[2].煤采区地表沉降是地下煤炭资源被持续采出后, 开采区域周围岩体的原始应力平衡受到破坏, 造成上覆岩层和地表失去支撑而发生强烈的非连续和非线性形变, 且随着时间推移应力重新分布达到新平衡的复杂过程, 是开采煤层几何形状、位置等几何状态参数在地表的综合直观体现.因此开展对煤田采空区的形变监测和煤层参数反演的研究对区域地质灾害防范与煤矿安全生产有重要的意义^[3-4].
传统的全球导航卫星系统(global navigation satellite system, GNSS)和水准等点式测量方法无法实现对沉陷区大范围的有效监测, 而雷达差分干涉测量技术(interferometric synthetic aperture radar, InSAR)是一种典型的面式对地观测手段, 具有高精度、高分辨率、空间连续覆盖的典型特征, 特别适合大范围、大梯度的地表沉降监测^[5-7].基于InSAR监测获得的覆盖整个煤采区的地表形变场, 结合地球物理非线性反演算法, 可快速准确反演得到煤层几何形状、位置以及产状等参数信息, 这对于分析近距离多煤层开采引起地表形变的机理和演变规律至关重要.
近年, 国内外****对煤层参数的反演研究逐渐增多, 但大多集中于两个方向.一种是基于FLAC^3D有限元或UDEC离散元等数值模型来对近距离多煤层共同采动引起的地表位移进行分析^[8-9], 这类方法对近距离多煤层开采深入阶段的实际状况进行了较为周密的考虑, 具有较高的实际应用指导价值, 但由于其模拟过程中需要已知详尽的多煤层位置、产状等先验参数信息, 并且缺乏实际监测的地表形变约束, 因此严重制约了在很多未知先验信息的煤层参数反演过程中的应用; 另一种是依托InSAR监测获得的地面沉降量信息, 基于Okada模型选取部分单个煤层巷道反演其几何状态参数^[10], 这类方法可以对单个煤层的几何位置等参数进行合理的反演解释, 但单煤层巷道不能反映多煤层开采的实际状况, 因此缺乏对近距离多煤层实际开采状态的合理表达.
为此, 本文提出了一种基于多源Okada模型的受地表形变约束的近距离多煤层参数反演方法, 有效避免FLAC^3D和UDEC等数值模型不受实际地表形变约束的缺陷.该方法依托小基线集技术(small baseline subset, SBAS)获取的研究区地面沉降信息作为观测数据, 实现多源Okada模型在近距离多煤层参数反演中的应用.本文的目的是通过研究近距离多煤层开采引发地表变形的规律, 以实测地表形变作为约束, 反演得到近距离多煤层的几何形状、位置以及产状等参数信息, 对于我国普遍分布的煤采区, 尤其是大量已关闭或废弃煤矿的地面沉陷范围圈定以及多煤层采空位置、产状等信息获取具有重要借鉴意义, 进而可为煤采区, 特别是已关闭或废弃工矿区地面沉降监测和综合治理提供依据.
1 近距离多煤层参数反演方法近距离多煤层共同开采是目前全球范围内各矿场在实际开采作业过程中所面临的一种复杂情况.研究其开采规律, 并根据实际监测的地表形变快速准确反演地下复杂的多煤层几何形状、位置、产状等参数对于煤矿采空区的综合治理极其重要.本文基于Okada模型^[11], 结合近距离多煤层分布的实际特征, 依托叠加原理, 提出了一种基于多源模型的近距离多煤层参数反演方法, 其示意图如图 1所示.
图 1(Fig. 1)

图 1 近距离多煤层的多源模型示意图Fig.1 The multi-source model for the close multi-coalseam (a)—煤层分布等效图；(b)—位错模型示意图；(c)—多源位错模型分布图.

图 1a是近距离多煤层开采过程中多个开采层和中间层的等效表示, 图 1b是Okada有限矩形源模型, 共有10个参数, 分别为开采面几何中心P在地表的投影坐标(X, Y)、开采面沿着走向的长度L、沿倾向的宽度W、煤层开采中心深度depth、走向角strike(开采掘进方向与正北方向的夹角)、倾向角dip(开采面与水平方向的夹角)、rake(断层破裂时相对于断层走向移动的角度)、破裂方向的位错分量slip、张性分量opening; 图 1c是依托多煤层的分布情况, 将Okada有限矩形源模型应用于多煤层的综合图.
基于多源模型的近距离多煤层参数反演方法首先使用基于梯度的自适应四叉树采样对InSAR监测数据进行降采样, 通过对较高位移梯度的区域进行更精细的划分来减少InSAR图像内部和边缘的数据间隙的影响, 在保证有效采样的基础上降低数据点位密度, 提升反演计算效率.同时通过计算较高位移梯度和数据边缘信息, 再结合InSAR形变图的空间分布特征和研究区多煤层的分布特征, 预估多源模型的数量和模型参数的初始值等.然后利用已有的先验信息和似然函数构建形变观测数据d和通过叠加理论计算得到的多源模型形变模拟数据G(t)之间的数学关系, 接下来采用非线性贝叶斯反演算法^[12-13]进行多煤层参数反演, 在反演过程中, 使用抽样计算模型参数的后验概率分布来判定多煤层参数拟合结果.多源模型形变模拟数据G(t)和后验概率分布的计算公式如下:

(1)

(2)

其中: G(t)_s和b_s分别表示多源模型中不同源模拟得到的形变和所占权重; s为多源模型的数量; d为M维观测数据; t为N维模型参数; p(d)为归一化常数; p(t)为多源模型参数的先验信息; p(t/d)为在已知先验信息下当前参数的后验概率; p(d/t)为多源模型的似然函数(likelihood function).
本文依据模型参数拟合值的分布特征, 选择多维高斯表达式L(t)来表示似然函数, 如式(3)所示, 其中C_d表示数据协方差矩阵.

(3)

在抽样计算多源模型后验概率时, 采样方法的选取对于后验概率的计算结果和效率至关重要.本文采用基于Metropolis-Hastings接受准则的马尔科夫链蒙特卡罗(Markov chain Monte Carlo, MCMC)算法^[14]在空间上进行采样, 该方法可以大大提高采样接受率, 弥补了传统蒙特卡罗方法的抽样分布不能随着模拟的进行而动态改变的缺陷.采样时, 首先根据先验信息和模拟退火(simulated annealing, SA)算法估计参数并作为初始状态参数, 然后在搜寻空间上建立一个平稳分布π(i)的各态遍历的马尔科夫过程, 如式(4)所示:

(4)

π(i)是t=1时刻处于状态q_i的概率.之后通过状态转移矩阵A生成大量样本, 如式(5)所示:

(5)

(6)

其中：a_ij是在t时刻处于状态q_i的条件下向t+1时刻转移到状态q_j的概率.其中i表示状态序列, q表示可能的状态.
最后计算所有样本多源模型反演参数的后验概率, 作为多煤层最优参数的选择指标.在寻找最优参数的过程中, 选择后验概率最大的一组多源模型参数即为最佳拟合参数.
2 研究区InSAR形变监测2.1 研究区与数据源康平煤田位于辽宁省沈阳市康平县, 属于地下开采的隐蔽型多煤层联合开采煤田(图 2).其中大平煤矿辖区面积28.61 km², 年生产能力330万t, 共划分为7个采区, 其中有5个采区处于三台子水库下; 小康煤矿辖区面积为30.08 km², 被划分为6个采区, 核定年开采量270万t.本文选用2018年1月至2019年8月间覆盖康平煤田的43景欧空局C波段Sentinel-1B降轨数据进行形变提取.SAR卫星的航向角为-166.19°, 入射角为37.13°.DEM数据选用德国TanDEM-X获取的空间分辨率为3 rad·s的DEM, 用于影像配准、平地效应去除、轨道精炼等.
图 2(Fig. 2)

图 2 研究区位置图Fig.2 Location of the studied area

2.2 InSAR形变场监测本文采用SBAS-InSAR技术对研究区域进行沉降监测.SBAS技术由Berardino等^[15]提出, 为了减弱时空失相干效应的影响, 小基线集技术从获取的SAR影像中选择时空基线均小于某一阈值的干涉像对生成差分干涉图, 同时根据基线情况将这些干涉像对划分成若干小基线集, 再基于最小二乘法解算每个小基线集的形变时间序列, 然后利用奇异值分解法(singular value decomposition, SVD)将多个小基线集联合求解, 得到整个时间段上雷达视线向(line of sight, LOS)的形变时间序列^[16-17].流程见图 3.
图 3(Fig. 3)

图 3 SAR数据处理流程图Fig.3 Flow chart of the SAR data processing

本文在获取研究区形变结果时, 基于研究区近距离多煤层开采的地表实际覆盖情况, 充分考虑时空基线和多普勒效应的影响, 设置空间基线阈值为临界基线的2%;时间基线阈值为72 d, 生成了216组干涉像对(图 4).
图 4(Fig. 4)

图 4 时间-空间基线图Fig.4 The spatio-temporal baseline

基于SBAS-InSAR技术, 获取了研究区在雷达视线向的累积沉降量, 结果如图 5所示.在监测范围内, 大平矿采区最大沉降速率为-99.7 mm/a, 最大累积沉降量为-136.6 mm; 小康矿近距离多煤层开采引起的最大沉降速率为-124.4 mm/a, 最大累积沉降量为-194.8 mm; 除此之外, 其他区域沉降现象不明显, 这表明该区域整体较为稳定, 地面沉降主要是由煤炭开采所致.
图 5(Fig. 5)

图 5 康平煤田累积沉降图Fig.5 Cumulative subsidence of the Kangping coalfield

为进一步分析两个沉降区不同位置的形变过程, 分别选取了Q1和Q2两个特征点, 刻画其沉降时间序列, 结果如图 6所示.蓝色实线表示大平矿采区特征点Q1的沉降时间序列, 红色实线表示小康矿采区特征点Q2的沉降时间序列.可以发现2018年1月~2019年8月间, 大平矿沉降区的Q1点累积沉降量为-111.7 mm; 小康矿沉降区Q2点的累积沉降量最大, 达到-178.4 mm.持续的地下煤炭开采已造成康平煤田区域产生了巨大的地表形变梯度和影响范围, 使当地居民住宅等设施被积水淹没或产生结构性损坏而被废弃，如图 5所示.
图 6(Fig. 6)

图 6 时间序列累积沉降图Fig.6 Diagram of the time series cumulative subsidences

3 小康矿多煤层参数反演与分析3.1 多源模型反演结果本文以时序SBAS-InSAR方法获取的小康矿形变监测数据为基础, 综合考虑形变反演过程中数据点位密度对反演计算效率的影响, 选用1.44×10^-4作为自适应四叉树方差阈值进行数据降采样, 然后通过较高位移梯度和数据边缘信息, 选取2和3作为多源模型的模型数, 分别使用非线性贝叶斯反演算法对近距离多煤层参数进行反演.Okada模型中slip表示断层破裂方向的位错分量, rake表示断层破裂时相对于断层走向移动的角度, opening表示断层垂直方向运动的张性分量.通过分析小康矿形变区的地表沉降特征, 发现小康矿开采过程中多煤层倾角很缓, 近乎水平开采, 且地表沉降主要集中在垂直方向, 因此本文简化了Okada模型, 假设断层倾角为0°, 只分析多煤层开采导致垂直方向的张性破裂引发的地表形变.在反演过程中使用独立的测量坐标系, 参考基准为[123.395°E, 42.65°N], 根据小康矿区域煤系地层岩性, 泊松系数确定为0.25.依托以上信息, 使用基于双源模型和三源模型的近距离多煤层参数反演方法分别反演小康矿近距离多煤层参数, 获得最佳拟合结果如表 1所示.
表 1(Table 1)

表 1 小康矿多源模型最佳拟合参数Table 1 Best fitted parameters of the multi-source model of the Xiaokang mine 参数 双源模型 三源模型 一组 二组 一组 二组 三组 X/m -511.14 -29.85 -905.76 348.67 -83.89 Y/m -355.36 241.14 -339.33 2.15 152.66 L/m 768.14 1415.11 834.02 1774.04 1560.53 W/m 401.74 400.01 400.70 577.76 447.77 depth /m 530.40 579.24 570.93 596.55 574.42 strike/(°) 3.56 178.94 176.78 141.98 177.36 opening/m -0.44 -0.47 -0.49 -0.08 -0.38

表 1 小康矿多源模型最佳拟合参数 Table 1 Best fitted parameters of the multi-source model of the Xiaokang mine

对比研究区近距离多煤层实际开采参数与多源模型反演最佳拟合参数, 评价多源模型参数反演方法的可靠性.本文以多煤层开采深度作为评价指标进行分析, 小康矿开采资料显示, 该形变区域对应S2N5和S2N7综采面的位置, 多煤层开采深度分别对应为427~593 m和534~614 m之间.由表 1可知, 双源模型反演深度为530.40和579.24 m; 三源模型反演得到多煤层深度分别为570.93, 596.55和574.42 m; 两种多源模型反演深度均与实际开采深度一致, 表明了多源模型在该区域多煤层参数反演中的可靠性.
3.2 对比分析基于双源和三源Okada模型反演得到的多煤层参数, 分别对地表形变进行正演, 模拟地表形变结果, 如图 7和图 8所示，分析模拟形变和观测形变之间的残差(图 7c和图 8c), 发现三源模型的残差控制更好, 符合多煤层开采引起的实际形变.为了进一步对比观测形变、双源模型和三源模型模拟形变之间的差异, 本文通过比较剖线A1A2A3上的形变信息, 评价不同反演模型的精度.图 9a中, 红色实点表示观测形变, 粉色虚线表示双源模型模拟形变, 蓝色虚线表示三源模型模拟形变.可以看出两种模型都反演得到两个漏斗式的沉降中心, 但三源模型拟合程度更好.图 9b显示了双源模型和三源模型的误差分布, 其中蓝色实线表示的三源模型误差在更靠近0的附近振荡, 相较于双源模型, 反演精度更高.
图 7(Fig. 7)

图 7 观测形变、双源模型模拟形变和残差结果对比图Fig.7 Comparison of the observed deformation, simulated deformation of dual-source model and residual results (a)—InSAR观测图; (b)—双源模型模拟图; (c)—双源模型残差图.

图 8(Fig. 8)

图 8 观测形变、三源模型模拟形变和残差结果对比图Fig.8 Comparison of the observed deformation, simulated deformation of three-source model and residual results (a)—InSAR观测图; (b)—三源模型模拟图；(c)—三源模型残差图.

图 9(Fig. 9)

图 9 沿剖线A1A2A3的观测形变和模拟形变对比分析图Fig.9 Comparison of observed and modelled deformation along profile A1A2A3 (a)—观测值与模拟值; (b)—模型误差.

为了定量评价不同反演模型的精度, 本文根据公式(7), (8)分别计算两种反演模型值与InSAR观测值在剖线A2A3上的均方根误差(root mean squared error, RMSE)和线性拟合优度(R²).

(7)

(8)

其中: n表示样本个数; d_i表示各样本所对应的InSAR观测值; m_i为多源模型模拟值; 表示模拟值的均值.
计算双源模型的RMSE为15.9 mm, R²为0.399 0;三源模型的RMSE为6.7 mm, R²为0.791 2.对比发现三源模型的均方根误差更小, 与InSAR观测结果的相关性更高, 证明三源Okada模型更符合该区域近距离多煤层共同采动的实际形态.
4 结论1) 研究区内主要分布有两个沉降漏斗, 其他区域无明显沉降现象.经实地调查验证, 发现沉降区对应于地下煤炭开采的位置, 说明该区域地面沉降主要是由煤炭采空所致.在2018年1月~2019年8月间, 小康矿采区最大累积沉降量为-194.8 mm, 最大沉降速率达到-124.4 mm/a.大平矿采区最大累积沉降量为-136.4 mm, 最大沉降速率为-99.7 mm/a.
2) 为了解决近距离多煤层参数的反演问题, 本文提出了一种基于多源模型的近距离多煤层参数反演方法.该方法首先依据数据降采样计算得到的较高位移梯度和数据边缘信息估计有限矩形源数目, 然后以叠加理论为原则, 将多个矩形源叠加作为反演模型, 最终通过非线性贝叶斯反演算法计算多源模型参数的后验分布概率, 选取后验概率最大的一组多煤层参数作为最佳拟合参数.
3) 将所提出的方法应用于近距离多煤层开采的小康矿进行验证.通过定量计算, 发现三源模型的RMSE相较于双源模型减少了9.2 mm, R²达到了0.791 2.结果表明三源模型反演精度更高, 更符合小康矿近距离多煤层共同采动的实际特征.
参考文献

[1]	吴立新, 高均海, 葛大庆, 等. 工矿区地表沉陷D-InSAR监测试验研究[J]. 东北大学学报(自然科学版), 2005, 26(8): 778-782. (Wu Li-xin, Gao Jun-hai, Ge Da-qing, et al. Experimental study on D-InSAR monitoring of surface subsidence in industrial and mining areas[J]. Journal of Northeastern University(Natural Science), 2005, 26(8): 778-782. DOI:10.3321/j.issn:1005-3026.2005.08.018)
[2]	Samsonov S, D'Oreye N, Smets B. Ground deformation associated with post-mining activity at the French-German border revealed by novel InSAR time series method[J]. International Journal of Applied Earth Observation and Geoinformation, 2013, 23: 142-154. DOI:10.1016/j.jag.2012.12.008
[3]	Zhang B, Wu S, Ding X, et al. Use of multiplatform SAR imagery in mining deformation monitoring with dense vegetation coverage: a case study in the Fengfeng Mining Area, China[J]. Remote Sensing, 2021, 13: 3091. DOI:10.3390/rs13163091
[4]	Chen B, Mei H, Li Z, et al. Retrieving three-dimensional large surface displacements in coal mining areas by combining SAR pixel offset measurements with an improved mining subsidence model[J]. Remote Sensing, 2021, 13: 2541. DOI:10.3390/rs13132541
[5]	Lu Z, Wicks Jr C. Characterizing 6 August 2007 Crandall Canyon mine collapse from ALOS PALSAR InSAR[J]. Geomatics, Natural Hazards & Risk, 2010, 1(1): 85-93.
[6]	Liao M, Zhang R, Lv J, et al. Subsidence monitoring of fill area in Yan'an new district based on Sentinel-1A time series imagery[J]. Remote Sensing, 2021, 13: 3044. DOI:10.3390/rs13153044
[7]	Guo L, Gong H, Ke Y, et al. Mechanism of land subsidence mutation in Beijing Plain under the background of urban expansion[J]. Remote Sensing, 2021, 13: 3086. DOI:10.3390/rs13163086
[8]	李全生, 张忠温, 南培珠. 多煤层开采相互采动的影响规律[J]. 煤炭学报, 2006, 31(4): 425-428. (Li Quan-sheng, Zhang Zhong-wen, Nan Pei-zhu. Mining impact pattern of mutil-coal seam[J]. Journal of China Coal Society, 2006, 31(4): 425-428. DOI:10.3321/j.issn:0253-9993.2006.04.004)
[9]	张立亚, 邓喀中. 多煤层条带开采地表移动规律[J]. 煤炭学报, 2008, 33(1): 28-32. (Zhang Li-ya, Deng Ka-zhong. The law of surface movement for multi-coal seam strip mining[J]. Journal of China Coal Society, 2008, 33(1): 28-32. DOI:10.3321/j.issn:0253-9993.2008.01.006)
[10]	任文静, 贾洪果, 闫斌. SBAS-InSAR方法支持下的矿区地表沉降监测及参数反演[J]. 测绘通报, 2021(3): 113-117, 155. (Ren Wen-jing, Jia Hong-guo, Yan Bin. Monitoring and parameter inversion of ground subsidence in mining area based on SBAS-InSAR method[J]. Bulletin of Surveying and Mapping, 2021(3): 113-117, 155.)
[11]	Okada Y. Surface deformation due to shear and tensile faults in a half-space[J]. International Journal of Rock Mechanics & Mining Sciences & Geomechanics Abstracts, 1986, 23(4): 128.
[12]	尹彬, 胡祥云. 非线性反演的贝叶斯方法研究综述[J]. 地球物理学进展, 2016, 31(3): 1027-1032. (Yin Bin, Hu Xiang-yun. Overview of nonlinear inversion using Bayesian method[J]. Progress in Geophysics, 2016, 31(3): 1027-1032.)
[13]	Bagnardi M, Hooper A. Inversion of surface deformation data for rapid estimates of source parameters and uncertainties: a Bayesian approach[J]. Geochemistry, Geophysics, Geosystems, 2018, 19(7): 2194-2211. DOI:10.1029/2018GC007585
[14]	Jiang Y, González P J. Bayesian inversion of wrapped satellite interferometric phase to estimate fault and volcano surface ground deformation models[J]. Journal of Geophysical Research: Solid Earth, 2020, 125(5): e2019JB018313.
[15]	Berardino P, Fornaro G, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline SAR interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2375-2383. DOI:10.1109/TGRS.2002.803792
[16]	Usai S, Klees R. SAR interferometry on a very long time scale: a study of the interferometric characteristics of man-mad features[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(4): 2118-2123. DOI:10.1109/36.774730
[17]	He L, Wu L, Liu S, et al. Mapping two-dimensional deformation field time-series of large slope by coupling DInSAR-SBAS with MAI-SBAS[J]. Remote Sensing, 2015, 7(9): 12440-12458. DOI:10.3390/rs70912440