融合D-InSAR与幂指数Knothe时间函数的老采空区地表残余变形预测

Prediction of Residual Surface Deformation in Old Goafs by Combining D-InSAR Technology and Power Exponential Knothe Function Model

摘要: 煤炭老采空区在我国广泛分布，随着老城改造和城镇化发展，越来越多基础设施不可避免需建造在老采空区，而老采空区地表残余变形预测既是这些工程建设的迫切需求，也是面临的重要难题。本文以河南平顶山某拟建变电站所在老采空区场地为研究区域，融合D-InSAR技术与幂指数Knothe时间函数模型对老采空区地表残余变形进行预测。首先利用合成孔径雷达差分干涉测量(D-InSAR)技术分析得到研究区域地表2017年6月~2022年6月的累计变形及逐年沉降量。在此基础上，根据地表累计变形与逐年沉降数据，确定了幂指数Knothe时间函数模型参数取值，建立了描述研究区域老采空区地表动态沉降规律的幂指数Knothe时间函数模型，并基于实测数据对该模型及参数进行验证。最后对变电站所在场地老采空区地表剩余残余变形及未来10年的逐年沉降量开展预测，研究成果可为老采空区工程建设提供依据。

Abstract: Old goafs are widely distributed in China. With the transformation of old cities and the development of urbanization, more and more infrastructure inevitably needs to be built in old goafs. However, the prediction of residual deformation of the surface in old goafs has always been an urgent demand and important challenge for these engineering constructions. In this study, an old goaf site in Pingdingshan, Henan Province is chosen as the research area, at which a substation is proposed to be built. The residual deformation of the old goaf surface is predicted by combining the D-InSAR technology and the power exponent Knothe function. Firstly, the surface deformation data of the research area from June 2017 to June 2022 were extracted using Synthetic Aperture Radar Differential Interferometry (D-InSAR) technology. As a result, the cumulative surface deformation and annual settlement of the old goaf during this period were obtained. Based on these cumulative surface deformation and annual settlement data, the parameters of the power exponent Knothe time function model were determined. Then, a power exponent Knothe function model was established to describe the dynamic settlement law of the old goaf surface in this research area. The power exponent Knothe function model and parameters were further verified based on these measured data. Finally, the residual deformation of the old goaf surface and the annual settlement within the next 10 years of the substation site are predicted. The results in this study can provide a basis for engineering construction.

1. 引言

我国当前煤炭资源开采主要还是采用井工方式，占比达96%左右，这种开采方式容易造成煤层顶板围岩失稳而产生位移、开裂、破碎跨落，导致上覆岩层整体下沉、弯曲并引起地表变形与破坏，进而形成大规模采空区和大面积地表沉陷 [1] [2] [3] 。

煤炭开采引起地表沉降变形在经初始期、活跃期和衰退期三个特征阶段后，就进入了煤矿终采后地表变形移动的残余下沉期 [4] 。煤矿终采后残余下沉期，采空区冒落破碎岩石、采动离层、裂隙等在上覆岩层载荷持续作用下被逐步压实压密 [5] ，地表缓慢下沉而引起的采空区地表残余变形一般会持续几年甚至几十年 [6] 。目前，我国对煤炭开采引起地表沉降变形监测预测已经形成了较为成熟的理论和方法，如可采用水准测量和GPS测量 [7] 、雷达干涉测量 [8] - [13] 、瞬变电磁法 [14] [15] 等技术与手段获取采空区沉降数据。此外，建立沉降模型对地表沉降全过程进行动态预测也是获得采空区地表沉降的一种重要方法 [16] 。常用的采空区地表沉降预测数学有双曲线模型 [17] 、Knothe时间函数模型 [18] [19] [20] [21] 、幂指数Knothe时间函数模型 [22] [23] [24] 等。

当前，随着我国经济社会发展，越来越多基础设施不可避免需建造在煤矿采空区，特别是煤矿老采空区 [25] [26] [27] [28] 。尽管这些老采空区一般已停止采掘3年及以上且由煤炭开采所引起的地表移动变形业已趋于稳定，但老采空区所建设的基础设施服务年限长，而采空区地表残余变形也需要经历一个漫长的时间过程。因此，预测老采空区地表残余变形量及其随时间的变化规律，是老采空区基础设施建设的迫切需求。

然而，多数情况下对老采空而言，由于开采终止时间较为久远，往往也缺乏采空区地表变形监测资料。同时，采空区煤炭开采的具体时间、矿层分布及其层数、厚度、深度、倾角等关键特征信息也往往不全甚至缺失，采空区场地的地质与水文条件、地层分布、岩土体性质等资料不完整，这些都给老采空区地表残余变形计算与预测带来了巨大困难。

本文以平顶山某拟建变电站所在场地老采空区为研究区域，利用合成孔径雷达差分干涉测量(D-InSAR)技术提取了研究区域2017年6月~2022年6月六景D-InSAR图像，分析获得了该时段研究区域的地表累计变形及逐年沉降量。在此基础上，根据地表变形数据求解得到了幂指数Knothe时间函数模型参数，从而建立描述研究区域老采空区地表动态沉降规律的幂指数Knothe时间函数模型，进一步基于实测数据对该模型及参数进行了验证，据此对该老采空区地表剩余残余变形及未来10年的逐年沉降量开展预测，研究成果可为变电站工程建设提供依据。

2. 研究区域概况及D-InSAR形变提取

2.1. 研究区域概况

本文选择河南平顶山地区某110千伏变电站所在煤矿老采空区为研究区域，该变电站场地老采空区位于平煤七矿，地势平坦，地层由新至老分别为第四系、二叠系、石炭系和寒武系，地质性质较差、环境复杂。开采煤层为己组，煤层开采深度约160 m，煤层厚度约3.5 m，开采时间为1987~1993年。

2.2. 研究区域D-InSAR形变提取与处理

本文采用的数据源为平顶山市区Sentinel-1A IW SLC图像(图1所示)，极化方式为VV极化，干涉宽幅模式为IW，波段类型为C波段，辅助12.5m分辨率DEM数据。

本研究过程中，选取了2017年6月~2022年6月的六景Sentinel-1A图像，将拟提取形变的两幅不同时间点图像作为像对，利用像对所生成的干涉图进行去平效应，去平后的干涉图可消除由于雷达侧视成像所造成地面上两个高度相同、与卫星距离远近不同的物体在水平方向上所产生的相位差而导致干涉条纹过密及增加相位解缠难度的问题。此后利用Goldstein滤波法去除干涉图的斑点噪声，经过解缠计算后，得到相位解缠结果及用于估算轨道精炼的修正参数的控制点位置。

对轨道进行精炼和重去平后，将得到的相位差转为形变数据及地理编码，即可得到研究区域地表形变。在此基础上，利用两形变对的成像几何关系及其与DEM数据结合模拟地形相位，用形变对相位减去地形对相位，即可得到该区域在该时间内地表沉降位。形变相位图转化为拟建变电站所在场地2017年6月~2022年6月间地表形变情况如图2所示。

进一步，可得到研究区域采空区2017年6月~2022年6月间地表逐年沉降量和累计变形如图3所示。结果表明，该时间段内研究区域老采空区地表沉降总体呈逐年减小趋势。

Figure 1. Cropping range of SAR image in this study

图1. 本文研究选取的SAR图像范围

Figure 2. The proposed substation site’s annual surface subsidence from June 2017 to June 2022

图2. 拟建变电站场地2017年6月~2022年6月间逐年地表形变量

Figure 3. Surface deformation from June 2017 to June 2022 extracted by D-InSAR

图3. 2017年6月~2022年6月间研究区域地表沉降情况

3. 基于D-InSAR变形的幂指数Knothe时间函数模型

波兰学者Knothe (1952) [17] 假设采空区地表上某点下沉速度与该点最终下沉量和该点在某时刻动态下沉量之差成正比，即：

$\frac{\text{d}W\left(t\right)}{\text{d}t}=c\left[W_{cm}-W\left(t\right)\right]$ (1)

式中：W(t)为自沉降开始起某时刻的地表累计沉降(mm)；t为自沉降开始的时间(a)；W_cm为最终沉降量(mm)；c为时间系数。

对式(1)求解得到t时刻累计沉降量表达式：

$W\left(t\right)=W_{cm}\left(1-{\text{e}}^{-ct}\right)$ (2)

当已知某两个时间点t₁和t₂沉降量差值 $\Delta W$ ，即可根据式(2)建立方程：

$\Delta W=W_{cm}\left[{\text{e}}^{-t_{2}c}-{\text{e}}^{-t_{1}c}\right]$ (3)

因此，利用本文前述研究区域采空区地表2017年6月~2022年6月间逐年沉降数据，基于式(3)可建立求解Knothe时间函数模型参数c和W_cm方程组，继而求解确定相应模型参数。

按照式(1)所示Knothe时间函数沉降模型，老采空区地表开始沉降时，地表变形速度最快。随着沉降时间的发展，地表沉降速率将不断减慢直至趋于稳定。但大量的已有监测资料表明 [4] [6] [28] ，在煤炭开采引起地表沉降过程中，地表下沉速率在初始阶段较为缓慢，其后不断增加直至达最大值，随后将不断减慢至趋于稳定。为更合理描述煤矿采空区地表沉降规律，刘玉成等(2009) [29] 在Knothe时间函数基础上增加幂指数常数k，提出了如式(4)所示的幂指数Knothe时间函数：

$W\left(t\right)=W_{cm}{\left(1-{\text{e}}^{-ct}\right)}^k$ (4)

本文研究区域的煤层开采时间为1987~1993年，由此假设式(1)和式(4)中初始时间起点为1987年，即t = 0时刻。如图3所示，2017年6月~2019年6月地表形变量10.5 mm，2019年6月~2022年6月地表形变量6.4 mm。由此可根据式(3)并利用Python + scipy.optimize科学计算库，基于最小二乘法拟合计算方法，求解得Knothe时间函数模型参数c = 0.23。

进一步地，如图3所示，研究区域采空区地表2017年6月~2022年6月间地表累计沉降量为16.9 mm，由此可由式(4)得到：

$W_{cm}{\left(1-{\text{e}}^{-0.23\times 35}\right)}^k-W_{cm}{\left(1-{\text{e}}^{-0.23\times 30}\right)}^k=16.9$ (5)

此外，对式(4)进行进一步求导，得到采空区地表的沉降速度表达式：

$v\left(t\right)=W_{cm}kc{\text{e}}^{-ct}{\left(1-{\text{e}}^{-ct}\right)}^{k-1}$ (6)

根据图3所示2017年6月~2022年6月间研究区域地表沉降情况，可计算得到研究区域在2019年的沉降速率为3.8 mm/a，即：

$\bar{v}\left(2019\right)=0.23W_{cm}k{\text{e}}^{-0.23\times 32}{\left(1-{\text{e}}^{-0.23\times 32}\right)}^{k-1}=3.8$ (7)

联立(5)式和(8)式，可求解得：k = 3.77，W_cm = 6508.7 mm。

至此，得到描述本文研究区域老采空区地表动态沉降变形幂指数Knothe时间函数模型表达式：

$W\left(t\right)=6508.7{\left(1-{\text{e}}^{-0.23t}\right)}^{3.77}$ (8)

根据式(8)，计算得到该采空区场地自开采以来地表累计沉降量随时间变化曲线如图4所示。

Figure 4. Settlement based on power exponential Knothe function model

图4. 基于幂指数Knothe时间函数地表沉降

根据文献 [30] 监测资料，本文研究区域采空区在2014年的累计沉降量为6.2 m。按本文式(8)所示幂指数Knothe地表沉降模型计算得到2014年的场地累计沉降为6.45 m，二者误差约为4.0%。结果表明，本文所建立的老采空区地表沉降模型能够较好反映研究区域沉降规律。

4. 采空区地表剩余残余变形及预测

根据前述计算结果，研究区域采空区地表预计最大沉降变形为6508.7 mm，截止2022年6月已经发生实际沉降量为6501.2 mm。因此，研究区域采空区地表累计剩余残余变形预计为7.5 mm。

进一步地，基于式(4)计算得到自2022年6月~2032年6月共10年内研究区域采空区地表逐年沉降量及累计变形量，结果如图5所示。预测结果表明，研究区域采空区地表沉降呈逐年减小，直至趋于完全稳定。采空区地表剩余残余变形对实际变电站工程建设影响较小。

Figure 5. Prediction of settlement within the future ten years from 2023 to 2033

图5. 研究区域场地未来10年(2023~2032)沉降预测

5. 结束语

1) 本文通过D-InSAR技术将卫星遥感图像的干涉相位转换为形变相位，确定了研究区域煤矿老采空区在过去5年内地表逐年沉降及累计变形增量，进一步基于该时间段内采空区沉降变形数据，采用参数反演方法建立了描述研究区域老采空区地表动态沉降变形规律的幂指数Knothe时间函数模型，并基于实测数据进行了验证。本文所提出的这一融合D-InSAR与幂指数Knothe时间函数的老采空区地表移动变形及剩余残余变形预测方法，可为老采空区地表残余变形及其随时间变化规律预测提供参考。

2) 本文融合D-InSAR与幂指数Knothe时间函数的老采空区地表残余变形预测结果表明，研究区域地表未来10年剩余残余变形总量为7.5 mm，且地表沉降量逐年减小，并趋于稳定。采空区地表剩余残余变形对实际工程建设影响较小。

基金项目

本项目研究得到平顶山电力设计院有限公司科技项目(PDSSJY-2022050501)和国网河南省电力公司科技项目(SGHAPD00JJJS2000744)资助。

NOTES

^*通讯作者。