基于时序InSAR技术的新井煤矿地表沉降监测

马 飞，李日印，刘小鹏，侯景怡，崔喆森，张垚杰

（1.长治学院计算机系，山西长治 046011；2.自然资源部第二地形测量队，陕西西安 710054）

大规模开发利用矿产资源会直接影响矿区安全[1]。例如：2023 年2 月22 日，内蒙古新井煤业有限公司露天煤矿发生大面积边坡失稳事故，造成人员财产损失。因此，对矿区地表变形进行监测具有重要的现实意义[2-4]。

差分干涉合成孔径雷达(DInSAR)、小基线集In-SAR(SBAS-InSAR)和永久散射体InSAR(PS-InSAR)技术都是监测采矿沉陷、冻土沉陷、建筑物沉陷和滑坡的常用方法[5-7]，具有空间覆盖面大、成本低、精度高、全天候观测能力等优点，而常规测量技术（如水准测量[8]、三维激光扫描[9]和全球导航卫星系统[10]等）无法做到这一点。为了满足开采沉陷长期监测的需要，削弱大气相位的影响，时序InSAR 提出了一种简单有效的方法，即通过对一系列展开的干涉图按时间跨度加权和平均来估计平均视线开采沉陷速度[11-12]。小基线集InSAR (SBAS-InSAR) 可以很好地估计变形结果[13-15]，在许多开采沉陷区得到了应用。现对内蒙古新井煤矿应用SBAS-InSAR技术，识别该矿区地表沉陷范围，对沉陷区下沉数值进行解算，分析灾前微小形变规律，以预防同类事故的再次发生。

1 SBAS-InSAR技术原理

SBAS-InSAR 通过设置特定的时间基线阈值来合并小的基线干扰对，根据最小参数准则执行最优的形变相位序列，以限制时间和空间消相干。假设已经采集了在时间序列(t0，t1，…，tn)相同区域中获得的（n+1）个SAR 图像。设置特定的时间基线阈值后，可以形成m个干涉对，且m满足：

SBAS-InSAR 在进行参数反演之前，需要将所有图像与超级主影像配准，并通过多视处理、地形差分和相位解缠来生成m幅差分干涉图图像。考虑在时刻ta和t（b假设ta

式中：φ(tb,x)和φ(ta,x)分别为相位形变值；λ为雷达波长；φn为随机噪声；d(tb,x)和d(ta,x)分别为卫星视线方向的形变量；v(x)为像元的形变速率；Ti为干涉对的时间基线长度。

这样可以组成m个方程组：

式中：A为m×n的系数矩阵，每行对应一幅干涉图，每列对应一个时间上的SAR 图像，主影像所在列+1，辅图像所在列-1，其余列为0。

如果m≥n，且A的秩是n，则最小二乘解得φ的估计值：

由式（4）可知，SBAS 技术的关键在于求时间序列干涉图间形变速率的最小范数最小二乘解，可以利用时序干涉图中解算出的形变相位值，作为时间域的约束条件解，简化形变速率解算的复杂度，达到反演形变速率的目的。

假设V为形变速率，P为形变参数，以下约束条件被满足：

将式（5）代入BV=δφ得：

在干涉图中任意像元（x,r）的解算模型转换为：

式中：LOS方向形变的平均加速度变化率∆、相位平均速率、平均加速度用下式即可解得：

2 研究区概况及实验结果

2.1 研究区概况

图1为内蒙古新井煤业露天矿区地理位置，该矿区位于贺兰山南段，矿区面积1.344 8 km2，地势北高南低，海拔高程1 159.03～1 490.70 m。矿区为低山丘陵区，属高原侵蚀性丘陵地貌，地形切割一般，基岩裸露，植被稀疏，为荒漠地区。气候属大陆性干旱荒漠气候，日温差大，夏季炎热，冬季严寒，气候干燥，雨量稀少，年最高气温38℃，最低-28℃，平均8℃～9℃，年降水量为146～ 198 mm，年蒸发量最高2 000 mm以上。地表植被属于荒漠化草原向草原化荒漠过渡地带。地表多沙质化、砾石化和有龟裂结皮。

图1 新井煤业露天煤矿地理位置

通过分析新井煤业露天矿区地表植被和地形地貌特征，该矿区地表属于荒漠地区，植被稀少，适用C波段数据，但地面建筑物、工矿构筑物、大型岩石等永久散射体物较少，适合采用SBAS-InSAR对矿区地表沉降情况进行监测。

2.2 SAR数据选取及预处理

表1为覆盖该矿区的Sentinel-1A数据列表，包括卫星拍摄时间、影像模式、极化方式和轨道方向等信息。根据欧空局(ESA)网站查询哨兵-1A 数据情况发现，该卫星2022年能够覆盖研究区的影像只有15景，观测时间是2022年1月到2022年11月。根据干涉处理实验需要，下载了该区域干涉宽幅模式数据，单幅影像幅宽为250 km，距离向分辨率为4 m，方位向分辨率为13 m，多视比例设置为3∶1。同时下载到所有影像的精密轨道数据，用于干涉数据处理过程中的图像配准、去平地相位等流程，提高数据处理精度。另外，为提高数据处理效率，对原始SAR影像进行了裁切。

表1 覆盖矿区的Sentinel-1A 数据列表

2.3 时序InSAR 技术数据处理流程

2.3.1 短基线集干涉图组合

图2 为干涉对时空基线图。将选取的15 景影像按照如下策略进行自由组合，形成干涉对：时间基线ΔT≤100 d，空间基线ΔB≤200 m，以20220401 为超级主影像。最终得到43个干涉对。

图2 SAR干涉对时空基线图

2.3.2 干涉图生成

图3 是在新井煤业露天矿区生成的以20220401为主影像的差分干涉图，该图像为SAR 斜距坐标系统。图中可以看到，矿区中心位置形变条纹非常清晰，随着时间的推移，矿区中心位置形变量在不断变化；时间基线越长，干涉条纹周期变化越多；在图(i)、(j)等干涉对的中心位置出现干涉条纹叠加和失相干。

图3 以20220401为主影像生成的的干涉图

2.3.3 SBAS时间序列形变和形变速率结果

图4 为新井矿业露天矿区形变速率图，红色代表形变方向与卫星视线方向一致，可解算为下沉值，蓝色代表形变方向与卫星视线方向相反，可解算为抬升值。图中可见，该矿区大部分区域相对稳定，只有在此次滑坡事故中的核心区发生剧烈形变，整体形变速率非常大，沉降中心位置最大形变速率达到-160～ -209 mm/a，说明此处滑坡2022 年开始已经出现不稳定状态。另外，在该区域的东南侧也出现一个小范围的形变区，该区域的形变速率在-60～ -110 mm/a，需要引起相关部门的高度关注，避免发生新的地质灾害。

图4 灾变核心区形变速率图

为进一步分析新井煤业露天矿区核心区灾前形变情况，将该形变区域数值提取出来并进行处理。该区域共解算得到2 475 个形变点，图4 显示该区域自边坡的坡顶向坡底发生形变量由大到小，2022年1月至11 月坡顶累计形变量超过170 mm，坡底累计形变量超过150 mm。

图5 为边坡在2022 年的形变速率剖面图，横坐标为距离边坡坡顶GCP1 号点的距离，纵坐标为形变速率值。为进一步分析该矿区边坡失稳前的形变规律，在边坡选取了四个点命名为GCP1～4进行分析，分别分析其剖面线和形变时序图。图中可知，边坡顶部下沉速率远远高于坡底位置，自坡顶向坡底形变速率逐渐下降，形变速率从-180～-140 mm/a。由于该区域边坡角度较大，坡顶区松散层不断下沉，经过一段时间的累积，需要采取任何进行边坡加固。

图5 形变速率剖面图

图6 为新井煤业露天矿区灾变核心区内边坡形变点的时序下沉变化情况。由图可知，在边坡顶部的GCP1 号点累积形变量最大，在底部的GCP4 号点累积形变量稍小，这与图5 中的剖面图一致；该边坡从2022 年2 月开始出现下沉，2022 年3 月至6 月累积形变趋势开始增大，这可能与该地区前期转暖冻融解冻有关，也可能与矿区露天开采活动有关。

3 结论

(1)利用SBAS-InSAR 方法可以有效识别出露天采矿区地表变形，这是一种快速、高效、低廉的技术，可以为矿区地表边坡失稳、地面塌陷、地面沉降、地裂缝、滑坡、崩塌等地质灾害提供预警。

(2)内蒙古新井煤业露天采矿区核心灾变区在2022 年1-11 月持续发生下沉形变。分析SBAS-In-SAR 方法获取的2 475 个形变点发现该区域最大形变量超过160 mm，沉降速率超过200 mm/a；在核心灾变区东南侧亦有一个缓慢形变区，累积形变量稍小，但需要有关部门关注，避免发生新的地质灾害。

(3)内蒙古新井煤业露天采矿区失稳边坡在灾前一年间持续发生下沉，且边坡坡顶形变速率高于坡底速率，坡顶形变速率超过180 mm/a，顶部持续微小的下沉形变是导致边坡应力失衡的重要原因。

(4)由于缺少地表监测点，故SBAS-InSAR 的形变监测值无法得到验证，但该技术获得的时间序列累积变形能更直观地反映沉陷盆地的时间序列变化趋势，在探测变形区域范围和形态方面具有优势，有助于定性确定形变类型。

致谢：本实验的Sentinel-1A 数据和SRTM 外部DEM数据均来自欧空局网站，在此表示感谢。