基于DS-InSAR的西部矿区地表时序沉降监测与分析

王 波 谭志祥,2 邓喀中,2

(1.中国矿业大学环境与测绘学院,江苏 徐州 221116;2.自然资源部国土环境与灾害监测重点实验室,江苏 徐州 221116)

多年来对东部矿区煤炭资源的高强度开发,使得 东部矿区煤炭产量不断下降。为满足国民经济发展对于煤炭资源的需求,我国加大了对西部矿区煤炭资源的开发利用力度。与东部矿区相比,西部矿区煤炭资源储量更大,而且具有煤层厚度大、倾角小等优点,因此有着更大的生产潜力[1]。然而由于历史原因,西部矿区煤矿开采时间比东部矿区要短,因此针对西部矿区开采引起的地表形变的监测工作开展较少,相应的研究成果也比较匮乏。而且由于东西部矿区之间地质采矿条件的巨大差异,导致以东部矿区为研究对象所得到的相关研究成果在西部矿区适用性不强[2],这给西部矿区地表形变预计、村庄保护煤柱留设以及各类建筑物保护带来了困难,不利于煤矿生产工作安全高效进行。因此有必要针对西部矿区特定的地质采矿条件,开展矿区地表形变监测工作,并以此为基础进行西部矿区地表形变规律和机理的研究,为后续煤矿生产工作提供技术依据。

水准测量等传统的地表形变监测方法在西部矿区的应用具有较大的局限性,主要体现在受西部矿区高原堆积型沙丘地貌的影响,地面水准点以及地表特征点的选取和观测困难。而且传统的测量方法往往费时费力,难以进行大范围面状地表形变监测[3]。合成孔径雷达差分干涉测量技术(Differential InSAR,D-InSAR)能够不依赖地面特征点来获取地表形变信息,因此被广泛应用于地表形变的监测与分析中[4]。然而该技术极易受到时空去相干因素的影响,难以获取矿区开采引起的地表形变范围。针对这一问题,科研工作者将注意力集中在了那些能够在长时间序列合成孔径雷达(Synthetic Aperture Radar,SAR)影像中保持高相干性的点目标上,提出了永久散射体干涉测量技术(Permanent Scatterer InSAR,PS-InSAR)[5],在此基础上进行了相干点目标识别[6]、三维相位解缠算法[7]等研究,将D-InSAR技术推进到了时间序列InSAR的新阶段。以PS-InSAR和小基线集方法(Small Baseline Subsets InSAR,SBAS-InSAR)[8]为代表的常规时序InSAR技术在干涉对优化的基础上,提取分析了影像中受时空去相干影响较小而表现出相对稳定散射特性的高相干像元在整个监测时段内的变化,有效克服了D-InSAR技术在应用过程中的失相干问题。然而矿区及周边地表的覆盖多为植被和沙地,导致常规的时序InSAR技术在矿区及其周围难以获得足够数量、分布均匀的高相干点,地表形变监测精度受到了很大限制。

近些年发展起来的融合分布式目标的时序InSAR技术[9]将研究的重点从散射特性相对稳定的高相干点目标扩大到了散射特性中等的分布式目标,可同时利用这两种目标进行地表形变监测,从而极大地提高了监测区域点目标的密度,有效克服了常规时序InSAR技术在低相干区域由于无法获取足够数量的点目标而无法对地表形变进行准确监测的不足。该技术先是在地震形变监测、山区滑坡灾害监测、冰川运动监测等方面进行应用并取得了良好的效果[10],证明了该方法与常规时序InSAR技术相比具有优越性。之后部分学者开始尝试将该技术应用到矿区地表形变监测中并取得了一定的成果[11-18]。但是上述研究的试验区域大多位于东部矿区,对于该技术在西部矿区地表形变监测方面的应用研究涉及较少,缺乏针对性的研究。为此,本研究在现有成果的基础上,利用一种针对西部矿区高原堆积型沙丘地貌的基于DS-InSAR的时序地表形变监测方法,获取了西部矿区内蒙古石拉乌素煤矿地表形变及其时空分布信息,并验证了该方法的可靠性,丰富了DS-InSAR技术在西部矿区地表形变监测应用方面的研究。

1 研究区概况与数据源

1.1 研究区概况

石拉乌素煤矿位于内蒙古自治区鄂尔多斯市伊金霍洛旗东胜煤田内,矿区面积约70 km2,设计生产能力为1 000万t/a。矿区位于鄂尔多斯高原的中南部,具有典型的高原堆积型沙丘地貌特征,地表大部分被第四系风积沙所覆盖,局部有白垩下统地层露出,植被稀疏,为沙漠—半沙漠地区。本文研究区域所在的221采区位于石拉乌素煤矿的中部(图1),采区地形总体趋势是东高西低、北高南低,地面标高+1 332～+1 411 m,采区内煤层埋深近700 m,上覆大厚度的中生代砂岩,上覆岩层性质具体如表1所示,其地质采矿条件与东部矿区相比存在明显差异。在监测时段内221采区有17、18、06A 3个工作面开采,采区北部的17、18工作面采用综合机械化一次采全高工艺开采,平均采厚约5.40 m,南部的06A工作面采用综合机械化放顶煤工艺开采,平均采厚约9.78 m。各工作面的参数取值见表2。

表1 煤层上覆岩层性质Table 1 Properties of overlying strata of coal seam

表2 221采区工作面参数Table 2 Parameters of working face in 221mining area

图1 研究区地理位置Fig.1 Geographical locationof the study area

1.2 数据源

Sentinel-1卫星由欧空局发射,并免费向全球提供卫星影像下载服务,为SAR相关研究提供了极大的便利。本研究选用了2016年9月—2018年6月共52景覆盖研究区域的IW模式下VV极化的Sentinel-1A升轨影像。Sentinel-1卫星影像的分辨率为5 m×20 m,卫星的重访周期为12 d。此外本研究采用30m分辨率的航天飞机雷达测图计划(ShuttleRadar Topography Mission,SRTM)数字高程模型(Digital Elevation Model,DEM)从干涉图中去除地形相位的影响和地理编码。

2 监测方法选择及其原理

2.1 高相干点目标选取

从时序SAR影像中准确地识别高相干点目标是时间序列InSAR分析的重要环节,也是实现联合高相干点目标与分布式目标精确提取地表形变信息的前提条件。受西部矿区高原堆积型沙丘地貌的影响,研究区域内能够在长时间序列中保持相位稳定的永久散射体目标的数量稀少,但存在数量较多在短时间内经滤波后相位失相关性变小而保持高相干性的点目标。本研究对所选SAR影像进行处理得到滤波后的差分干涉图集,通过幅度差分离差法[19]选取出了此类点目标作为高相干点目标。幅度差分离差法原理可用下式表示

式中,DΔA为幅度差分离差值;σΔA为主、辅影像幅度差的标准差;μA为差分干涉图幅度的均值。通过设定幅度差分离差阈值,便可实现高相干点目标的选取。

2.2 分布式目标识别与相位优化

2.2.1 FaSHPS同质像元识别

受矿区特定地质条件的影响,研究区域内高相干点目标的空间分布极不均匀,大部分高相干点目标集中分布在地表裸露岩石较多的区域,而在地表大范围的积沙区域内此类点目标数量十分稀少。为实现对研究区完整地表形变信息的提取,必须在高相干点目标的基础上充分利用广大积沙区域内的分布式目标。参数假设检验方法中的快速同质点选取方法(Fast Statistically Homogeneous Pixel Selection,FaSHPS)[20]具有同质像元选取效率高且提取出的同质像元之间的异质性小的优点,为此本研究选用该方法对研究区域内的同质像元进行识别,并通过设置时间相干性阈值实现对分布式目标的提取。FaSHPS方法的思路为,假设时序SAR影像数据集内任一像元时间维上振幅的均值服从高斯分布,从而可以将传统的假设检验方法转化为置信区间的估计问题,只要利用简单的逻辑计算来判断两个像元服从的分布函数是否相同,便可实现同质像元的识别。

设由N景时序SAR影像组成数据集,其中任一像元L在时间维上的振幅均值为(L),由中心极限定律可知,随着样本数目N加大,振幅均值(L)趋于高斯分布,此时其区间估计可表示为

式中,P{·}表示概率;μ(L)为像元L的数学期望;z1-α/2为标准正态分布中置信度为α时对应的分位点;Var(A(L))为像元L时间维上振幅的真实方差。

当样本的数量足够大时,单视复数SAR影像的时序振幅信息近似服从瑞利分布,此时变异系数CV为常量,其值约为0.52。若像元在时间维上的后向散射特性稳定,则式(2)可表示为

当给定置信度α时,式(3)便给定了一个区间,若待估像元在时间维上的平均强度值在目标像元所对应的置信区间内,则待估像元与目标像元为同质像元。

2.2.2 特征值分解相位优化

分布式目标像元内地物匀质分布[21],其回波信号由像元内各子散射体相干叠加而成,易受时空失相干影响,相位稳定性较差。因此需要基于分布式目标像元散射特性的整体统计特征来去除噪声相位,实现对相位的优化。特征值分解法在相位优化的过程中无需迭代计算,因此其优化效率较高,同时该方法的相位优化质量也较好。故本研究选用该方法对分布式目标像元相位进行优化。特征值分解法的核心思想为从相干矩阵中分离,估计出多元散射机制并得到其对应的相位分量,通过分离出主导散射机制的相位分量来达到优化相位的目的[22-23]。

式中,Λ为N阶对角矩阵,对角线上的元素为相干矩阵的N个特征值,用λi表示;每个不同的特征值分别对应着不同的特征向量,用μi表示;特征值不同,其对应的散射机制也不同。特征值越大,在像元整体散射信号中其对应的子散射体所占比重就越大。因此将最大特征值对应的特征向量作为主散射体,其余的特征向量作为噪声,就可以实现同质像元相位的优化,此时相干矩阵可以表示为

在完成同质像元识别以及相位优化之后,通过计算相位优化前后干涉相位之间的拟合度便可对相位优化的质量进行评价。相位优化质量评估完成之后,设置时间相干性阈值,从分布式目标中筛选出相干性较好、信噪比较高的部分。将其与筛选出的高相干点目标一并融入时序InSAR处理流程中,便可对整个研究区域地表形变进行高精度监测。针对西部矿区高原堆积型沙丘地貌的基于DS-InSAR的时序地表形变监测方法的主要流程如2所示。

3 监测结果与可靠性分析

3.1 SBAS-InSAR与DS-InSAR监测结果对比

本研究分别利用传统的SBAS-InSAR方法和图2所示的基于DS-InSAR的时序地表形变监测方法提取了2016年9月—2018年6月研究区地表沉降的时空分布信息。两种方法获取的研究区地表卫星视线向的年平均沉降速率如图3所示。传统的SBASInSAR方法共提取到了57 803个监测点,监测点数目较少,且空间分布极不均匀。SBAS-InSAR方法监测到研究区地表存在一处较为明显的下沉盆地,最大年平均沉降速率出现在该下沉盆地的中心位置,为139 mm/a,在该下沉盆地的南部同样监测到一处较为明显的沉降区域,但由于在该沉降区域内获取的监测点数量稀少,只能大致反映该区域的外围轮廓,无法得到其整体的沉降信息。DS-InSAR方法共提取到了239 750个监测点,与传统的SBAS-InSAR方法相比该方法提取到的监测点数量显著增加,约为前者的4倍,且监测点空间分布均匀。DS-InSAR方法监测到研究区地表存在南北方向相邻的2个明显的下沉盆地,且在2个下沉盆地中北部下沉盆地的沉降更为明显,沉降范围也更大。在北部下沉盆地的中心位置监测到了研究区域内的最大年平均沉降速率,其值为179 mm/a。

图2 基于DS-InSAR的时序地表形变监测方法数据处理流程Fig.2 Data processing flow of time series surface deformation monitoring method based on DS-InSAR

图3 研究区地表卫星视线向年平均沉降速率(单位:mm/a)Fig.3 Annual average subsidence rate of LOS direction in study area

结合卫星影像底图对比分析两种方法的监测结果可知,受西部矿区高原堆积型沙丘地貌的影响,传统的SBAS-InSAR方法只能在矿区地表裸露岩石较多的区域获取足够数量的高相干点目标,在积沙区域由于地物的后向散射特性较弱获取的观测点数量稀少,监测结果难以准确反映矿区地表整体的沉降信息,不能对沉降区域范围进行准确探测。本研究基于DS-InSAR的时序地表形变监测方法有效克服了传统SBAS-InSAR方法在西部矿区地表沉降监测中的不足,在充分利用裸露岩石较多区域的高相干点目标的基础上,通过提取研究区域内的分布式目标显著增加了积沙区域观测点的数量,从而能够对矿区地表形变进行更高精度地监测,获取矿区地表整体丰富的沉降细节信息。

3.2 DS-InSAR时序监测结果

以SAR影像数据集中首期影像的日期为开始时间,计算其余影像对应时间相对于开始时间的累计沉降值,将观测结果从卫星的视线方向转换为垂直方向,即可得到矿区地表的时序累积沉降。为更加直观地显示整个监测时段内研究区地表沉降区域的形成与地下工作面开采之间的关系,将地下工作面的位置叠加到时序累积沉降图上,并提取出研究区域内受采动影响地表发生沉降的区域,即可得到沉降区域的累计沉降及沉降细节信息。研究区地表沉降区域范围如图4所示,沉降区域内的累计沉降及沉降细节信息如图5所示。

图4 研究区地表沉降区域范围Fig.4 Scope of ground subsidence in the study area

图5 研究区地表累计沉降及沉降细节信息Fig.5 Surface cumulative subsidence and sunsidence details in the study area

由图5可知:地表发生沉降的区域与地下工作面的位置一致。17、18两个工作面的开采导致研究区地表北部下沉盆地形成,本次试验监测到的最大累计沉降出现在该下沉盆地的中心位置,为402 mm。受06A工作面开采的影响,研究区地表形成了南部下沉盆地,该盆地内监测到的累计沉降较小,其中心位置的最大累计沉降约220 mm。监测结果表明,针对西部矿区高原堆积型沙丘地貌的基于DS-InSAR的时序地表形变监测方法在研究区域内获取了更多数量、空间分布均匀的监测点目标,能够对矿区地表开采引起的沉降区域的形成和发展过程进行直观地描述,并准确地显示地表沉降区域的范围。

3.3 DS-InSAR监测结果可靠性分析

本研究选取了17、18工作面上方2条观测线在4个不同时段内的水准实测数据,通过将观测线上DS-InSAR监测结果与水准实测数据进行对比分析来验证其可靠性。观测线与工作面的位置关系如图6所示。因水准实测数据的日期与SAR影像日期并非完全一致,故将DS-InSAR监测结果按时间进行插值使其与水准实测数据的日期保持一致。本研究在2条观测线上共选取了93个监测点目标,观测线上DS-InSAR沉降监测结果与水准实测数据对比如图7所示。

图6 观测线与工作面位置关系Fig.6 Position relationship between observation lines and working faces

由图7可知:DS-InSAR监测结果与水准实测数据拟合较好,除极个别点该方法监测结果与水准实测数据偏差略大外,各观测线上二者得到的累计沉降值及沿观测线的沉降趋势基本一致。为更加直观地体现两种方法监测结果间的相关性与误差分布,绘制了如图8所示的相关性及误差分布图,并通过计算每条观测线上的误差均值、绝对误差均值及标准差,对监测结果的精度进行评价,结果见表3。

图7 DS-InSAR监测结果与水准实测数据对比Fig.7 Comparison of DS-InSAR monitoring results and leveling measured data

图8 DS-InSAR监测结果与水准实测数据相关性及误差分布Fig.8 Correlation and error distribution between DS-InSAR monitoring results and leveling measured data

由图8和表3可知:两种方法监测结果间的相关性较高,皮尔逊相关系数为0.97,二者间的误差较小且整体分布合理,证明针对西部矿区高原堆积型沙丘地貌的基于DS-InSAR的时序地表沉降监测方法在研究区地表非充分采动情况下的监测结果具有较高的可靠性。

表3 DS-InSAR监测结果精度评价Table 3 Accuracy evaluation of DS-InSAR monitoring results

4 沉降监测结果分析

4.1 地表沉降与地下开采关联性分析

结合采工图和DS-InSAR时序监测结果可知,2016年9月17工作面开始掘进,但此后较长一段时间内地表并未监测到沉降,这是由于工作面的采深大且上方覆有大厚度的中生代砂岩,工作面开采影响尚未传播到地表。地表开始发生沉降时,工作面推进距离已达采动距,之后随着工作面继续推进,地表沉降区域的范围和累计沉降值不断增大。DS-InSAR监测结果显示,2016年11月份左右研究区地表开始出现沉降,12月份地表已经形成一处较为明显的沉降区域,2017年5月份17工作面开采结束后,沉降区域发展成为一个明显的下沉盆地,其形状近似为圆形,盆地内各点的下沉值随着与盆地中心距离的增加而减小。2017年6月,17工作面西邻的18工作面开始回采,由于两个工作面之间相距很近,18工作面的开采并没有在矿区地表形成新的下沉盆地,而是使原有下沉盆地的范围和累计沉降值进一步增大,但受两个工作面之间留设的120 m宽保护煤柱以及上覆岩层应力状态的影响,地表仍未达到充分采动且下沉极不充分。2017年12月份18工作面回采结束,其平均采厚达5.47 m,但直至2018年6月监测结束该下沉盆地内监测到的最大沉降也仅为402 mm。

06A工作面从2017年9月开始掘进,因此根据2017年12月地表沉降监测结果可知,除了在开采已经结束的17、18工作面上方存在一处下沉盆地外,在其南部也出现了一处较为明显的沉降区域。2018年4月06A工作面推进结束,到2018年6月该工作面上方已经形成了一处明显的下沉盆地,且其沉降范围已经与北部下沉盆地相当。通过比较南北两个下沉盆地的形成和发展过程不难发现,从地表开始发生沉降到形成明显的下沉盆地过程中南部下沉盆地所用时间要明显短于北部下沉盆地。这是因为虽然3个工作面的采深和上覆岩层的性质十分接近,但06A工作面采用综合机械化放顶煤工艺开采,其采厚达9.78 m,远大于17、18两个工作面,因此其开采影响传播到地表所用的时间更短,但受该矿区特殊地质采矿条件的影响,直至监测时段结束该工作面上方地表同样未达到充分采动。

4.2 南部沉降区地表动态变形分析

为进一步验证基于DS-InSAR的时序地表沉降监测方法的性能,本研究利用该方法获取了图9所示06A工作面开采期间南部沉降区地表沿该工作面走向和倾向主断面的动态沉降曲线,对该区域缺少的水准实测数据起到了很好的补充作用。为充分反映工作面推进过程中地表的动态变形情况,将沉降监测结果进一步处理得到了沿06A工作面走向主断面的动态倾斜与动态曲率曲线,分别如图10和图11所示。由图10、图11可知:在06A工作面推进过程中,沿该工作面走向的倾斜和曲率变形不断增大,至2018年4月底开始趋于稳定,倾斜在20号点附近(采空区中心)近似为0,在8号点和26号点附近(采空区边界)达到最大,为0.41 mm/m;地表曲率变形较小,均小于0.02 mm/m2。

图9 06A工作面走向与倾向主断面动态沉降曲线Fig.9 Dynamic subsidence curves of main section of strike and incline of 06A working face

图10 06A工作面走向主断面动态倾斜曲线Fig.10 Dynamic inclination curves of main section of strike of 06A working face

图11 06A工作面走向主断面动态曲率曲线Fig.11 Dynamic curvature curves of main section of strike of 06A working face

5 结 论

本研究针对西部矿区高原堆积型沙丘地貌,以内蒙古石拉乌素煤矿为例,采用基于DS-InSAR的时序地表形变监测方法对2016—2018年间覆盖研究区域的52景Sentinel-1A影像数据进行处理,获取了研究区域的地表时序沉降信息,并将监测结果结合井下开采资料进行了分析。所得结论如下:

(1)与传统的SBAS-InSAR方法相比,本研究方法在充分利用高相干点目标的基础上,依靠FaSHPS方法对研究区域内的分布式目标进行提取,并通过特征值分解对其相位进行优化,之后将二者联合进行时序处理来获取地表形变信息,可有效提升研究区地表监测点密度和形变监测精度。

(2)在监测时段内,该方法监测到研究区地表形成了南北相邻的2处明显的沉降区域,最大沉降约400 mm。与井下开采资料对比分析表明,该方法监测到的地表沉降区域的位置和范围较为准确,且沉降区域的形成及发展能与地下工作面开采进程相对应,可为开采沉陷影响边界确定及其机理研究提供参考。

(3)验证了基于DS-InSAR的时序地表形变监测方法的可靠性,并利用该方法获取了矿区06A工作面开采期间南部沉降区域沿06A工作面走向和倾向主断面的动态沉降曲线,对该区域缺少的水准实测数据起到了很好的补充作用,并将沉降监测结果进一步处理,得到了沿06A工作面走向主断面的动态倾斜和曲率变形曲线,可为地表动态变形规律研究及形变控制提供依据。