基于SBAS InSAR的新疆哈密砂墩子煤田开采沉陷监测与反演

沙永莲，王晓文，刘国祥，张 瑞，张 波

(1.西南交通大学地球科学与环境工程学院,成都 611756；2.北京理工大学重庆创新中心,重庆 401120;3.高速铁路运营安全空间信息技术国家地方联合工程实验室,成都 611756)

0 引言

煤炭是我国能源结构的主体，对煤炭资源开采活动进行合理的规划和管理对我国经济建设和生态环境保障具有重要意义[1]。我国煤炭资源多是通过地下开采，而这种开采模式往往会破坏地层内部原有的力学平衡状态，导致采空区上方地表发生移动和变形[2]。煤田矿区地表不均匀变形可能导致地面形成裂缝或者塌陷，破坏原有地质环境条件并引发地质灾害，从而严重影响地表相关建筑物等基础设施的安全使用，给矿区当地经济可持续发展造成隐患[3]。

针对矿区不均匀地表形变的时空演化特征进行监测，可为当地地表沉陷隐患治理，安全生产保障和当地生态环境保护等提供关键支撑信息。目前，水准测量、全球导航卫星系统(GNSS)和合成孔径雷达干涉(InSAR)等多种大地测量方法已在矿区地表沉陷监测中得到广泛应用[4-5]。相比于传统大地测量手段，InSAR具有全天候观测、覆盖范围广和观测效率高等优势[6-7]。常规基于两景SAR影像进行干涉处理测量地表形变的两轨法DInSAR (Differential InSAR)，由于无法对形变观测模型中的误差项(如地形误差、轨道误差和大气延迟误差等)进行校正，常被用于瞬时、大幅度地表形变(如地震和火山)测量[8-9]。近十年来，众多学者基于多景SAR影像提出多种时序InSAR分析方法以克服常规DInSAR相位失相干、大气延迟等影响，并实现长时间序列地表形变监测[10-11]。其中，小基线集(small baseline subset，SBAS)时序InSAR方法能够有效利用冗余干涉对进行干涉相位建模，可有效抑制形变恢复模型中的误差项，已在城市地表沉降、活动断层滑移和滑坡运动监测等领域中得到广泛应用。

SBAS时序InSAR技术已在煤田矿区地表形变特征监测、采空区位置和边界探测等方面展现了独特的应用潜力[12]。赵伟颖等[13]应用SBAS InSAR对9景ALOS PALSAR数据进行处理，获取了徐州市某煤矿采动区中村庄在2007—2011年的累计沉降量，表明了利用SBAS InSAR监测矿区地表移动边界及建筑物区域地表移动临界区域的可靠性；李达[14]基于13景TerraSAR-X 数据分析了陕西省榆林市某煤矿区地表形变的时空演化特征，发现不同工作面开采时间与地表沉陷速率有一定相关性；阎跃观等[15]基于Sentinel-1A 影像采用SBAS InSAR 测量了河北某煤矿的地表形变，并根据形变结果计算了该煤矿综采面走向边界角。有学者将煤田地下综采面简化为一个矩形位错模型，建立了综采面参数(开采深度，采煤厚度)与位错模型张裂分量之间的关系，王亚男[16]根据该模型反演了SBAS InSAR测量的内蒙古上湾煤矿地表形变，探讨了采煤深度和厚度与地表形变量特征间的联系。朱煜峰[17]反演了SBAS InSAR观测的江西省丰城市某煤矿地表累积沉陷量，得到了InSAR观测期间位错模型张裂分量；然后利用反演的矩形位错张裂量模拟了矿区地表沉陷，发现与水准测量数据一致，证明了简化位错模型在煤矿地下开采参数反演中的可行性。综上可见，由InSAR观测获取的高分辨率形变场能够反演采煤作业面深度等关键参数，但是基于InSAR地表形变为约束反演煤矿开采量的研究目前还比较少见。在反演得到综采面参数基础上，能否结合煤层视密度构建综采面参数与开采量之间的关系仍值得进一步研究。

本文以新疆哈密砂墩子煤田为例研究了基于时序SBAS InSAR观测形变反演地下煤矿综采面参数及煤矿开采量的可行性。砂墩子煤田是国家重点煤矿，设计产能5.00 Mt/a，自2012年开始试生产以来，至今仍在开采中；该煤田处于半干旱地区，矿区表面及周边植被覆盖较少，非常有利于采用时序InSAR准确恢复其地表形变特征。本研究选取2018年9月—2019年10月共30景Sentinel-1A升轨影像，首先利用SBAS方法恢复了该矿区一年时间内的地表形变；其次，以SBAS InSAR测量的地表形变速率为约束，基于Okada矩形位错模型反演了该煤矿综采面参数，进而根据煤层视密度、综采面参数和开采量的关系估算了砂墩子矿开采量，并与已有矿区资料进行了对比验证。

1 研究区概况及数据源

砂墩子煤矿位于新疆哈密市三道岭矿区西部，地理坐标为E92°26′09″～92°32′10″，N43°08′33″～43°12′26″，自西向东构成“L”形，北距312国道5 km，南距兰新线6 km (图1)。图上，白色虚线为矿区边界线，红框表示InSAR分析区域，光学影像为研究区Sentinel-2光学影像，白色三角为主井位置。井田东西走向长约7 km，南北最宽约8 km，矿井分布三个井筒，主井位置为E92°28′05″，N43°10′04″。从煤炭资源赋存来看，该矿区资源储量877.44 Mt，设计产能5.00 Mt/a，一期设计产能3.00 Mt/a[18]。该矿区水文地质条件简单，属简单型矿床，主采煤层顶底板以粉砂岩-泥岩为主[19]。整个矿区气候干旱，昼夜温差大，属于典型温带大陆性干旱气候；降雨量少但集中，雨期多在夏季，易形成短时洪流。

图1 砂墩子煤矿地理位置及其周边地形图Fig.1 Map of Shadunzi coal site

本文选取30景Sentinel-1A升轨SAR影像对研究区域进行SBAS InSAR形变监测分析，影像获取时间为2018年9月7日—2019年10月8日。经过裁剪后监测区域面积约为66.30 km2，SAR影像空间分辨率约为13.95 m(方位向)×2.33 m(距离向)，SAR影像中心入射角为39.34°。

2 算法原理

2.1 SBAS InSAR 时序地表形变提取

SBAS InSAR是一种经典的时序InSAR方法，它通过选择具有短时空基线的SAR影像对进行差分干涉处理，并针对相干目标点的相位进行建模和解算[20]，在一定程度上能够克服DInSAR常常面临的干涉相位失相干影响。此外，由于SBAS方法应用了尽可能多SAR干涉对参与解算，其能够改正DInSAR测量中的地形残差、轨道误差和大气延迟误差等多项误差，进而得到高分辨率和空间上连续的地表形变场。根据地表形变模型、干涉图组合方法和选点方法等不同，众多学者提出了不同的SBAS 时序InSAR处理策略。本研究采用了Hopper等[21-22]提出的StaMPS SBAS方法进行砂墩子煤矿地表形变测量。StaMPS结合振幅离差指数法和相位稳定性分析法选取滤波后相位微失相干 (slowly decorrelating filtered phase,SDFP)像素，能够在容易失相干地区保障高密度的高相干目标点，以获取连续可靠的形变测量结果[23]。同时，由于StaMPS SBAS方法没有预先设定地表形变的时序变化特征(如假设地表形变随时间发生一阶线性变化)，因而更适宜于矿区等可能在时间序列上具有非线性高阶变化特征的地表形变场提取。同时，StaMPS采用三维相位解缠方法恢复缠绕干涉图的绝对相位，相比于经典二维解缠方法能够有效提升相位解缠的精度和效率[24]。

在SAR数据处理过程中，首先根据SAR影像垂直基线和时间基线进行干涉对选取，其中最大垂直基线和时间基线分别是150 m和109 m，共组成100个干涉对。图2展示了所有干涉图基线组合情况。干涉图生成过程中，设置相干系数阈值为0.4，对低相干区域像素进行掩模。设置幅度差分离差阈值0.4，利用StaMPS相位稳定性分析方法在整个研究区域获取了SDFP点6 663个，平均监测密度为101点/km2，表明基于高分辨率Sentinel-1 SAR影像利用StaMPS SBAS方法能够获取高密度的监测点。

图2 SBAS InSAR时空基线分布图 Fig.2 The spatial-temporal baseline for SBAS InSAR analysis

2.2 基于矩形位错模型的开采参数反演

(a)俯视图 (b)正视图

本研究采用贝叶斯后验反演算法进行砂墩子煤矿综采面参数反演。贝叶斯反演算法是基于统计理论的随机反演算法，关键步骤是在反演过程中进行大量采样，以获得模型后验概率分布，这里常用到马尔科夫链条蒙特卡洛采样方法(Markov chain Monte Carlo method，MCMC)，该方法可进行有效采样并刻画后验概率分布[29]。在反演中主要以贝叶斯框架为理论基础，结合MCMC方法采样，通过获取模型参数后验概率密度函数，进而得到模型参数最优解。在贝叶斯算法中，后验概率密度函数p(m|d)是模型参数m对形变观测量数据d的条件概率函数，公式为:

(1)

式中，p(m)为模型参数的先验概率密度函数。其中，先验信息可来自资料、经验等，常用假设先验概率密度分布有均匀分布、高斯分布等。p(d)是与m无关的归一化常量，p(d|m) 是给定形变数据d情况下模型参数m的似然函数，它能够反映观测数据与反演模型之间的匹配程度。似然函数可定量描述模型参数的不确定性，一般假设多为高斯分布，可表示为：

(2)

式中：N为点数；∑d为方差-协方差阵；G(m)为模型参数m的非线性模型函数。获取先验信息和似然函数后，利用贝叶斯公式获取后验概率密度函数p(m|d)，假设模型参数为n个，第i个参数的边缘概率密度为：

(3)

本研究反演综采面关键参数的具体流程如图4所示,其中输入的InSAR观测形变数据d可为地表形变速率或累积形变量。在反演计算时为了减少不必要的数据冗余，一般首先需要对输入数据进行降采样处理，主要方法有均匀重采样和四叉树重采样，同时构建观测区InSAR形变方差-协方差矩阵∑d，用于计算似然函数；其次，根据先验信息获取先验概率密度函数p(m)，并利用确定的初始模型参数mi(其中mi∈p(m))正演获取模拟形变 (这里正演是指利用模型参数m和非线性模型函数G(m)计算得到地表形变结果)，然后基于以上结果利用式(2)计算似然函数；最后，根据 MCMC采样方法和Metropolis-Hastings法则，自动选择步长来获取新的模型参数，从而计算新似然函数，并根据 MCMC和Metropolis-Hastings法则对两个似然函数进行比较且保留对应的模型参数，以此来对模型参数进行大量有效的抽样，最终获得模型参数后验概率密度函数，并取最大后验概率解为参数最优值，即输出模型参数的最优估值。

图4 基于矩形位错模型的贝叶斯参数反演流程图Fig.4 Flow chart of the Bayesian parameter inversion based on rectangular dislocation model

3 实验结果

3.1 砂墩子煤田地表形变场

由于SAR卫星的侧视成像方式，利用SBAS InSAR测量得到的是沿卫星视线方向(line of sight，LOS)地表形变。考虑煤田开采沉陷主要表现在垂直向发生位移，为更直观地体现矿区沉陷特征，根据雷达波入射角将LOS向形变转换至垂直向。图5(a)展示了利用SBAS InSAR测量的砂墩子矿在2018年9月7日—2019年10月8日期间的年平均沉陷速率图，图中灰色虚线为过沉陷漏斗的剖面线，P1，P2，P3，P4为沉陷区选取的特征点，黑色虚线为矿区边界。图5(a)显示煤田主井西北部存在一个显著的沉陷漏斗，面积约为1.5 km2，其中红色最深的区域表示沉陷中心，此处最大年均沉陷速率约为150 mm/a。通过提取穿越沉陷漏斗的两条剖面线A1A2和B1B2上年均沉陷速率(图5(b))，发现接近沉陷中心的两条剖面线上的沉陷速率变化趋势一致，均表现为从边缘至中心急速持续性下沉趋势，两剖面线处最大年均沉陷速率约为140 mm/a。

(a)基于SBAS InSAR的砂墩子矿年均沉陷速率 (b)沿A1A2和B1B2剖面线地表沉陷速度变化

为直观地表现砂墩子矿地表形变的时序变化特征，图6显示了在沉陷中心区域选取的特征点P1，P2，P3，P4在观测时间段内的时序累计沉陷量变化，图中红色线条为形变序列的拟合曲线。可以看出，4个特征点在2018年9月—2019年6月期间均呈现持续性下沉趋势；而在之后时间段里，4个特征点累计沉陷量逐渐趋于稳定。4个特征点时序形变特征较为一致，交叉验证了SBAS InSAR恢复砂墩子煤田地表形变结果的可靠性。根据《砂墩子煤矿安全生产标准化煤矿申报表》显示，砂墩子矿分多个采区并分阶段开采，2019年6月后砂墩子矿累计沉陷量趋于稳定，可能是该沉陷区对应综采面逐渐减小开采量或者停止开采的结果。

图6 SBAS InSAR获取的4个特征点时序累积沉降量Fig.6 Time series accumulative subsidence of the 4 points from SBAS InSAR analysis

此外，4个特征点在2019年6月—2019年8月发生相对抬升。根据水文地质调查，砂墩子矿地质条件富含较强的含水带或导水带，利于地表降雨和融雪水渗入地下。当地气象资料及《潞新煤化工有限公司矿井改扩建工程报告》显示，矿区降雨多集中于夏季(夏季月均降雨量约30 mm)，融雪水通常也在6至7月汇流而下，易形成短时洪流。另外，砂墩子主采煤层顶底板以粉砂岩-泥岩为主，遇水易膨胀。因此，SBAS InSAR观测到的抬升现象可能是当地的地质条件与气候环境因素综合导致。值得注意的是，在矿区当地发生暴雨和山洪时期，井矿可能面临突发性充水而更易发生事故，应采取积极有效预防措施，以免造成不必要的人员和财产损失。

3.2 Okada位错模型反演结果

煤矿区地下开釆活动引发的地表形变幅度和形态取决于多方面因素，如开采深度、工作面尺寸等。本文基于Okada矩形位错模型，以SBAS InSAR测量的2018—2019年地表年均形变速率为约束，反演了沉陷漏斗区对应综采面的关键参数。注意由于输入反演模型的观测量是地表年均形变速率，故最终反演估计得到的开采量也代表年均值。反演所得参数结果会影响后续开采量的计算，因此验证反演所得综采面参数的准确性十分重要。本文通过以下3个方面来说明反演所得综采面参数的准确性：①对比反演得到的综采面参数与搜集的矿区真实资料；② 分析InSAR观测形变与模拟形变之间的残差大小和分布特征；③对比由模型计算得到的开采量和已有资料公布的矿井产能。

采用后验贝叶斯反演方法获取的砂墩子矿综采面最优参数如表1所示。该综采面中心投影至地表的平面坐标相对于坐标原点(E 92.44°，N 43.20°，即图7左上角顶点)为 968.64 m 东，1 739.45 m 南，开采深度约349.89 m；综采面走向与北方向夹角即走向角约177.41°，这与该煤田自北向南开采的实际情况相符；综采面走向长(沿着煤层走向布置的工作面顺槽长度)约1 001.27 m，倾向宽约211.80 m；综采面与水平面的夹角即综采面倾角约6.01°，根据矿区已有资料显示，砂墩子矿主采4号煤层，煤层倾角5°～25°，平均倾角8°，该矿区采用综采采煤工艺，综采面大致沿煤层层面布置，也就是综采面倾角基本反映该区域煤层倾角，即本研究反演得到的倾角基本符合实际情况。此外，本文将反演的综采面参数结果与二维形变图叠加展示(图7)，可以看出反演的综采面中心投影至地表后基本位于地表形变场中心位置，综采面长度与地表形变场的空间尺度吻合，进一步表明了反演所得最优综采面参数的可靠性。

表1 砂墩子矿综采面最优拟合参数Tab.1 Optimal fitting parameters of fully mechanized mining face in Shadunzi Mine

图7 反演参数与InSAR形变结果的叠加示意图Fig.7 InSAR surface deformation map with the superposition of the inverted optimal parameters

图8中(a)—(c)分别展示了SBAS InSAR观测形变、模拟形变，以及由二者做差得到的残差分布。从图8(a)，(b)形变空间分布来看，InSAR观测形变与模拟形变的分布特征基本一致；煤矿开采主要形变区域(图8白色矩形框)对应的残差值基本小于20 mm，较大残差值主要分布于沉陷漏斗的北缘(图8(c))；整个观测区的残差均值为-2.4 mm，说明残差值相对InSAR观测形变较小。残差分布特征和大小表明基于Okada位错模型能够完好地模拟砂墩子煤矿地表形变场，反演得到的煤矿综采面参数较为可靠。

(a)InSAR观测形变 (b)模拟形变 (c)InSAR观测形变与模拟形变差异

开拓煤量即开采量指在矿井可采储量范围内已完成设计规定的主井等开拓掘进工程构成的煤储量，并减去开拓区内水文地质损失量、设计损失量和在可采期内不能回采的临时煤柱及其他开采量。由于损失煤量通常占比较小可以忽略，本文根据煤田开采设计规范[30]将开拓煤量计算公式简化为：

Q=L×W×M×D，

(4)

式中：L为煤层两翼已开拓的走向长度；W为采区平均倾斜长；M为采区煤层平均厚度；D为煤的视密度，指单个煤粒的质量与外观体积(包括煤的空隙)之比。

经反演获得砂墩子矿综采面走向长约1 001.27 m，倾向宽约211.80 m，根据《潞新煤化工有限公司矿井改扩建工程报告》显示，砂墩子矿主采的4号煤层平均厚度约10.00 m，最大厚度32.48 m，煤层视密度为1.50 t/m3，将上述参数代入式(4)，可大致估算砂墩子矿年均开采量约3.18 Mt。根据新疆维吾尔自治区生态环境保护产业协会发布的《砂墩子矿环境影响评价》公告，砂墩子煤矿自2017年起实际生产规模达已到3.00 Mt/a，表明本文所估算的开采量与该矿实际产能基本相符。

4 结论

本文利用2018年9月7日—2019年10月8日期间覆盖砂墩子矿区的30景Sentinel-1A升轨影像，基于时序SBAS InSAR方法对该矿区进行了地表沉陷监测。SBAS InSAR监测结果显示砂墩子矿区主井西北侧存在一个明显的沉陷漏斗，对应该矿目前投产的一个综采面，沉陷面积约1.5 km2，最大年均沉陷速率约150 mm/a。综采面上方地表累积形变时序曲线显示在2018年9月—2019年6月，沉陷漏斗区地表表现持续性下沉趋势；2019年6月后累积沉陷量趋于稳定，可能是开采作业放缓或者停止的结果。

基于Okada矩形位错模型，本文采用贝叶斯后验反演方法利用SBAS InSAR观测形变，反演得到砂墩子矿区综采面最优参数解，其中采深约349.89 m，综采面走向与北方向夹角约177.41°，走向长约1 001.27 m，倾向宽约211.80 m，倾角约6.01°，与已有矿区资料基本符合。将由最优参数解算的模拟形变与SBAS InSAR观测形变结果对比，发现它们空间分布特征一致，两者残差较小。根据反演所得综采面走向长与倾向宽参数，估算得到砂墩子矿在InSAR观测期间的年均开采量约3.18 Mt，与已有矿区资料报道的矿区年产能3.00 Mt/a一致。本文研究结果表明了基于时序InSAR观测形变反演和估算煤田综采面关键参数的可行性。

志谢：感谢欧空局提供Sentinel-1A SAR影像及美国地质调查局(UGSG)提供的SRTM-1地形数据。