基于SBAS 技术的煤矿采空区探查与稳定性评价

2022-04-29刁鑫鹏谭秀全李爱军

煤炭科学技术 2022年3期

张丰，刁鑫鹏，谭秀全，李爱军，杨隽

（1.山东省鲁南地质工程勘察院（山东省地质矿产勘查开发局第二地质大队），山东济宁 272100；2.自然资源部采煤沉陷区综合治理工程技术创新中心，山东济宁 272100；3.中国矿业大学环境与测绘学院，江苏徐州 221116；4.江苏省资源环境信息工程重点实验室，江苏徐州 221116）

0 引言

煤炭作为我国的主体能源，在国民经济和社会发展中起着至关重要的作用。济宁市煤炭资源丰富，开采历史悠久，辖区内查明的煤炭资源储量达140 亿t；经过长期的开采，地下已形成面积庞大的采空区。近年来，城市建设快速发展，土地资源日益紧张，煤矿沉陷区土地成为可持续发展的重要资源，故而导致地下采空与地面建设间的矛盾日益尖锐［1］。因此，针对该地区进行高精度和大覆盖面积的采煤沉陷地调查与采空区稳定性评价工作，有利于缓解两者间的矛盾，对城市的可持续发展亦具有重要的现实意义。

传统地表形变监测技术（水准测量、GPS、全站仪等）的测量精度高，但却存在监测周期长、成果不连续、空间覆盖率低的缺陷［2－3］，不易于实现大范围采煤沉陷地的调查。而当前三维激光扫描与无人机低空摄影测量的形变监测精度相对较低，难以满足采空区稳定性评价的要求［4］。随着高分辨率、短重访周期SAR 卫星的连续升空以及数据处理手段的不断发展，InSAR 获取地表形变的精度已达毫米级［5］，使该技术在煤矿区地表变形监测领域得到广泛应用［6－13］。

笔者选用短基线集技术（small baseline subset，SBAS），通过2018 年9 月至2019 年3 月的14 景Sentinel－1 数据开展济宁城市规划区的地面形变解译，获取了研究区地表形变的时间演变过程和空间分布特征；并结合相应时间段内的地表实测数据进行了解译精度的评价。研究表明，SBAS 能够实现采煤沉陷地毫米级精度的形变监测，且能够显著区分不同开采条件下地表残余变形的特征。其成果有利于提高煤炭资源采空区地质调查的有效性，能够为采空区稳定性评价提供数据参考。

1 SBAS 基本原理

SBAS 技术的提出是为削弱时空失相干影响，提高InSAR 地表形变解译精度。 BERARDINO、LA⁃NARI 等［14－15］先后探讨了该技术在地表形变监测中的应用，而后HOOPER 等［16］又对处理过程中高相干点的选取方法进行了优化。

SBAS 技术的基本思想是通过设置时、空基线条件，将同一研究区的多期SAR 影像分组成若干个短基线集合；利用最小二乘准则，获取每个短基线集合的地表形变时间序列；通过奇异值分解方法再将分组后的短基线集进行联合求解；进而得到研究区整个观测时间段内时序的地表形变信息［17］。具体可表述［18－19］为：假定同一研究区成像时间为t0，t1，t2，…，tn－1的N幅SAR 影像，在短基线距的条件下可形成M幅差分干涉图，且满足N／2≤M≤N（N－1）／2。对于干涉图i，在去除平地相位后，干涉图中任意位置（x，y）的干涉相位可表示为：

其中，tA、tB为干涉图i所对应SAR 影像的获取时间；δφdef为tA～tB时间段内地表视线向形变相位；δφε为目标区地形相位；δφα、δφn分别为大气相位和噪声相位。其中前3 项可以表示为

式中，λ为雷达波长；R为斜距；B⊥为垂直基线；Δz为DEM 高程差；θ为入射角。

若不同干涉图间的形变速率为vk，k＋1，则tA～tB的累积形变量可表达为

对M幅干涉图进行三维时空相位解缠和地理编码后，便可求出不同SAR 影像获取时间的形变速率。

2 研究区域与数据来源

2.1 研究区概况

研究区域为济宁城市规划区，面积约为950 km2，包括任城区、兖州区，邹城市、曲阜市、嘉祥县的街道和部分乡镇。研究区内矿井众多，井田范围与城市规划区重叠面积约559.36 km2。研究区的空间位置、范围和各矿井分布情况如图1所示。

图1 济宁城市规划区空间位置、范围及矿井分布情况Fig.1 Spatial location，scope and coal minesdistribution of Jining urban planning area

据统计，仅2018 年上述各矿井的煤炭资源产出量即超过3 500 万t。煤炭资源的长期大范围开采，导致地下已形成面积庞大的采空区。为缓解城市建设与地下采空之间的矛盾，开展采煤沉陷地的调查与老采空区的稳定性评价工作迫在眉睫。

2.2 数据来源

用于研究区形变信息提取的SAR 数据为2018年9 月至2019 年3 月的14 期Sentinel－1 影像，C 波段，地面分辨率为5 m×20 m，相关参数见表1。 Sen⁃tinel－1 影像的轨道参数需根据影像的成像时间下载，欧空局提供了哨兵影像的Precise Orbit Ephem⁃erides（POD 精密定轨星历数据），此外，为消除地形相位，外部DEM 采用的是美国宇航局SRTM3 数据。

表1 济宁Sentinel－1 影像数据相关参数Table 1 Relevant parameters of sentinel－1 images in Jining

3 数据处理与形变解译过程

SBAS 的数据处理主要包括：差分干涉对的生成、点目标选取、差分干涉图的解缠，以及时间形变序列的获取等几个步骤，详细的处理流程如图2所示。

图2 SBAS 地表形变解译流程Fig.2 Flowchart of SBAS method

本次SBAS 地表形变解译采用的是SARscape软件。数据处理过程中参数设置如下：①连接图生成时，为充分利用影像数据，时间基线设置为200 d，空间基线设置为临界基线的10%；共连接成91 个像对，其中时间基线最长为180 d，空间基线最长为178.49 m（平均62.32 m）。 ②干涉处理过程中距离向与方位向的视数比设置为4 ∶1，影像配准利用精密轨道数据采用强度互相关算法，滤波方法选用自适应滤波法，解缠方法采用最小费用流法（相干性阈值设置为0.3）。 ③轨道精炼和重去平处理过程中，地面控制点选择相干性高、相位好的点，并且所选控制点位置的形变需为0。④第一次解译反演选用较为稳定的线性模型。其余参数均选用SARscape 处理Sentinel－1 数据的推荐参数。

最终解译的地表形变速率分布和各时间段内的地表垂直位移情况分别如图3 和图4 所示。

图3 研究区地表形变速率分布Fig.3 Distribution of surface deformation rate in study area

4 监测结果分析与精度评价

4.1 监测结果分析

图3 和图4 分别展示了沉陷速率大于10 mm／a的区域和各沉陷位置在2018－09－12—2019－03－11的演化过程。可以发现：在SAR 影像监测期间济宁城市规划区范围内产生了10 余处明显的沉陷盆地；通过叠加相应时间段内各矿井地下开采情况（图3），确定所监测的地表变形区域与开采工作面具有很好的一致性。验证了SBAS 技术在监测地下煤层开采所引起地表形变的可行性，其监测结果可以有效提高采煤沉陷地调查的针对性。通过对图3 和图4 进一步分析可以得到以下结论：

1）地表形变速率超过100 mm／a 的位置主要集中在存在地下开采活动的各矿井，如：唐口煤矿、新河煤矿、济宁二号煤矿、何岗煤矿、许厂煤矿、岱庄煤矿、运河煤矿等。

2）各矿井范围内开采位置以外，仍监测到大面积形变速率大于20 mm／a 的区域（主要为唐口煤矿和济宁城市建成区）；对于唐口煤矿，监测所得的地表形变速率约为20 mm／a，其原因主要是该区域为深部开采，前期开采影响即老采空区的残余变形仍较为显著；对于济宁城市建成区，监测所得的地表形变速率约为15 mm／a，因该区域不存在地下煤炭资源的开采活动，分析地表变形原因除形变解译误差外，主要与城市建设和地下水开采有关。

3）研究区域以外的其他位置，地表沉陷速率大都低于10 mm／a，表明相应位置的老采空区已基本趋于稳定。

4.2 精度评价

Sentinel－1 数据SBAS 时序处理结果的精度是通过与相近时间段内的水准测量结果的比较进行评价。项目组为评价SBAS 解译结果的可靠性，在研究区内进行了水准测量工作。本次精度评定，以唐口煤矿和岱庄煤矿的3 条观测线（TKL4 测线、TKL9测线、DZL15 测线）的95 个观测点的外业水准结果作为参考，以10m 范围内的高相干点的平均值作为InSAR 形变监测结果。其中，唐口煤矿TKL4 测线、TKL9 测线和岱庄煤矿DZL15 测线水准实测时间段分别为2019－01－06—2019－03－08、2018－11－25—2019－01－07、2018－12－26—2019－03－01；InSAR 监测时间段分别为2019－01－10—2019－03－11、2018－11－23—2019－01－10、2018－12－29—2019－03－11。

图5 和表2 为各测线的InSAR 解译与水准测量的对比结果与统计情况。

图5 各时间段地表垂直沉陷图Fig.5 Surface subsidence curve in different periods

表2 SBAS 地表形变解译精度统计Table 2 Statistical surface deformation accuracy

从图5 中各测线SBAS 解算结果与水准测量成果的对比情况可以看出：当地表形变在InSAR 可探测的量级范围内时，两者在沉陷的趋势上基本保持一致，具有较高的相关性。表2 中的统计数据显示，TKL4、TKL9 和DZL15 三条测线上两者最大差值分别为12.1、3.0 和11.6 mm，平均差值小于3.3 mm；监测精度高，但所能监测的最大沉陷量有限。

综合以上分析可以得到，SBAS 技术所监测的研究区沉陷盆地的状态、位置和沉陷量，具有良好的准确性和可靠性，监测结果有利于指导矿山地质环境调查工作高效、有针对性地开展。

4.3 老采空区形变监测可行性分析

前述SBAS 地表形变的监测结果和精度分析表明，InSAR 时序处理方法可以获得研究区地表形变的时空演变过程，且其形变监测精度能够达到毫米级。因此，该技术完全可以应用于老采空区地表残余形变的监测，而且对于变形缓慢的老采空区地表，SBAS 应更具有优越性。

另一方面，为验证SBAS 技术鉴别不同停采时间和不同煤层采厚老采空区的能力，做了进一步分析。以唐口煤矿为例，其中图6 为不同停采时间工作面上方地表点的时间演变曲线，图7 为不同煤层采厚条件下地表点的时间演变曲线。

从图6 和图7 中可以得出：①老采空区上方不同位置的地面点，在SAR 影像监测期间，仍然存在着持续的地表下沉现象。 ②随着工作面停采时间的增长，地表下沉量逐渐减小；其中，停采时间超过2、5、10 a 的采空区上方在180 d 内发生的地表残余形变量分别小于20、15 和10 mm；满足监测点连续6个月累积下沉值小于30 mm 的指标［20］，表明停采时间超过2 年的采空区上方地表已趋于稳定。 ③老采空区上方地表的残余变形量，随煤层采厚的增加而增加，且两者近似成比例关系。

上述结论表明，基于SBAS 的地表形变解译结果，不仅能够区分不同停采时间的采空区，对于不同采厚的采空区亦能有所区别。同时，由于SBAS 技术毫米级的形变监测精度，故其解译结果可以作为老采空区稳定性评价的依据（图6、图7 中图注，例如：2008（3.50）含义为2008 年开采，开采厚度为3.50 m）。