利用DS-InSAR技术考虑地形效应的山区煤田三维表面位移重建 李鹏宇、 魏少军、 夏晓明、 傅耀坤

《Processes》：3D Surface Displacement Reconstruction of Mountainous Coalfields Considering Topographic Effects Using DS-InSAR Pengyu Li, Shaojun Wei, Xiaoming Xia and Yaokun Fu

　　摘要：为了解决山区采矿区严重地表起伏、InSAR几何畸变严重以及无法直接从单轨视线（LOS）数据计算三维（3D）位移的问题，本文提出了一种3D变形重建方法，该方法将分布式散射体干涉合成孔径雷达（DS-InSAR）与改进的概率积分模型（PIM）相结合，同时考虑了地形滑动效应。传统的概率积分模型（PIM）忽略了地形坡度引起的额外滑动，因此在山区应用时会产生显著偏差。本研究引入了一个非线性滑动影响函数，并构建了一个地形校正模型，该模型结合了滑动强度、位置偏移和形态衰减参数，以定量描述采矿和地形耦合作用下的地表运动模式。基于此，建立了一个模型参数驱动的单轨InSAR观测方程，并采用自适应遗传算法（AGA）利用从DS-InSAR获得的高密度LOS变形作为约束条件，从而反演出整个盆地的3D位移场。来自中国太行山地区一个典型煤矿的实验结果表明，该方法有效克服了从单轨InSAR数据反演3D位移的病态问题。最大垂直沉降量为1050毫米，最大水平位移量为370毫米，与水准测量结果一致（垂直RMSE：75.1毫米；水平RMSE：27.2毫米）。与未进行地形校正的传统PIM方法相比，所提出的模型将3D位移的RMSE降低了约35%，显著提高了地形起伏严重山区的计算精度。通过对沿走向和倾向方向分布的水准测量点的验证，确认了重建的3D位移场的可靠性。该方法不仅恢复了地形滑动的物理特性，还为复杂地形下的采矿损伤监测提供了一种低成本、高精度的解决方案。

1. 引言
地下煤炭开采不可避免地破坏了岩层的原始应力平衡，导致地表沉降和移动，对采矿区的生态环境和基础设施构成严重威胁[1,2,3]。在平坦地形区域，地表运动模式已通过经典的概率积分模型（PIM）得到了很好的描述[4,5]。然而，中国的主要煤炭能源基地（如山西、陕西、贵州等）大多位于沟壑纵横、地形破碎的山区。在这些地区，地形起伏不仅改变了岩体运动的传播路径，还由于自重作用在斜坡上引起了显著的额外滑动变形，使得基于平坦地面假设的传统监测和预测理论失效。干涉合成孔径雷达（InSAR）技术凭借其全天候、实时和高精度的优势，已成为采矿区沉降监测的主流方法[6,7,8,9]。然而，InSAR固有的侧视成像几何特性决定了它只能获得雷达视线（LOS）方向的一维变形。尽管理论上可以通过整合上升和下降的多几何数据来计算3D位移，但由于卫星采集的限制，实际应用中通常使用单轨数据，这使得从一维观测反演3D位移成为一个典型的病态问题[10,11]。为了解决这些问题，一种有效的策略是将先验物理模型作为约束条件引入。近年来，一些学者尝试将InSAR观测与传统的PIM相结合[12,13,14,15]，通过参数反演间接获得3D位移。然而，大多数现有研究没有考虑山区特有的地形滑动效应，导致斜坡区域的水平运动方向与实际重力方向偏离，物理意义不明确。关于山区地表运动机制，一些研究已将滑动影响函数模型应用于山区沉降研究[16]。该函数通过引入非线性参数显著提高了对复杂斜坡滑动形态的拟合能力，但目前该模型主要基于稀疏的水准点数据，并未与高分辨率InSAR技术深入整合。鉴于此，本文提出了一个从山区单轨InSAR数据计算3D位移的新框架。本研究的主要贡献包括：（1）首次提出了一种结合非线性滑动影响函数的地形校正模型，用于山区采矿区的单轨InSAR 3D位移重建；（2）开发了一种结合DS-InSAR高密度观测数据和物理约束PIM参数的联合反演框架，实现了从单轨数据可靠地计算3D位移。本研究利用DS-InSAR技术克服了山区植被去相关性的限制；它结合了改进的PIM和滑动影响函数，并通过自适应遗传算法（AGA）解决了高维非线性参数反演问题。本研究的具体目标是：（1）开发一个包含非线性滑动影响函数的地形校正PIM模型，定量描述采矿-地形耦合作用下的地表运动模式；（2）建立模型参数驱动的单轨InSAR观测方程，并实施基于AGA的采矿和滑动参数的联合反演；（3）重建整个盆地的3D位移场，并在典型的山区煤矿对方法进行验证；（4）量化与传统方法相比的精度改进，并分析地形滑动参数的敏感性和适用性。

2. 研究区域和数据
2.1. 研究区域
研究区域位于吕梁山脉的东侧，地形复杂。西南部海拔较高，东北部较低，Tunlan河谷的最大海拔为1521.2米，最小海拔为1065.5米，最大相对高度差为455.7米，属于中山区。在研究区域为工作面建立了两条观测线（走向和倾向）。观测站区域的地面海拔范围大约在1280米到1400米之间，最大相对高度差为120米，平均坡度约为21%。该工作面的采矿方法是全机械化长壁后退开采，采用全层开采方式。工作面的地面海拔范围为1297米到1411米，工作面海拔范围为882.1米到943.2米，覆层厚度范围为402米到511米，平均值为463米。走向长度为1709米，倾向长度为176米，带状和轨道入口的总宽度为9米。煤层倾向角范围为1°到8°，平均值为2°。煤层厚度范围为1.66米到2.69米，平均值为2.25米。

2.2. SAR数据和辅助数据
本研究使用了欧洲航天局（ESA）提供的Sentinel-1A卫星SAR影像数据，时间跨度为2021年11月至2023年6月，共46个上升轨道数据集。卫星数据覆盖范围和研究区域位置如图1所示。使用分辨率为10米的数字高程模型（DEM）去除了地形相位影响。

3. 方法论
为了解决山区采矿位移中单轨InSAR 3D计算的地形滑动、去相关性和病态特性问题，本文提出了一种参数反演策略，该方法将分布式散射体InSAR（DS-InSAR）技术与考虑地形滑动效应的改进概率积分模型相结合。技术方法主要包括三个部分：基于DS-InSAR的高密度变形场获取、考虑地形滑动因素的改进PIM构建，以及基于AGA的模型参数反演和3D位移场重建。DS-InSAR在山区采矿区具有三个优势：首先，植被覆盖的山区斜坡缺乏稳定的人造结构，导致用于PS-InSAR的持久散射体（PS）数量不足。DS-InSAR利用自然地面表面的分布式散射体，显著增加了相干点密度。其次，DS-InSAR的统计均匀像素选择有效过滤掉了异质山区地形中的去相关像素，提高了相位可靠性。第三，DS-InSAR提供了空间连续的覆盖范围，这对于捕捉沟壑破碎地形中的完整采矿沉降盆地至关重要。

3.1. 基于DS-InSAR的山区地表变形获取
山区植被茂密，传统持久散射体InSAR（PS-InSAR）的稀疏点目标难以充分反映沉降盆地的特征。本文采用分布式散射体（DS）干涉技术，通过统计均匀像素过滤和相位优化算法来提高山区的相干点密度。

3.1.1. 信号模型和统计均匀像素选择
假设有两组共配准的SAR图像。每个DS像素的时间序列复杂观测向量可以表示为，表示第n幅图像的复杂观测向量。每个DS像素的协方差矩阵的最大似然估计（MLE）可以表示为[17,18,19]：
在上述方程中，是时间段内的观测数据标准差向量，是标准化观测向量，包含所有干涉对的相干性和干涉相位信息，的模是目标像素的相干系数矩阵。为了识别具有统计均匀性的邻近像素（统计均匀像素，SHP），使用双样本Kolmogorov–Smirnov（KS）检验或Anderson-Darling（AD）检验[20,21,22]。对于像素和，它们的累积分布函数分别为和，KS检验统计量为：
如果小于给定置信水平的阈值，则确定为的统计均匀像素。

3.1.2. 相位优化和变形计算
使用选定的SHP集合，对中心像素的协方差矩阵进行最大似然估计（MLE）。为了从噪声干涉相位中恢复最优相位序列，使用相位链接和特征值分解算法进行求解[19,21]。

3.2. 考虑地形滑动因素的改进概率积分模型
传统的PIM主要适用于水平地面条件下的采矿沉降预测，在山区应用时会产生显著偏差。根据山区地表运动的修正理论，山区的地表运动是采矿引起的理论水平条件下的运动与地形引起的滑动的矢量叠加。

3.2.1. 传统概率积分方法
假设煤层倾向角较小，根据传统PIM，在平坦地面条件下，地表点的沉降和水平运动计算如下：
其中是主要影响半径，是水平运动系数；和分别是走向和倾向方向上的等效采矿长度，由实际采矿长度、偏移距离、采矿影响传播角和煤层倾向角确定。

3.2.2. 适用于山区的改进PIM
在任何点的地形坡度和下倾角下，引入滑动影响函数来校正水平运动和沉降。校正后的3D位移场表示为：
其中和是由地形引起的滑动增量：
其中是地表特征系数（凹坡为负，凸坡为正）；是滑动影响函数值，是待预测的水平运动方向。
在上述方程中，是滑动影响函数的参数，可以用代替进行计算。

3.3. 基于遗传算法的参数反演和3D计算
由于单轨InSAR无法直接分解3D位移，本文构建了一个目标函数，通过反演上述模型参数来实现山区地表的三维场重建。

3.3.1. 单轨投影观测方程
建立待反演的模型参数集与InSAR观测值之间的映射关系：
3.3.2.**fitness函数构建**
为了找到最优参数，构建了一个基于某种范数的fitness函数，其目标是最小化模拟的LOS变形与DS-InSAR观测值之间的残差：(10) 反演过程如下：
(1) 初始化：在参数空间内随机生成初始种群。
(2) 正向计算：对于每个个体，将其代入方程(6)计算3D位移，然后投影到LOS方向。
(3) 遗传操作：计算适应度，执行选择、自适应交叉和变异操作。
(4) 收敛输出：获得全局最优参数。
(5) 3D位移重建：将其代入方程(6)直接计算地表沉降和位移，从而解决从单轨数据中进行3D计算的问题。
整体流程图如图2所示。

**4. 结果与分析**
**4.1. DS-InSAR LOS变形获取**
与传统PS-InSAR技术相比，DS-InSAR通过统计上均匀的像素过滤显著提高了山区点的密度。在研究区域内，共提取了75,212个有效相干点，并成功获取了植被覆盖的山坡和山谷区域的监测点。InSAR计算出的累积LOS沉降图显示，最大沉降中心位于工作面中间，最大LOS位移超过900毫米。需要注意的是，图3中与工作面相连的东北侧的变形是由另一个工作面引起的，而非所关注的工作面。因此，在后续的参数反演中仅计算该工作面影响范围内的位移。

**4.2. 参数反演**
11个参数的范围设定如下：(沉降系数)，(水平移动系数)，(主要影响角正切)，(拐点偏移距离)，(采矿影响传播角)，等。采用自适应遗传算法来反演山区的PIM参数。算法参数设置如下：种群规模150，最大迭代次数150，收敛阈值设定为连续20代且适应度变化小于某个值，交叉概率0.9，变异概率0.01。反演得到的主要参数如下表1所示。表1. 参数反演结果。AGA在所有独立运行中均在120-135代内一致收敛。最终适应度值（RMSE）稳定在约50.3毫米，表明模型预测与DS-InSAR观测结果吻合良好。
反演得到的参数与采矿工程的实测值一致性良好。沉降系数=0.78，表明覆盖层条件为中硬，这与研究区域的地质剖面相符；=1.81处于中国煤矿的典型范围内；采矿影响传播角（=89.2°）与浅层采矿深度与厚度比（H/m≈206）一致。

**4.3. 3D位移计算结果**
基于反演得到的参数，计算了研究区域的3D位移场，如图4所示。结果显示，最大垂直沉降为1050毫米，位于工作面上方；最大水平位移为370毫米，朝向沉降盆地中心。为了验证这些值的物理合理性，我们将其与理论预测进行了比较。对于全层开采（平均厚度m=2.25米）且沉降系数q=0.78的情况，理论预测的最大沉降为1755毫米。但由于该工作面的宽深比较大，根据不完全开采校正，实际最大沉降约为950-1050毫米，因此实测的最大沉降1050毫米在预期范围内。最大水平位移370毫米对应的水平移动系数b约为0.35，符合中国煤矿的典型范围（0.20-0.35）。由于山区地形的影响，水平移动系数往往较高。斜坡上的水平移动向量会向坡下方向偏移，偏离理论方向（朝向盆地中心）。在坡度超过25°的情况下，滑动位移分量ΔU_slip的最大值为85-120毫米，占总水平位移的40-55%。这表明忽略地形滑动会导致对水平移动方向的严重误解，尤其是在陡峭地形中。

**5. 讨论**
**5.1. 地形滑动参数的敏感性分析**
本研究中引入的滑动参数对反演结果有显著调节作用。分析表明，滑动强度系数与平均地表坡度呈正相关。当坡度小于10°时，反演值接近1，此时可以忽略滑动效应；当平均坡度超过30°时，如果使用传统模型，即，从InSAR反演得到的沉降系数会因试图通过增加沉降来补偿观测到的较大LOS变形而产生误估。本文提出的方法通过将采矿沉降与地形滑动解耦，恢复了真实的物理过程。

**5.2. 滑动影响函数的优点**
与简单的线性坡度校正相比，非线性滑动函数通过引入位置偏移参数，能够准确描述滑动峰相对于采空边界的不对准现象。实验数据的残差分布显示，线性校正模型在采空边界处仍有较大残差，而本文采用的非线性滑动函数将边界处的拟合残差减少了约40%，证明了该数学模型的优越性。

**5.3. 反演结果的唯一性和稳定性分析**
反演结果的唯一性和稳定性是非线性参数估计中的关键问题。为确保获得可靠且唯一的解，本研究采取了以下策略：
(1) 参数搜索空间严格受限于来自采矿工程经验的物理意义范围内的值（表1）。沉降系数q根据覆盖层岩性限制在[0.4, 0.9]之间；主要影响角的正切tanβ根据岩体性质限制在[1.5, 3.0]之间；水平移动系数b根据与采矿深度的经验关系限制在[0.1, 0.4]之间。
(2) DS-InSAR提供了75,212个相干测量点，形成了高度过定的系统，大大降低了多解的可能性。
(3) 自适应遗传算法（AGA）采用自适应交叉和变异概率，能够高效探索整个参数空间，并降低收敛到局部最优解的风险。

**5.4. 错误来源分析**
所提出方法中的残差误差来源于多个方面：
(1) InSAR数据中的大气相位残差，估计对RMSE的贡献为10-20毫米；
(2) DEM误差导致地形相位去除不准确，尤其是在陡峭斜坡上；
(3) 滑动影响函数的简化可能无法完全捕捉复杂的三维坡度力学。

**6. 结论**
为了解决山区单轨InSAR监测的局限性，本文提出了一种基于DS-InSAR和改进的滑动影响函数的3D位移计算方法。通过引入地形滑动机制，该方法成功解决了从单轨数据中分离3D变形的病态问题。使用滑动影响函数构建了一个考虑地形因素的正向模型，并通过自适应遗传算法实现了滑动参数和概率积分参数的联合反演。结果表明，地形滑动显著改变了地表运动场的分布模式，最大滑动位置沿斜坡向下移动。在典型山区矿场的验证表明，该方法重建的3D位移场的准确性明显优于传统方法，尤其是在高陡坡区域，3D位移的计算误差减少了约35%。本研究为缺乏多轨SAR数据的山区矿山提供了科学、经济且高精度的采矿沉降监测技术。
尽管本研究证明了所提出方法的有效性，但仍存在一些局限性。影响函数中嵌入的平行于坡度的滑动假设虽然简化了计算复杂性，但可能过于简化了控制高度不稳定斜坡的复杂三维破坏机制。此外，当前的静态框架没有明确捕捉变形的时间演变过程，可能忽略了对于早期预警应用至关重要的瞬态变形阶段。未来的工作将集中在结合时间依赖的流变参数来模拟采矿引起的变形的运动学演变，从而实现动态危险评估而非静态分析。