基于InSAR技术的缓倾煤层开采诱发顺层岩体地表变形模式研究

2020-06-02姚佳明李凌婧任开瑀刘星洪

水文地质工程地质 2020年3期

姚佳明，姚鑫，陈剑，李凌婧，任开瑀，刘星洪

(1.活动构造与地壳稳定性评价重点实验室/中国地质科学院地质力学研究所，北京 100081；2.中国地质大学(北京)工程技术学院，北京 100083)

合成孔径雷达干涉测量法(InSAR)是一种测量天然或人为造成的地表形变的方法，在国内外均有较多应用。Schmidt D A等[1]在2005年利用13景ERS SAR影像干涉叠加得到了美国Hayward断层1992—2000年卫星距离向变形速率，并结合地面GPS数据对断层滑移速率进行了反演；Beavan J等[2]利用合成孔径干涉雷达观测资料估算了新西兰7.8级地震同震场和震后早起滑移分布；Chen F L等[3]结合C波段和L波段SBAS技术，分析了青藏铁路沿线冻土变形与基础设施的相互作用关系；Raucoules D等[4]使用ERS干涉雷达监测法国Vauvert盐矿地表缓慢变形，计算得到了可靠的沉降盆地范围，也证明InSAR技术可以应用于大范围地表沉降测量；Strozzi T等[5]对瑞士阿尔卑斯地区变质基岩中引水隧洞上方地表沉降和隆起进行了探测[5]；Jiang L M等[6]以中国北方的乌达为研究对象，分析SBAS技术在监测地下煤矿火灾引起地表变形的应用前景；赵超英等[7]利用ALOS PALSAR数据对大面积滑坡进行探测和监测。D-InSAR利用两幅相关合成孔径雷达(SAR)图像的相位差，能够在数公里范围内探测到精度在几毫米以内的地表位移[8]。近年来，这些方法被广泛应用于地面位移测量和三维地表信息的有效检索[9]，其原理是通过自主发射和接收微波来探测地表，通过干涉相位变化测量地表变形[10-11]。目前，D-InSAR利用遥感卫星多时相的复视雷达图像相干信息,进行地表垂直形变量的提取，其精度已经达到了毫米级[12-13]。

SAR应用于地球观测的首批研究之一是马森内特等人的研究，其中南加州上空的一幅差分干涉图记录了1992年里氏7.3级兰德斯地震的同震位移。几年后，Sandwell和Price引入了一种叠加技术(Stacking technique)，该技术将干涉对组合成一个具有更好信噪比的产品。几年后，Ferretti等人引入了永久性散射体(PS)方法，该方法将SAR的适用性扩展到原本在干涉测量上低相干的植被密集区域。最后，Berardino和Usai引入了小基线子集(SBAS)方法[14-15]，该方法通过将所有符合时间、空间基线的SAR数据集合利用最小二乘法求解每个小集合内的地表形变时间序列。Hooper将PS和SBAS技术结合到一个方法中计算最小二乘解来产生地面变形的时间序列，经过国内外学者的对PS-InSAR技术的研究，目前利用PS-InSAR计算得到滑坡的时序变形曲线对滑坡垮塌预警起到一定作用[16]，此外还提出了各种修改，其中包括用于研究三维变形的开发[17-18]。

SAR是一种使用微波探测地表目标的主动式成像传感器，与传统的地面水准测量技术、GPS技术相比，InSAR技术具有全天时、全天候、观测范围大、灵敏性高、高性价比、观测时间可回溯等优势，并能穿透某些地物表面[19]，同时InSAR与地面水准数据的结合使InSAR计算结果提高了可靠性[20]。InSAR在大范围测量地面沉降、地质灾害隐患的变形有大量应用[21-23]，而且对于大范围的采矿采空区地质灾害快速识别地质灾害区具有重要意义[24]。

近年来，贵州省矿山引发的地质环境问题被重新分析与研究[25-26]，InSAR技术也被越来越多地应用到矿区地表变形监测中[27-30]。本文以贵州省贞丰县某煤矿作为研究区，利用15期3 m空间分辨率L波段升、降轨PALSAR-2 SAR为数据源，在多期D-InSAR测量地表变形的基础上，分解计算了采空三维变形，进而通过与地下开采过程的相关性分析，揭示了采矿滑坡变形模式。

1 研究区概况

研究区位于贵州省黔西南布依族苗族自治州贞丰县境内西南方向，处于贞丰县与安龙县之间，距贞丰县城区9 km(图1)。

图1 研究区位置、SAR数据覆盖范围与岩性图

研究区地层从石炭系至第四系均有出露，其中以三叠系地层最为发育(图2)。矿区内地层发育特征如下：石炭系(C)主要分布于煤矿区北部，岩性为浅海相碳酸盐岩、碎屑岩及硅质岩。二叠系(P)分布于矿区北部，下二叠统岩性为开阔台地相质地较纯的碳酸盐岩，厚度变化大，百余米至千余米，与下伏石炭系地层呈整合或假整合接触；上二叠统岩性主要为浅海相碳酸盐岩及碎屑岩，与下伏地层呈整合接触，厚200～300 m。三叠系(T)为台地型沉积岩，下统为碎屑岩及碳酸盐岩，中统以碳酸盐岩为主，上统由碎屑岩组成，自下而上由海相向陆相逐步过渡。第四系(Q)主要分布于沟谷及低洼地带，由残坡积物、冲积物等组成，厚0～25 m不等。

贞丰县主要含煤地层为上三叠统火把冲组(T3h)，由海陆交互相和冲积相含煤沉积组合形成，也是本研究区煤矿主要开采煤层。T3h地层由石英砂岩、粉砂岩、黏土岩、碳质黏土岩及数十层煤层(线)组成多个沉积旋回，与下伏地层上三叠统把南组(T3b)整合接触。矿区内仅出露至火把冲组第三段(T3h3)：浅灰、灰、黄灰绿色等厚层至块状粗至细粒石英砂岩、细砂岩、黏土岩组成的不等厚韵律层，砂岩中具有斜层理。在煤矿采空区中，含有可采煤层(主采)二层，编号为JK、K2。JK煤层厚0.70～0.80 m，在矿区稳定开采；K2煤层在西南部厚1.06～1.10 m，平均厚1.08 m，其余地区厚1.20～1.30 m，煤层在区内展布稳定[31](图3)。

研究区内构造特征表现为NW向、NWW向与NE向、近SN向交切，以NW向、NE向为主。主要褶皱构造为龙头山向斜，矿区位于龙头山向斜北段的东扬起端，轴长55 km，轴向285°，地层倾角10°～50°。其中含煤岩系组成向斜轴长约28 km，区域内地质活动较弱。矿区内构造比较简单，岩层为单斜构造，倾向290°～350°，倾角6°～15°，平均为12°。

图2 研究区遥感影像鸟瞰图煤矿回采区

图3 研究区地质剖面图(位置见图1剖面线I-I′)

矿区断层较发育，有柏枝树断层、大榜断层(图1红实线)。柏枝树断层分布于矿区中部，柏枝树至大石堡一线，长4 km，走向NE—SW向，倾向NW，倾角65°～75°，断层破碎带宽5～10 m，水平断距50 m左右，为逆断层；大榜断层分布于矿区中部大榜一线，长3 km，走向NE—SW向，在柏枝树南西1.3 km附近与柏枝树断层相交，断层倾向SE，倾角60°～65°，断层破碎带宽3～8 m，水平断距400～500 m，地层断距40 m，为逆断层。煤矿工作区与断层无接触关系，故采煤工作面未受到断层影响。

经过实地考察，研究区煤矿矿井设计采用斜井开拓方式，全矿井划分为三个采区开拓全井田(图5右)，采区内各区段为下行式开采，煤层开采顺序为K2、JK联合开采，平均层间距8～10 m，煤层倾角6°～15°，平均倾角12°。煤矿东区井田K2、JK联合开采，JK层较薄，采深160 m，采厚1.6 m，工作面平均走向长度220 m，倾向长度110 m，开采面积约24 000 m2；煤矿西北区井田K2、JK联合开采，采深200 m，采厚1.8 m，工作面走向长度280 m，倾向长度100 m，开采面积约28 000 m2；煤矿西南区井田K2、JK联合开采，采深250 m，采厚1.8 m，工作面走向长度420 m，倾向长度100 m，开采面积约42 000 m2。

2 数据与方法

2.1 SAR数据的选取

本文数据选取了日本JAXA宇宙航空研究中心研发的L波段ALOS-PALSAR 2卫星条带模式数据，其优势在于重访周期短(14 d左右)、波长较长、分辨率高，能较好地克服植被失相干、适应大变形，适合采矿地质[32]。工作区石漠化较严重，SAR信号反射强度高，InSAR相干性好。因为雷达卫星是侧视成像，因此会产生顶底倒置、阴影、叠衍等问题，研究区最大高差达1600 m，坡度陡，在SAR成像方面存在以上问题，而升降轨两种观测角度可以很好地互补，所以选用选取升轨8景、降轨7景共15景SAR数据(表1)。使用美国30 m分辨率的SRTM作为高程数据去除干涉相位中的地形相位。

表1 SAR数据基本参数

2.2 处理方法与过程

图4 PALSAR-2升、降轨基线图

本文使用GAMMA商业软件处理，D-InSAR方法第一步是根据时间和空间基线选取主影像与其他像对进行配准计算；第二步是DEM与SAR影像的配准与裁剪，采样出DEM与主影像位置一致的范围，用以去除地形相位、解缠干涉图；第三步利用配准好的影像对计算干涉相位，利用配准后的DEM计算地形相位，从干涉相位中去除地形相位、平地相位等，生成差分干涉图。最后对干涉图滤波、解缠得到解缠相位。采用8景升轨、7景降轨数据分别进行了干涉处理。升、降轨分别选择20171224景和20180228景为主影像进行配准，地表变形量大时间基线过长会影响解缠效果，故干涉对的选择使用多参考以150 d时间基线和120 m空间基线为阈值计算出来，另额外填加20170514—20170611、20171224—20180415、20170514—20180415、20170806—20180415、20170426—20171206、20170426—20180228、20170426—20180509、20170524—20171206、20170524—20180228、20170524—20180509、20170715—20180228、20170715—20180509干涉对保持时间分辨率连续(基线图如图4)。综合单视影像幅长、幅宽与方位、距离向分辨率并尽可能地保持高分辨率，确定了2∶2的多视比，居民区较为稳定，故选择煤矿东北侧居民区为解缠起算点，使用最小费用流法进行解缠。影像覆盖区域相对于金沙江河谷高差较小，且研究区范围内河流较少，水汽含量均一，大气误差较小，误差源主要为轨道误差。故本文在处理D-InSAR过程中，在两次移除轨道误差后，高程带来的大气误差的去除和适应性滤波去除湍流大气均采用大窗口滤波，减少对矿区实际变形量的影响，最后投影到UTM坐标，最终共得到25景干涉图。

上述方法存在一些限制：1)只能处理一个SAR数据集；2)只产生视距变形时间序列，在复杂信号情况下难以解释；3)产生的时间序列具有有限的时间覆盖和粗糙的时间分辨率等。

如果有多个轨道的SAR数据，可以产生地面变形的多维时间序列，包括东西向水平分量和垂直分量，同时利用SAR数据集中的所有数据，可实现不间断的时间覆盖，提高时间分辨率。通常的情况来估计地面变形是将InSAR测量值除以余弦值，即假设变形区无水平方向位移，只在竖直方向的位移量条件下获得的。然而，这种不良的假设可能会导致对实际可能发生在垂直和水平方向上运动的错误解释[33]。此外，如果表面位移处在垂直于LOS方向的方位方向上，那么表面位移就会被丢失。通过利用从至少三个成像几何形状获得的多个InSAR测量，我们理论上可以将位移矢量扩展到3D[19,21]。

3 计算结果与实际对比

3.1 时序计算结果

本文收集了煤矿2017—2018年的地下采空过程和范围(图5)，与InSAR变形图对比可以看到在2017年以前煤矿地表几乎无变形出现；2017年下半年以后，随着开采的进行，采区地表开始出现明显沉降变形，局部变形量很大并造成影像干涉失相干现象(图6)。

由2017年4月26日至2018年8月19日获得的差分干涉图中(图6a—i)，在研究区内分解出多个形变区。

3.2 地下采煤诱发的地表变形的滞后时间

图6中黄色框区域a对应2017年下半年回采区，从图6(c)—(g)影像中得出变形范围主要处在采空区上方及西南侧附近200 m内地区，变形时间在2017年7—12月，与地下采矿时间2017年6—12月存在约1个月的滞后时间。结合图6(c)—(d)，影像5—7月为出现变形，6—8月开始出现少量变形，判断出变形开始时间为2017年7月下旬；从图6(e)中看出，8—11月为变形最大阶段覆盖采区上方及附近区域，干涉图中局部区域出现由于变形量很大造成的失相干现象(遥感影像呈黑色)；在图6(f)、(g)中看出，仅11—12月期间回采a区西侧有局部变形，2018年时间段该区几乎无变形出现。

图5 煤矿回采区分布图

图6 D-InSAR观测结果与地下煤开采布置对照图(a、b、c表示煤矿回采区)

b区域为煤矿2018年1—6月回采区，在布置图中的时间信息看出该回采区自西向东掘进，从图6(g)—(h)干涉图中看出，变形区处于各段回采区中心自西向东移动，与布置图中掘进方向一致，在遥感影像中表现为以采空区为中心的沉降漏斗式变形。根据影像中变形位置与时间对照推测变形时间与实际采空时间存在1个月左右的滞后时间。图6(g)为2018年初至4月的累计变形，而影像中变形区范围仅对应于煤矿2018年1—3月回采区范围；图6(h)为2018年3—5月累计变形值，变形区对应于煤矿2018年4月采区；图6(i)为2018年7—8月中旬累计变形，变形区对应于煤矿2018年6月采区。从变形时间与地下采煤时间对照发现，地表变形时间滞后地下采煤时间约30 d。c区域开采时间在7月份之后，图6(i)中未出现变形，也证明地下采空反应到地表存在一定的滞后时间。

根据b区域工作面较长，地下开采时间与SAR数据干涉对时间分辨率较高的特点，沿b区域走向方向作出时间滞后曲线(图7)。横轴代表时间，纵轴为地面投影距离。观察到红蓝线横向偏差即滞后时间。

图7 b区域横向地下采空与地上沉陷时域对比图

3.3 野外现场验证

围绕采空区及附近山体进行了野外现场调查，发现采空区上方有大面积塌陷坑出现，规模不等；采空区上方山体坡度较陡处，地面有拉张和挤压变形破坏(图8)。

4 地下采空与地表变形关系

4.1 InSAR观测结果三维分析

国内缓倾地层采空区滑坡地表变形研究较少，仅有李宏杰[34]、李滨[35]等人做过研究。塌陷区与沉降区分布范围大，地面实测变形较难开展，且D-InSAR获得的干涉影像为LOS视线向变形值不是实际地表形变位移场，故本文使用了P2升、降轨时间基线较近的两对干涉影像20170514—20171224、20170524—20171206做三维分解处理来获取地表三维形变值，分解区域对应于煤矿2017年下半年采区，三维分解示意图见图9。因为只有升、降两个方向的SAR影像，对于模糊度较高的南北方向变形量不易准确计算，所以本文没有北方向变形量分解。得到2017年5—12月地表的垂直向、东西向的变形量。计算公式如下：

(1)

(2)

ai=cosϑinc,i(i=1,2)

(3)

(4)

式中：du、de——垂直向、东西向变形量；

Γ——转换系数；

dlos,1、dlos,2——升、降轨两干涉对变形值；

ϑinc,i、αaz,i——入射角和方位角。

4.2 顺层采煤地表变形规律

图10为对2017年5月至12月的升、降轨影像三维分解得到该地区垂直向和东西向变形图。在垂直方向上，受地下采空影响采空区上方地表整体为以采空区为中心的地表下沉为主，最大下沉值30 cm位于采空区西侧角点，采区变形量大出现失相干现象；回采区四周100 m范围内均有沉降产生，尤其在西、南两侧沉降较大，最大沉降量10 cm；在东西方向上，采空区地表西侧、西北及西南侧附近岩土体以向东运动为主最大位移值22 cm，东侧岩土体表现出较小位移值且运动方向不一。

推测形成上述地表变形特点的原因有：

(1)该区域为向斜缓倾地层倾向北方，且该处地表坡度方向为东，煤层采空后，地下岩层受重力作用有沿产状及坡向运动趋势而发生塌陷，其上覆岩体的移动便可能涉及到地表，引起地表下沉，造成上方地表沉陷区向西侧发展。

图8 地表变形破坏的野外验证

图9 三维分解示意图[33]

(2)在岩体倾斜成层条件下，自重方向不与岩层层面垂直，在自重作用下，岩体除发生垂直于层面方向的弯曲外，还产生沿层理方向的顺层滑移，外加地表坡度因素造成采空区地表在东西向以东方向运动为主。而东侧地表的运动形成原因可能是受顺层滑移和地下采空塌陷使外侧岩土体向中心水平移动的综合结果。其结果使采空区上山部分岩石受拉，下山部分受压。采空区两侧岩土体在侧向压力下有水平运动趋势，造成采空区上方附近地表形成与运动方向垂直的拉裂缝。

4.3 边界角的计算

根据我国开采沉陷研究的现有规范，本文根据地表移动盆地主断面与10 mm下沉等值线的连线与水平线在煤柱一侧的夹角来计算出上、下山移动角。煤矿2017年下半年回采区内有对应与该开采时间段的垂直向变形值，故根据遥感影像变形等值线计算出了该地区缓倾岩层采空区的上山边界角和下山边界角，如图11～12所示数据，计算得出上山边界角70.23°，下山边界角58.69°，为类似地区判断地下采空范围提供依据。

图10 基于升、降轨观测的三维形变量分解

5 讨论

图11 2017年下半年回采区引发的地表垂直沉降等值线图(白色范围为失相干区)

图12 缓倾斜煤层采矿地面沉降边界角计算示意图(图11中A—A′为计算等值剖面线)

本文虽然通过D-InSAR方法计算出地下采空与地表变形的滞后时间与东西向、垂直向变形值，但计算与分析过程存在误差和一些不确定性因素，包括存在失相干区的最大形变量的确定，滞后时间的确定，InSAR解算的精度问题。在4.2节中，在有失相干的情况下最大沉降量、位移值等为已知数据得出的最大值不包括失相干区。在煤矿采空区地表大变形区域，很容易在差分干涉解缠过程中造成解缠错误，形成失相干现象。在InSAR技术发展中Offset Tracking技术利用强度信息来提取变形值可以解决大变形带来的失相干问题[36]，但偏移量追踪窗口与步长选取以及后处理过程中的滤波方法会引入较大误差，计算精度很难掌握，且Offset Tracking结果不易与D-InSAR结果相融合进行三维分解计算，故本文未进行相应计算。本文根据煤层地下开采资料与地表变形时间确定的变形滞后时间为该区域地质条件下的时间滞后，采深、采厚、岩层、岩性、开采工艺等因素未考虑在内。分析结果存在一定的误差，误差主要来源有：地下开采资料与实际施工存在时间误差；地表变形距离的确定为10 mm沉降边界，滤波处理中会造成沉降量计算误差；地表变形时间段根据升、降轨各时间段是否变形的互补性确定，较低的时间分辨率会在分辨率方面影响实际变形区段的确定；本文使用的三维分解方法为升、降轨D-InSAR反演东西向与垂直向变形，缓倾地层煤矿塌陷地表位移复杂，无法通过假设性原则估算出南北向运动，但D-InSAR计算出的LOS变形与Offset Tracking计算出的距离向和方位向变形可以很好地反演出南北向运动，对煤矿采空区地表运动有较大应用价值，也是本文作者以后的研究方向。