摘要
为了分析内蒙古白音华一号矿区地面沉降的时空变化特征及地表形变规律，本文利用Sentinel-1B卫星2 条相邻降轨数据，结合相位叠加技术（Stacking）和小基线集（Small Baseline Subset，SBAS）技术，提取了矿区的沉降形变信息。通过对2条轨道InSAR监测结果的对比分析发现，矿区沉降主要集中在南、北排土场区域。其中，南排土场的累计形变量达527 mm，年均形变速率为148 mm/a；北排土场的最大累计形变量为 400 mm，年均形变速率为 112 mm/a。为进一步分析矿区的时序形变特征，本文选取了南北排土场的 8个特征形变点进行研究，结果表明两区域的形变规律较为一致，主要经历了“缓慢形变—加速形变—缓慢形变”3个阶段。其中，加速形变阶段集中在每年的1—6月，其变化规律与矿区工期规划、冻融影响有关。研究表明，矿产开采与排弃计划、冻融等气候条件是导致矿区地表沉降的主要原因；研究成果可为矿区地质灾害防治、工期优化和环境保护提供重要的参考依据。
关键词
Abstract
To analyze the spatiotemporal characteristics and surface deformation patterns of ground subsidence in the Baiyinhua No.1 mining area in Inner Mongolia, this study utilized Sentinel-1B satellite descending orbit data from two adjacent tracks. Phase stacking (Stacking) and Small Baseline Subset (SBAS) techniques were employed to extract subsidence and deformation information of the mining area. Comparative analysis of InSAR results from the two tracks revealed that subsidence is primarily concentrated in the south and north dumping sites. The cumulative deformation in the south dumping site reached 527 mm, with an average annual deformation rate of 148 mm/year, while the maximum cumulative deformation in the north dumping site reached 400 mm, with an average annual deformation rate of 112 mm/year. To further explore the temporal deformation characteristics of the mining area, eight characteristic deformation points in the two dumping sites were selected for analysis. The results indicate that both sites exhibit similar deformation patterns, characterized by three stages: "slow deformation—accelerated deformation—slow deformation". Among them, the accelerated deformation stage is concentrated from January to June each year, and its variation law is related to the mining area construction period planning and freeze-thaw influence. Studies have shown that mineral mining and dumping plans, freeze-thaw climate conditions, etc., are the main causes of surface subsidence in mining areas. The research results provide important reference basis for geological disaster prevention, construction period optimization, and environmental protection in mining areas.
0 引言
随着中国对能源的需求量不断增大，矿产资源的开采强度逐渐增加，这对采矿安全和矿区周围建筑物造成了巨大威胁（何满潮等，2005）。矿区地表沉降监测是采煤沉陷区可持续发展的重要保障（张庆斌等，2021）。近些年，随着白音华一号矿区开采活动的进行，该矿区采空区的范围也逐渐扩大，引起了较为明显的地表形变，加大了地质灾害发生的风险（樊昱初等，2024）。随着该矿区开采活动的不断进行，开采边坡的高度也在逐年增加，特别是软岩边坡的滑坡现象数量众多、规模巨大，且防治难度高，给矿业生产和环境带来了巨大的危害（孟繁和耿谏，2019）。传统对地观测方法有水准测量、三角高程测量、全站仪测量、全球导航卫星系统 （global navigation satellite system，GNSS）等（杨红磊等，2021）。使用 GPS 这类散点式观测手段虽然有较高的精确度（陶延林和张明，2019），但人力物力耗费严重且监测范围有限，难以达到大区域范围内高效监测的目的（李启亮等，2024）。
近些年发展起来一种全新的对地观测技术：合成孔径雷达干涉测量（Interferometric Synthetic Aper‐ ture Radar，InSAR），该技术具有全天时、全天候的特点，可以从空间直接获取大范围、高精度的地表高程和形变信息。InSAR 技术已经成为目前发展迅速、极具潜力的新型对地观测及测绘技术，成为各国地学界研究的热点之一（王刘宇，2022）。汤志刚等（2024）利用InSAR技术分析地下煤炭资源开采引起的地面沉陷监测，充分发挥了技术优势，快速获取雷达影像覆盖范围的地表形变信息，已成为大范围矿区开采沉陷监测必不可少的监测方法。以往研究获取的露天矿边坡形变信息主要是以地基检测设备为主（王天宇等，2019），但使用地基监测设备成本代价都较高而且不能进行长时间、大尺度的进行监测。利用时间序列InSAR技术可以对矿区开采沉陷、滑坡等地表形变进行有效的监测（朱建军等，2017）。时间序列InSAR技术在整体上对一定区域的地质灾害分布以及其严重程度有着更加深刻的认识，并且提高了地质灾害监测者对监测工作的管理水平，为地质灾害的防治工作提供基础性的资料（刘雨鑫，2019）。白音华煤矿为全掩盖煤矿，位于大兴安岭西坡南端北侧，二连盆地东端乌尼特坳陷中，构造总体方向为北东向，呈长条状展布，地层倾角较为平缓（张建强等，2016）。在白音华矿区范围内，学者们对矿区开采和地表稳定性方面开展了较多研究。尚涛等（2007）分析了白音华矿区的开采规模和开采工作推进的模式；王浩等（2009）利用 RTK技术对白音华露天矿边坡位移进行了监测，发现边坡变形与调查结果基本吻合；刘永杰和王兆刚 （2014）介绍了白音华露天煤矿西端帮滑坡情况，阐述了滑坡治理措施及滑坡治理后的稳定性分析；白润才等（2017）对白音华露天矿北帮边坡稳定性进行了分析；王立文（2022）利用 GNSS 技术对白音华矿区的边坡稳定性进行了监测。本文利用欧洲航天局发射的 Sentinel-1B SAR 影像对白音华一号露天矿进行监测，利用 Stacking 技术得到近似形变速率进行参考，最终使用 SBAS-InSAR 技术获取该矿区地表的形变特征信息，本文采用了 2 个相邻轨道的哨兵-1B 数据，在时间段为 2018 年 1 月—2021 年 12 月分别获取了 90 景和 87 景数据进行长时序分析，利用 2 个轨道数据获取的地表形变信息进行对比分析，在前人研究基础上重点研究2018—2021年这 4 年期间矿区地表沉降的形变特征，以期为矿区防灾减灾提供技术支撑。
图1矿区附近高程分布图（a）及研究区概况图（b）
表1Sentinel-1B影像参数
1 研究区及数据概况
1.1 研究区概况
白音华一号矿区位于内蒙古自治区锡林郭勒盟西乌珠穆沁旗巴彦花镇，是内蒙古自治区中部特大型煤田之一，面积形状呈东西长，南北较窄型。地理坐标位于北纬 44°46'～44°50'，东经 118°23'～118°29'。白音华一号露天煤矿位于白音华煤田西南部，目前正在开采首采区。如图1a所示，矿区内地形呈南高北低、东高西低状，海拔高度为 1000～1050 m，主要地貌特征为高原低山缓坡丘陵地形。按照地层赋存与边坡产状之间的空间关系，可将矿区边坡划分为 6 个工程地质区域，即采场东帮、西帮、北帮，内排土场，南排土场及北排土场（图1b）。边坡稳定性主要受排弃物堆积形态和基底的影响 （贾兰等，2023）。
1.2 数据概况
本文选择欧洲空间局（ESA）发射的Sentinel-1B 卫星，获取了覆盖研究区域分别为 90 景和 87 景单视复数影像（SLC）作为实验数据，数据范围如图2所示。表1显示的是所用 SAR 影像的基本参数情况，其中获取了 2018 年 1 月—2021 年 12 月期间的 2 个轨道的影像数据。数据极化方式为 VV，轨道方向为降轨，分辨率为 5 m×20 m。此外，本文使用的数字高程模型（DEM）数据为美国宇航局（NASA）公开的SRTM数据，其分辨率为30 m×30 m。
2 时序InSAR技术与方法
本文利用 Sentinel-1B 数据、高精度 DEM 数据，采用 Stacking 技术和 SBAS-InSAR 技术获得高精度差分干涉形变图、时序形变速率图、累计形变图，然后结合强度影像图、阴影叠掩分布图、相干系数图、光学影像图等，按一定阈值消除误差图斑，探测识别地表形变聚集区，作为识别地质灾害隐患的前提。得到研究区速率图和累计形变图，同时针对重点形变区域生成其累计形变曲线，开展进一步分析。
图2白音华一号矿区周边影像及位置
2.1 Stacking技术
Stacking（Sandwell and Price，1998）技术也称相位叠加技术。该技术假设：大气误差随机存在于每个独立的干涉图中，而形变在时间域上表现出线性变化特征。该技术在一定程度上可以削弱InSAR技术中存在的随机的轨道误差、大气延迟相位误差及地形误差等影响。Stacking技术的数学模型可以表述为：
(1)
式（1）中，v 为线性形变速率（cm/a）；λ 为波长 （cm）；ϕi为第 i个解缠后的差分干涉相位图；ti为第 i 个差分干涉相位的时间间隔。相应的误差传播公式为：
(2)
式（2）中，E 为假定的单幅干涉图相位误差；M 为干涉图数量。
在 Stacking 数据处理过程中，要使用去除大气误差的解缠干涉图，并对解缠干涉图进行人为评判剔除掉一些误差较大的干涉图，从而获取较为准确的形变速率（熊国华，2022）。
2.2 SBAS-InSAR技术
SBAS（Berardino et al.，2002）技术，又称小基线集技术，是目前有代表性的时序 InSAR 方法之一。该方法通过设定合理时间和空间阈值，选取时空基线较短且相干性满足解算要求的干涉对，可以有效减小时空失相关对形变结果的影响。采用 SBAS-InSAR时序形变解算方法可以有效抑制地形误差和大气误差（杨超等，2023）并提高形变监测精度。 SBAS 技术假设获得了同一地区不同时间的 N+1 景数据，根据所选的时间基线和空间基线为约束条件，并在这些影像中筛选一景，把这景作为主影像并将其余影像分别主影像进行配准，并对所选择的干涉对进行差分干涉处理，得到差分M幅干涉图，M 满足：
(3)
若所获得的 SAR 影像干涉对的时间分别为 Ta、 Tb，且 Tb>Ta，生成的差分干涉图为第 j 景，则差分图的干涉相位可以用如下公式表示：
(4)
式（4）中，1≤j≤M；、分别是 Tb、Ta时间相对于 T0时间 LOS 方向的累积形变量（cm）。ϕj、、分别表示 j 时刻、Tb 时刻和 Ta 时刻的相位（rad）。 Δϕjtop、Δϕjarm、Δϕjnoise分别为地形相位、大气相位、噪声相位（rad）。去除 Δϕjtop、Δϕjarm、Δϕjnoise 后，将式（4）改写为：
(5)
式（5）中，ϕj、、分别表示 j 时刻、Tb时刻和 Ta时刻的相位（rad）。、分别是 Tb、Ta时间相对于T0时间LOS方向的累积形变量（cm）。
假定相邻间隔时间内，将地表形变速率认为是线性的，第j景干涉图的形变相位信息为：
(6)
式（6）中，Δϕj为相位差（rad）。vk，v 为线性形变速率（cm/a）。A 是 M×N 的矩阵，由于在干涉对生成时对时空基线阈值进行了限制，因此可能造成不连续干涉图子集，此时矩阵 A 为秩亏阵。为了解决秩亏问题，采用奇异值分解（SVD）法对多个小基线集联合求解，再用最小二乘法求得最小范数下各高质量相干点的形变量，进而得到时间序列形变量（梁涛，2014）。
2.3 时序InSAR技术处理流程
本文获取了 2 个轨道的 Sentinel-1B 数据，选择了 2019.08.13和 2019.08.08日期分别为轨道 47-441 和 149-443 的配准主影像，其他影像依次与主影像配准。根据数据的临界基线和时间失相干条件，本文设置垂直最大垂直基线 230 m，最大时间基线在 150 d，共得到286对干涉组合（图3）。
图3空间基线图
a—轨道47-441；b—轨道149-443
（1）通过导入外部DEM数据来模拟干涉图中的地形相位，而后将其从差分干涉图中减去得到不含地形相位的干涉图。并选用 Goldstein 方法选择窗口大小对得到的差分干涉图进行滤波处理，最后采用基于 Delaunay 三角网的最小费用流法（minimum cost flow，MCF）对滤波后的差分干涉图进再次进行差分、滤波和解缠。如图4所示，进行精化基线可以有效减小基线轨道误差，进行 Goldstein 滤波可以有效减少相干斑噪声误差可以避免错误点影响计算结果（李陶等，2007）。
图4使用Goldstein滤波前后的差分干涉图和解缠干涉图
a1—轨道号47-441的20180207—20180303初始差分干涉图；b1—轨道号47-441的20180207—20180303初始解缠干涉图；c1—轨道号47-441 的20180207—20180303改正后差分干涉图；d1—轨道号47-441的20180207—20180303改正后解缠干涉图；a2—轨道号149-443的 20180214—20180310初始差分干涉图；b2—轨道号149-443的20180214—20180310初始解缠干涉图；c2—轨道号149-443的20180214— 20180310改正后差分干涉图；d2—轨道号149-443的20180214—20180310改正后解缠干涉图
（2）利用 Stacking结果生成的平均速率形变图，识别出形变区与稳定区并进掩膜操作掩盖掉形变范围，对非形变区进行大窗口空间滤波消弱干涉图中的大气影响和解缠过程中的一些错误点。而后对掩膜掉的形变区进行插值处理，最终将滤波结果从干涉图中减去。如图5可知，空间滤波后解缠干涉图相位标准差有着明显的降低，从而达到去噪处理。使用空间滤波技术对解缠相位进行分大小窗口滤波使干涉图的平均相位标准差降低了 33.9% （轨道号 149-443）；从 1.136 降低到 0.629，提高了 44.64%（轨道号 47-441），可以看到 2个轨道的数据在进空间滤波后的相位标准差都明显降低。
（3）利用奇异值分解（SVD）法估算平均形变速率和残余地形相位，并对干涉图去平，重新进行相位解缠和精炼，生成更为精确的结果。通过进一步反演，进行定制的大气滤波，估算去除大气相位，得到更精确的时间序列位移结果。
（4）将生成的结果进行地理编码，得到研究区雷达视线向（line of sight，LOS）的平均形变速率和累积形变量。由于本研究所用的2个轨道的航向角非常接近，因此可以直接对比 2 个轨道的结果不用将结果转换到垂直向。数据处理技术流程图如图6所示。
图5空间滤波前后相位标准差折线图
a—轨道47-441；b—轨道149-443
图6数据处理技术流程图
3 矿区时序监测结果分析
3.1 矿区地表形变结果分析
在研究时段 2018 年 1 月—2021 年 12 月近 4 年的时间里，本文通过 SBAS-InSAR 技术处理得到研究区的累计形变量（图7）和年平均形变速率（图8），其中负值表示远离卫星方向（下沉），正值表示靠近卫星方向（抬升）。从累积形变图可以看出，研究区整体形变趋势以下沉为主，采矿工程影响范围内的 2个轨道的累积最大沉降量可达524 mm和527 mm，最大平均年沉降速度达 148 mm/a 和 119 mm/a。综合2个轨道的数据可以得到南排土场、内排土场、北排土场以及采场东西北帮的形变速率。其中南排土场的形变速率为 53~95 mm/a；内排土场形变速率为 102~148 mm/a；采场东西北帮的形变速率较小，为 10~26 mm/a。由此可知矿区的主要形变位置在南北排土场和内排土场，其余地方沉降不明显。结合相关资料分析可知，研究区的主要变形区域与区内开采和排弃工作面位置空间分布一致。该区域内分布有采煤工作区和排土工作区，地面沉降主要受开采活动和排土活动的影响，快速沉降区主要集中在北排土场、东端帮南部和南排土场周围，沉降量从沉降中心向周围递减，矿区地面普遍受到沉降灾害的影响；而研究区东端帮东北部和西端帮北部地区在监测期内有较为明显的抬升运动，最大累积抬升量可达 301 mm 和 211 mm，最大平均年沉降速围递减。由此可见，白音华一号矿区现有的采区和排土区域都分布有快速沉降区域，并随着开采和排弃强度增大，沉降速度有可能继续加速。此外，为验证SBAS-InSAR技术的可靠性，本文采用Stacking 技术进行了对比验证。如图9所示，Stacking技术提取的形变规律与 SBAS-InSAR 结果基本一致，表明 SBAS-InSAR 技术在形变监测中的可行性。然而，由于Stacking技术在叠加处理过程中可能导致形变信号的细节被平滑化，其对短期快速变化或小尺度特征的捕捉能力有限。因此，本文主要以 SBAS-InSAR 技术为核心，后续的特征点时序分析也基于 SBAS-InSAR技术展开。
3.2 特征点形变结果分析
为了进一步研究地面沉降速度及其对矿区周围的影响，分别在南、北排土场的大形变区域设置了 8 个监测站，并在每个点的周围选择 3 个相邻的点来共同表示该特征点的时序形变曲线图。其中 A、B、C、D为一组，位于白音华一号矿区南排土场的区域，E、F、G、H 为一组，位于白音华一号矿区北排土场的区域，具体特征点位分布如图10所示。
图7矿区地表累计形变量
a—轨道47-441；b—轨道149-443
图8矿区平均形变速率
a—轨道47-441；b—轨道149-443
图9利用Stacking技术获取的形变速率图
a—轨道47-441；b—轨道149-443
图10特征点位分布图
如图11所示为特征点 A、B、C、D 的时序形变图，从图中可以看出在南排土场的北边的 A 点的沉降变化范围最大，可以看到 2 个轨道的最大累计沉降值分别可达到 530 mm 和 490 mm；从时间上来看 A 点在 2018 年 1 月—2018 年 7 月形变速率超过 100 mm/a；从 2018 年 7 月—2018 年 9 月 A 点有着短暂的稳定期；从2018年9月—2019年12月期间A点形变持续增加，平均形变速率可达 130 mm/a；2019 年 12 月—2020年8月，A点形变又处于平缓阶段，形变速度非常缓慢；从2020年8月—2021年12月A点又处于加速形变的过程，形变速率可达 134 mm/a。A 点的形变规律可以总结为：形变—短暂稳定—加速形变—稳定—加速形变。B、C、D 点的最大累计沉降值分别可达 357 mm、350 mm、398 mm，且形变规律与A点较为一致。从站点位置上来看A点和D点属于内排土场的工作重点区，排弃工作强度较大，而B 点和 C 点位于南排土场边缘地带，受开采工作影响相对较小，因此形变量级也大约为A、D点形变量级的2/3大小。比较2个轨道的同名点的形变曲线，二者的时序形变曲线极其相似，滑动趋势基本一致，形变速率也基本吻合。
如图12所示为特征点 E、F、G、H 的时序形变图。从图中可以看出在北排土场的E点和H点的累积沉降量分别为350 mm和400 mm；F点和G点的累积沉降量分别为340 mm和350 mm。在图上可以观察到这 4 点的形变规律较为一致；从 2018 年 1 月— 2018 年 8 月形变速率为 110 mm/a，处于加速沉降阶段；从 2018年 8月—2021年 8月形变速率为 95 mm/ a，地表处于匀速沉降阶段；从 2021年 8月—2021年 12月形变速率为10 mm/a，此时地表逐渐趋于稳定。从点位分布上来看，点 E 和点 G 位形变规律与量级更为相似，这是由于 E、G 位于北排土场的外侧，受开采与排弃活动影响较弱；而点 F 和点 H 则位于北排土场的靠内一侧，受到开采与排弃活动比较剧烈，则形变量级比 E、F 两点要稍大一些。对比 2 个轨道的同名点时序形变曲线可以看到形变规律基本一致。
结合矿区采矿计划发现在该时间段内处于持续开采和排土工作中，这与本文分析的时序形变规律基本一致。总结来说，白音华一号矿区南排土场区域的沉降情况远比北排土场区域的沉降情况严重，南侧区域的排弃活动强度更大、车流更多且更为频繁，应该重点监测该区域地表形变情况，制定合理的开采与排弃计划，避免出现煤矿坍塌的情况 （钱鸣高等，2018）。
图11特征点A、B、C、D时序形变图
图12特征点E、F、G、H时序形变图
4 结论
本文采用 SBAS-InSAR 方法，对 2018 年 1 月— 2021 年 12 月间白音华一号矿区的地表形变情况进行了监测，基于 Sentinel-1B卫星影像的分析得出了以下主要结论：
（1）研究发现，南、北排土场为地表沉降的主要区域，沉降规模较大，沉降速率最高可达 148 mm/a，累计沉降量达 525 mm。南排土场的形变量和形变速率均较大，应重点关注。
（2）在南北排土场的高形变区域选取了 8 个特征点，绘制并分析了时序沉降曲线。其中，特征点A 的形变规律表现为“缓慢形变—加速形变—缓慢形变”，加速阶段集中在1—6月，其沉降趋势与矿区开采、排弃计划、自然沉降和冻融等气候因素相关。其他特征点的形变规律与 A 点相似，说明矿区形变主要与矿区开采、排弃计划、自然沉降和冻融等气候因素密切相关。
（3）通过时序 InSAR 技术对白音华一号矿区地表形变研究，能够发现矿区变形较大的区域，进而对变形较大区域采取有针对性的监测和治理措施，确保矿山采场和排土场的边坡安全。
