基于InSAR的矿区地表三维形变场提取方法应用*

杨 隽 汪 娟 王 瑞 郭金城 冯 瀚 向淇文

（1.贵州省第一测绘院；2.贵州省北斗导航位置服务中心；3.赣南科技学院）

InSAR 是一种新兴的雷达遥感技术，近20 a 来得到了迅速发展，能够高精度地获取地表形变信息［1］。相较于传统方法，该方法具有空间分辨率高、覆盖范围广、全天候监测和不受天气影响等优势，可以有效地揭示区域性地表宏观趋势［2］。许多国内外学者已经对InSAR在矿区监测方面进行了大量研究，比如杨泽发等［3］利用水准测量和Levenberg-Marquard 算法模拟矿区动态沉降过程，陈磊［4］提出了基于Konthe模型的矿区大尺度形变信息的预计方法。

然而，在矿山开采过程中地表沉陷是三维的，传统的InSAR只对视线向的形变敏感，对于东西向和南北向变形监测精度不高。为获取地表真实变形，本文提出了一种单视线向解算方法，结合开采沉陷梯度关系来提高对东西向和南北向变形监测的精度。该方法利用开采引起的水平移动与倾斜之间的数学关系，对反演三维形变场提供约束条件，并构建解算模型，通过该模型可以对开采矿区地表进行三维形变提取，不仅能够探索矿山开采过程中地表形变的规律，而且为矿山的安全生产和地表灾害的预测预防提供了可靠的技术支撑。

1 InSAR视线向模糊问题分析

合成孔径雷达自身侧视成像的工作模式存在着视线向模糊的问题，即地表东西（E-W）、南北（S-N）、竖直（U）三维方向上的真实的地表形变分量为视线向（LOS）上的投影，这使得观测到的视线向形变并不是地表的真实形变［5］。将地表点的变形看成在垂直、东西、南北（U，E，N）3个方向分量的组合，并且3个方向的分量对于LOS 形变的贡献各不相同，图1 为In-SAR视线向形变量与地表三维形变场之间的关系，αh表示卫星飞行方向与正北方向的夹角，θ表示雷达的入射方位角。

假设地面目标视线向形变远离雷达为正，靠近雷达为负时，则dLOS与东西dE、南北dN和垂直dU向形变量之间的关系可以有：

由式（1）可知，利用一维形变求解三维形变会使得计算矩阵产生秩亏，因此，需要有不同视角的雷达观测量进行联合解算，但是在实际应用中存在许多制约和不足，如需要3个不同SAR传感器的观测值比较难以满足；此外该方法还假定观测目标在多平台观测的时间间隔内是不变或者线性变化的，这对于矿区的非线性形变来说，受到了一定的限制。

地下开采引起的岩层及地表移动过程是一个复杂的时间-空间过程，在移动盆地内地表各点的移动方向和移动大小都不相同，岩层的移动会逐渐传递直至地表，使得地表发生下沉和水平移动，形成一个三维形变场。这导致了InSAR 在矿区地表形变监测中无法准确地获得地面形变特征和规律。

2 基于开采沉陷知识的三维形变提取方法

2.1 开采引起的水平移动与倾斜的关系

本方法将结合采矿区开采沉陷规律，利用开采引起的水平位移与垂直位移梯度之间的函数关系作为约束条件，更好地将差分干涉图反演为地表真实的三维形变场。

随机介质模型是开采沉陷工程中预测采动位移最常用的方法之一［6］，在这个模型中，一个地下采掘面可以分成许多无穷小的采掘单元，岩体内部（x，z）点受单元开采影响产生的水平移动量Ue（x，z）和下沉量We（x，z）如图2所示。

根据Knothe，地表单元下沉函数服从正态分布：

式中，z为开采深度，m；x为开采边界投影到地表的距离，m；rz表示在z处开采影响的球体半径，m。

根据随机介质模型可知，无穷小的采掘单元引起的地表位移等于其位移之和，假设工作面板上的岩体不可压缩，则在二维平面上x、z方向上的法向应变之和为零，通过数学推导可得到

式中，b为水平位移常数；β水平位移主要影响角，（°）；T(x,H)为沿x方向的倾斜。

式（3）表明，在x-z平面上，地表点的水平位移Ue(x,H)与垂直下沉的水平梯度成比例，同理，该方法在y-z平面也同样适用。

2.2 SAR三维形变场解算模型

对于矿区而言，开采所影响的地表点在水平位移与垂直下沉之间存在一定的函数关系，可以使用单对干涉差分图反演地表三维变形，但是由于合成孔径雷达失相干导致视线向形变场的缺失，将导致解算模型失去连续的梯度从而给结果引入误差，因此需要通过克里金插值来补充形变场的空值区域［7］。

假设地理编码后的视线向形变场维数为n×m，则点（i，j）（i=1，…n；j=1，…m）的东西TE（i，j）和南北向TN（i，j）的倾斜可以表示为

式中，W（i，j）为地表点的下沉；ΔE、ΔN分别为差分干涉图在东西向和南北向的空间分辨率。

由于当前SAR 卫星大多都是侧视成像，所以SAR影像的像元阵列如图3所示。

联立式（3）、式（4），可以进一步推导EW、SN向的水平位移为

式中，H（i，j）为开采深度。

在解算地表三维变形的过程中，首先利用dLOS求解W（n，m），再将其代入到前面的方程直到求解出W（1，1），然后可以求解出东西向和南北向的水平位移。单轨干涉合成孔径雷达反演矿区地表三维位移场的处理流程如图4所示。

3 应用实例

3.1 研究区概况

万年矿位于峰峰矿区北部，该矿区横跨武安市伯延、磁山、午汲3个乡镇，北距武安市10 km，东距邯郸35 km。矿井的设计生产能力为180万t/a。在井田范围内，井田南北长约9 km，东西宽1～4 km，面积约21.167 km2，矿井的开采标高为325～-750 m。研究区在万年矿范围内，属于伯延镇辖区内，如图5所示。

3.2 试验处理

本次试验使用的是Sentinel-1A 卫星于2020 年2—4 月间获取的峰峰地区的两景数据，SAR 数据的基本信息见表1，由于VV 极化方式对于地表的穿透更强，因此试验采用了该极化方式成像的数据。

试验中选择2019 年11 月的数据作为主影像，2020 年2 月的作为从影像，2 次影像获取的时间基线为48 d，经过主辅影像配准、多视、干涉图生成、去平地效应、噪声滤波、相位解缠、地理编码［8-10］等处理步骤，最终得到研究区地表视线向形变，由于受地表不同地物相干性的影响，导致相干性低的区域出现形变场的缺失，这会使得图像像素失去连续性，从而影响模型的解算精度。因此，在反演研究区的三维形变前需要对LOS 向形变缺失区域进行插值，本文选用克里金插值法，利用方差图作为权重函数，基于最小二乘法对数据点进行空间建模和插值，如图6 所示。

3.3 三维变形反演

研究区位于万年矿东边，靠近井田边界，含煤地层主要为太原组、二叠系的山西组和石炭系的本溪组，根据万年矿以往的开采实测资料，可以获得该地质条件下的开采沉陷参数见表2。

通过使用基于开采沉陷知识的三维形变提取算法，将单视线向LOS 形变代入模型计算，并对其进行迭代求解得到三维形变场（图7）。

图7（a）为南北向水平移动，向北方向移动为负值，最大移动量为31 mm，向南方方向移动为正值，最大移动量为42 mm；图7（b）为东西向水平移动，向东方向为负值，最大移动量为35 mm，向西方向为正值，最大移动量为27 mm；图7（c）为垂直向下沉，最大下沉为61 mm。可以看出，利用该模型提取的三维形变场基本符合矿区的开采沉陷规律［11］，在沉陷盆地中，东西向和南北向的变形是向盆地中心进行移动，而垂直向下沉为一个明显的沉降漏斗。

在垂直向下沉中，出现了较为清晰的沉降阻隔线，在东西向水平移动变形图中，该区域在远离盆地的情况下出现大幅向西方向的移动，将在监测期间获取的三维形变场进行矢量合成（图8）。结合地质资料合理推断是由于工作面在采动过程中岩层的静态结构发生变化，影响到地质弱面断层的稳定性所导致的，断层作为岩层中的弱面，力学强度远远低于周围的岩体［12-13］，当开采的影响作用于断层时，其活化会使得岩体内部的形变传递变得异常。

4 总结

本文采用“理论分析—模型构建—实例分析”的主体研究思路，针对传统单轨InSAR只对视线向形变敏感，而对于东西向和南北向变形监测精度不高的不足，分析了InSAR 技术在矿区监测的局限，并结合矿区开采沉陷中岩层内部的梯度关系构建地表的三维解算模型。得到以下结论：

（1）通过SAR 侧视成像的几何关系，分析其在矿区地表形变监测的局限性，导致InSAR无法获得准确的矿区地表三维形变场。

（2）提出一种联合开采沉陷梯度关系的单视线向解算方法，通过开采引起的水平移动与倾斜的数学关系为反演三维形变场提供很好的约束条件，构建地表形变的三维解算模型。

（3）基于开采沉陷知识的三维形变提取方法，完成万年矿区的地表形变监测，提取了沉陷盆地的三维形变场，有效地揭示了该区域真实的地表形变特征，分析了地质弱面断层受采动影响活化使得岩层内部形变传递出现异常。