基于InSAR三维分解技术的贵州省贞丰某煤矿地下采空区反演*

2020-09-10姚佳明吴作启蒋永亨刘星洪

工程地质学报 2020年4期

关键词：反演采空区煤层

姚佳明姚鑫吴作启陈剑蒋永亨刘星洪

（①活动构造与地壳稳定性评价重点实验室，中国地质科学院地质力学研究所，北京100081，中国）

（②中国地质大学（北京）工程技术学院，北京100083，中国）

（③煤炭科学技术研究院有限公司，北京100013，中国）

0 引言

贵州贞丰县煤矿采空区众多，小煤窑的开采易造成地表塌陷、滑坡等一系列地质灾害。而且在一些煤炭资源丰富的地区非法采矿事件时有发生，不仅严重破坏矿产资源和生态环境，易诱发斜坡变形开裂甚至失稳（史文兵等，2016），造成矿山灾害和巨大的经济损失，威胁安全生产和社会稳定。由于矿山分布广泛，地下非法采矿活动隐蔽性强，仅依靠矿山执法部门的调查，很难及时、准确地发现这些行为，需要耗费大量的时间、精力和财力（Xia et al.，2018）。了解地下开采情况可以对将来可能由采矿引起的地表破坏进行防护，减少生命财产损失，并且对煤矿越界开采起到一定监测作用。我国煤矿采空区探测的技术方法较多，比如地球物理方法、钻探方法、三维激光扫描等。地球物理方法利用波在不同地层和岩石的弹性、速度和运动方式差异来确定地下采空区（吴成平等，2007），然而物探类方法受环境噪音、适用性、分辨率、多解性、不确定性等多种条件的限制，无法做到准确定量预报（王升阳等，2015）；钻探在工程施工、探测方面应用较多，但施工难度大、成本高，探测成果受钻孔数量、钻孔密度、配套的监测、测试方法等限制，分辨率较低，无法查清周边采空区准确边界（薛国强等，2018）；三维激光扫描技术可以高精度、高效率地测量出地下采空区，但性价比较低（薛国强等，2018）。对于已关闭或已废弃煤矿，上述方法虽然可以在一定精度上测量出地下采空区范围，但需要消耗大量人力和财力且有一定的作业危险性。

近年来很多学者将InSAR技术应用于监测与地球或人类活动有关的地表变形（Amelung et al.，1999；Du et al.，2017；Feng et al.，2015；Milillo et al.，2016a，2016b；Xu et al.，2016a，2016b；Zhao et al.，2016；Yang et al.，2017a，2017b；Blachowski et al.，2019），与传统的地面水准测量技术、GPS技术相比，InSAR技术具有全天时、全天候、观测范围大、灵敏性高、高性价比、观测时间可回溯等优势，并能穿透某些地物表面（吴立新等，2004）。煤矿开采产生的大区域变形可以使用InSAR技术获取，在国内外已有大量学者将InSAR技术应用于矿区地面塌陷监测，并通过实地GPS等地面测量验证了InSAR计算结果的准确性。朱建军等（2011）介绍了InSAR技术在矿区地表形变监测中的应用现状及进展，分析了D-InSAR技术相比于传统测量手段的优势；Xia et al.（2018）采用D-InSAR采矿沉陷监测方法，建立了地表变形与以沉陷为特征的地下开采的时空关系模型；Fan et al.（2015）采用相叠加和亚像素偏移跟踪相结合的方法，对采动引起的沉陷机理进行了研究；Wang et al.（2018）结合D-InSAR和概率积分法提出了一种基于单视距D-InSAR的概率积分参数模拟沉降的计算方法。

煤矿地下采空与地表破坏存在一定变形规律，其变形量、变形范围与地下开采方式、煤层埋深、开采煤层厚度、地层倾角、岩层组合、地表地形、坡面与地层的组合关系以及开采工作面的布置方向有直接关系。通过实地调查发现，研究区内地质构造简单，地下煤层展布特征基本一致，顺层缓倾岩层以砂泥岩为主，通过周边煤矿调查得知该区域开采方式一致皆采用下行式开采，且地下开采方向较为一致（自西向东开采），故在以上条件限制下可以根据该地区开采沉陷参数进行地下采空范围反演，所选参数参照于《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》（以下简称《三下采煤规范》）中的边界角概念，建立InSAR短期动态监测煤矿的影响角计算方法，本文目的是通过该方法计算地下采空地表变形趋于稳定之后，在开采时间与稳定时间段内的短期地表变形与地下采空变形间的影响角规律，而后通过计算获取研究区已知煤矿地下开采与地表沉降变形关系，来推算反演已废弃煤矿的地下开采范围及时间。

1 研究区概况

研究区位于贵州省黔西南布依族苗族自治州贞丰县境内，坐标范围为北纬25°20′～25°21′，东经105°31′～105°33′，变形区的坡度10°左右，坡向近北偏西30°（图1）。

研究区地层从石炭系至第四系均有出露，其中以三叠系地层最为发育。贞丰县主要含煤地层为三叠系上统火把冲组（T3h），由海陆交互相和冲积相含煤沉积组合形成，也是本研究区煤矿主要开采煤层（图2）。T3h地层由石英砂岩、粉砂岩、黏土岩、碳质黏土岩及数十层煤层（线）组成多个沉积旋回，与下伏地层把南组整合接触，矿区内仅出露至火把冲组第三段（T3h3）：浅灰、灰、黄灰绿色等色厚层，块状粗至细粒石英砂岩、细砂岩、黏土岩组成的不等厚韵律层，砂岩中具有斜层理。该区含有可采煤层两层，编号为JK、K2。JK煤层厚0.70～0.80 m；K2煤层平均厚1.08 m，煤层倾角平均为10°，倾向290°～350°，矿区内构造比较简单，岩层为单斜构造。矿区断层较发育，有柏枝树断层、大榜断层（图1红实线），但变形区与断层无接触关系且相距较远，故排除受到断层影响。

经过实地考察，受地下煤层和当地条件限制，研究区内矿井设计采用下行式开采，煤层开采顺序为K2、JK联合开采，平均层间距8～10 m，煤层倾角6°～15°，平均倾角12°，采深约250 m，采厚1.8 m，工作面走向长度420 m，倾向长度100 m，开采面积约42 000 m2。

2 研究与计算方法

2.1 开采沉陷参数计算与地下采空区反演

图1 研究区位置与地质图Fig.1 Location and geological map of the study area

图2 研究区地质剖面图（位置见图1剖面线I－I′）Fig.2 Geological profile of the study area

本文研究区变形为已知煤矿开采和未知煤矿开采引起的两部分变形。根据对本文研究区附近煤矿的调查，由于开采煤层与地形因素的影响，该地区煤矿开采方式皆为下行式开采，地层产状、开采煤层、厚度等和邻近煤矿开采条件较为一致。可以根据已知煤矿区计算的沉陷参数与变形时间滞后关系，反演未知区煤矿地下开采时间和范围。

根据各时段地表LOS形变图判断地下掘进方向与地表变形时间以及在变形时段地表的累计变形，进行三维分解计算来获取地表垂向变形信息。研究从两个方面进行（图3）：（1）在获取已知煤矿地下开采资料和地质条件的情况下，根据已知煤矿信息与地表变形信息计算出该区域开采沉陷参数（上、下山及走向影响角）和变形滞后时间；（2）根据已知煤矿开采计算获得的参数和未知废弃煤矿开采引起的地表变形反演其地下采空范围和地下开采时间。

2.2 计算方法

本文使用D-InSAR方法的计算结果，对煤矿采空区三维地表变形进行反演，结合数值模拟技术对反演结果进行核查验证。本文选取了日本JAXA宇宙航空研究中心研发的L波段ALOS-PALSAR2（以下简称P2）卫星条带模式数据，升、降轨各8景SAR影像。数据的斜距向和方位向分辨率分别为1.430 422 m和2.217 833 m，影像参数如表1。对比其他数据，本文所选的数据波长较长、分辨率较高、穿透力强，更有能力识别出矿区地表变形而不出现失相干等现象，适用于矿区地表变形监测与分析。

使用D-InSAR方法分别对升、降轨影像做了差分干涉处理，为了排除地表大梯度变化带来的失相干，时间基线就仅选取36对相邻的影像对，来减少时间去相关造成的影响。处理过程中为了保持较高分辨率和信噪比，我们选择了2×2的多视比，采用偏移多项式对影像进行配准，使用美国SRTM的30 m分辨率DEM去除地形相位。对干涉图像使用最小费用流法（mcf）进行解缠，居民区较为稳定，故选择煤矿东北侧居民区为解缠起算点，干涉时间基线图如图4所示。影像覆盖区域相对于河谷区高差较小，且研究区范围内河流较少，水汽含量均一，大气误差较小，误差源主要为轨道误差。故本文在处理D-InSAR过程中，在两次移除轨道误差后，高程带来的大气误差的去除和适应性滤波去除湍流大气均采用大窗口滤波，减少对矿区实际变形量的影响，最后投影到UTM坐标。

煤矿区地表变形为水平向与垂直向的三维变形，而D-InSAR计算出的结果为视线向（LOS）变形，若忽略水平运动而直接将LOS变形值进行转化，求得的垂直向变形会存在很大误差。国内外很多学者将InSAR三维分解技术应用于形变监测与地表变形（王晨兴等，2011；Li et al.，2015；Yang et al.，2017a，2017b；Huang et al.，2017；李凌婧等，2017；Liu et al.，2018；Yang et al.，2018；Zheng et al.，2019）。Hu et al.（2014）详细介绍了目前国内外学者在InSAR三维分解技术方面的发展情况，分别分析了各种方法的计算原理及其优势与劣势。通过利用从至少3个成像几何形状获得的多个InSAR测量，我们理论上可以将位移矢量扩展到3D（吴立新等，2004）。本文使用的方法是基于升、降两轨计算方法，利用两轨D-SAR影像干涉结果，根据升、降轨雷达成像几何原理进行三维变形求解，因为只有升、降两个方向的SAR影像，对于模糊度较高的N－S方向变形量不易准确计算，所以本文没有N方向变形量分解，计算公式如下：

式中：du、de为垂直向、东西向变形量；Г为转换系数；dlos，1、dlos，2为升、降轨两干涉对变形值；ϑinc，i、αaz，i分别为入射角和方位角。

3 结果与分析

图4 升、降轨基线图（升－左、降－右）Fig.4 Baseline diagram of ascending and descending（ascend-left，descend-right）

通过两轨干涉变形信息判断，研究区变形可以分为两部分，由两个煤矿开采诱发：一部分为南侧已知煤矿（已获取地下开采资料）开采诱发的地表变形；另一部分为北侧（已废弃，无地下开采信息）未知煤矿开采带来的地表塌陷变形（图5）。通过南侧煤矿开采信息（地下开采时间、开采范围）和地表形变信息来计算开采沉陷参数，再根据未知煤矿地表变形信息反演地下开采时间和范围。

3.1 已知区域开采沉陷参数和变形滞后时间计算

图6为升、降轨在已知煤矿开采时间地表各时段地表变形信息，地下煤矿为自西向东开采，开采时间在2018年1月开始，2018年6月结束。两个轨道的干涉影像时间不同，有助于形变区变形时间确定，在升、降轨共6个变形时间段内，地表变形区随时间的推移由西移动至东侧，与地下开采方向一致，地表变形的时间与地下开采的时间存在一定延迟，滞后时间规律影响因素较多，在不同地区无法进行判断与分析，但在地质条件、煤层条件、开采方式与方向等相同的邻近区域，滞后时间较为相近。

3.1.1 D-InSAR三维分解

根据煤矿开采及前期钻探资料，已知煤矿区开采时间段为2018年1月至6月，为覆盖变形时间段，本文的升、降轨数据分别选用了20171112～20180819和20171206～20180620共四景影像用于两两干涉处理。通过两个干涉对进行三维分解解算，计算出地面垂向变形量。如图7所示，在煤层开采的走向和倾向方向分别对垂向变形做变形剖面图，开采面中心300 m长度范围内沉降量在20 cm以上，最大沉降30 cm，周边区域变形量递减。

图6 煤矿开采过程多时间节点地表变形Fig.6 Known surface deformation diagram in each period of coal mining time

图7 三维分解的垂向变形（a为已知煤矿地表变形垂向变形）Fig.7 Vertical deformation of three-dimensional decomposition（a shows the vertical deformation of known coal mine surface deformation）

3.1.2 开采沉陷参数计算

参考《三下采煤规范》中的边界角概念，建立以地表InSAR监测10 mm变形边界与地下开采边界连线与水平线的角度为InSAR监测煤矿开采变形影响角，本文根据三维变形中的垂向变形量，结合地下煤矿开采边界对已知煤矿开采区的走向，上、下山影响角做了计算（图8），计算得出上山影响角75°，下山影响角80°，走向影响角83°。

3.1.3 变形时间滞后关系

地下采空反映到地表形成塌陷需要一定的时间，开采条件、地质条件、开采方向等因素皆会影响变形滞后时间，而在本研究区上述条件较为类似，根据SAR影像干涉对时间节点的地表变形和地下开采时间，在本研究区做如图9地表变形时间与地下采空时间的关系曲线，两条曲线在横轴（时间轴）方向的差值或偏移即为地表变形滞后时间，根据曲线走向趋势可以判断出在研究区的地质和开采条件下，地下采空至地表变形处1个月的时间滞后关系。

3.2 未知区地下采空区反演

图10为升、降轨不同时间段的干涉影像，地表变形移动与煤矿采空引起的变形一致，两轨影像在近乎相同时间内的地表形变量、形变区大小基本一致也验证D-InSAR结果的准确性。通过两轨干涉图地表变形信息分析得出，该地区煤矿自西向东开采，变形时间在2017年4月26日至2017年10月25日之间。随时间的推移，地表变形区也随之向东移动，该地区属干旱地区，年均降雨量少、植被稀疏，地下水位变化不大；在历史SAR干涉影像中未开采时段该地区无变形产生，且该地区地表坡度、地下岩层较缓，不易发生滑坡，故也可以排除因地下水位变化、滑坡等引起的地表变形。

从图10中8景干涉对中得出变形时间大致可以确定为：2017年4月末至2017年10月末。在考虑影像获取时间与误差允许范围内，我们对该区升轨20170514～20171112和降轨20170426～20171025两个相近时间段做三维分解解算。

3.2.1 D-InSAR三维分解

图8 走向、上山、下山影响角计算示意图（横纵轴1:1）Fig.8 Schematic diagram for calculating the ascending and descending influence angles（horizontal and vertical axis 1:1）

图9 地下采空与地上沉陷时域对比图Fig.9 Time-domain comparison of underground mining and overland subsidence

图10 a～d为升轨干涉图，e～h为降轨干涉图Fig.10 a～d shows the ascending orbital interference pattern，e～h shows the descending orbital interference pattern

获取的升轨影像入射角 ϑinc，i、方位角 αaz，i分别为39°、－10°，降轨影像入射角 ϑinc，i、方位角 αaz，i分别为43°、－170°，经过2×2重采样并地理编码后的P2 SAR影像距离向和方位向分辨率分别为4.5 m和4.4 m，故三维分解的垂直向和东西向分辨率也为4.5 m和4.4 m。对升、降轨两个干涉对进行等尺寸裁剪获得1290×1290的干涉对，根据如上参数和公式计算得出了垂直和东西向变形量。

在图11a、图11b以及5条剖面线结果中看出，图11a为垂直向形变结果，整体变形较为集中，最大形变量在区域中心为0.5 m，周围大部分区域变形为0.2～0.4 m左右；图11b为东西向形变结果，红色值和蓝色值分别代表运动向东和向西，从EE’剖面线中可见明显分界线。中心左侧整体为向东运动，符合采矿引起的地表塌陷变形特征。最大位移值0.5 m，中心右侧整体向西运动，最大位移值0.36 m。

3.2.2 采空区边界反演及开采时间推算

本文在废弃煤矿采空沉降中心区域建立5条剖面线来反演地下开采范围。在开采倾向（近自西向东）方向的4条剖面线最大沉降、10 mm沉降边界位置以及沉降趋势大致一致，故取4条剖面线平均值来反演开采倾向边界；主断面上地表移动盆地的范围最大，移动最充分、移动形变量大，故在开采走向方向最大区域变形处拉取一条剖面线作为走向边界的反演，根据已知影响角、地下开采深度、煤层产状以及地表高程和形变信息通过三角函数几何计算出地下开采边界（图12，图13）。

综合走向、倾向确定的开采范围，可以确定开采面位置如图14。

根据地表变形时间和滞后关系，本文对地下开采时间进行了推测（图15）。在邻近已知矿区的地下开采与地表破坏30 d左右的时间滞后关系判断，地表破坏对应的地下采空时间应为地表变形时间的1个月之前，在两轨D-InSAR各时段节点变形结果中做出地表变形在煤层开采走向上的时间－距离曲线（蓝、黄虚线），两曲线变形时间－距离曲线基本吻合，再根据滞后时间推算出地下煤矿的开采时间曲线（红色虚线）。

3.2.3 采掘范围数值模拟验证

根据该区地质岩层剖面与煤层信息，在MIDAS中建立简化FLAC3D计算模型（图16），计算模型范围包含了未知采空区，尺寸为680 m×400 m×300 m，模型上表面为该地区30 m分辨率DEM双线性重采样至5 m分辨率构成，依据当地煤矿钻孔勘察报告该地区煤层上覆地层主要为泥质砂岩与泥质粉砂岩，为方便计算地层划分过程中忽略了较薄的泥岩、粉砂等地层；在煤层位置设置东西向长380 m、南北向长150 m、倾向330°、倾角10°的开采煤层，具体分层见图16。

图11 基于D-InSAR的升、降轨三维分解Fig.11 D-InSAR based three-dimensional decomposition of ascending and descending orbits，where Fig.11a represents vertical deformation and red integrity represents settlement；Fig.11b is the east-west deformation，the red negative value represents the east，and the blue represents the west；A－A′equilinear lines are shaped variable profile lines；the NaN white areas are incoherent and without data areas

图12 开采煤层倾向方向的地下开采区反演Fig.12 Inversion of underground mining area in mining seam inclination direction

图13 开采煤层走向方向的地下开采区反演Fig.13 Inversion of underground mining area in strike direction of coal seam

计算模型左、右边界x＝0im和x＝680 m两个垂直面仅约束边界面法向位移，平面内无约束；前、后边界y＝0 m和y＝400 m两个垂直面仅约束边界面法向位移，平面内无约束；模型底部z＝－300 m水平边界采用固定约束；水平地表面为自由面，收敛标准为体积最大不平衡力与典型内里的比率小于10－5。岩体采用六面体八节点单元模拟，本构模型采用基于弹塑性理论分析的莫尔－库仑模型。研究区降雨较少，开采面上无地下水赋存，未考虑孔隙水、静水压力，岩体物理力学参数选取见表2。

图14 开采区位置图Fig.14 Location of mining area

图15 地下开采时间推测Fig.15 Estimation of underground mining time

表2 岩体物理力学参数表Table 2 Physical and mechanical parameters of rock mass

计算过程如下：

（1）依据弹性求解法，利用莫尔－库仑模型计算初始地应力，以不平衡力与典型应力的比值小于1×10－5为收敛条件判断初始应力平衡，以模型角点（5，5，－250）与除顶面所有面为固定边界，并固定底面x、y方向速度。

（2）在走向、倾向两个主要变形线设置34个监测点用于变形量验证。

（3）煤层开挖，计算步长设置为2000步，保证不平衡力收敛时足够小，使模型内部应力状态充分调整。

计算结果如下，数值模拟结果显示变形区主要集中在开采区上方地表附近，其他地区明显变形，开采中心上方形变最大，垂向沉降达30 cm左右，变形特征与煤矿开采沉陷特征相似，中心沉降最大，四周呈梯度减小（图17）。

在数值模拟过程中对走向、倾向剖面线设置了34个监测点用于变形量的检验。与D-InSAR结果对比变形区内南北、东西主要两条主要剖面如图18所示，变形量总体变形趋势一致，但有较小差别，在倾向、走向沉降剖面的差值均方根都在6 cm左右（图19），在矿区大变形的条件下6 cm误差可以接受，根据数值模拟对未知采矿区变形的计算结果来看，两者大致趋势保持一致，变形值相差不大，标准方差较小，计算结果较为可靠。

3.2.4 野外验证

图16 FLAC3D计算模型Fig.16 FLAC3D calculation model

图17 数值模拟垂向形变量Fig.17 Numerical simulation of vertical shape variables

我们围绕采空区及附近山体进行了野外现场调查，调查发现采空区上方地表存在多处挤压和拉张裂缝以及大规模的岩土体变形塌陷现象，塌陷变形区主要在开采区上部，与InSAR计算沉降区位置基本吻合；在开采区四周，有多条地面水平运动产生的拉张裂缝（图20），由地面不均匀水平运动产生，证实了计算结果的合理性和可靠性。

4 讨论

本文只针对研究区InSAR观测周期内进行分析讨论，煤矿采空造成的地表塌陷影响因素很多，包括地下开采方式、煤层埋深、开采煤层厚度、地层倾角、岩层组合、地表地形、坡面与地层的组合关系以及开采工作面的布置方向等，且开采沉陷参数在各个地区各不相同，故本文只针对贵州贞丰县地区对煤层采空区反演。反演不仅需要该区地表垂向形变，还需要地层信息、煤层深度、煤层产状等信息，且开采沉陷参数获取来源为该区附近煤矿，仅对研究区煤矿有一定应用意义。由于该采区煤矿已经废弃，无法获得实际地下采空资料，该区地表也无水准设备做GPS测量，故本文应用数值模拟技术进行了验证，验证结果与反演结果大致相同，一定程度验证了反演的准确性。另外，由于研究区植被覆盖、地物反射特性、季节因素等，不同SAR数据的波段选择、分辨率、极化方式等皆会对计算结果产生影响。本文研究区地表以裸露砂岩或草植被覆盖砂岩为主，而采用的L波段3 m分辨率数据对植被穿透力较强，且地物反射信号较高，计算数据结果较为可靠。

本文误差来源主要有：D-InSAR计算过程中对较大变形区的解缠可能会丢失相位信息造成变形量结果计算的误差；两轨D-InSAR三维分解使用的两对干涉对的时间有一定差别，造成煤矿采空变形周期中某轨形变量的缺失，对三维变形量分解的结果产生影响；本文采用的开采沉陷参数为定义的InSAR监测地表10imm沉降变形曲线计算出的影响角，未使用以0 mm沉降量边界计算的移动角，而开采沉陷边界参数的选取可能对地下采空区的反演结果有一定差异；因煤矿变形一般2 a左右才会达到沉降的最大值，从SAR数据量积累和InSAR变形计算所允许时间基线的角度来看，都无法实现矿区最大沉降量的观测计算，直接使用《三下采煤规范》中边界角概念可能会影响计算准确度，故本文仅能在矿区地表短期沉降稳定后获取该时段沉降，并计算新定义的煤矿开采沉陷影响角；反演结果的验证使用数值模拟技术，该区地层条件比较复杂，故在数值模拟中进行了模型简化，不同地层的物理力学性质不同，可能造成一定的误差。

图18 数值模拟与D-InSAR三维分解对比图Fig.18 Numerical simulation and D-InSAR three-dimensional decomposition comparison

图19 数值模拟与D-InSAR三维分解误差柱状图Fig.19 Numerical simulation and D-InSAR 3D decomposition error histogram

图20 野外验证地表变形破坏Fig.20 Field verification of surface deformation and failure

Offset-Tracking技术可以应用于矿区变形测量中，其应用于变形量大、变形范围大的地区有较大优势，其计算结果与窗口选取、滤波技术的使用有较大关系，误差较大，但是本文研究区变形较小，变形范围较小，而且本文为L波段3 m高分辨率数据测量精度较大D-InSAR计算结果会更加精确，故未使用Offset-Tracking技术进行变形量计算。