基于单轨道InSAR的矿区地表三维形变监测研究

姜志伟 常英豪 孔晓岚 李高峰

（内蒙古工业大学，内蒙古 呼和浩特 010052）

0 引言

长时间、大面积的地下煤层开采使矿区地表出现大量沉降，造成道路塌陷、房屋受损和水土流失等问题。对矿区地表沉降进行监测具有重要意义。D-InSAR 技术能够获取高精度地表微小形变信息，具有全天时、全天候、高空间分辨率等优势［1］。煤矿开采周期长，受时空失相干及大气延迟等因素的影响严重，因此发展出了时序InSAR 技术。SBASInSAR 技术通过设定时间基线和空间基线的阈值，从而获取研究区域的时序累积沉降和形变速率［2］。

InSAR 技术在地表形变监测领域应用广泛，但矿区形变梯度大、速率快的特点极易造成形变超出相位解缠的临界阈值［3］。本研究以时序D-InSAR技术和概率积分模型融合的方法对红庆河煤矿进行长时序地表三维形变监测，并对开采沉陷规律进行分析总结。

1 研究方法

1.1 时序D-InSAR融合概率积分模型方法

Baran 等［4］对影像相干性γ、InSAR 监测地表形变Dmax和InSAR 技术所能监测到的最大形变dmax的关系进行分析，得到三者间的关系式为式（1）。

假设有N组干涉对，统计影像的平均相干性，根据式（1）可得该条件下的最大形变梯度，基于形变梯度理论求得N组干涉对在InSAR 监测中的最大累积形变量dmax/InSAR为式（2）。

根据上述公式得到的dmax/InSAR作为时序DInSAR 结果与概率积分模型融合的阈值，从而得到完整的矿区地表沉降［5］，为式（3）。

将概率积分模型结合观测线上的实测数据进行加权处理，使概率积分法的沉降结果更趋近于对应时间段内的实测数据。时序D-InSAR 融合概率积分方法的处理流程如图1所示。

图1 时序D-InSAR融合概率积分方法的处理流程

1.2 单轨道矿区三维形变解算模型

对于SAR影像的像元（i，j），东西向和南北向的梯度可以表示为式（4）［6］。

式中：ΔE和ΔN为SAR影像中单个像元在东西向和南北向的长度，由式（4）得到东西向形变和南北向形变为式（5）。

式中：μE=b·H（i,j）/tanβ·ΔE，μN=b·H（i,j）/tanβ·ΔN。μ中ΔN、ΔE和H与煤矿开采位置及SAR数据分辨率有关，b和tanβ为概率积分模型中的参数。

地面三维形变与SAR 卫星LOS 向的几何关系为式（6）。

式中：α为雷达运行轨道的方位角，θ为雷达入射角。将式（5）代入式（6）得：

式中：未知量为n×m个，而方程组的个数为（n-1）×（m-1）。基于矿区边缘形变量较小的经验，假设图像行和列的结尾不存在水平向形变，对于最后一行和最后一列表示为式（8）。

通过式（8）可得到垂直向形变W(i,j)，代入式（7）可得到东西和南北向形变。

2 研究区域概况

红庆河煤矿位于鄂尔多斯高原的东部，隶属于神东煤炭基地。区内植被较少，属沙漠～半沙漠地区。研究区域为3-1301 工作面（见图2），其走向长3 549 m，倾斜长300 m，工作面面积1 064 733 m2，赋存标高664～685.3 m，平均674 m，煤层倾角1°～5°。该区域属于高原侵蚀性低中山地貌特征，井田中南部存在多条冲沟，沟深5～20 m。

图2 红庆河煤矿3-1301工作面示意

3 试验结果及分析

3.1 时序处理结果

试验中数据为Sentinel-1A 卫星影像，其重访周期为12 天，空间分辨率为5*20 m。利用时序DInSAR 技术处理21 景Sentinel-1A 影像，共生成20组干涉对。卫星距离向和方位向的视数比设置为5：1；制图分辨率为20 m；采用Goldstein 滤波；相干性阈值设为0.3；通过最小费用流法进行相位解缠；最终经过地理编码得到研究区域LOS 向时序形变结果。除此之外，试验发现2020 年2 月25 日—2020年3月8日影像相干性较差，选择跳过副影像，将2020年2年25日和2020年3月20日的SAR影像组成干涉对，提高平均相干性，进一步增加形变监测精度。

3.2 时序D-InSAR和概率积分模型融合结果

时序D-InSAR 和概率积分模型融合方法如图3 所示，3-1301 工作面的部分参数见表3。根据上述公式并结合表中开采参数得到研究区域在2019年8 月15 日至2020 年4 月30 日间的开采沉陷模型。如图3（a）所示，黑色框线为3-1301 工作面位置，与开采沉陷模型中沉降区域一致。

表3 红庆河煤矿3-1301工作面开采参数

图3 时序D-InSAR和概率积分模型融合方法

根据上节所述方法对研究区域时序D-InSAR结果进一步处理得到完整的矿区时序形。如图3（b）、图3（c）所示，融合结果具有良好的连续性，弥补了时序D-InSAR无法监测到的空值区域；沉降区域与工作面位置一致，沉降范围受地形及周边工作面的影响向西南方向扩散，与时序D-InSAR结果表现一致。LOS 向形变无法反映地表的真实形变，还需解算三维形变。

3.3 三维形变解算

通过时序D-InSAR 和概率积分模型融合方法获取了研究区域在2019 年8 月17 日至2020 年4 月25 日间的LOS 向时序形变结果。在此基础上建立基于单轨道InSAR 的三维形变解算模型，提取垂直向、东西向和南北向的形变信息，并结合实测数据对比验证试验精度。

时序D-InSAR 融合概率积分模型的三维形变解算结果如图4 所示，垂直向形变范围为-548.0～21.7 mm，东西、南北向形变范围为-81.2～54.4 mm和-88.4～52.1 mm；图中融合结果表明在垂直向形变中沉降影响范围沿工作面开采方向逐渐减小；东西、南北向形变在工作面起始部分呈对称分布，但随着开采向东南方向推进，研究区域内大部分表现为南北向形变。结合该区域地形情况进行分析，工作面开采推进过程中地表存在大量东西走向的沟壑，造成地表形变沿南北向进行。

图4 时序D-InSAR融合概率积分模型的三维形变解算结果

3.4 试验结果分析

根据试验结果的形变范围选取走向观测线AB 上A1～A79 观测点的三维形变结果及实测数据绘制其沉降曲线，如图5 所示。实验结果在垂直向形变上与实测数据之间误差较大，在观测线前端属于植被覆盖区域，SAR 影像的相干性受到一定的影响，实测数据中在23 号观测点左右出现了阶梯式沉降，使沉降过程中多次在短时间、小范围内的出现大量级沉降，从而导致垂直向形变误差。

图5 走向观测线AB上的沉降曲线

4 结论

本研究提出时序D-InSAR 和概率积分模型融合的方法对红庆河3-1301 工作面进行监测，结果表明该方法能有效提高大量级形变区域的监测效果，通过对比融合结果和SBAS-InSAR 结果的三维形变，结合实测数据进行分析，实验表明融合结果和实测数据的均方根误差为：垂直向0.100 m、东西向0.018 m、南北向0.025 m，而SBAS-InSAR 结果的均方根误差为：0.173 m，0.016 m，0.021 m。三维形变解算效果良好，融合结果地表三维方向的精度均有提升，垂直向精度更高，水平向与SBAS-InSAR结果相接近。本研究监测结果与实测数据的沉降趋势基本一致，为煤矿地表沉降的监测提供了参考。