刘翠芝, 王兴杰, 贺黎明, 唐永亮

(1. 东北大学 资源与土木工程学院, 辽宁 沈阳 110819; 2. 辽宁省地质环境监测总站, 辽宁 沈阳 110032)

煤矿开采带来的采空沉陷是我国诸多资源开采型城市面临的严重的地质灾害之一[1],由于煤采区开采的持续性和区域多煤层地质结构的复杂性,地表沉降主要表现为非线性和大梯度的特征,给当地的经济发展和人民生活带来了极大的危害,严重制约资源型城市的可持续性发展[2].煤采区地表沉降是地下煤炭资源被持续采出后,开采区域周围岩体的原始应力平衡受到破坏,造成上覆岩层和地表失去支撑而发生强烈的非连续和非线性形变,且随着时间推移应力重新分布达到新平衡的复杂过程,是开采煤层几何形状、位置等几何状态参数在地表的综合直观体现.因此开展对煤田采空区的形变监测和煤层参数反演的研究对区域地质灾害防范与煤矿安全生产有重要的意义[3-4].

传统的全球导航卫星系统(global navigation satellite system,GNSS)和水准等点式测量方法无法实现对沉陷区大范围的有效监测,而雷达差分干涉测量技术(interferometric synthetic aperture radar, InSAR)是一种典型的面式对地观测手段,具有高精度、高分辨率、空间连续覆盖的典型特征,特别适合大范围、大梯度的地表沉降监测[5-7].基于InSAR监测获得的覆盖整个煤采区的地表形变场,结合地球物理非线性反演算法,可快速准确反演得到煤层几何形状、位置以及产状等参数信息,这对于分析近距离多煤层开采引起地表形变的机理和演变规律至关重要.

近年,国内外学者对煤层参数的反演研究逐渐增多,但大多集中于两个方向.一种是基于FLAC3D有限元或UDEC离散元等数值模型来对近距离多煤层共同采动引起的地表位移进行分析[8-9],这类方法对近距离多煤层开采深入阶段的实际状况进行了较为周密的考虑,具有较高的实际应用指导价值,但由于其模拟过程中需要已知详尽的多煤层位置、产状等先验参数信息,并且缺乏实际监测的地表形变约束,因此严重制约了在很多未知先验信息的煤层参数反演过程中的应用;另一种是依托InSAR监测获得的地面沉降量信息,基于Okada模型选取部分单个煤层巷道反演其几何状态参数[10],这类方法可以对单个煤层的几何位置等参数进行合理的反演解释,但单煤层巷道不能反映多煤层开采的实际状况,因此缺乏对近距离多煤层实际开采状态的合理表达.

为此,本文提出了一种基于多源Okada模型的受地表形变约束的近距离多煤层参数反演方法,有效避免FLAC3D和UDEC等数值模型不受实际地表形变约束的缺陷.该方法依托小基线集技术(small baseline subset,SBAS)获取的研究区地面沉降信息作为观测数据,实现多源Okada模型在近距离多煤层参数反演中的应用.本文的目的是通过研究近距离多煤层开采引发地表变形的规律,以实测地表形变作为约束,反演得到近距离多煤层的几何形状、位置以及产状等参数信息,对于我国普遍分布的煤采区,尤其是大量已关闭或废弃煤矿的地面沉陷范围圈定以及多煤层采空位置、产状等信息获取具有重要借鉴意义,进而可为煤采区,特别是已关闭或废弃工矿区地面沉降监测和综合治理提供依据.

1 近距离多煤层参数反演方法

近距离多煤层共同开采是目前全球范围内各矿场在实际开采作业过程中所面临的一种复杂情况.研究其开采规律,并根据实际监测的地表形变快速准确反演地下复杂的多煤层几何形状、位置、产状等参数对于煤矿采空区的综合治理极其重要.本文基于Okada模型[11],结合近距离多煤层分布的实际特征,依托叠加原理,提出了一种基于多源模型的近距离多煤层参数反演方法,其示意图如图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)和后验概率分布的计算公式如下:

G(t)=b1G(t)1+b2G(t)2+…+bsG(t)s,

(1)

(2)

其中:G(t)s和bs分别表示多源模型中不同源模拟得到的形变和所占权重;s为多源模型的数量;d为M维观测数据;t为N维模型参数;p(d)为归一化常数;p(t)为多源模型参数的先验信息;p(t/d)为在已知先验信息下当前参数的后验概率;p(d/t)为多源模型的似然函数(likelihood function).

本文依据模型参数拟合值的分布特征,选择多维高斯表达式L(t)来表示似然函数,如式(3)所示,其中Cd表示数据协方差矩阵.

L(t)=

(3)

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

π(i)=P(i1=qi),i=1,2,…,N.

(4)

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

A=[aij]N×N.

(5)

aij=P(it+1=qi|it=qi),i,j=1,2,…,N.

(6)

其中：aij是在t时刻处于状态qi的条件下向t+1时刻转移到状态qj的概率.其中i表示状态序列,q表示可能的状态.

最后计算所有样本多源模型反演参数的后验概率,作为多煤层最优参数的选择指标.在寻找最优参数的过程中,选择后验概率最大的一组多源模型参数即为最佳拟合参数.

2 研究区InSAR形变监测

2.1 研究区与数据源

康平煤田位于辽宁省沈阳市康平县,属于地下开采的隐蔽型多煤层联合开采煤田(图2).其中大平煤矿辖区面积28.61 km2,年生产能力330万t,共划分为7个采区,其中有5个采区处于三台子水库下;小康煤矿辖区面积为30.08 km2,被划分为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 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 SAR数据处理流程图Fig.3 Flow chart of the SAR data processing

本文在获取研究区形变结果时,基于研究区近距离多煤层开采的地表实际覆盖情况,充分考虑时空基线和多普勒效应的影响,设置空间基线阈值为临界基线的2%;时间基线阈值为72 d,生成了216组干涉像对(图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;除此之外,其他区域沉降现象不明显,这表明该区域整体较为稳定,地面沉降主要是由煤炭开采所致.

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

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

图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 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 Comparison of the observed deformation,simulated deformation of dual-source model and residual results(a)—InSAR观测图; (b)—双源模型模拟图; (c)—双源模型残差图.

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

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

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

(7)

(8)

计算双源模型的RMSE为15.9 mm,R2为0.399 0;三源模型的RMSE为6.7 mm,R2为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,R2达到了0.791 2.结果表明三源模型反演精度更高,更符合小康矿近距离多煤层共同采动的实际特征.