基于SBAS-InSAR技术和Logistic模型的矿区沉降动态预测方法

doi:10.6046/zrzyyg.2021354

基于SBAS-InSAR技术和Logistic模型的矿区沉降动态预测方法

徐子兴^,¹, 季民^,¹, 张过², 陈振炜²

1.山东科技大学测绘与空间信息学院,青岛 266590

2.武汉大学测绘遥感信息工程国家重点实验室,武汉 430079

Method for dynamic prediction of mining subsidence based on the SBAS-InSAR technology and the logistic model

XU Zixing^,¹, JI Min^,¹, ZHANG Guo², CHEN Zhenwei²

1. College of Geodesy and Geomatics, Shandong University of Science and Technology, Qingdao 266590, China

2. National Key Laboratory of Surveying and Remote Sensing Information Engineering, Wuhan University, Wuhan 430079, China

通讯作者: 季民(1970-),男,教授,研究方向为合成孔径雷达干涉测量与应用等。Email:jamesjimin@126.com。

责任编辑: 李瑜

收稿日期: 2021-10-25 修回日期: 2022-03-12

基金资助:

国家自然科学基金项目“X波段双天线SAR卫星干涉应用方法与设计指标研究”(418013973)
山东省重大科技创新工程项目“省级自然资源监测监管大数据应用服务平台建设”(2019JZZY020103)

Received: 2021-10-25 Revised: 2022-03-12

作者简介 About authors

徐子兴(1997-),男,硕士,研究方向为合成孔径雷达干涉测量与应用。Email: zixingxu@whu.edu.cn。

摘要

根据矿区开采沉降规律对后续沉降进行预测是评估矿山开采风险、调整开采规划的关键。对使用Logistic时间函数模型进行矿区沉降预测中的可用条件进行了分析和模拟实验,并提出了一种基于小基线集合成孔径雷达干涉测量(SBAS-InSAR)技术和Logistic模型的矿区沉降动态预测方法。首先,通过SBAS-InSAR获得矿区的时序沉降数据; 然后,以时序沉降数据作为拟合数据,采用信赖域算法逐像元计算其Logistic模型参数,根据Logistic模型可用条件,确定出可以对后续沉降进行预测的像元范围; 最后,根据Logistic模型对可预测范围内的后续沉降进行预测。以内蒙古鄂尔多斯市某矿区为研究区对上述预测方法进行了实验,采用对应日期的InSAR监测结果对预测结果进行了验证,结果表明: 36 d后和108 d后预测结果的均方根误差分别为0.010 1 m和0.023 6 m,预测误差小于0.03 m的比例分别达到98.9%和89.3%,表明该动态预测模型预测精度较高。

关键词： SBAS-InSAR; Logistic模型; 煤矿开采沉降; 动态预测

Abstract

Predicting the subsequent subsidence in mining areas according to the law of mining subsidence is the key to assessing mining risks and adjusting mining planning. This study determined the available conditions of the logistic model for mining subsidence prediction through analysis and simulation experiments and proposed a method for the dynamic prediction of mining subsidence based on small baseline subset (SBAS)-interferometric synthetic aperture radar (InSAR) technology and the logistic model. Firstly, the time-series subsidence data of a mining area was obtained using the SBAS-InSAR technology. Then, taking the time series subsidence data as the data for fitting, the parameters of the logistic model were calculated pixel by pixel by using the trust region algorithm. Then, the pixel range in which the subsequent subsidence can be predicted was determined according to the available conditions of the logistic model. Finally, according to the Logistic model, the subsequent subsidence within the predictable range was predicted. This method was applied to a certain mining area in Erdos City, Inner Mongolia for tests, and the prediction results were verified using the InSAR monitoring results of corresponding dates. The predicted results after 36 days and 108 days of ming had the root mean square error (RMSE) of 0.010 1 m and 0.023 6 m, respectively, and their proportion with prediction errors of less than 0.03 m reached 98.9% and 89.3%, respectively. These results indicate that the method for dynamic prediction proposed in this study has high prediction accuracy.

Keywords： SBAS-InSAR; Logistic model; coal mining subsidence; dynamic prediction

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

本文引用格式

徐子兴, 季民, 张过, 陈振炜. 基于SBAS-InSAR技术和Logistic模型的矿区沉降动态预测方法[J]. 自然资源遥感, 2022, 34(2): 20-29 doi:10.6046/zrzyyg.2021354

XU Zixing, JI Min, ZHANG Guo, CHEN Zhenwei. Method for dynamic prediction of mining subsidence based on the SBAS-InSAR technology and the logistic model[J]. Remote Sensing for Land & Resources, 2022, 34(2): 20-29 doi:10.6046/zrzyyg.2021354

0 引言

地下采矿活动会破坏开采工作面上覆岩层的原始应力平衡,从而导致矿区地表发生水平移动和下沉,不仅会对地面建筑和生产设施造成破坏,还可能诱发地质灾害,给矿区居民的生产和生活留下安全隐患^[1]。地下开采导致的地表沉降不仅与煤层的水文地质条件有关,而且也与开采方法、开采速度等密切相关。利用矿山开采沉降规律对地下开采活动引起的后续地表沉降进行有效预测,是评估矿山开采风险并进一步采取防护措施的关键。

以往基于水准或GPS测量的沉降仅可对少量离散点进行沉降预测,而合成孔径雷达干涉测量(synthetic aperture Radar interferometry,InSAR)技术的发展和应用为大范围的沉降预测提供了数据基础。InSAR技术凭借其大范围监测、高空间分辨率、穿云透雾等优点被广泛应用于地表形变的监测中^[2]。矿区地表通常有植被覆盖,并且容易发生大梯度形变,这使得矿区容易发生失相干现象^[3]。而小基线集干涉测量技术(small baseline subset-InSAR,SBAS-InSAR)作为一种常用的InSAR形变监测技术,可以有效减弱时空失相干和大气效应的影响,在矿区沉降监测中得到了广泛的研究与应用^[4⇓⇓-7]。

矿区开采沉降过程大致可分为3个阶段: 初始沉降期、主要沉降期和残余沉降期,沉降增长呈“S”型^[8]。描述该过程的时间函数模型主要有Knothe模型^[9⇓-11]、Richards模型^[12-13]、正态分布时间函数^[14-15]及Logistic模型^[16-17]等。其中,Logistic模型作为典型的“S”型增长曲线,与矿山开采沉降的3个阶段比较符合。张文志等^[16]用Logistic模型曲线拟合开采沉陷-时间关系曲线,并进行了单点开采沉陷预测,结果表明该模型应用于矿区沉降拟合和预测具有较高的精度。杨泽发等^[17]通过将合成孔径雷达差分干涉测量(differential interferometry synthetic aperture Radar,D-InSAR)监测到的矿区沉降累加,结合Logistic模型对矿区全盆地沉降规律进行了分析,并对全盆地的后续沉降进行了预测。但仅当整个沉降盆地的沉降过程接近结束时,全盆地的所有点才满足可对后续沉降进行预测的条件。在煤炭开采过程中,整个下沉盆地会有部分点满足预测条件、部分点不满足预测条件,上述不考虑可用条件的预测方法并不适用。

本文对使用Logistic时间函数模型进行矿区沉降预测中的可用条件进行了分析和模拟实验,并提出了一种基于小基线集InSAR(SBAS-InSAR)技术和Logistic时间函数模型的矿区沉降动态预测方法。首先,采用SBAS-InSAR技术对矿山开采引起的沉降进行监测,获取矿区沉降时间序列; 然后,以时序沉降数据作为拟合数据,采用信赖域算法逐像元计算Logistic模型参数,并根据Logistic模型的可用条件判断出可对后续沉降进行预测的像元范围; 最后,根据Logistic模型计算对可预测范围内的所有像元的后续沉降进行预测。

1 方法与原理

1.1 Logistics模型

选择Logistic模型以确定矿区沉降动态预测中的可预测条件和预测矿区沉降,这是因为其对于矿区沉降具有良好的拟合精度,且具有最大加速度点、最大速度点和最小加速度点3个特征点,利于沉降阶段的划分。Logistic模型是经典的“S”型增长曲线,其公式为:

(1)W(t)=

\frac{W_{0}}{1 + a e^{- b t}}

式中: W(t)为某一点t时刻的下沉值; W₀为该点最大下沉值; a和b为形状参数; 待求参数为W₀,a和b。对式(1)分别求一阶导数和二阶导数,可得到Logistic模型下沉速度和下沉加速度的表达式,分别为:

(2)v(t)=

\frac{d W (t)}{d t}

\frac{a b W_{0}}{e^{b t} (a e^{- b t} {+ 1)}^{2}}

(3)a(t)=

\frac{d^{2} W (t)}{d t^{2}}

\frac{2 a^{2} b^{2} W_{0}}{e^{2 b t} (a e^{- b t} {+ 1)}^{3}}

\frac{a b^{2} W_{0}}{e^{b t} (a e^{- b t} {+ 1)}^{2}}

。

Logistic模型具有3个特征点,分别为最大加速度点、最大速度点和最小加速度点。根据这3个特征点设置了拟合数据对沉降过程的不同覆盖程度,并探究其对Logistic模型拟合精度的影响。假定一个W₀,a和b参数分别为-1 m,120和0.04的Logistic模型,时间范围设置为0~240 d,每12 d获取一个沉降值得到该Logistic模型下的沉降时间序列。为了使模拟的数据符合实际情况,在沉降时间序列中加入了0.02 m的正态分布随机误差。根据拟合数据对沉降过程的不同覆盖程度设置了4组模拟实验,采用信赖域算法进行拟合得到Logistic曲线,如图1所示。图上3个红色三角形自左到右分别代表最大加速度点、最大速度点和最小加速度点。

图1

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

图1 参与拟合数据对沉降过程的不同覆盖范围对拟合正确性的影响

Fig.1 The influence of different coverage of the fitting data on the subsidence process on the fitting accuracy

从图1中可以看出,随着拟合数据覆盖范围的增加,拟合曲线越来越逼近正确曲线。直到拟合数据覆盖范围超过最右边特征点(最小加速度点)时,拟合曲线才非常接近正确曲线。因此,将拟合数据覆盖范围是否超过最小加速度点作为Logistic模型的可用条件,并根据该可用条件逐像元进行判断得到可对后续沉降进行预测的范围。

1.2 SBAS-InSAR矿区沉降监测

采用SBAS-InSAR对矿山开采引起的沉降进行监测,SBAS-InSAR可以有效减弱常规InSAR方法中易受大气延迟误差干扰、时空失相干等问题。假设有N+1幅覆盖研究区域的SAR影像时,选取其中一幅作为超级主影像,将其余N幅SAR影像都配准到超级主影像上,根据小基线原则和矿区形变速度快、形变量级大的特点设定时空基线阈值进行基线组合,得到M个干涉对。M满足:

(4)(N+1)/2<M<N(N+1)/2,

其中任一干涉对主影像和从影像获取时间分别设为t_A和t_B(t_B> t_A),对干涉对进行干涉、去平地效应和去地形相位后,得到的第j幅干涉图上高相干点(x,y)上的相位值可以表示为^[18]:

(5)

\begin{array}{l} δ φ_{j} (x, y) = φ_{B} (x, y) - φ_{A} (x, y) \approx \frac{4 π}{λ} [d (t_{B}, x, y) - d (t_{A}, x, y)] \\ + Δ φ_{j}^{t o p o} (x, y) + Δ φ_{j}^{o r b} (x, y) + Δ φ_{j}^{a t m} (x, y) + Δ φ_{j}^{n o i s e} (x, y) \end{array}

式中: λ为波长; d(t_B,x,y)和d(t_A,x,y)为相对于d(t₀,x,y)=0时的雷达视线方向的累积形变量; Δ ${φ_{j}}^{t o p o}$ (x,y),Δ ${φ_{j}}^{a t m}$ (x,y),Δ ${φ_{j}}^{o r b}$ (x,y)和Δ ${φ_{j}}^{n o i s e}$ (x,y)分别表示残余地形、轨道信息不精确、大气延迟和噪声导致的相位误差,去除这些误差相位后,式(5)可简化为:

(6)δφ_j(x,y)≈

\frac{4 π}{λ}

[d(t_B,x,y)-d(t_A,x,y)],

(7)d(t_B,x,y)-d(t_A,x,y)=v_i(t_B-t_A),

式中: v_i为t_A至t_B时间段内雷达视线方向的平均形变速率。将式(6)和式(7)转化为矩阵形式为:

(8)BV=δφ,

式中: B为参数矩阵; V为平均形变速率矩阵; δφ为相位矩阵。采用最小二乘法和奇异值分解法对式(8)求解,得到各时间段内平均形变速率的最小范数解,并对其进行积分得到视线(line of sight,LOS)方向的时间序列形变量^[19]。

利用式(9)将LOS方向形变转换为沉降值,公式为:

(9)W=d_LOS/cosθ,

式中: W为沉降值; d_LOS为LOS方向形变; θ为入射角。

1.3 信赖域算法求解Logistic模型参数

采用信赖域算法逐像元将SBAS-InSAR监测的时序沉降进行Logistic模型拟合,计算其最优Logistic模型参数,即求解式(10):

(10)minf(x)=

\overset{N}{\sum_{i = 1}}

(W_InSAR-

\frac{W_{0}}{1 + a e x p (- b t)}

设x_k为第k次迭代点, $f_{k}_{}$ = f(x_k), g_k = g(x_k) = ∇ f (x_k),B_k为∇ ²f (x_k)的第k次近似,则第k次迭代步信赖域二次模型子问题如式(11)所示^[20]:

(11)

\begin{array}{l} m i n φ_{k} (d) = f_{k} + {g^{T}}_{k} d + \frac{1}{2} d^{T} B_{k} d \\ s . t . | | d | | \leq Δ_{k} \end{array}

式中: Δ_k >0为信赖域半径。设子问题的解为d_k,通过式(12)反映二次函数模型与目标函数之间的近似程度,公式为:

(12)r_k=

\frac{f (x_{k}) - f (x_{k} + d_{k})}{f (x_{k}) - φ_{k} (d_{k})}

若 ${r_{k}}_{}$ < 0,则 ${x_{k}}_{}$ + d_k不能作为下一个迭代点,此时需要缩小信赖域半径; 若 ${r_{k}}_{}$ 接近于1,可以认为模型函数与目标函数在信赖域范围内具有良好的一致性,可以扩大信赖域半径; 反之,若r_k离1 较远,则需要缩小信赖域半径。不断迭代直至得到收敛解,流程图如图2所示。

图2

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

图2 信赖域算法计算Logistic模型参数流程图

Fig.2 Flow chart of calculating Logistic model parameters by trust region algorithm

1.4 矿区沉降动态预测方法

本预测方法首先采用SBAS-InSAR对开采沉降进行监测,获取沉降时间序列; 然后以沉降时间序列作为拟合数据,采用信赖域算法逐像元进行Logistic模型求参,并根据可预测条件判断可预测范围; 最后利用Logistic模型对可预测范围内的像元进行沉降预测。流程图如图3所示。

图3

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

图3 矿区沉降动态预测方法流程图

Fig.3 Flow chart of subsidence dynamic prediction method in mining area

2 矿区沉降预测实验与结果

2.1 研究区与数据源

所选研究矿区为内蒙古某矿区,位于鄂尔多斯市东南部,地处毛乌素沙漠的东北边缘,与陕西省榆林市毗邻,如图4所示。研究矿区表面植被覆盖率较低,受时间去相干影响较小。所处位置地势较为平坦,平均海拔为1 310 m。利用40景IW模式、VV极化的Sentinel-1A影像进行SBAS-InSAR处理,时间跨度为2017年9月20日—2019年2月18日,影像参数如表1所示。其中2017年9月20日—2018年11月02日的时序沉降用于作为拟合数据计算Logistic模型参数,并分别对2018年12月08日和2019年02月18日的累积沉降进行预测,采用对应日期的InSAR监测累积沉降对预测结果进行了验证。相邻影像时间间隔大部分是12 d,少量为24 d,小的时间间隔可以最小化大部分去相干源的影响。采用30 m分辨率的航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission,SRTM)数据用于模拟和去除地形相位。

图4

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

图4 研究矿区位置

Fig.4 Location of the study mining area

表1 Sentinel-1A影像参数

Tab.1 Parameters of Sentinel-1A images

序号	日期	入射角/(°)	序号	日期	入射角/(°)	序号	日期	入射角/(°)
1	2017/9/20	33.914	15	2018/3/19	33.910	29	2018/9/27	33.915
2	2017/10/2	33.913	16	2018/3/31	33.912	30	2018/10/9	33.912
3	2017/10/14	33.912	17	2018/4/12	33.917	31	2018/10/21	33.912
4	2017/10/26	33.910	18	2018/5/6	33.915	32	2018/11/2	33.910
5	2017/11/7	33.909	19	2018/5/18	33.915	33	2018/11/14	33.908
6	2017/11/19	33.909	20	2018/5/30	33.916	34	2018/11/26	33.908
7	2017/12/1	33.908	21	2018/6/11	33.915	35	2018/12/8	33.909
8	2017/12/13	33.908	22	2018/6/23	33.916	36	2018/12/20	33.906
9	2017/12/25	33.907	23	2018/7/5	33.915	37	2019/1/1	33.907
10	2018/1/6	33.907	24	2018/7/17	33.917	38	2019/1/25	33.908
11	2018/1/30	33.909	25	2018/7/29	33.915	39	2019/2/6	33.908
12	2018/2/11	33.909	26	2018/8/22	33.915	40	2019/2/18	33.909
13	2018/2/23	33.909	27	2018/9/3	33.915	轨道号	Path: 11 Frame: 121
14	2018/3/7	33.910	28	2018/9/15	33.917	极化方式	VV

新窗口打开| 下载CSV

2.2 SBAS-InSAR监测

根据小基线集的原则和矿区大梯度形变的特点,设置空间基线和时间基线的阈值分别为150 m和36 d。经过配准、差分干涉、滤波、相位解缠、基线精化、估算线性形变速率、去除大气延迟误差、估算最终形变速率、相位转形变和地理编码,并将形变转到垂直方向后,得到了该矿区2017年9月20日—2019年2月18日的39幅时序沉降图,部分如图5所示。

图5-1

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

图5-1 时序沉降图

Fig.5-1 Time series subsidence

图5-2

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

图5-2 时序沉降图

Fig.5-2 Time series subsidence

图5(h)中的A—D这4个点于下节中展示Logistic模型拟合情况。

2.3 Logistic模型拟合

以SBAS-InSAR处理得到的2017年9月20日—2018年11月02日的时序沉降数据作为拟合数据,对最终沉降大于0.05 m的像元逐像元采用信赖域算法计算Logistic模型参数。图6展示图5(h)中的A,B,C,D这4个点的Logistic模型拟合情况。

图6

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

图6 Logistic模型拟合情况

Fig.6 Logistic model fitting situation

2.4 沉降预测结果

根据拟合得到的Logistic模型的参数,以参与拟合数据的覆盖范围是否超过最小加速度点作为判断依据逐像元进行判断,得到可预测范围。分别对2017年9月20日—2018年12月08日和2019年02月18日的累积沉降进行预测,如图7所示。

图7

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

图7 预测沉降图

Fig.7 Prediction subsidence diagram

2.5 预测精度验证

采用SBAS-InSAR获得的2017年9月20日—2018年12月08日和2019年02月18日的累积沉降对预测结果进行验证。将InSAR监测值和预测值作一致性验证图,如图8所示,可以看出大部分点分布在y=x附近,两者一致性较高。可预测范围内全部像元的预测沉降值与InSAR监测沉降值的RMSE分别为0.010 1 m和0.023 6 m。统计预测误差作柱状图如图9所示,预测误差小于0.03 m的比例分别达到98.9%和89.3%,表明该动态预测模型预测精度较高。

图8

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

图8 一致性验证图

Fig.8 Consistency verification diagram

图9

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

图9 预测误差统计柱状图

Fig.9 Statistical bar graph of prediction error

2.6 预测精度影响因素分析

预测误差的空间分布如图10所示,从图中可以看出: 垂直于工作面开采方向,预测误差较大的像元基本处于沉降盆地的中间位置,因为沉降盆地中间位置下沉量大、下沉速度快,容易导致大的预测误差; 沿着工作面开采方向,预测误差较大的像元大部分处于开采方向末端,分析原因是开采方向末端像元拟合数据对最小加速度点后下沉过程的覆盖程度较低。将最小加速度点出现时间作为x轴,预测误差作为y轴,绘制散点图,如图11所示。从图11可以看出,最小加速度点出现时间越晚,出现较大预测误差的概率也随之增加。因为拟合数据的时间范围是一个定值,最小加速度点出现时间越晚,最小加速度点之后的拟合数据就越少,对最小加速度点后下沉过程的覆盖程度就会越低,会对后续沉降的预测精度产生影响。

图10

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

图10 预测误差空间分布

Fig.10 Spatial distribution of prediction error

图11

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

图11 最小加速度点出现时间与预测误差之间的关系

Fig.11 Relationship between the occurrence time of minimum acceleration point and prediction error

3 结论

本文对Logistic模型在矿区沉降动态预测中的可用条件进行分析与模拟实验,并基于SBAS-InSAR技术和Logistic模型提出了一种新的矿区沉降动态预测方法。以内蒙古某矿区为实验区,利用2017年9月20日—2019年2月18日的40景Sentinel-1A数据进行了矿区开采沉降预测实验,结果显示36 d后和108 d后预测结果的RMSE分别为0.010 1 m和0.023 6 m,预测误差小于0.03 m的比例分别达到98.9%和89.3%,表明该预测方法精度较高。该方法具有一定的实用性,对评估矿山开采风险、指导矿区开采规划具有一定的理论借鉴价值和实际应用意义。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

何国清, 杨伦, 凌赓娣, 等.

矿山开采沉陷学[M]. 徐州: 中国矿业大学出版社, 1991.