面向复杂矿区的Stacking技术辅助DS-InSAR地表形变监测方法

doi:10.6046/zrzyyg.2024104

面向复杂矿区的Stacking技术辅助DS-InSAR地表形变监测方法

李志^,¹, 张书毕^,¹, 李鸣庚², 陈强³, 卞和方¹, 李世金¹, 高延东¹, 张艳锁¹, 张帝¹

1.中国矿业大学环境与测绘学院,徐州 221116

2.中铁上海设计院集团有限公司,上海 200040

3.兖矿能源集团股份有限公司济宁三号煤矿,济宁 272000

Stacking-assisted DS-InSAR method for monitoring surface deformations in complex mining areas

LI Zhi^,¹, ZHANG Shubi^,¹, LI Minggeng², CHEN Qiang³, BIAN Hefang¹, LI Shijin¹, GAO Yandong¹, ZHANG Yansuo¹, ZHANG Di¹

1. School of Environment and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China

2. China Railway Shanghai Design Institute Group Corporation Limited, Shanghai 200040, China

3. Yankuang energy Group Company Limited Jining No.3 coal mine, Jining 272000, China

通讯作者: 张书毕(1965-),男,博士,教授,主要研究方向为GNSS数据处理与大气遥感。E-mail:zhangsbi@cumt.edu.cn。

责任编辑: 陈昊旻

收稿日期: 2024-03-21 修回日期: 2024-06-13

基金资助:

国家自然科学基金“高寒关键生态走廊滑坡机理与天空地协同智能监测预警方法”(U22A20569)
“融合GNSS和InSAR数据的动态节点基水汽层析理论与同化方法研究”(42271460)
“基于机器学习与多源数据融合的多基线相位解缠方法研究”(42001409)
中国博士后科学基金资助项目“基于机器学习与数据压缩的矿区DSInSAR数据处理方法研究”(2022M723376)
江苏省卓越博士后计划“‘学习型’理论与时间概率积分模型驱动下DSInSAR矿区监测研究”(2023ZB277)

Received: 2024-03-21 Revised: 2024-06-13

作者简介 About authors

李志(1999-),男,硕士研究生,主要研究方向为时序InSAR技术形变灾害监测。Email: liz@cumt.edu.cn。

摘要

合成孔径雷达干涉测量(interferometric synthetic aperture Radar, InSAR)技术在矿区植被覆盖密集且存在大梯度地表形变复杂环境下进行监测时,存在监测点数量不足、监测精度不高等问题。针对这些问题,该文提出一种Stacking技术辅助下的改进分布式目标InSAR(distributed scatterer InSAR, DS-InSAR)方法。该方法采用置信区间假设检验算法识别同质像元并基于相位三角剖分算法完成相位优化,随后去除先期利用Stacking技术模拟的线性形变相位获取残余相位,进而稀疏形变相位条纹,提高后续DS-InSAR处理框架中时空滤波与三维解缠结果的精确性,最终补偿模拟相位获得完整形变场。通过处理2015年10月—2016年3月期间覆盖新巨龙煤矿的Sentinel-1A合成孔径雷达(synthetic aperture Radar,SAR)影像,解译了该时段内矿区时序地表形变特征,得到以下结论: 监测期间,矿区存在3处显著形变,雷达视线向最大累积形变量达到-313 mm;所提方法相较常规短基线集干涉测量(small baseline subset InSAR,SBAS-InSAR)技术能够反演出分布更加均匀的监测点位信息,其密度约是SBAS-InSAR的12.9倍;对比水准数据的均方根误差(root mean squared error,RMSE)约为6.82 mm,精度较SBAS-InSAR提高了约3.0 mm。

关键词： Stacking; DS-InSAR; 地表形变; 残余相位; 矿区监测

Abstract

Interferometric Synthetic Aperture Radar (InSAR) faces the challenges of the insufficient number of monitoring points and low monitoring accuracy when applied to complex environments with dense vegetation and large-gradient surface deformation in a mining area. To address these challenges, this study proposed an improved distributed scatterer InSAR (DS-InSAR) method assisted by stacking technology. This method identified statistically homogenous pixels using a confidence interval hypothesis test and achieved phase optimization utilizing a phase triangulation algorithm. Subsequently, the residual phases were derived by removing the linear deformation phases determined via stacking-based simulation. This step contributed to sparse deformation phase fringes, thereby enhancing the accuracy of spatiotemporal filtering and three-dimensional phase unwrapping within the subsequent DS-InSAR processing framework. Finally, the simulated phases were compensated, and thus complete deformation fields were determined. By processing Sentinel-1A SAR images covering the Xinjulong Coal Mine from October 2015 to March 2016, this study interpreted the time-series surface deformation characteristics in the mining area during this period. The findings revealed three significant deformation sites in the mining area, with a maximum cumulative radar line-of-sight (LOS) deformation of up to -313 mm. Compared to conventional small Baseline Subset (SBAS) InSAR, the proposed method yielded more uniformly distributed monitoring points via inversion, with a density approximately 12.9 times higher. The root mean squared error (RMSE) of the inversion was approximately 6.82 mm relative to leveling data, representing an accuracy improvement of about 3.0 mm compared to the SBAS results.

Keywords： Stacking; DS-InSAR; surface deformation; residual phase; mining area deformation monitoring

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

本文引用格式

李志, 张书毕, 李鸣庚, 陈强, 卞和方, 李世金, 高延东, 张艳锁, 张帝. 面向复杂矿区的Stacking技术辅助DS-InSAR地表形变监测方法[J]. 自然资源遥感, 2025, 37(4): 12-20 doi:10.6046/zrzyyg.2024104

LI Zhi, ZHANG Shubi, LI Minggeng, CHEN Qiang, BIAN Hefang, LI Shijin, GAO Yandong, ZHANG Yansuo, ZHANG Di. Stacking-assisted DS-InSAR method for monitoring surface deformations in complex mining areas[J]. Remote Sensing for Land & Resources, 2025, 37(4): 12-20 doi:10.6046/zrzyyg.2024104

0 引言

煤炭资源一直占据我国能源结构主体地位^[1-2]。伴随着煤矿开采,地表塌陷、煤层滑坡和周边建筑与设施破坏等众多安全隐患不断增多^[3-4],对矿区及周边地区开展长时序地表形变监测与分析具有重要意义^[5]。合成孔径雷达干涉测量(interferometric synthetic aperture Radar, InSAR)技术因具备全天时、低成本等优势,有效弥补了传统形变监测方法(水准、全球定位系统测量等)的缺陷,已被广泛应用于大范围的对地观测场景。特别是近年来,星载合成孔径雷达(synthetic aperture Radar, SAR)系统持续蓬勃发展,填补了大量数据空缺^[6],极大地推动了InSAR技术的进步。

在形变监测领域,为克服差分InSAR(differential InSAR, D-InSAR)技术^[7]易受失相干及大气延迟等因素影响,Sandwell等^[8]提出了干涉图堆叠(Stacking InSAR)技术并获得了广泛应用; Xu等^[9]借助Stacking技术提取了大同煤矿的地表形变速率信息,并结合实验验证了该技术可以有效识别矿区形变量级与范围;吴琼等^[10]基于Stacking技术获取相位数据,协助深度学习方法识别出处于变形期的滑坡隐患位置,在广域地灾隐患识别中取得较好的应用成效。此外,各种基于时间序列分析的InSAR技术也陆续涌现,如永久散射体干涉测量(persistent scatterer InSAR, PSInSAR)技术^[11]、短基线集干涉测量(small baseline subset InSAR, SBAS-InSAR)技术^[12]和分布式目标InSAR(distributed scatterer InSAR, DS-InSAR)技术^[13]等。特别是DS-InSAR技术,在缺乏典型地物的裸地、荒原等低相干区域,能够有效提高监测点密度,进而提供更为丰富、可靠的形变细节。Zhao等^[14]将分布式目标(distvibuted scatterer,DS)与永久散射目标(persistent scatterer,PS)进行联合,并基于此开发了一种新的非线性DS相位估计方法,利用该方法对美国图森市西部地区进行长时序监测,所获监测点密度与精度较PSInSAR结果明显提升;蒋弥等^[15-16]针对DS-InSAR技术中的同质样本选取问题,提出了一种改进的最优参数统计同质样本选择方法,该方法在快速统计同质像元识别(fast statistically homogenous pixel selection, FaSHPS)算法基础上建立了更加精确的置信区间,进一步减少了异质像素;Wang等^[17]结合DS对湖北省大冶市西南矿群开展沉降监测,在大规模和长时序沉降监测后,监测数据与实际结果匹配良好,证明融合DS可以有效增强对矿区沉降情况的探测能力。

而在实际形变观测场景中,时序InSAR监测精度与相位解缠的效果密切相关。当目标区域存在大梯度变形时,干涉条纹会更加密集,从而降低相位解缠的精确性^[18]。高延东等^[19]针对地形起伏变化大、强噪声地区相位解缠问题,提出一种基于势函数改进权值的最大流/最小割相位解缠算法,为大梯度相位解缠提供了有效的思路; Feng等^[20]在D-InSAR技术基础上,通过移除偏移量追踪技术在估算出粗略形变场后进行解缠,从而达到降低相位梯度的目的,该方法可更准确地估算出物体表面位移速度。

综上所述,时序InSAR技术在矿区形变监测中已展现出较为理想的应用成效,然而其在矿区复杂环境下变形梯度过大时的监测能力尚需进一步深入研究。因此,本文改进了常规DS-InSAR处理路线,提出一种面向复杂矿区的Stacking技术辅助下的DS-InSAR方法,将山东省新巨龙煤矿作为研究区,获取该区域2015年10月至2016年3月时序地表形变信息并展开后续分析论证。

1 研究区概况及数据源

1.1 研究区概况

山东省菏泽市巨野县新巨龙煤矿(图1)位于巨野县西北部、巨野煤田中南部,地处村镇交界地带,交通便利。区域内为黄河流域冲洪积平原,平均地面高程约为43.3 m,地形相对平坦。地质构造以单斜构造为主,大致呈南北走向、向东分布,并伴有一定数量的断层。矿井类属于华北型石炭-二叠纪煤田,该类煤田是我国主要的原煤产地。煤系以中、下奥陶系石灰岩作为基底,上覆地层主要为新近系和第四系,含煤地层主要为石炭-二叠系山西组和太原组,总共含煤24层。现主采山西组3煤层,煤层走向为343°~18°,倾向为73°~108°,该煤层平均厚度8.8 m,现已被鉴定为弱冲击倾向性煤层,其地面标高43.9 m,井下标高-830.8~-797.5 m,结构环境复杂。矿区及其周边地区覆盖耕地、湿地和人造地表等7种地物,其中耕地面积占比达85 %。该矿井于2009年正式启动生产。

图1

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

图1 研究区地理位置

Fig.1 Geographical location of study area

1.2 数据源

本文基于2014年欧空局发射的Sentinel-1A SAR卫星影像数据开展实验分析,卫星工作波段为C波段,重访周期12 d。实验选取卫星在IW成像模式、VV极化方式下共10景升轨单视复数影像(single look complex, SLC)时序数据集,时间跨度: 2015年10月3日—2016年3月31日(表1),覆盖范围: 35.272° ~35.365° N,115.854° ~115.921° E。

表1 干涉对时空基线

Tab.1 Interference pair space-time baseline

序号	主影像日期	辅影像日期	时间基线/d	垂直基线/m
1	20151003	20151015	12	77.60
2	20151015	20151027	12	76.98
3	20151027	20151120	24	-122.34
4	20151120	20151202	12	7.22
5	20151120	20151214	24	32.01
6	20151120	20151226	36	92.14
7	20151202	20151214	12	24.79
8	20151202	20151226	24	84.92
9	20151202	20160107	36	108.82
10	20151214	20151226	12	60.13
11	20151214	20160107	24	84.02
12	20151226	20160107	12	23.90
13	20160107	20160307	60	-77.03
14	20160307	20160331	24	-62.22

新窗口打开| 下载CSV

除此之外,实验还使用了90 m分辨率的航天飞机雷达测图计划(shuttle Radar topography mission, SRTM)数字高程模型(digital elevation model, DEM)、Sentinel-1A精密定轨星历数据以及InSAR通用大气校正在线服务(generic atmospheric correction online service for InSAR, GACOS)产品进行辅助数据处理。

2 研究方法

2.1 Stacking技术

Stacking技术是一种利用加权平均方法计算地表形变速率的时序InSAR技术,数据处理流程相对简单高效,可以在短时间内获取大范围地表形变特征。其方法原理基于一个先前假设: 对于生成的每一幅独立差分干涉对,相对应的大气延迟相位都是相等且随机的,地表形变可以近似看成线性形变。同时,为减缓失相干因素的影响,将时空基线加以精细化约束进而筛选出M幅具有高相干性的差分干涉影像,利用筛选影像进行线性叠加,估算出年平均形变速率,其公式为:

(1)v_{disp_rate}=

\overset{M}{\sum_{m = 1}}

T_m·ψ_m/

\overset{M}{\sum_{m = 1}} T_{m}^{2}

式中: m为差分干涉对序列;v_{disp_rate}为干涉叠加得到的形变速率;ψ_m为第m幅干涉图的解缠相位;T_m为第m幅干涉图的时间基线。

由于只凭借线性叠加估算来获取形变速率,Stacking技术提取到的仅仅是线性形变信息,该技术在线性形变趋势明显的区域具有良好适用性,但在复杂环境下还需要结合其他手段进行更深入的研究。利用已获得的Stacking监测结果,可以模拟出各差分干涉对线性形变相位,模拟相位公式如下:

(2)ψ_lindef,_m=-

\frac{4 π}{λ}

·v_{disp_rate}·T_m,

式中: ψ_lindef,_m为依据Stacking监测结果模拟的第m幅干涉影像线性形变相位; λ为卫星传感器波长。

2.2 Stacking辅助下的DS-InSAR技术

前已论述,Stacking技术可以高效快速地获取矿区线性时序形变特征。而在实际矿区复杂环境背景下,完整的形变场通常由众多相位分量联合贡献,包括线性形变相位、高程误差、非线性形变相位、大气延迟相位和噪声相位等。时序InSAR技术可以将各分量分离并进行精化处理,最终反演高精度形变场,但以上过程需要基于准确解缠结果,当矿区存在大梯度形变时,解缠精度无法得到保证,各相位分量的精化与去除存在很大误差。

因此,所提技术在前期依据Stacking监测结果模拟出各差分干涉对线性形变相位后,同步将其引入DS-InSAR处理流程,辅助DS-InSAR降低相位梯度,进而达到利于解缠的目的。DS-InSAR技术中的2个关键环节分别为同质像元识别和相位优化。在同质像元识别领域存在许多经典方法,本文采用置信区间假设检验算法(hypothesis test of confidence interval, HTCI)进行同质像元识别。HTCI算法是对FaSHPS参数统计方法的一种修改,在FsSHPS基础上,该方法通过结合似然比检验估计无偏性可控制类型I误差的特点识别初始同质像元集合,然后依据SAR影像强度服从指数分布的特点,其累加和服从伽马分布,去除高斯假设限制,进而获得更精确的置信区间,增强同质像素选取的准确性。HTCI算法所构筑的区间为:

(3)p

\{g_{α / 2; N} \cdot \frac{{\hat{μ}}_{r e f}}{N} < \bar{A} (X) < g_{(1 - α) / 2; N} \cdot \frac{{\hat{μ}}_{r e f}}{N}\}

=1-α,

式中: p $\{\cdot\}$ 为概率; ${\hat{μ}}_{r e f}$ 为参考像素时间均值; $\bar{A}$ (X)为任意邻域像素时间均值;α为置信度;g_α/_2;_N和g_(1-_α₎_/_2;_N分别为标准伽马分布的α和1-α分位点; N为SAR影像数量。

相位优化本质是从干涉图集中估计出一组质量最优的干涉相位结果,并且优化后满足相位一致性准则。相位三角剖分算法(phase triangulation algorithm, PTA)基于最大似然估计量来建立相位优化模型,被证实具有良好的应用效果。定位到时间序列影像集N中的一个确定像素P,其样本的协方差矩阵公式为:

(4)

\hat{C}

=E[y(P)y(P)^H]≈

\frac{1}{|Ω|} \sum_{P \in Ω}

y(P)y(P)^H,

式中: (·)^H为共轭转置; E[·]为期望运算; Ω为识别出的同质像素集合; y(P)为该样本像素P的时间序列观测值向量。当其幅度数据被归一化后(E[|y_j|²]=1,j=1,2,…,N),可以得到一个相应的最大似然估计复相干矩阵 $\hat{Γ}$ ,其公式为:

(5)

\hat{Γ}

{\hat{γ}}_{a b}

e^{i ψ_{a b}}

式中: ${\hat{γ}}_{a b}$ 为a和b时刻2景SLC影像组成干涉对的相干系数;ψ_ab为干涉对的干涉相位; i为虚数单位。进一步地,假设样本真实复相干矩阵Γ可以表示为:

(6)Γ=ΦrΦ^H,

式中: r为N×N矩阵,矩阵元素对应所有干涉对的真实相干值。Φ=diag(e^iφ),是一个包含像素P真实相位值的N×N对角阵; φ为真实的SLC相位值。故同质像素的联合概率密度函数(probability density function, PDF)可表示为:

(7)p(

{\hat{Γ}}_{Ω}

|Φ)∝

e^{- t r a c e (Φ r^{- 1} Φ^{H} \hat{Γ})}

式中: p( ${\hat{Γ}}_{Ω}$ |Φ)为像素的概率密度函数; ∝为正比子; trace()为矩阵逐运算。

通过最大化上述PDF即可得到相位序列 $\hat{φ}$ 的最大似然估计量其公式为:

(8)

\hat{φ}

=argmax[θ^H(-r^-¹°

\hat{Γ}

)θ],

式中: θ=e^jφ,表示N维真实相位向量;°为Hadamard乘积; arcmax()为求极大值点集合函数。由于真实相干性矩阵r未知,算法采用r估计量的绝对值进行代替,公式为:

(9)

\hat{φ}

=argmax{θ^H(-

| \hat{Γ} |^{- 1}

\hat{Γ}

)θ} 。

获取到研究区线性形变相位与干涉优化相位后,将二者在复数域进行共轭相乘,即可在一定程度上去除线性形变相位,实现相位梯度降低,该步骤完成后可得到各干涉对优化后的残余相位值,其计算公式为:

(10)ψ_res,_m=

e^{ψ_{o p, m} i}

·conj(

e^{ψ_{l i n d e f, m} i}

式中: ψ_res,_m为第m幅干涉图的残余相位;ψ_op,_m为第m幅干涉图的优化相位;ψ_lindef,_m为第m幅干涉图的模拟线性形变相位; conj(·)为共轭运算。

通过设置相干性阈值,继续利用优化后的残余相位筛选高相干性DS点,并融合PS点进行后续时序InSAR处理。直到相位解缠完毕后,结合外部GACOS大气改正数据分离大气相位等其他相位分量,同时将先期步骤中去除的模拟线性形变相位加回,最终实现对研究区完整、高精度的形变监测。以上所提出的融合方法数据处理流程如图2所示。

图2

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

图2 技术流程图

Fig.2 Technical flow chart

3 监测结果与分析

3.1 时序监测结果

按照图2所示的数据处理流程,首先利用融合方法获取了新巨龙煤矿时序地表沉降信息,同时利用常规的StaMPS-SBAS技术采用相同时间基线组合方式对该地区进行对照处理。通过2种方法获取的山东省菏泽市新巨龙煤矿2015年10月—2016年3月雷达视线方向(line of sight, LOS)地表形变速率如图3所示。由图可知,2种手段监测到的形变区域分布特征具有较好的一致性。研究时段内,该区域共存在3处显著形变区,即北部、东部和南部形变区(图中对应标记为1、2、3),分别位于于家庄、毕海村以及张楼村附近,周边其他村落与积水区环绕。形变区挨靠大量住宅地带,为当地居民的生产生活带来了一定的影响。

图3

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

图3 研究区时序监测结果

Fig.3 Time series monitoring results of the study area

利用SBAS技术共计获取到36 071个监测点,约占总像素个数的3.7%,密度仅约795像素/km²,其中识别出的最大沉降速率仅约300 mm/a,位于东部形变区,对应累积最大沉降量仅约112 mm。进一步分析可知,SBAS技术识别出的形变场较为模糊,特别是北部形变区,其覆盖范围无法有效判别,形变量级也与前期Stacking实验结果存在出入。此外,该技术识别出的监测点大多集中分布于村落建筑物附近,在植被区、耕地区的监测点密度非常稀疏,而通过查阅该试验区相关地质资料可知,新巨龙煤矿所存在的3处主要的形变区均属于耕地覆盖区域,因此SBAS监测结果无法准确反映开采沉陷区的变形状况。相比之下,基于融合方法识别出的监测点密度与覆盖度都有非常明显的提升,3个主要形变区得到更为完整的呈现,可以更加细致地刻画形变的空间分布。利用融合方法获取的最大沉降速率点同样位于东部形变区,约为555 mm/a,对应累积最大形变量达到了约-313 mm,监测结果为后续的地表形变分析提供了丰富的解译信息。结合整个实验评估,所提改进方法共识别出了464 715个有效监测点,像素占比达到47.4 %,监测点密度约为10 236像素/km²,2项指标分别较SBAS技术提升了约43.7 %和12.9倍。可以看出,所提方法获取的监测结果能够较好地反映出新巨龙矿区该时段内所发生的地表沉降情况。

为分析融合方法监测结果的可靠性,以SBAS监测结果作为参照,对2种技术公共监测点的形变速率进行相关性分析。由图4可知,通过分析2种监测结果的同名点信息,其二者之间的Pearson相关系数(R)约为0.74,说明所提方法的监测精度具有较好可靠性。

图4

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

图4 同名点形变速率相关性分析

Fig.4 Correlation analysis of deformation rate of the same name points

3.2 水准数据对比分析

为对2种方法的差异性进行准确评价,本研究结合高精度水准数据对监测结果进行分析。选取了研究时段内(2015年11月15日—2016年3月15日)26个高精度水准点作为基准开展分析对比,水准点布设在南部形变区,近似垂直分布。由于数据采集的6个时间段内个别水准点数据缺失,因此实验过程将剔除部分水准点。其中沿走向方向选取12个,由北到南分别标注为C19-C1(C16,C13-C10,C8,C4除外),倾向方向选取14个,由东到西分别标注为T1-T15(T3除外),具体布设情况如图1所示。

考虑到SAR数据的特殊性,一维LOS向的形变结果实际上是三维(垂直向、东西向及南北向)位移分量在LOS向上的投影之和,并且InSAR结果对垂直向位移的敏感度最高。因此,通过假设地表只存在垂直向形变,可以将监测结果从LOS向转型为垂直向。同时,顾及影像日期与水准采集日期的差异性,统一SAR数据集与水准数据集时空基准。本次实验在InSAR垂直向结果基础上,通过插值算法将其获取的累积沉降量进行修正,最终得到对应水准采集日期的形变监测结果,并将T1点作为实验公共时空基准。

3种数据的累积最大形变量结果如图5所示。

图5

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

图5 融合技术、SBAS及水准结果比较

Fig.5 Comparison of fusion technology, SBAS and level

由图5分析可知,在走向方向上,2种方法监测结果与水准数据都有着较好的吻合度,融合方法结果一致性更好,累积最大形变量点位出现在C7点处,与水准结果相同。而倾向方向上融合方法监测结果则有更为明显的优势,没有出现异常误差,整体趋势更吻合水准数据。实验还针对SBAS技术和融合方法监测结果在相同时间段内与水准实测结果之间的差异性,从定量角度出发对二者误差做了具体的统计分析,如表2所示。SBAS结果最大误差为28.07 mm,平均误差约7.25 mm,均方根误差(root mean squared error, RMSE)为9.79 mm;相比之下,融合方法结果最大误差仅为22.55 mm,平均误差为4.90 mm,RMSE为6.82 mm。融合方法结果精度较传统SBAS技术提升了约3.0 mm。

表2 时序技术精度对比分析

Tab.2 Comparative analysis of timing technical accuracy

监测方法	最大误差/mm	平均误差/mm	RMSE/mm
SBAS	28.07	7.25	9.79
融合方法	22.55	4.90	6.82

新窗口打开| 下载CSV

为保证研究过程的严密性与合理性,实验同样对矿区进行了DS-InSAR对照处理。DS-InSAR所获取的监测点个数与密度较本文所提方法监测结果略低,其对比水准数据的RMSE为7.2 mm。对于所应用SAR卫星监测能力范围之内的地表形变,新巨龙煤矿虽然形变量级较大,但未达到本文所提方法的适用量级,因此其监测精度较DS-InSAR技术提升有限,但相较SBAS技术则有着全面的提升。后续可结合长波段卫星及更大梯度形变区进一步验证所提方法相较DS-InSAR技术的优越性。

综合上述分析可知,无论是定性还是定量方面,本文所提改进方法都表现出更优异的效果,因而在一定程度上有助于实现对复杂环境下矿区长时序、高精度的地表形变动态监测。

3.3 重点区域地表形变监测分析

结合前文所述,为对新巨龙煤矿地表形变特征有更加深入直观地了解,实验将基于融合方法监测结果探究其具体的变化现象。考虑到东部和北部形变场沉降量级与范围都较大,故本节选取二者作为重点区域展开详细分析。图6(a)首先展示的是东部形变区的沉降速率结果,该区域为新巨龙煤矿最大的形变场,区域内最大沉降速率达到了约555 mm/a,对应为图中的Z点。形变范围近似椭圆形,且西北-东南方向形变距离较长,西南-东北方向形变距离较短。结合地质开采资料分析,该区域存在着3个开采工作面,对应图中分别标记为L1,L2和L3,工作面同样呈现西北-东南方向分布,从而证实了东部形变区的地表沉降主要是由充填开采活动以及采空区的存在导致的。北部形变区(图6(b))在融合方法的处理下,其形变范围与量级较SBAS结果已经得到了更为良好的呈现。研究时段内存在L4,L5和L6这3个开采工作面,同样沿西北至东南走向分布。整体观察,北部形变场也大体呈现椭圆状范围,从监测点的分布特征来看,该区域已初步显现出一个沉降漏斗,可识别出的监测点最大形变速率约为-403 mm/a,并且形变区域几乎涵盖了下方开采工作面,进一步证实了该区域的地表沉降现象主要是煤矿长期的采掘活动造成的。

图6

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

图6 重点形变监测区

Fig.6 Key deformation monitoring area

除此之外,形变区沉降漏斗中心区域监测点出现数量不足的现象,其原因主要是该区域沉降变化量级过大,一定程度上超出C波段Sentinel-1A SAR卫星可监测范围,致使漏斗中心在干涉图中的相位信息太过模糊而无法提取。但相较于SBAS的监测结果,其提取到的监测点个数、展现出来的形变场信息的全面度和翔实度都得到了有效提升。

为充分挖掘东部形变区的时序形变演化规律,沿L1工作面设置了一条走向剖面线A-A',沿剖面线提取其时序形变量。由图7(a)可以看出,A-A'沿线沉降量分布呈明显的“抛物线”形状,累积最大沉降量约为250.35 mm,对应纬度约为35.308° 。分析其时序形变数据可知,2015年10月27日以前,剖面沿线基本处于一个整体下沉趋势,27日后开始出现突出沉降区,纬度范围约为35.308°~35.312° (剖面线中心区域),形变速率快,形变量级大。此后剖面沿线都是呈现出以突出沉降区为引领的规律下沉状态,沉降折线逐渐呈对称形状,越靠近中心区域,沉降量越大。同时,对特征点Z的时序形变量展开单独分析,特征点的时序形变结果在一定程度上也能反映出研究时段内该区域的形变特质。从图7(b)看出,特征点Z展现出较为显著的线性沉降趋势,整个观测期内,沉降量逐渐增大,并没有出现明显的减缓迹象,直至2016年3月31日,沉降量达到最大,约为312.7 mm。对此现象展开深入探究可知,观测时段内,L3工作面仍处于开采状态,开采工作面的持续影响导致了特征点的形变过程并没有明显减缓。

图7

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

图7 东部形变区时序数据

Fig.7 Time series data of eastern deformation region

为细致分析形变趋势,建立了B-B'和C-C'这2条剖面线做定量分析,如图8所示。剖面线B-B'沿走向设置,其完全贯穿整个L5工作面并经过该区域的沉降漏斗中心,剖面线C-C'沿倾向设置。从图中可以看出,虽然形变区已初步显现出一个沉降漏斗,但由于东南部的监测点稀少,其真实形变面貌仍无法准确获知。B-B'沿线上,最大累积沉降量出现在纬度约35.334 6° 处,大小约为179.65 mm,然而结合矿区开采沉陷规律,从该处沿B-B'继续推进,可能还会继续呈现沉降增大趋势,真实下沉曲线会比图8(a)拟合度更高。通过时序数据分析,在2015年10月3日—27日期间,B-B'沿线处于广泛形变期,没有突出的大沉降区域。10月27日之后,该区域沉降重心逐渐凸显,主要集中于纬度范围约35.332° ~35.340°,即L5工作面中心沉降地区,并逐渐向东南部偏移,其中2016年1月7日后发生了明显的偏移,并逐渐达到最大沉降。如图8(b)所示,在C-C'沿线上,自2015年11月20日起,经度约115.877°~115.879°处出现了明显的沉降中心,最大沉降量在经度约115.878° 处出现,约为218.45 mm。

图8

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

图8 北部形变区时序数据

Fig.8 Time series data of the northern deformation region

4 总结

针对时序InSAR技术在大梯度复杂矿区进行监测时,监测精度会受到监测点密度不足、干涉条纹密集导致的相位解缠困难等问题限制,本文提出了一种面向复杂矿区的Stacking技术辅助DS-InSAR地表形变监测方法。通过实验获取了山东省菏泽市巨野县新巨龙煤矿2015年10月—2016年3月时序动态地表形变特征,并结合相关煤矿开采和地质资料进行分析,总结如下:

1)提出的改进方法相较SBAS技术能够反演出更高覆盖度、更均匀的有效监测点,其密度达到10 236 像素/km²,数量约为后者的12.9倍。利用所提方法能够获取到试验区域更加精确的地表形变信息。

2)对比同期水准数据,所提方法的监测精度达到6.82 mm,在SBAS技术的基础上提高了约3.0 mm。方法在一定程度上解决了常规时序InSAR技术应用于大梯度矿区监测时存在的监测点密度不足、监测精度不高的问题。

3)实验监测时段内,新巨龙煤矿共存在3处显著形变场,区域存在的开采工作面以及充填开采活动是诱发地表形变的主要原因。研究区最大沉降点位于东部形变区内,最大形变速率达到-555 mm/a。提出的改进方法较为准确地探测出了形变场位置与形变范围,并且监测结果较为可靠,能够为大梯度复杂矿区煤炭开采活动提供一定的形变监测技术参考。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

王国法, 任世华, 庞义辉,

等.

煤炭工业“十三五” 发展成效与“双碳” 目标实施路径

[J]. 煤炭科学技术, 2021, 49(9):1-8.