蔚县矿区地面沉陷InSAR 多维形变监测

杜建涛，闫 丽，赵超英

(1.长安大学地质工程与测绘学院，陕西 西安 710054；2.中煤科工集团西安研究院有限公司，陕西 西安 710077)

我国是一个煤炭能源大国，已探明的煤储量位居世界前列，而矿产资源的过度开采容易破坏矿区地质结构，进而引发地表开裂、塌陷甚至滑坡等地质灾害[1]。因此，矿区沉陷监测是煤矿安全开采和可持续发展的重要工作，对预防潜在地质灾害具有重要意义[2]。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Rader，InSAR)因其高精度、大覆盖范围及可重复观测等技术优势在采矿沉陷监测中得到广泛应用，极大地克服了传统测量手段空间分辨率低、周期长、成本高等局限。Zhao Chaoying等[3]、Yang Chengsheng 等[4]对山西大同矿区大量地表塌陷进行InSAR 监测；Zheng Meinan 等[5]、Du Sen等[6]基于多种先进InSAR 技术分析了河北峰峰矿区的沉陷机理；Yang Zefa 等[7]、Ma Chao 等[8]、贺卫中等[9]针对陕西榆林矿区地面塌陷采用InSAR 技术进行大范围、多维形变监测。然而，目前针对采矿沉陷的InSAR 时空演化分析和精细化监测研究不足，且该技术在河北蔚县矿区沉陷监测的应用也较少。为了揭示该矿区的地表形变演化特征，本文利用Sentinel-1A/B 数据对采矿沉陷开展InSAR 监测，获取了矿区升降轨数据沿视线向(LOS)形变监测结果，并在此基础上重点分析了西细庄矿近1 a 内地表二维形变的时空演化特征。

1 研究区概况

蔚县位于河北省张家口市最南边，境内南部为高山地区；中部地势平坦，壶流河蜿蜒而过；北部为丘陵区，沟壑密布，地下蕴藏着丰富的煤炭资源。蔚县矿区总面积264 km2，距离蔚县县城约10 km，是河北省最大的地下采煤区之一。矿区井田分布如图1 所示，崔家寨井田、单侯井田、南留庄井田和西细庄井田为目前主要开采区域；郑沟湾矿、兴源矿及周边部分矿井已开采殆尽。受小煤窑破坏，崔家寨井田西翼四采区目前处于停采状态。矿区含煤地层为侏罗系下花园组，分为上、下两段，上段为杂色含煤地层，可靠性较差，均不可采；下段含煤8 层，位于底部的1 号煤层(西细庄矿主采层，平均厚度3.16 m)和中部的5 号、6 号煤层(崔家寨矿主采层，平均厚度分别为0.9 m 和2.04 m)，层位稳定，厚度和分布面积均较大，为矿区主采煤层，埋藏深度为143.6～336.5 m。煤层顶板岩性为泥岩、粉砂岩、细砂岩，底板岩性为粉砂岩、细砂岩及泥岩。地表第四系松散层厚140～245 m，基岩厚度30～55 m，松散层下段砂砾石含水层厚7.5～34.3 m，岩性为粗砂、砾石及砾卵石[10]。

图1 河北蔚县地理位置及SAR 影像覆盖范围Fig.1 Location of Yuxian County and SAR image coverage

2 数据采集与处理

本文采用61景Sentinel-1A/B干涉宽幅(Interferometric Wide swath，IW)模式数据进行矿区形变监测，包括33 景降轨数据(轨道号T120，时间为2017 年7 月至2018 年9 月)和28 景升轨数据(轨道号T142，时间为2017 年8 月至2018 年8 月)，影像具体参数如表1 所示。此外，本文还利用精密轨道数据AUX_POEORB文件来精化影像初始轨道状态矢量，利用30 m 分辨率的SRTM DEM 数据来模拟和去除地形相位。

表1 Sentinel-1A/B SAR 数据参数Table 1 SAR data and their imaging parameters

对获取的61 景升降轨Sentinel-1A/B 数据分别进行小基线集(Small Baseline Subsets，SBAS)干涉处理[11]，然后采用融合升降轨数据的多维形变时序估计方法，生成垂向和东西向的二维形变时间序列[12]。关键处理步骤如下：

a.Sentinel-1 预处理 由于Sentinel-1 数据的特殊成像模式，首先需对原始影像进行预处理，将其转化成可以干涉的单视复数(SLC)格式，具体包括采用AUX_POEORB 轨道文件进行基线精化，对主、从影像进行去斜和配准，使方位向配准精度达到千分之一像元。

b.SLC 影像干涉 设置时间和空间基线阈值生成小基线差分干涉图。为提高干涉图的相干性，并抑制噪声，对差分干涉图进行自适应频率域滤波处理，采用基于Delaunay 三角网的最小费用流算法进行相位解缠。为抑制时空去相干的影像，本文时间基线阈值取36 d，空间基线阈值取200 m。干涉组合的时空基线分布如图2 所示。

图2 干涉组合时空基线分布Fig.2 Baseline of InSAR combination

c.生成形变图 挑选高质量干涉组合，估计并去除趋势项和大气延迟相位，构建相干像元的形变相位解算模型，采用最小二乘法进行参数求解，得到SAR 影像获取时刻间的平均形变速率。

d.计算二维形变时间序列 融合多轨道SAR 数据的多维形变时序估计技术，利用奇异值分解法来计算多源SAR 数据在重叠时间段内的二维形变速率和形变时间序列。为了避免观测方程的病态，需要对设计矩阵进行正则化处理[13]。观测方程表示为式(1)。

式中A=(-cosαsinθΔt,cosθΔt)为设计矩阵；α为卫星飞行方向的方位角；θ为雷达波入射角；Δt为相邻SAR影像时间间隔；ξ为正则化参数；I为单位矩阵；V=(VE,VU)T，VE和VU为待求参数，分别表示相邻SAR 影像获取时刻的东西向和垂向形变速率；D为形变相位矩阵。待求参数估计值可表示为：

最后，对相邻SAR 获取时刻间的形变速率在时间域积分，即可得到地表东西向和垂向的形变时间序列：

式中d为累积形变时间序列；n为不同轨道影像时间跨度重合的总数量。

3 结果分析

3.1 矿区沉陷监测数据

将获取的升降轨SAR 数据按上述流程处理，生成蔚县矿区LOS 向年平均形变速率图，重采样到地理坐标系下，如图3 所示，其中，图3a 为降轨形变监测结果，图3b 为升轨形变监测结果。

图3 河北蔚县矿区2017—2018 年地表形变场Fig.3 Surface deformation site of Yuxian mining area of Hebei Province from 2017 to 2018

从图3 可以看出，升、降轨数据获取的矿区形变场基本一致，整个矿区在监测期间均出现不同程度的地表形变现象，且沉陷区普遍呈近椭圆形的漏斗状分布。矿区崔家寨井田、单侯井田和西细庄井田内采煤造成的地面沉陷较为严重，南留庄井田在监测期间形变相对不明显。其中，崔家寨井田沉陷面积最大，且沉降速率达-30 cm/a。资料显示，矿区在2015 年调整井田边界和优化资源配置，提高了各矿井生产能力，其中，崔家寨煤矿生产能力180 万t/a，多个立井开拓；单侯煤矿生产能力150 万t/a，3 个立井开拓；西细庄矿生产能力60 万t/a，一对立井开拓；南留庄矿生产能力30 万t/a，一对立井开拓[14]。地下煤炭资源的长期高强度开采破坏开采区周围原始应力平衡，从而导致覆岩移动并向上波及地表，使地表发生形变和沉陷。

为了定量分析矿区地面沉陷灾害情况，本文对不同量级年沉陷速率区域的面积进行统计分析。图4 为矿区沉陷等值线图，等值线由内向外依次为-20、-15、-10、-5 cm/a，从图4 可以看出，崔家寨、麦子破、回回墓等多个村庄受到采矿沉陷影响。表2 为不同量级年沉陷速率区域面积统计结果，其中年沉陷速率超过-10 cm/a 的区域达到2.16 km2，年沉陷速率超过-20 cm/a 的区域达到0.10 km2。大面积的地表沉陷对矿区的基础建设及地表生态环境造成了严重的影响[15]。

图4 河北蔚县矿区沉陷等值线(单位：cm/a)Fig.4 Subsidence contours in Yuxian mining area in Hebei Province

表2 河北蔚县矿区沉陷面积统计Table 2 Subsidence area statistics in Yuxian mining area of Hebei Province

对比图3a、图3b 可以发现，升、降轨数据获取的监测结果中部分区域在形变范围和量级上存在一些差异。为此，本文提取跨越西细庄矿塌陷区的剖线AA′，对其形变空间特征进行分析(图5)。从图5 可以看出，升、降轨数据获取的西细庄矿最大年形变速率相差约5 cm/a，且沉陷中心位置也具有一定偏差。究其原因，可能是由于升、降轨数据具有不同的轨道方向和入射角度，其对地表真实形变的敏感度不同。另外，采矿沉陷的LOS 向形变中可能包含水平形变的影响。

3.2 二维形变时序

图5 蔚县西细庄煤矿年形变速率剖线Fig.5 Annual deformation rate along profile A—A′ of Xixizhuang mine in Yuxian County

由于西细庄井田的干涉图整体相干性较高，便于时序分析，因此，在升、降轨一维形变分析基础上，仍以该矿为试验区，采用融合多轨道SAR 数据的多维形变估算法，研究矿区二维形变时空演化特征。将统一到相同地理坐标系下的升降轨形变图，按式(1)构建方程组，求解不同时间段内的二维形变速率，然后依据式(3)，获取西细庄矿2017 年8月至2018 年8 月期间垂向与东西向的形变时间序列，如图6 所示。分析图6 可以发现，西细庄矿地表在监测期间以垂向形变为主，也具有明显的东西向水平形变，且东、西水平形变方向相反，量级也存在一定差异。垂向形变在空间上呈规则的椭圆形，东西向水平形变则表现为两个半圆形，二者形变范围均随时间逐渐扩大，形变量级也相应增加。

为定量分析矿区地表形变的时空演化特征，共选取4 个特征点P1、P2、P3、P4(图6)，进行二维形变时间序列分析(图7)，其中P1、P2 点分别位于沉陷中心区域及采矿工作面倾向边缘，P3、P4 点位于采矿工作面走向边缘。

由图7 可以看出，在2017 年8 月至2018 年8月期间，P1 点的垂向累积形变最大，达–30 cm，且具有持续下沉的趋势；P3、P4 点位于东、西向水平形变最大的区域，形变量分别达到–10 cm 和7 cm，其形变方向相反，其中正值表示向东移动，负值表示向西移动，且形变都指向沉陷漏斗中心，同时二者均具有较大的垂向形变，形变量分别达–22 cm 和–13 cm。可以看出，P3 点的垂向和东西向形变量都较点P4 大，这与工作面的开采方向有关，P2 点的垂向形变达–20 cm，但和P1 点都具有水平方向保持不变的特点。此外，不同于城市地面沉降的线性特征，采矿沉降的形变过程与时间呈现出非线性关系，从图中P1 点的垂向形变时间序列曲线(图7a)可以看出，该点处的地表形变具有先缓慢下沉后加速沉降最后趋于稳定的特征，这一形变演化过程对应于矿区地表单点沉降的初始沉降、主要沉降和残余沉降3 个时期，与“S”型增长的Logistic 时间函数模型能够很好地吻合[16-17]，更能反映采矿沉陷的真实规律，因此，利用InSAR 数据进行多维形变估算获得的地表沉陷值，与实际情况较为吻合，可为沉陷的分析、治理提供参考数据。

图6 蔚县西细庄矿二维形变时间序列Fig.6 Time sequence of 2D deformation of Xixizhuang mine in Yuxian County

图7 蔚县西细庄矿特征点二维形变时间序列Fig.7 Time sequence of 2D deformation of feature points of Xixizhuang mine in Yuxian County

4 结论

a.采用InSAR 技术获取了河北蔚县采矿区大范围地面沉陷的空间分布特征，且部分区域年形变速率达到-20 cm/a。该结果显示出InSAR 技术用于矿区地表形变监测的优越性，为矿区安全开采和生态恢复提供科学依据。

b.通过对研究区西细庄井田地表形变的精细化监测与分析，获得采矿沉陷的时空演化特征，即矿区地表动态沉降过程符合Logistic 时间函数模型，该特征揭示了真实的采矿沉陷规律。

c.基于时间函数模型参数的采矿沉陷预测和矿区形变三维分解将是下一步研究计划。