地球物理学报  2021, Vol. 64 Issue (10): 3554-3566   PDF    
InSAR与概率积分法联合的矿区地表沉降精细化监测方法
陈洋, 陶秋香, 刘国林, 王路遥, 王凤云, 王珂     
山东科技大学 测绘与空间信息学院, 山东青岛 266590
摘要:结合InSAR与概率积分法的优势,提出一种InSAR和概率积分法联合进行矿区地表沉降的精细化监测方法.该方法首先计算InSAR时序累积沉降盆地,进而建立判别大梯度形变的约束条件,区分沉降边缘与沉降中心.对于形变较小的沉降边缘,保留InSAR结果,而对于大梯度沉降中心,则结合InSAR与概率积分模型建立矿区工作面的沉降盆地,并通过空间插值,获取地理坐标系下连续的地表沉降信息,最终得到完整的矿区地表沉降结果.论文以山东某矿区为研究区域,采用2016年10月16日—2018年3月4日期间的21景SAR影像和工作面水准实测数据对该方法的可行性和精度进行了验证.结果表明,该方法能够在减少水准监测工作量的前提下,获得与实际情况相吻合的沉降结果,其监测能力明显优于常规InSAR和概率积分法,可有效弥补两种技术单独在矿区地表沉降监测中的不足,获取更为准确、可靠的矿区地表沉降信息.
关键词: InSAR      概率积分法      矿区地表沉降      形变梯度     
Detailed mining subsidence monitoring combined with InSAR and probability integral method
CHEN Yang, TAO QiuXiang, LIU GuoLin, WANG LuYao, WANG FengYun, WANG Ke     
College of Geodesy and Geomatics, Shandong University of Science and Technology, Qingdao Shandong 266590, China
Abstract: By combining the advantages of interferometric synthetic aperture radar (InSAR) and the Probability Integration Method, this paper proposes a refined method for monitoring mining subsidence. In this method, we first calculate the time-series cumulative basin subsidence monitored by InSAR. Then, we establish the constraint conditions for assessing the large gradient subsidence. Subsequently, the subsidence edges and centre can be distinguished. For the edges, with small subsidence, the InSAR results are retained; for the centre, with large gradient subsidence, the InSAR results and the probability integration method are combined to establish the subsidence basin. Finally, a complete surface subsidence result for the mining area is obtained. This study took a mining area in Shandong Province, China, as the study area and used 21 SAR images acquired from 16 October 2016 to 4 March 2018 and levelling data of the working face level to verify the feasibility and accuracy of the method. The results showed that this method can obtain subsidence results consistent with actual conditions on the premise of reducing the workload of levelling monitoring. In addition, its monitoring capability is significantly better than that of conventional InSAR or the probability integration method. It can effectively compensate for the deficiencies of the two technologies that exist when they are used alone and obtain more accurate and reliable mining subsidence information.
Keywords: InSAR    Probability Integration Method    Mining subsidence    Deformation gradient    
0 引言

煤炭作为重要的能源基础,是工业现代化和社会发展的有力支撑.但随着能源需求的日益提升,煤炭开采在为世界经济创造巨大价值的同时,其引发的灾害问题也不容忽视:矿区大规模持续开采造成的沉降、塌陷,不仅损害地面基础设施,还容易诱发各种地质灾害,对人民财产安全造成巨大威胁(Yang et al., 2019; Saedpanah and Amanollahi, 2019; Sahoo and Khaoash, 2020).因此,为预防灾害发生,保证安全生产,必须对矿区进行长期有效的沉降监测与防治.

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR)是一种地面沉降监测新技术,可监测到厘米级甚至毫米级的地面变化,是探测地震、火山、矿区等形变的有效工具(Klein et al., 2017; 韩宇飞等, 2010; 申文豪等, 2019).20世纪90年代,该技术开始用于矿区地表沉降监测,所获结果在空间、时间上均连续,并且可以很好地反映矿区沉降的结构特征(Raymond and Rudant, 1997).相较于精密水准、GPS等传统测量手段,InSAR表现出全天候、低成本等的独特优势,为矿区地面沉降监测提供了新的思路与技术支持(Yang et al., 2018; Liu et al., 2013).但矿区开采引起的地面沉降往往呈现沉降速度快、量值大,空间不连续等特征,InSAR技术受限于SAR影像的固定特性,在下沉较为严重的沉降盆地中心会出现失相干现象,无法得到准确的形变信息,很大程度上限制了该技术在实际工程中的推广应用(Luo et al., 2019; Carlà et al., 2018; Zhang et al., 2020).

为解决InSAR监测矿区大梯度沉降时出现的失相干问题,国内外学者首先从技术本身出发,通过减少对流层延迟、选取长波段数据等方法提升InSAR的大梯度沉降监测能力,后又结合实际地质模型展开进一步的研究与探索,其中概率积分法是应用最为广泛的一种矿区沉降预计模型(Xu et al., 2006; Wang et al., 2020; Ge et al., 2007; 陶秋香等, 2012; Fan et al., 2014).该方法是由我国学者刘宝琛、廖国华等基于随机介质理论提出的一种沉降预计方法,利用高精度监测点,通过选定适当的数学模型,完成参数优化反演,获取沉降中心的下沉信息,有助于解决InSAR在矿区大梯度沉降监测中的失相干问题(刘宝琛和廖国华, 1965; 刘宝琛和戴华阳, 2016).但该方法由于受到模型自身特性的限制,在沉降边缘处收敛过快,尤其模型通常依赖成本高、耗时多的水准数据,导致其在矿区沉降监测中仍存在一定的局限性.

综合InSAR与概率积分法在矿区地表沉降监测中的优缺点,本文联合InSAR与概率积分两种监测方法,获取准确的矿区沉降盆地.该方法首先选取覆盖矿区工作面的SAR影像集,利用InSAR技术,按照时间序列两两差分,并对其进行时序累积,获取矿区多时段累积沉降监测结果;然后基于相干性条件,选取InSAR沉降盆地边缘的高相干点,与沉降中心附近少量插值后的水准数据作为特征点,通过曲面拟合法对概率积分模型参数进行优化反演,得到概率积分预计的矿区地表沉降盆地;进而以相干性为基准,通过形变梯度理论确定两种监测结果的融合阈值,获取完整的沉降盆地.论文以山东某矿区为研究区,利用该方法获取矿区地表沉降信息,结果表明,该方法结合了InSAR与概率积分法的优势,既能反映矿区开采对周边环境的影响,又能提取沉降中心准确的地表下沉信息,实现矿区沉降的高效、准确提取.

1 基本原理 1.1 InSAR形变监测及最大形变梯度

实际观测中,InSAR获取的相位与雷达参数、天线位置、天线入射角及地面目标高程紧密相关,干涉过程中不可避免会受到各种误差因素的影响,导致干涉相位信息中不仅包含高程和形变信息,还存在其他干扰项(Massonnet and Feigl, 1998),即

(1)

由式(1)可以看出,InSAR获取的相位信息主要包括:平地相位φflat,地形相位φtopo,地表形变相位φdef,电离层延迟和对流层延迟引起的大气相位φatm以及噪声相位φnoise.显然,要想获取地表形变相位,必须将干扰相位从干涉图中分离.

当去除各类干扰后,其形变相位在干涉图中表现为一系列干涉条纹,且干涉条纹所包含的相位信息并不是绝对相位,而是各个像元之间的相对相位.为此,Massonnet和Feigl(1998)提出形变梯度的概念,用以代表任意两点间的相对位移,并给出了InSAR可监测的最大形变梯度为

(2)

式中,dmax表示InSAR理论上能够监测到的形变梯度最大值,μ表示单个SAR影像像元的尺寸,λ表示雷达波长.由式(2)可见,dmax理论上与雷达本身参数有关,雷达波长越长、单位像元越小,能够监测到的形变梯度越大.以多视处理后的C波段Sentinel(分辨率为20 m)和L波段ALOS(分辨率为10 m)数据为例,经过1∶4多视处理的Sentinel影像,理论上能监测到1.4×10-3的形变梯度,而经过1∶3多视处理的ALOS影像,理论上能监测到11.5×10-3的形变梯度.

式(2)是假定理想条件下,InSAR能监测到的最大形变梯度.然而在实际情况中,InSAR不可避免地会受到轨道误差、时空失相干、大气延迟效应等的影响,最大形变梯度监测值往往比理论值要小.针对这一现象,Baran等(2005)分析研究了相干性对InSAR技术形变梯度的影响,给出了以影像相干性为基准的形变梯度阈值表达式,即

(3)

式中,Dmax为实际InSAR最大可监测形变梯度,γ为影像相干性.由式(3)可以看出,InSAR可监测最大形变梯度与影像相干性呈线性正相关,相干性越高,实际可监测最大形变梯度越大.对于SAR影像,一般存在某一相干阈值,使得Dmax=0.例如,对Sentinel影像,相干系数γ=0.3时,Dmax=0,即对于Sentinel数据来说,当影像相干性小于0.3时,就无法监测到地表形变.通常情况下,矿区常常具有植被覆盖率高、短时间内形变量值大的特点,为了更为有效地提取形变信息,在矿区利用InSAR技术进行地表沉降监测往往对相干性有着更高的要求.

1.2 概率积分模型

开采沉降研究中的岩体往往成因复杂,由于各类地质作用,使得岩体生成过程中被各种断层、节理等不均匀介质切割,产生一系列次生裂隙及非连续面,而矿山开采会再次扰动岩体的原生结构,呈现出更加明显的非连续性,因此非连续介质模型更适合开采沉降的应用研究.20世纪50年代,Litwiniszyn(1956)基于非连续介质模型,首次提出随机介质理论,并将其用于岩层移动规律的研究,认为随机介质的颗粒介质模型可用来描述由于矿山开采引发的地下岩层及地表移动规律.20世纪60年代,我国学者刘宝琛和廖国华(1965)在此基础上提出概率积分法,该方法将矿山岩层移动看作一个服从正态分布的随机现象,结合概率密度函数模拟地表沉降盆地的形态.对于走向长为D3、倾向长为D1的工作面,其开采引发的地面任意点(x, y)的沉降值可表示为

(4)

其中:

(5)

式(4)、(5)均在工作面坐标系下,其中W0表示最大下沉值,W0(x)、W0(y)分别表示待求点在走向、倾向主断面上投影点的下沉值,u表示概率积分的开采单元,m为煤层厚度,H0为平均开采深度,H1H2分别为下山、上山开采深度,q为下沉系数,tanβ为主要影响角正切,tanβ1、tanβ2分别表示下山、上山主要影响角正切,S3S4分别表示左右边界的拐点偏移距,S1S2分别表示下山,上山拐点偏移距,θ为开采影响传播角.

1.3 InSAR与概率积分法联合的矿区地表沉降监测方法

InSAR技术是获取矿区沉降的高效手段,能够准确获取沉降的边缘信息,但由于SAR影像自身特性的限制及各种噪声误差的干扰,该技术在矿区沉降量大的中心区无法探测到准确的沉降信息,监测结果与实际情况相差较大.概率积分法在沉降盆地中心具有较高的监测精度,但开采煤层容易受外力影响产生裂缝,加大煤柱上覆岩层的变形,导致概率积分法预计的沉降盆地边缘收敛过快,在沉降边缘的监测精度降低(王正帅和邓喀中, 2012).综合InSAR与概率积分法两种方法的沉降监测优势,可以获取准确完整的矿区地表沉降信息.其具体步骤如下:

(1) 选取覆盖矿区的N+1景SAR影像,按照时间序列对其进行两两差分处理,得到N个干涉对,并通过时序累积,获取矿区时序形变信息dInSAR.

(2) 收集矿区实际资料及工作面水准监测数据,考虑到水准点高程值在各个监测时段内近似呈线性变化,对其进行时间域上的分段线性插值,实现InSAR与水准在时间上的统一,预测出SAR成像日期下,各水准点的实际沉降信息,即

(6)

式中,dlevling表示水准点在SAR成像日期下的沉降信息,Tlevling表示水准点观测日期,TInSARTInSAR分别表示距水准点观测日期最近的前后两个SAR成像日期,dInSARdInSAR为这两个SAR成像日期下该水准点的InSAR沉降监测结果.

(3) 选取InSAR沉降盆地边缘的高相干点,与沉降中心附近少量插值后的水准数据作为特征点,并将其投影至工作面坐标系.假设坐标为(x, y)的特征点的沉降为D(x, y),则按照曲面拟合法可将其表示为一个与特征点坐标(x, y)及概率积分预计参数B有关的函数,即

(7)

DZk(k=1, 2, …, n)表示各特征点对应的实际沉降,并使曲面拟合值D(x, y)与实际监测结果的偏差满足平方和最小,即

(8)

式中,V为各特征点实际监测值相对于最小二乘拟合值的偏差.从而即可得到矿区概率积分参数的反演结果为

(9)

(4) 根据反演的概率积分参数,按照式(4)、(5),建立矿区概率积分沉降模型,并将沉降结果转换至地理坐标系下.进一步利用克里金插值法,对区域化沉降进行最优估计,获取连续的概率积分沉降结果dPIM.

(5) 统计N个干涉对的相干情况,计算其平均相干系数γavg,根据式(3)得到该相干条件下InSAR可监测到的最大形变梯度Dmax.基于形变梯度理论,考虑到Sentinel-1A影像的分辨率是20×20 m,可以计算得到N个干涉对能监测到的最大形变量dmax|InSAR(Massonnet and Feigl, 1998),即

(10)

dmax|InSAR作为InSAR与概率积分法沉降结果融合的阈值,从而得到完整的矿区地表沉降盆地,即

(11)

其数据处理流程如图 1所示.

图 1 InSAR与概率积分法联合的矿区沉降监测数据处理流程 Fig. 1 Data processing of mining subsidence monitoring using the method combined InSAR and Probability Integral Method
2 研究区及数据源 2.1 研究区概况

试验选取山东某矿区为研究区域,其地理位置及范围如图 2a所示.矿区地属黄河冲积平原,地形平坦,地面标高+40.53~+45.89 m,自然地形坡度为0.01%.矿区内交通便利,与邻近各县、乡均有公路相通,且有兖新铁路、220国道在其南部通过.区内河渠众多,主要有宋金河、鄄郓新河、郓巨河等,水源丰富,地表以植被为主.区内村庄众多,人口稠密,因此对于矿区来说,一方面要保证开采安全,另一方面,又要考虑到周边建筑、道路等地物的保护问题.同时兼顾矿区安全开采和地面维护,完成高效率、低损耗的生产,迫切需要有效的矿区地面沉降监测.

图 2 研究区地理范围 (a) 矿区地理位置及范围;(b) 工作面及水准点位置. Fig. 2 Geographical scope of thestudy area (a) Location and scope of the mining area; (b) Positions of the working surface and leveling points.

图 2b中白色矩形表示该矿区1301工作面的开采范围,该工作面上覆地表主要以植被为主,西北部为已开采的1300工作面老采空区,且有宋金河、郓城新河自其上方穿过.工作面平均煤厚6.8 m,自2016年10月开始开采,倾斜宽223.4 m,计划开采长度2265.8 m,至2018年3月,工作面推进1385.6 m,平均开采深度为938 m,地面标高为+41.8~+44 m.图 2(b)中绿色标记点表示工作面上方的水准点,沿倾向线自西向东布设有H1~H69,共69个水准点,沿走向线自南向北布设有Z1~Z92,共92个水准点,测量工作自2016年10月16日开始,至2018年3月24日结束,共19期.

2.2 数据源

哨兵1号(Sentinel-1)卫星于2014年4月成功发射,载有C波段合成孔径雷达.该卫星在轨运行高度为693 km,数据更新周期为12天,是一个可实现时间、空间连续观测的雷达成像系统,实现了单、双极化等若干种不同的极化方式.可提供连续影像(白天、夜晚和各种天气),以其大范围、多模式、多应用的特点为更多用户提供了数据服务.为获取矿区的地面沉降情况,本试验选取覆盖矿区的21景C波段、VV极化的Sentinel-1A升轨数据(成像中心入射角约为38.92°),时间跨度为2016年10月16日—2018年3月4日.具体参数如表 1所示.

表 1 SAR影像参数 Table 1 Specific parameters of SAR images

InSAR技术需要通过外部数字高程模型(DEM)模拟地形相位,并将其从SAR干涉相位中移除,以提取地表沉降相位和沉降值.由于研究区域地形平坦,实验选取3*3弧秒的SRTM3 DEM,可以满足矿区地面沉降监测的精度需求.SRTM3 DEM数据的采集间隔为3弧秒,地面分辨率为90 m,平均精度为16 m.

3 联合法提取的矿区地表沉降结果及分析 3.1 InSAR矿区沉降监测

表 1中的20个干涉对按照时间序列进行两两差分处理,考虑到干涉相位不可避免会受到多种噪声的影响,利用自适应滤波方法进行噪声去除,得到滤波后增强的差分干涉相位及相应的相干系数图,如图 3图 4所示.

图 3 滤波后的差分干涉图示例 Fig. 3 Examples of filtered differentialinterferograms
图 4 相干系数图示例 Fig. 4 Examples of coherence images

图 3可见,矿区在监测初期,干涉条纹较少,之后随着时间推移,矩形区域内干涉条纹变得密集,因而判定该区域在研究时段内可能发生地面沉降.结合图 4可知,第18、19个干涉对在监测时段内,相干性较好,相应的差分干涉图中显现出清晰完整的干涉条纹;相比而言,第9、12个干涉对噪声较大,出现严重的失相干,这是由于夏季植被覆盖率较高或沉降梯度过大,超出影像的实际监测能力导致的,对比此时的差分干涉图,在沉降区域,干涉条纹杂乱,甚至未形成完整的条纹形状.

将滤波增强后的差分干涉相位进一步进行解缠处理,并将解缠后的真实相位转换成雷达视线向上的地表形变值,并通过投影、地理编码得到地理坐标系下垂直向上的真实形变量.为了实现研究区域时序沉降盆地的监测,假设起始监测日期2016年10月16日无形变发生,将InSAR监测结果进行时序累加,得到2016年10月16日—2018年3月4日期间各成像时刻的累积沉降结果,如图 5所示.

图 5 InSAR累积沉降结果 Fig. 5 Cumulative subsidence monitored by InSAR

图 5分析可以得到,1301工作面自2016年10月开始开采,至2016年12月,InSAR监测到工作面起采线附近开始出现漏斗状的沉降盆地,说明此时工作面的开采活动已对上覆地表产生影响.但由于开采时间较短,沉降量较小,截至2016年12月3日,InSAR监测到的最大累积沉降为22 mm.随后,由于开采工作面的持续推进,矿区沉降范围及沉降量相对于监测初期而言,明显增大,至2017年6月1日,InSAR监测到的累积沉降达到225 mm.随着工作面持续向北推进,监测结果显示沉降中心开始沿工作面走向逐步向北移动,至2017年11月28日,监测到的沉降中心距离起采线约160 m,最大累积沉降为420 mm.由此可见,InSAR监测到的沉降位置、范围及空间变化趋势与矿区实际开采工作面相符,但在沉降中心处的大梯度沉降监测能力明显不足.随着工作面的持续开采,矿区沉降继续加大,至2018年3月4日,水准监测到工作面最大累积沉降达到2813 mm,而InSAR监测到的最大沉降仅为548 mm,与矿区实际存在较大差异.

3.2 特征点选取

为了构建概率积分模型,进一步提取InSAR沉降结果中的高相干特征点.由InSAR监测最大形变梯度式(3)可知,对于Sentinel-1A SAR影像来说,当相干性系数小于0.3时,InSAR监测结果是不可信的.在实际情况中,由于矿区短时间内沉降过大,其监测精度对相干性有更高的要求.因此,根据图 4矿区的相干情况,这里将0.5作为相干性阈值,认为相干性大于0.5的特征点为高相干点,提取InSAR沉降盆地边缘的13个高相干点,并根据概率积分原理,选取临近矿区最大沉降点和拐点的5个水准点,共18个特征点,其叠加在工作面的分布情况如图 6所示.表 2给出了各点的相干性系数及沉降信息.

图 6 特征点分布图 Fig. 6 Distribution of feature points
表 2 特征点相干系数及沉降信息 Table 2 Coherence coefficient and subsidence information of the feature points
3.3 概率积分沉降预计

结合选取的18个特征点和概率积分模型,利用曲面拟合法,反演得到该矿区的概率积分参数为

将这些参数代入式(4),计算得到矿区1301工作面的概率积分沉降预计盆地,如图 7所示.为获取连续的概率积分地面沉降,考虑到Sentinel-1A影像分辨率为20×20 m,在矿区1301工作面沉降盆地处,以20 m为间隔共选取18396个特征点,计算其沉降值,并将结果转换至地理坐标系下,进一步利用克里金插值法,对区域化沉降进行最优估计.图 8给出了插值后空间连续的概率积分矿区沉降结果.

图 7 概率积分矿区沉降预计盆地立体图 Fig. 7 Three-dimensional map of the subsidence basin of the mining area predicted by the Probability Integral Method
图 8 连续的概率积分沉降结果 Fig. 8 Continuous subsidence monitored by the Probability Integral Method

图 7图 8可以看出,研究时段内,概率积分法探测到1301工作面上方呈现明显的漏斗状沉降盆地,最大沉降为2528 mm.但该方法测得工作面周边几乎无沉降发生,沉降盆地在边缘处收敛过快,其原因是在实际情况下,矿区沉降往往存在复杂地形、地下水流失及其他采空区开采等外部因素的影响,而概率积分模型仅考虑开采活动引发的地面沉降,导致其监测到的地面沉降集中在工作面上方.

3.4 联合法提取完整的矿区地表沉降盆地

图 4所示矩形区域内的相干系数结果进行均值滤波处理,并统计20个干涉对中沉降区域的相干情况,统计结果如表 3所示.分析表 3可知,2016年10月16日—2018年3月4日期间,矿区沉降区域的相干性为0.39~0.64.其中夏季4—8月,为农作物生长时节,植被覆盖率高,获取的SAR影像相干性普遍较差,秋冬季节SAR影像质量相对较好.

表 3 矿区20个干涉对的相干系数统计 Table 3 Coherence coefficient statistics of 20 interferometric pairs in the coal mine

根据表 3各干涉对的相干结果,计算得到研究时段内20个干涉对的平均相干系数γavg为0.502.结合式(3)、式(10),考虑到Sentinel-1A的影像分辨率为20×20 m,进一步计算得到此时InSAR能监测到的最大形变量dmax|InSAR为162 mm.以此为阈值,对InSAR沉降边缘和概率积分沉降中心结果进行融合,得到矿区完整的沉降盆地,其结果如图 9所示.

图 9 联合法沉降监测结果 Fig. 9 Subsidence monitored by the Combined Method

图 9可以看出,InSAR与概率积分法通过融合,可以得到矿区完整的沉降结果,但在融合的边界却出现结果的不连续现象.为解决这一问题,本文在融合边界的两侧寻找一定区域内的InSAR与概率积分法预计的沉降数据,对其进行加权平均.

由InSAR监测结果可知,20个干涉对的平均相干系数γavg为0.502,但处于夏季的干涉对干涉效果差,相干性低,其相干系数低于平均相干系数,监测精度较低.因此,在加权平均中,选择相干性大于平均相干系数γavg的第1、2、3、4、5、17、18、19个干涉对,计算其平均相干系数为0.58.在该相干条件下,结合式(3)、式(10),考虑到Sentinel-1A的影像分辨率为20×20 m,计算得到此时InSAR能监测到的最大形变量为90 mm.以此作为加权平均的下临界阈值,当沉降量小于90 mm时,以InSAR监测值作为最终沉降结果.

对于概率积分法而言,拐点通常位于沉降盆地主断面上沉降曲线曲率为零的分界点,根据概率积分预计的沉降结果,以工作面走向起采线与终采线处概率积分沉降结果的平均值917 mm作为加权平均的上临界阈值.当沉降量大于917 mm时,以概率积分预计的沉降值作为最终沉降结果.

当沉降量介于90 mm到917 mm时,通过对InSAR与概率积分预计的沉降值进行加权平均,确定最终沉降结果.首先选择沉降量介于90~917 mm的水准点,共48个,计算其InSAR、概率积分法沉降结果与水准的绝对差值,如图 10所示.然后计算48个水准点InSAR监测结果、概率积分法预计结果的累积均方差,分别为186 mm、208 mm.根据

(12)

图 10 部分水准点的误差情况 Fig. 10 Error of partial levelling points

计算得到PPIM=0.8PInSAR.这里取PPIM=0.44,PInSAR=0.56,则由

(13)

得到加权平均后完整的矿区沉降盆地,如图 11所示.

图 11 加权融合后的联合法沉降监测结果 Fig. 11 Subsidence monitored by the Combined Method after weighted fusion

图 11可知,2016年10月16日—2018年3月4日期间,联合法监测到矿区出现地面沉降,且沉降盆地位置与矿区1301工作面高度吻合.研究时段内,该方法监测到工作面周边存在小量级沉降,这是因为工作面附近存在一系列小断层,开采工作的持续推进,扰动了其上覆岩层,产生微小形变,对周边环境造成不同程度的影响,与矿区实际开采进度及地质情况相符,说明该方法在沉降盆地边缘的监测效果优于概率积分法.由联合法在沉降盆地中心的监测结果分析可以得到,截至2018年3月4日,该方法监测到1301工作面的沉降中心位于距离起采线约600 m处,沉降中心向东部偏移,这主要是因为1301工作面存在约10°的地面倾斜,沉降中心向下山方向偏移,与矿区实际情况相符,其监测到的最大累积沉降达到2528 mm,在沉降中心的监测能力明显优于InSAR手段.

4 联合法监测结果精度评估

为进一步验证联合法监测结果的精度,提取矿区1301工作面倾向线上H1~H69,走向线上Z1~Z92,共161个水准点的联合法监测结果、InSAR累积沉降、概率积分法预计沉降,并收集实际水准测量数据,对三种监测方法的测量结果进行定量对比分析,图 12给出了对比结果图.以InSAR能监测到的最大形变量dmax|InSAR=162 mm为区分沉降盆地边缘和中心的阈值,分别统计三种监测方法在沉降盆地边缘和中心的误差,表 4给出了三种监测方法的误差对比情况.

表 4 误差统计 Table 4 Error statistics

图 12表 4三种监测方法的对比结果,分析可以得到:

图 12 监测结果对比图 Fig. 12 Comparison of monitored results

(1) 2016年10月16日至2018年3月4日期间,联合法、InSAR和概率积分法均探测到1301工作面上方出现地面沉降,监测得到的沉降盆地位置、沉降曲线形态及变化趋势基本一致,且与水准监测结果相符.

(2) 通过与概率积分法结果对比,联合法在沉降边缘处与水准吻合程度更高.分析图 12(c)(e)中二者在倾向线上的监测结果可以得到,在沉降边缘处,由于1301工作面西部存在小断层,受其影响,导致沉降范围向西部扩大,联合法监测到这一沉降趋势,在H1~H23范围内,其测得的沉降曲线与水准吻合较好.对比此时概率积分法的监测结果,发现该方法探测到的沉降范围略小,沉降曲线在边缘处收敛较快,分析是由于概率积分法在反演过程中并未考虑外部条件对沉降的影响,造成该方法对于受到断层影响而引起的沉降不敏感.分析图 12(d)(f)中走向线上二者的沉降曲线同样发现,在沉降盆地的边缘区域,联合法可以获得与实际更相符的监测结果.在沉降边缘处,联合法监测结果的平均绝对误差为30 mm,均方根误差为40 mm,最大绝对误差为104 mm,相比概率积分法而言,三种误差分别提升53%、69%和84%.

(3) 通过与InSAR结果对比,联合法在沉降中心处的大梯度监测能力明显提升.图 12 (a, e)二者的监测结果显示,水准监测到倾向线上32点处出现2813 mm的最大沉降.InSAR虽也探测到此处为沉降中心,但监测到的最大沉降点与水准略有差别,出现在29点处.监测结果与水准相比,最大绝对误差达到2445 mm,这主要是由于矿区沉降中心沉降梯度过大,超出Sentinel-1A影像的实际监测范围,加上地表植被、噪声等降低了SAR影像质量,导致其在沉降中心处无法探测出真实的沉降信息.对比联合法获得的沉降结果可以看出,在沉降中心处,该方法监测结果与水准吻合较好,其监测得到的最大沉降点也与水准结果一致,监测到的最大累积沉降为2510 mm,与水准相差303 mm.图 12(b, f)为走向线上两种技术在沉降盆地中心的监测结果,表现出与倾向线上相同的特征,InSAR与联合法监测到的最大沉降分别为469 mm和2510 mm,与水准的最大绝对误差分别为2396 mm和706 mm.在沉降中心处,联合法监测结果的平均绝对误差为198 mm,均方根误差为258 mm,最大绝对误差为706 mm,最大沉降点处的绝对误差为303 mm,相比InSAR而言,四种误差分别提升81%、80%、71%和88%.

(4) 根据三种监测方法的误差统计分析可知,矿区倾向线、走向线上所有水准点的联合法监测结果的平均绝对误差为121 mm,均方根误差为192 mm,最大绝对误差为706 mm,最大沉降点处的绝对误差为303 mm.相比InSAR和概率积分法而言,整体监测精度有所提升,四种误差相较于InSAR分别提升79%、80%、71%和88%,平均绝对误差与均方根误差相较于概率积分法分别提升27%和17%.

5 结论

本文针对InSAR和概率积分法目前各自监测矿区地表沉降存在的不足,基于影像相干性及InSAR形变梯度理论,提出了两种技术联合进行矿区沉降监测的方法.以山东某矿区1301工作面为研究区,利用该方法进行沉降监测,获取了研究时段内矿区沉降盆地.通过与同期水准实测数据、InSAR累积沉降和概率积分监测结果进行对比分析,得出如下主要结论:

(1) 联合法相比概率积分法而言,在沉降盆地边缘的监测精度有所提升,可以在减少水准观测的基础上,很好地避免模型本身的缺陷.该方法在沉降边缘处与水准的平均绝对误差、均方根误差和最大绝对误差相较于概率积分法分别提升53%、69%和84%,可以解决概率积分法在矿区地表沉降监测中收敛过快的问题,监测结果与水准更为吻合.

(2) 联合法相比InSAR技术而言,在沉降盆地中心的监测能力明显提升.该方法在沉降中心处的平均绝对误差、均方根误差、最大绝对误差和最大沉降点处的绝对误差相较于InSAR分别提升81%、80%、71%和88%,能够有效弥补InSAR在监测矿区大梯度地表沉降时的不足.

(3) 联合法结合InSAR与概率积分法各自的优点,可以准确探测出矿区沉降盆地的位置、范围及沉降变化趋势,反映矿区开采对周边环境的影响,且监测得到的沉降中心符合矿区开采的地表移动规律,能够提取沉降中心准确的下沉信息.相比概率积分法和InSAR技术而言,该方法的整体监测精度有所提升,可以获得与实际情况更为吻合的地表沉降.

致谢  感谢欧洲航天局(European Space Agency, ESA)提供免费的Sentinel-1A SAR数据和美国国家航空航天局(National Aeronautics and Space Administration, NASA)提供SRTM3 DEM数据.
References
Baran I, Stewart M, Claessens S. 2005. A new functional model for determining minimum and maximum detectable deformation gradient resolved by satellite radar interferometry. IEEE Transactions on Geoscience and Remote Sensing, 43(4): 675-682. DOI:10.1109/TGRS.2004.843187
Carlà T, Farina P, Intrieri E, et al. 2018. Integration of ground-based radar and satellite InSAR data for the analysis of an unexpected slope failure in an open-pit mine. Engineering Geology, 235: 39-52. DOI:10.1016/j.enggeo.2018.01.021
Fan H D, Gu W, Qin Y, et al. 2014. A model for extracting large deformation mining subsidence using D-InSAR technique and probability integral method. Transactions of Nonferrous Metals Society of China, 24(4): 1242-1247. DOI:10.1016/S1003-6326(14)63185-X
Ge L L, Chang H C, Rizos C. 2007. Mine subsidence monitoring using multi-source satellite SAR images. Photogrammetric Engineering & Remote Sensing, 73(3): 259-266. DOI:10.14358/PERS.73.3.259
Han Y F, Song X G, Shan X J, et al. 2010. Deformation monitoring of Changbaishan Tianchi volcano using D-InSAR technique and error analysis. Chinese Journal of Geophysics (in Chinese), 53(7): 1571-1579. DOI:10.3969/j.issn.0001-5733.2010.07.008
Klein E, Vigny C, Fleitout L, et al. 2017. A comprehensive analysis of the Illapel 2015 MW8.3 earthquake from GPS and InSAR data. Earth and Planetary Science Letters, 469: 123-134. DOI:10.1016/j.epsl.2017.04.010
Litwiniszyn J. 1956. Application of the equation of stochastic processes to mechanics of loose bodies. Archives of Applied Mechanics, 8(4): 393-411.
Liu B C, Liao G H. 1965. Basic Law of Surface Movement of Coal Mine (in Chinese). Beijing: China Industry Press: 50-56.
Liu B C, Dai H Y. 2016. Research development and origin of probability integral method. Coal Mining Technology (in Chinese), 21(2): 1-3.
Liu Z G, Bian Z F, Lv F X, et al. 2013. Monitoring on subsidence due to repeated excavation with DInSAR technology. International Journal of Mining Science and Technology, 23(2): 173-178. DOI:10.1016/j.ijmst.2013.04.024
Luo H B, Li Z H, Chen J J, et al. 2019. Integration of range split spectrum interferometry and conventional InSAR to monitor large gradient surface displacements. International Journal of Applied Earth Observation and Geoinformation, 74: 130-137. DOI:10.1016/j.jag.2018.09.004
Massonnet D, Feigl K L. 1998. Radar interferometry and its application to changes in the Earth's surface. Reviews of Geophysics, 36(4): 441-500. DOI:10.1029/97RG03139
Raymond D, Rudant J P. 1997. ERS-SAR interferometry: Potential and limits for mining subsidence detection. //Guyenne T D, Danesy D eds. Third ERS Symposium on Space at the Service of Our Environment. Florence, Italy: European Space Agency, 414: 541.
Saedpanah S, Amanollahi J. 2019. Environmental pollution and geo-ecological risk assessment of the Qhorveh mining area in western Iran. Environmental Pollution, 253: 811-820. DOI:10.1016/j.envpol.2019.07.049
Sahoo S, Khaoash S. 2020. Impact assessment of coal mining on groundwater chemistry and its quality from Brajrajnagar coal mining area using indexing models. Journal of Geochemical Exploration, 215: 106559. DOI:10.1016/j.gexplo.2020.106559
Shen W H, Li Y S, Jiao Q S, et al. 2019. Joint inversion of strong motion and InSAR/GPS data for fault slip distribution of the Jiuzhaigou 7.0 earthquake and its application in seismology. Chinese Journal of Geophysics (in Chinese), 62(1): 115-129. DOI:10.6038/cjg2019L0541
Tao Q X, Liu G L, Liu W K. 2012. Analysis of capabilities of L and C-band SAR data to monitor mining-induced subsidence. Chinese Journal of Geophysics (in Chinese), 55(11): 3681-3689. DOI:10.6038/j.issn.0001-5733.2012.11.015
Wang LY, Deng K Z, Zheng M N. 2020. Research on ground deformation monitoring method in mining areas using the probability integral model fusion D-InSAR, sub-band InSAR and offset-tracking. International Journal of Applied Earth Observation and Geoinformation, 85: 101981. DOI:10.1016/j.jag.2019.101981
Wang Z S, Deng K Z. 2012. Edge-amended model of probability-integral method for subsidence prediction. Journal of Xi'an University of Science and Technology (in Chinese), 32(4): 495-499.
Xu C J, Wang H, Ge L L, et al. 2006. InSAR tropospheric delay mitigation by GPS observations: A case study in Tokyo area. Journal of Atmospheric and Solar-Terrestrial Physics, 68(6): 629-638. DOI:10.1016/j.jastp.2005.11.010
Yang Z, Li W P, Li X Q, et al. 2019. Assessment of eco-geo-environment quality using multivariate data: A case study in a coal mining area of Western China. Ecological Indicators, 107: 105651. DOI:10.1016/j.ecolind.2019.105651
Yang Z F, Li Z W, Zhu J J, et al. 2018. Locating and defining underground goaf caused by coal mining from space-borne SAR interferometry. ISPRS Journal of Photogrammetry and Remote Sensing, 135: 112-126. DOI:10.1016/j.isprsjprs.2017.11.020
Zhang L L, Cai X X, Wang Y, et al. 2020. Long-term ground multi-level deformation fusion and analysis based on a combination of deformation prior fusion model and OTD-InSAR for longwall mining activity. Measurement, 161: 107911. DOI:10.1016/j.measurement.2020.107911
韩宇飞, 宋小刚, 单新建, 等. 2010. D-InSAR技术在长白山天池火山形变监测中的误差分析与应用. 地球物理学报, 53(7): 1571-1579. DOI:10.3969/j.issn.0001-5733.2010.07.008
刘宝琛, 廖国华. 1965. 煤矿地表移动的基本规律. 北京: 中国工业出版社: 50-56.
刘宝琛, 戴华阳. 2016. 概率积分法的由来与研究进展. 煤矿开采, 21(2): 1-3.
申文豪, 李永生, 焦其松, 等. 2019. 联合强震记录和InSAR/GPS结果的四川九寨沟7.0级地震震源滑动分布反演及其地震学应用. 地球物理学报, 62(1): 115-129. DOI:10.6038/cjg2019L0541
陶秋香, 刘国林, 刘伟科. 2012. L和C波段雷达干涉数据矿区地面沉降监测能力分析. 地球物理学报, 55(11): 3681-3689. DOI:10.6038/j.issn.0001-5733.2012.11.015
王正帅, 邓喀中. 2012. 概率积分法沉陷预计的边缘修正模型. 西安科技大学学报, 32(4): 495-499. DOI:10.3969/j.issn.1672-9315.2012.04.016