一种基于对称性的镜像重建方法，用于填补采矿沉降TLS（Tilt-Stacking）调查中的大型数据缺口

《Frontiers in Earth Science》：A symmetry-guided mirror reconstruction method for filling large gaps in TLS surveys of mining subsidence

　　摘要
煤炭开采引起的地表沉降是矿区最显著的人为地质灾害之一。在沉降盆地的地面激光扫描（TLS）调查中，由于不可避免的地表条件（包括复杂的地形、陡峭的斜坡和人造结构的遮挡），点云数据往往不完整。为了解决这个问题，提出了一种镜像栅格重建（MRR）方法，该方法利用近水平煤层中沉降盆地发展的固有对称性来重建缺失的DEM数据。为了提高在实际工程条件下的适用性，引入了一个系统偏差角（SDA）来考虑与理想对称性的偏差。选择208工作面作为研究区域，在多个时间段进行了协调的TLS和UAV LiDAR观测。首先对3月18日和4月14日获得的TLS和UAV LiDAR DEM进行差异处理，以得到相应的沉降盆地DEM。然后使用MRR方法重建4月14日TLS DEM中的缺失区域，并将重建结果与同步的UAV LiDAR观测结果进行验证。比较结果显示RMSE为0.23米，RRMSE为7.18%，证明了所提方法的可靠性。进一步使用6月15日在实际工程条件下获得的TLS数据评估了MRR方法的适用性。基于重建的DEM，反演了沉降盆地的关键角度参数，并与同时期的GNSS测量结果进行了比较。重建后，沿H线的拉伸角度与GNSS结果相差约7°，而沿C线的边界角度和拉伸角度分别相差约10°和5°。与传统的空间插值方法（如逆距离加权（IDW）和克里金法）相比，所提出的方法用受沉降力学约束的确定性重建过程替代了随机插值。该方法有效缓解了由于现场遮挡造成的TLS数据缺失问题，并为矿区常规的大规模沉降监测提供了可扩展的解决方案。它还为解决岩土工程遥感应用中的数据不完整性提供了一个基于物理的框架。

1 引言
地下煤炭开采破坏了上覆岩层的力学平衡，导致地表沉降盆地的形成，这对地表基础设施、生态系统和公共安全构成了严重威胁。目前，煤炭约占全球一次能源消耗的27%，而每年有超过150万公顷的土地因煤炭开采活动而退化（Arif等人，2025；Behera和Rawat，2023）。因此，准确监测沉降盆地对于土地复垦、保护地表结构和优化采矿计划至关重要（Jiang等人，2024）。采矿引起的沉降盆地的形成遵循既定原理。最大垂直沉降通常达到采矿深度的60%-90%，而水平位移可能超出工作面边界0.2-0.4倍。沉降盆地通常呈现碗状形态，边缘坡度陡峭，中心坡度相对平缓（Cai等人，2023；Zhu等人，2023）。然而，在山区矿区，地表地形与采矿引起的变形相互作用，导致更复杂的三维位移场，大大增加了监测难度（Lian等人，2024）。

沉降监测技术已从单一传感器观测发展到多源集成方法。全球导航卫星系统（GNSS）测量提供了高精度的基于点的观测数据，典型的水平和垂直精度分别为±1-2厘米和±2-3厘米（Tonda?等人，2023a；Tao等人，2019；Tonda?等人，2023b；Bo等人，2023）。干涉合成孔径雷达（InSAR）能够进行大尺度变形测绘，最近结合SBAS-InSAR和深度学习模型的研究报道预测误差控制在约1.44毫米以内（Yang等人，2020；Zhu等人，2024a；Wang等人，2022；Chen等人，2021）。然而，在盆地中心附近的大变形区域，InSAR测量常常遭受相干性损失，导致这些关键区域的可靠性降低（Wang等人，2022；Chen等人，2021）。地面激光扫描（TLS）因其高空间分辨率和密集的三维点云测量而成为详细沉降盆地表征的强大工具，实现了厘米级的变形精度（Matwij等人，2021；Gu等人，2020）。尽管如此，实际的TLS应用经常受到严重数据缺失的影响。Matwij等人（2021）报告称，尽管TLS点云配准精度很高，但数据完整性常常受到影响。同样，Gu等人（2020）证明，即使在王家塔煤矿采用无目标扫描策略，扫描盲点和大量数据缺失仍然不可避免。

TLS数据缺失由多种因素引起。首先，建筑物、植被和岩石露头等物理障碍物阻挡了激光传播，形成了阴影区。其次，茂密的植被限制了激光对地表的穿透，导致地形信息丢失（Blachowski等人，2024；Petras等人，2023）。第三，扫描几何形状的限制限制了对陡峭斜坡和凹陷的有效覆盖。第四，矿区广阔的空间范围加上有限的扫描站点数量，使得全面覆盖不切实际。最后，崎岖的山地地形加剧了遮挡效应，增加了盲区的普遍性（Lian等人，2024）。这些因素的综合影响导致TLS衍生的DEM中存在广泛的数据缺失，严重限制了沉降盆地形态的准确表示和关键变形参数的可靠提取。虽然UAV LiDAR可以从空中角度部分补充缺失的信息，其报告的精度约为10毫米（Zheng等人，2022；Yang等人，2023），但其运营成本和时间灵活性有限，限制了其常规应用。

为了缓解单一来源观测的局限性，多源数据融合已成为一个活跃的研究领域。Zhu等人（2024b）证明，整合D-InSAR、SBAS-InSAR和UAV观测可以获得比任何单一技术更可靠的结果。然而，大多数现有研究仍依赖传统的空间插值方法来填补TLS数据缺失。例如，Wang等人（2024）使用局部多项式插值整合了InSAR、UAV和LiDAR数据，获得了0.040米的RMSE。在大变形区域，GNSS-InSAR融合方法的误差范围为13至35毫米（Tonda?等人，202a；Zhou等人，2022；Tao等人，2025）。Chen等人（2013）结合D-InSAR和3D激光扫描，并采用逆距离加权填充不连贯区域，获得了低于5%的相对误差。其他常用的技术包括克里金插值（Chen J.等人，2024；Dong等人，2022）、薄板样条（Chen等人，2022；Li Y.等人，2022）和径向基函数（Meng等人，2019）。近年来，深度学习方法（如生成对抗网络（GANs）和基于Transformer的模型）也被用于DEM间隙填充（Qiu等人，2019；Gavriil等人，2019；Li S.等人，2022；Li等人，2023）。尽管方法多样，但这些方法都有一个根本的局限性：它们将采矿沉降盆地的数据重建视为纯粹的空间插值问题，而没有明确考虑控制盆地形成的地质力学过程。传统插值方法依赖于空间自相关假设，实际上构成了统计推断。虽然深度学习方法可以学习复杂的地形特征，但它们同样缺乏与沉降力学相关的明确物理约束。这种数据驱动而非物理驱动的重建范式导致了两个主要问题。首先，在变形梯度陡峭或拐点明显的区域，插值结果往往表现出过度平滑或虚假振荡，与真实的沉降行为明显偏离（Dong等人，2022；Meng等人，2019；An等人，2024）。其次，重建的表面可能违反采矿沉降的基本物理定律，削弱了关键参数（如拉伸角度和边界角度）提取的可靠性（Zhu等人，2023；Liu等人，2023）。

更根本的是，现有方法忽略了一个被采矿沉降理论反复验证的核心物理原理：在近水平煤层开采条件下，具有相对均匀地质性质的沉降盆地关于工作面中心轴表现出明显的两侧对称性（Zhou等人，2014；Liu等人，2023；Tang等人，2019）。这种对称性不是经验现象，而是岩石质量力学的直接结果。对称的采空边界产生对称的应力重分布模式，进而产生对称的地表变形场（Cai等人，2023；Zhu等人，2023）。作为中国主要的沉降预测框架，概率积分法（PIM）明确将盆地对称性作为核心假设（Zhou等人，2014；Liu等人，2023）。从这个角度来看，利用对称性作为物理约束——而不是忽视它——代表了TLS数据重建的逻辑和基于物理的策略。

基于对称性的约束的有效性已在多个地球科学应用中得到证明。Tang等人（2019）开发了一个三维位移反演模型，将沉降对称性纳入InSAR监测中，并表明对称性约束显著提高了反演精度。在冰川动力学研究中，对称性特征被用于数据稀疏区域的运动场重建（Milenkovi?等人，2019），而在滑坡监测中，几何对称性被用于约束变形反演（Armesto等人，2009）。这些研究表明，当变形系统表现出固有的物理对称性时，将这一属性作为先验约束可以提高重建可靠性，同时减少观测冗余。然而，将对称性约束应用于基于TLS的采矿沉降重建面临两个关键挑战。首先，由于地质异质性（如断层和褶皱）、非均匀的采矿进展和复杂的地表地形，实际沉降盆地会偏离理想对称性（Zhu等人，2023；Lian等人，2024；Jiang等人，2023；Pawluszek-Filipiak等人，2020），直接镜像重建可能会引入较大误差。其次，目前还没有成熟的算法框架能够在栅格DEM背景下高效且稳健地实现基于对称性的重建。这些挑战限制了对称性约束在TLS沉降监测中的实际应用。

为了解决这些问题，本研究提出了一种镜像栅格重建（MRR）算法，建立了一种将物理约束与工程校正相结合的新范式。MRR的方法论贡献有两个方面。首先，从理论角度来看，采矿沉降的对称性原则被明确引入作为TLS数据重建的物理约束。与基于空间接近性的传统方法不同，MRR基于这样一个力学前提：对称轴两侧对应位置的变形幅度应该是相等的。这种表述将一个病态的空间插值问题转化为一个受约束的确定性重建过程，确保重建的表面符合沉降力学而非统计假设。其次，从工程角度来看，引入了一个系统偏差角（α），将理想化的对称性约束扩展到实际地质条件。该参数通过全局优化确定，适应性地捕捉由地质结构、采矿进展和地形变化引起的沉降盆地的几何偏差。这种结合物理约束与工程校正的混合框架避免了纯物理模型的过度理想化，同时克服了数据驱动方法的非物理限制。它为基于TLS的沉降监测提供了一种稳健且可扩展的技术途径，对提高变形监测精度、支持土地复垦规划和增强矿区地表结构的保护具有重要意义。

2 研究区域和数据来源
2.1 研究区域
研究区域位于中国陕西省榆林市神东煤田的霍家屯矿第二盘区（图1）。该区域开发了多个煤层，包括1-2上煤层、1-2下煤层和2-2煤层。1-2上煤层和1-2下煤层已经完全开采，而2-2煤层的208工作面目前正在积极开采中，代表了该区域的第三次循环开采。2-2煤层在研究区域内横向连续且可完全开采。其地质结构从简单到中等复杂不等，有1-6层局部发育的夹层，厚度在0.02至0.89米之间。这些夹层主要由泥岩和粉砂岩组成。直接顶板主要由细粉砂岩和泥岩组成，局部夹有碳质泥岩，而底板主要由粉砂岩、细砂岩和泥岩组成，偶尔含有粘土岩。在工作面，底板海拔高度在1,085米至1,090米之间。工作面的走向长度和倾角长度分别约为2,700米和325米。平均开采厚度为4.2米，沿走向方向的埋藏深度约为140米。煤层以2°–3°的缓倾角倾斜，因此被归类为近水平煤层。208工作面的开采始于2023年4月，并于2023年11月30日完成，平均推进速度约为每天15米。该开采作业采用了全机械化的长壁后退式开采方法，并实施了全面的顶板管理。图1 研究区域。(a) 沙姆市的地理位置 (b) 研究区域的地理位置 (c) 208工作面GNSS点的详细布局。

2.2 数据来源
2.2.1 208工作面的GNSS数据采集
由于其高精度、全天候操作能力和短观测周期，GNSS技术非常适合用于地表条件复杂且变形迅速的矿区沉降监测。在本研究中，使用实时动态（RTK）GNSS定位技术获取了地表高程数据。通过在已知坐标的基站和多个流动站之间形成载波相位差分观测，实时校正了由卫星时钟偏差、电离层延迟和对流层延迟引起的误差，实现了厘米级的定位精度。为了确保整个监测期间的位置稳定性，使用了混凝土标志建立了永久性控制点。在沉降盆地内布置了两条地表观测线，分别标记为H线和C线。C线沿着煤层倾角方向布置，包含22个观测站。2023年3月31日至6月15日期间，在208工作面共进行了44次现场观测，观测间隔根据气象条件介于1至3天之间。

2.2.2 TLS点云数据的采集
鉴于研究区域地形起伏较大且地貌复杂，使用RIEGL VZ-1000扫描仪采集了地面激光扫描（TLS）数据。该仪器具有宽广的视野、高操作灵活性以及最大1,400米的测量范围。其360°全景扫描功能显著提高了在崎岖地形中的数据采集效率。RIEGL VZ-1000采用红外激光，并通过全波形数字化技术支持高速扫描，从而能够快速获取高密度、高精度的地表点云数据。由于其对环境干扰的强抵抗力，该系统在能见度低和恶劣的现场条件下仍能保持稳定的目标检测能力，特别适用于山区采矿环境。所采集的点云具有均匀的空间分布，并能有效穿透植被覆盖（如灌木和草地），有助于可靠地重建地表形态。此外，扫描仪配备了高分辨率数码相机，可以生成同时包含几何精度和表面纹理信息的点云数据。TLS测量在三个观测阶段进行，以捕捉不同采矿阶段的地表形态。数据采集分别于2023年3月19日、4月14日和6月15日进行。3月18日的测量代表了采矿前的地表状况，而后续测量则是在逐步开采煤炭后进行的，以记录采矿引起的地表沉降（图2）。

2.2.3 UAV LiDAR点云数据的采集
UAV LiDAR数据是使用配备D-LiDAR2000激光扫描模块的FEIMA D200无人机获取的。该集成系统结合了惯性导航单元、激光扫描仪和卫星定位传感器，能够同步获取点云数据以及精确的定位和姿态信息。飞行作业在相对高度100米处进行，相邻飞行线之间的横向重叠率为40%。208工作面的LiDAR测量覆盖面积约为3平方公里，产生了约1.9 × 10^8个点，平均点密度为每平方米64个点。两次UAV LiDAR测量分别于2023年3月18日和4月14日进行，累计飞行时间约为3小时。为了确保DEM的准确性，在数据采集过程中部署了实时动态检查点以进行动态质量控制。从UAV LiDAR数据中提取沉降盆地信息的整体工作流程如图3所示。

3 研究方法
3.1 基于镜像栅格重建（MRR）的栅格数据重建原理
为了重建由于LiDAR覆盖范围有限而产生的DEM间隙，在Python中实现了一种基于对称性的栅格镜像算法。该方法定义了一个穿过工作面的对称轴，并应用双线性插值和条件替换，将数据丰富的区域的地形和变形特征映射到被遮挡的区域。缺失（NoData）的像素仅由其镜像对应像素替换，从而保留了原始地形结构并最小化了重建过程中引入的人为不连续性。
根据经典的采矿沉降理论，在近水平煤层中单条长壁开采时，地表变形在走向方向和倾角方向上表现出明显的两侧对称性，这两个方向分别以采空区为中心，并垂直或平行于长壁推进方向。沉降盆地内的垂直沉降和水平位移都遵循相对于这些主轴的对称分布模式（Tang等人，2019年）。因此，沉降监测方案通常沿走向和倾角方向布置观测线，从而有效捕捉盆地的变形特征，同时降低现场成本。本研究采用的观测布局遵循这一原则。基于208工作面的矩形矢量表示，通过连接矩形两个短边的中点来确定倾斜主轴的对称轴。提取轴端点的坐标作为镜像过程的输入参数。然后将栅格值对称地映射到轴的对面，以恢复地表形态的连续性和完整性（图4）。

3.2 基于MRR的DEM数据重建工作流程
(a) 地表沉降盆地的对称示意图 (b) TLS数据缺失区域的示意图 (c) 基于MRR的缺失栅格数据重建示意图 (d) 通过与UAV生成的DEM数据比较进行精度验证
根据两点确定一条直线的数学原理，给定直线LL的端点(x1,y1)和(x2,y2)，其一般方程可以表示为：
Ax + By + C = 0 (1)
在方程中，系数A = (y2 - y1)，B = (x2 - x1)，C = (x1y2 - x2y1)，这两个端点决定了直线的方向。
在使用MRR方法计算所有栅格单元的“镜像”值之前，需要确定每个单元相对于对称轴的相对位置。这是通过将向量叉积应用于每个栅格单元来实现的，以确定该点位于对称轴的“左侧”还是“右侧”，从而确定哪一侧包含原始数据，哪一侧需要镜像重建。设直线LL上某点的方向向量为：
→v = (x2 - x1, y2 - y1) (2)
设从直线起点A(x1,y1)到点P(x,y)的向量为：
→w = (x - x1, y - y1) (3)
使用二维向量叉积确定点的相对位置的标准如下（方程4）：
cross = vx * wy - wy * vx
如果cross > 0，则点P位于向量→的“左侧”，表示原始数据；如果cross < 0，则点P位于“右侧”，表示需要补充数据的区域；如果cross = 0，则点P位于直线上。
由于对称轴需要“填充”的一侧也包含原始观测值，这些值不应被MRR方法重建，因此针对每个栅格单元(i,j)构建了以下镜像重建逻辑：
如果DEM[right(i,j)] = NoData，则DEM[left(i,j)] = DEM[left(i′,j′] (5)
如果栅格单元位于对称轴的右侧并且包含有效值，则保留其原始值；如果右侧单元是NoData点，则检查对称轴左侧对应的镜像位置(i′,j′)是否包含有效值，如果存在有效值，则将其赋给原始单元(i,j)，从而通过镜像映射过程重建缺失数据。
MRR方法的核心计算原理是通过向量投影来映射点位置。对于直线LL，其法向量定义为→n = (A, B)。为了计算方便，将该向量归一化为单位长度（方程6）：
→n = (A, B) / √(A^2 + B^2) (6)
点P(x,y)相对于直线的镜像点P′(x′,y′)可以表示为：
P′ = P - 2 * [(P · n) / √(A^2 + B^2)] (7)
由于镜像位置不一定与栅格单元中心重合，因此使用双线性插值从四个最近的有效邻居估计像素值。尽管基于对称性的填充有效恢复了缺失区域，但重建的DEM依赖于严格的镜像对应关系，并没有明确考虑沉降盆地的内在空间连续性。因此，在原始数据和镜像数据之间的交界处可能会出现人为不连续性。为了减轻这些影响并确保盆地形态的物理一致性，对过渡区域应用了带屏障的扩散插值（DIB）技术。该过程选择性地平滑了重建的接缝，同时保持边界约束，从而得到与实际地面变形模式更符合的连续且真实的沉降盆地表面。DIB通过求解稳态扩散方程（拉普拉斯方程）建立了一个连续值场。通过数值模拟空间扩散过程，该方法产生了在整个域内平滑且连续的插值结果，包括一阶导数。数学上，DIB可以表示为以下二维偏微分方程：
?2Z(x,y) = 0 (8)
其中u(x,y)表示目标变量（例如地面高程或沉降）的空间分布，?2表示拉普拉斯算子。观测样本点ui作为边界条件施加（方程9）：
?2Z = ?2Z/?x2 + ?2Z/?y2 (9)
插值计算使用观测样本值zi作为边界条件（方程10）：
Z(xi,y) = ∑_{i=1}^{n} wi * zi (10)
方程在网格上离散化，并使用有限差分或有限元方法求解，得到所有未知网格点上的平滑连续估计Z(x,y)（方程11）：
Z(x,y) = ∑_{i=1}^{n} wi * zi
权重wi(x,y)由扩散过程自动确定，满足∑_{i=1}^{n} wi = 1，并且随着相邻点之间距离的增加而逐渐减小。

3.3 概率积分法（PIM）
概率积分法是预测采矿引起的沉降最广泛使用和最成熟的分析方法之一。其理论基础最初由Knothe（1953年）提出，基于渐进式覆盖层移动的概念，后来由Whittaker和Reddish（1989年）、Peng（1992年）等人进行了形式化和工程应用。PIM的核心思想是将地下采矿对覆盖层-地表系统的影响建模为“传播-衰减”过程。任何点的地表沉降可以视为工作面内所有开采单元贡献的概率叠加，遵循正态分布。因此，由采矿引起的给定点(x,y)处的沉降可以表示为：
W(x,y) = 1/2 * ∫_{x=0}^{x} W_0 * e^(-u^2) * ∫_{y=0}^{y} W_0 * tan(β) * e^(-u^2) (13)
其中W_0(x)和W_0(y)分别表示沿走向和倾角方向的目标点的最大沉降值；mm是开采厚度；qq是沉降系数；H0是平均开采深度；H1和H2分别表示向下和向上开采的深度；tanβ表示主要影响角的正切值；θ是采矿影响传播角；S3和S4分别表示左右边界拐点的水平偏移距离；S1和S2分别表示向下和向上开采的拐点偏移距离。

3.3 沉降角的计算
地下煤炭开采会导致地表移动盆地，控制这些盆地的边界通常需要确定沉降角。这些参数主要用于描述采空区与沉降盆地之间的空间和时间关系，包括边界角、拉伸角和完全开采角。由于煤层的倾角约为2°–3°，属于近乎水平的煤层，因此在确定沉降角时不会区分上倾和下倾方向。本研究采用的方法的整体工作流程如下图（图5）所示。图5 沉降盆地重建方法的工作流程。需要注意的是，最初的无人机LiDAR和TLS测量分别于3月18日和3月19日进行，而煤层开采直到4月初才开始。因此，两次观测日期之间相差1天的时间不会影响实验结果。为了在后续分析中保持一致性，这两个数据集被视为在同一天获取的。

4 结果与分析
4.1 TLS和无人机数据的差异结果及准确性评估
4.1.1 TLS和无人机数据的差异结果
图6展示了截至4月14日，通过TLS和无人机LiDAR观测得到的208工作面的1米分辨率DEM。如图6a所示，TLS测量受到穿越南侧的高速公路的部分遮挡，导致沉降盆地覆盖不完整，并且在倾斜方向上缺少沉降信息。相比之下，无人机LiDAR观测（图6b）通过空中数据采集克服了这一限制，几乎完全覆盖了沉降盆地。由此产生的沉降模式呈现出近似轴对称的几何形状，为近水平煤层开采引起的沉降盆地的对称特性提供了实证支持。由于无人机LiDAR数据是在同一观测期间获取的，并且完全包含了TLS测量区域，因此它们被用作验证基于TLS的重建结果的参考数据。
图6 2023年3月18日至4月14日的沉降观测：(a) TLS DEM (b) 无人机LiDAR DEM。
4.1.2 TLS和无人机LiDAR数据的可靠性分析
原则上，TLS和无人机LiDAR是应用于同一研究区域和同一时间段的两种不同观测技术；因此，衍生DEM中相应像素的高程值应该是相同的。然而，实际上由于传感器精度、激光波长和覆盖范围、入射角、飞行姿态以及两个平台之间的点密度差异，不可避免地会出现差异。在应用本研究提出的重建算法之前，有必要使用独立的地面GNSS观测来评估这两个数据集的可靠性。为此，在沉降盆地内选择了十个GNSS观测点，并提取了这些位置处TLS和无人机衍生DEM的相应高程值。比较结果总结在表1中。
表1
点编号 GNSS(m) TLS(m) UAV(m) GNSS/TLS(m) GNSS/UAV(m)
H0 9 -0.443 -0.562 -0.516 0.118 0.07
H10 -0.993 -1.195 -1.170 0.201 0.176
H11 -1.579 -1.595 -1.729 0.015 0.149
H12 -2.32 -2.388 -2.541 0.068 0.221
H13 -2.718 -2.801 -2.944 0.083 0.226
H14 -2.675 -2.823 -2.648 0.147 -0.026
H15 -2.079 -2.228 -2.064 0.204 -0.014
H16 -1.392 -1.372 -1.398 -0.02 0.005
H17 -0.681 -0.782 -0.679 0.1 -0.002
H18 -0.325 -0.318 -0.306 -0.006 -0.018
RMSE 0.118 0.126
RRMSE 3.6% 3.9%
TLS-UAV-GNSS数据的准确性分析。
可以看出，TLS和无人机LiDAR在沉降盆地上的测量结果相对于GNSS观测的RMSE分别为0.118米和0.126米，分别占最大沉降量的3.6%和3.9%。这些结果表明，这两种数据类型都能可靠地捕捉地表变形，适合用于构建高精度的地表沉降模型。
4.2 TLS数据重建性能及准确性评估
4.2.1 TLS数据重建结果
4月14日的TLS沉降数据使用MRR方法进行处理，包括镜像映射、重采样和使用DIB进行平滑处理。重建后的TLS DEM如图7b所示，有效地填补了缺失区域并恢复了沉降盆地的连续性。
图7 MRR重建结果。(a) 4月14日的原始TLS DEM数据 (b) 使用MRR方法重建后的TLS DEM。
如图7b所示，基于MRR的重建产生了良好的几何结果，沉降盆地表现出高度的对称性。尽管如此，在原始区域和重建区域之间的交界处仍存在局部不连续性。为了更直观地评估重建性能，沿着最大沉降方向（图7中的黑色虚线所示）以1米间隔采样了高程值。下面展示了MRR重建的TLS数据与无人机LiDAR参考数据沿此轴向截面的对比：
如图8所示，从TLS重建得到的沉降曲线与无人机LiDAR剖面基本一致，最大沉降接近3.2米。该剖面可以近似认为关于点200对称，该点作为有效的对称轴。在轴的左侧，比较了测量的TLS和无人机LiDAR值，发现在拐点附近存在较大的差异。先前的研究表明，尽管TLS提供了高密度和高分辨率的点云，并且在快速和大幅度变形区域表现良好，但其灵敏度在沉降缓慢或幅度较小的区域会降低。在沉降梯度较大的区域，变形信号足够强，可以区分背景噪声；而在沉降较小的区域，低信噪比阻碍了微小变形信号的可靠识别（Soudarissanane等人，2011年）。在剖面的右侧，即重建区域，预测值与左侧的测量值表现出近似对称的行为，尽管仍存在局部偏差。总体而言，重建和参考剖面在沉降盆地的中心部分表现出很强的一致性，差异主要局限于孤立点，而不是系统性趋势。
图8 沿主要倾斜方向的剖面值提取。
4.2.2 错误分析
为了评估基于TLS的预测准确性，通过将重建的高程与相应位置的无人机LiDAR测量值进行比较来进行定量错误分析。第一步是定义验证的空间范围。由于重建的沉降盆地在视觉上很清晰，因此通过手动数字化一个包含重建区域的矢量边界来划定验证区域。在这个边界内，同时提取了TLS重建DEM和无人机LiDAR DEM中具有有效高程值的像素对，共得到1,515对像素对。为了表征重建错误的空间分布，进行了多种统计分析。使用Moran’s I来评估空间自相关性，而方差分析用于评估误差分散。Kolmogorov–Smirnov检验用于检查误差是否遵循均匀分布。最后，应用DBSCAN聚类来识别空间聚集的高误差区域。误差统计结果总结在表2中。
表2
分析方法 统计量 P值 结果
Moran’s I = 0.642 N/AS 强正相关
误差显著聚集
CVCV = 0.677 N/A 高变异性
空间分布极不均匀
KS检验 KS = 0.54 <0.0001 拒绝H0 不遵循均匀分布
DBSCAN 4个簇 N/A 4个明显簇
高误差系统性聚集
误差分析。
如表2所示，Moran’s I和变异系数都表明误差分布具有强正空间自相关性，表明重建误差是空间聚集的而不是随机分布的。Kolmogorov–Smirnov检验在5%的显著性水平上拒绝了零假设；在这种情况下，p < 0.0001强烈反对均匀误差分布的假设。这些结果证实了重建误差表现出明显的空间异质性。误差的空间模式及DBSCAN聚类得到的相应分类结果如下所示。
如图9a,b所示，误差主要集中在数据集的上缘和下缘，而中心区域的误差相对较小。图9c中的聚类结果显示，簇1–4（c1–c4）的平均误差范围为0.295至0.47米，其中c4的平均误差最大（0.47米），c2的最小（0.295米）。通过将四个高误差簇与高分辨率卫星图像进行比较，可以看出整个采样区域位于人工林种植区内。在TLS数据采集过程中，激光回波经常被多层植被遮挡或反射。因此，接收到的信号通常代表来自多个散射体的能量叠加或衰减，而不是真实的地面反射，导致测量误差增加。就误差点的数量而言，簇c2–c4包含的误差点比c1多。就平均误差大小而言，c3和c4的误差显著大于c1和c2。从空间上看，簇c2–c4主要分布在沉降盆地的边缘，而c1相对远离盆地边界。如4.2.1节所讨论的，沉降盆地边缘的变形监测准确性通常会降低。因此，c2–c4中观察到的较大误差可以归因于植被干扰和边界相关变形特征的综合影响。尽管c2的平均误差略小于c1（仅0.008米），但c2中的误差点数量明显更多，表明累积误差的潜力更大。值得注意的是，用于图8中剖面分析的采样点位于c4簇内，这进一步解释了在沉降曲线拐点附近观察到的较大误差。鉴于4月14日盆地的最大沉降值约为3.2米，使用MRR方法得到的预测结果的RMSE约为0.23米，相对RMSE约为7.18%，与无人机观测结果相比。由于该值低于通常接受的10%阈值，因此所提出的MRR方法被认为适用于重建缺失的TLS沉降数据。
图9 (a) 误差空间分布图 (b) 误差空间分类图 (c) DBSCAN聚类分析图。
5 讨论
5.1 工程应用
基于2023年4月14日208工作面沉降盆地重建的准确性评估，MRR方法在解决由于不可避免的现场条件导致的TLS数据缺失方面显示出明显的潜力。因此，该方法进一步应用于重建2023年6月15日获取的TLS数据集中的缺失区域，并推导出相应沉降盆地的沉降角。随后将重建结果与独立从GNSS观测获得的沉降角进行了比较。需要注意的是，208工作面正在进行第三轮重复开采。如图10a所示，第一和第二轮开采工作面的布局与第三轮开采工作面（208工作面）不垂直，而是呈现出大约11°的系统性偏差角α。因此，沉降盆地的发育沿着这个偏差方向延伸。如果在不考虑这种系统性偏差的情况下应用MRR方法，将会在重建的盆地中引入一致的位置偏差，如图10b所示。
图10 偏差角示意图。(a) 208工作面与第一和第二轮开采工作面的相对位置 (b) 考虑系统性偏差角α的MRR。
在图10b中，MM表示盆地的对称轴。PP表示原始像素位置，而PmPm表示对称轴右侧PP的真实位置，考虑了系统性偏差角αα。p′p′表示如果忽略偏差角αα时MRR方法会重建像素的位置。如图所示，偏移角的存在导致MRR算法将值分配给p′p′而不是正确的位置PmPm，从而在重建的沉降盆地中引入了系统性偏差。为了纠正这种效应并尽可能准确地恢复真实的盆地几何形状，通过计算偏移量ΔPΔP来调整每个像素的坐标，并相应地重新分配相应的沉降值，如下所示（公式14）：
Pm=P′+(2?d?tan(α))?→v
Pm=P′+2·d·?tanα·v→
(14)
其中，PmPm表示点的校正坐标，dd是从PP到对称轴的垂直距离，αα是偏差角，→vv→是对称轴的方向向量，决定了平移方向。通过将这种角度偏差校正集成到MRR方法中，重建了2023年6月15日的208工作面TLS数据集。重建后的沉降盆地如图11所示。
图11 (a) 6月15日的原始TLS DEM数据 (b) 使用MRR方法重建后的TLS DEM。
如图11所示，MRR方法有效地重建了缺失的TLS数据段。研究结果表明，该方法在保持对称性方面表现出色，相较于原始数据集具有更优的填充性能，并能够近乎完整地重建工作面沉降盆地。尽管采用了DIB技术来平滑数据“接缝”处的过渡，但仍存在一些微小的不连续性，这反映了该方法固有的局限性。这些差异的产生是因为原始沉降盆地在各个方向上并非完全对称。虽然对称性假设为数据预测提供了基础，但实际的盆地几何形状受到地质条件变化的影响——例如层理差异、节理和含水层的变化、煤层厚度的微小变化或倾角、不均匀的采矿进度以及上覆地层的异质性。三维可视化进一步便于全面检查，使得重建的沉降盆地可以在三维空间中查看（图12）。

图12显示了208工作面的沉降盆地，其轮廓大致呈“W”形，由两个明显的低洼区和三个高峰组成。第一个低洼区的最大沉降量为5,139毫米，第二个低洼区的最大沉降量为5,347毫米。这种波浪状的模式主要是由于该地区反复采矿的历史造成的。在2-2煤层上方的上覆地层中，12号上煤层和12号下煤层都已被开采过，为了适应连续的一级和二级开采阶段，保留了三个煤柱。当工作面形成采空区时，上覆地层失去支撑而下沉，而保留的煤柱则继续提供结构支撑，从而限制了局部沉降。采空区和煤柱的交替排列形成了“工作面-煤柱-工作面”的模式。因此，每个采空区上方的地层经历了明显的沉降，对应于地表沉降曲线中的低洼区；而煤柱上方的区域沉降较小，对应于高峰。采空区无支撑与煤柱支撑的交替作用共同形成了工作面上观察到的特征性“W”形地表沉降模式。

为了确定沿走向和倾角方向的沉降角度，分别以1米的间隔提取了走向（ab）和倾角（cd）剖面上的像素值。由于TLS观测的覆盖范围有限，以及沉降盆地的广阔范围，在走向或倾角接近-10毫米的位置缺乏数据，导致无法直接识别这些点，从而无法计算边界角度。为了解决这个问题，使用PIM对观测数据进行了拟合以估计沉降盆地的边缘点。理论上，在没有煤柱约束的情况下，走向和倾角主剖面应能完全反映采矿活动的影响（He, 1991）。确定标准如下（公式15）：
{Lstrike > 1.2H0 ～ 1.4H0, Ldip > 1.2H0 ～ 1.4H0}
其中H0表示平均采矿深度，在本例中为140米，Lstrike和Ldip分别表示采空区的走向和倾角长度。截至6月15日，走向方向的采矿长度为900米，倾角方向的采矿长度为325米。计算表明1.4H0等于196米，而Lstrike和Ldip均超过了这个值，表明地表已达到超充分采矿影响的状态。因此，使用PIM在有限采矿条件下拟合走向和倾角主剖面，从而确定了-10毫米沉降点。图13和图14分别展示了沿走向和倾角主剖面的二维沉降剖面，突出了超过最大沉降值的区域以及使用PIM获得的拟合结果。

图13显示，走向剖面呈现出明显的“W”形。PIM在沉降梯度较大的区域拟合效果有限，主要是因为PIM是为单煤层采矿条件下的地表沉降盆地设计的。尽管工作面内存在煤柱，但复杂的地形无法仅用一个沉降低洼区来充分表示。然而，沿盆地侧翼的拟合仍然令人满意。显然，在没有煤柱影响的情况下，走向主剖面将符合PIM模型的假设。因此，PIM提供了近似的拟合，并能可靠地识别-10毫米沉降边界。

图14表明，PIM在倾角主剖面上几乎实现了完美拟合，尤其是在盆地侧翼。在盆地中心观察到轻微的抬升，这主要是由于MRR预测在“接缝”过渡处引入的不连续性所致。尽管本研究涉及第三次重复采矿周期，但倾角主剖面仍然有效地反映了单煤层采矿情景，其中所有三个开采周期都产生了整体沉降，且没有保留煤柱的影响。与单煤层采矿相比，沉降系数可能略有增加；尽管如此，整体沉降模式仍与PIM假设一致，允许使用该方法进行可靠的盆地建模。

在确定了走向和倾角方向上对应-10毫米沉降的点之后，可以计算TLS预测的沉降角度。由于GNSS数据由间隔5-30米的离散点组成，且不与最大沉降的走向或倾角主剖面相交，仅从GNSS观测得到的沉降角度可能与实际移动沉降盆地的边界有显著偏差。因此，使用TLS预测的数据来优化沉降角度。总结的沉降角度如下：

表3显示，与GNSS数据相比，使用MRR方法预测的一些参数显示出显著差异。沿H剖面，GNSS和TLS-MRR数据得出的边界角度和完全开采角度大致相似，而抽采角度相差近7°。沿C剖面，虽然完全开采角度差异不大，但边界角度和抽采角度分别相差约10°和5°（表3）。因此，通过所提出的方法获得的沉降角度可以有效校正GNSS测量得到的值，为准确描述地表沉降模式提供了有价值的指导。

尽管所提出的方法表现出有效性，但仍存在一些局限性。MRR算法的核心步骤依赖于局部对称性的假设。实际上，沉降受到断层、关键地层和异质扰动等因素的影响，可能导致复杂的非对称变形模式。因此，该方法最适合地质条件相对均匀、地表坡度平缓、煤层水平或接近水平的地质环境。如果在原始观测数据可用的对称轴一侧存在少量缺失数据，可以先构建一个符合沉降盆地力学的物理约束趋势面（例如PIM）来填补这些空白。然后可以使用完成的数据集作为基于MRR的重建的输入。然而，如果在原始数据一侧发生大规模数据丢失，MRR方法的应用将不可避免地将这种缺陷传播到另一侧，导致大面积数据缺失。在这种情况下，MRR方法的基本前提被违反，该方法失去了实际意义。

对于非矩形或不规则的工作面，或者相邻工作面之间存在相对位置偏移的情况，仅从工作面矢量几何形状得出的对称轴可能与地表沉降盆地的实际对称轴不一致。但是，只要满足局部对称性的基本假设（条件（1）），就可以使用基于SAR的监测来识别大规模地表变形特征（例如沉降边界）。然后可以利用这些变形模式来判断沉降盆地是否具有近似对称性，之后可以通过视觉解释或适当的形态学算法估计对称轴。

该方法在单层、水平或接近水平的煤层采矿条件下表现最佳。在本研究考虑的第三次重复采矿周期中，尽管RRMSE结果符合标准精度要求，但仍可能出现潜在误差。在理想化条件下——仅考虑垂直位移、忽略水平移动、煤层倾角较小以及地质结构复杂的情况下，保留煤柱上方的地表沉降模式通常呈现“W”形。在本研究中，第一次和第二次采矿周期中保留煤柱的位置和宽度大致相同。因此，没有煤柱的区域沉降加剧，而有煤柱的区域则表现出更大的峰谷差异，导致第二次采矿周期的沉降较为明显。第三次采矿周期在没有保留煤柱的情况下，导致了整体沉降的增加。理论上，三次采矿周期的对称轴不应有显著差异。然而，实际情况更为复杂，包括其他未提及的因素，这些因素在应用单一对称轴时可能会引入偏差。因此，在重复采矿情景中，当数据覆盖有限且不需要超高精度时，该方法最适合快速近似地表沉降盆地。

为了解决由于客观现场条件导致的部分TLS数据丢失问题，本研究提出了一种基于采矿引起的沉降盆地对称性原理的镜像栅格重建（MRR）方法。该方法结合了镜像映射、区域匹配和边界约束优化，无需依赖传统的插值假设即可重建缺失区域，实现了形态学上一致的预测。它可以有效地通过预测由单层或多层接近水平煤层开采引起的地表变形导致的缺失TLS观测值来重建三维沉降盆地。

主要发现如下：（1）使用MRR重建的TLS数据与同时期的UAV测量结果相比，RMSE为0.23米，RRMSE为7.18%，证明了MRR对缺失像素的预测价值。（2）误差的空间分布和聚类分析显示，四个高误差簇主要集中在重建区域的边缘，主要是由于森林地形的干扰。对于第4个簇，最大的误差来自森林冠层遮挡和靠近沉降盆地边界的综合效应。（3）将MRR应用于6月15日的TLS数据集，并引入系统偏差角度后，沉降角度的准确性得到了提高：沿H线，抽采角度相对于GNSS观测结果提高了约7°；沿C线，边界角度和抽采角度分别提高了约10°和5°。（4）该方法为涉及单层或接近水平煤层采矿的沉降盆地重建提供了一种快速且经济有效的方法。尽管如此，该方法主要依赖于单一局部对称轴的假设，在复杂现场条件下（如数据缺失严重、地质环境异质或重复采矿操作可能导致对称轴偏移），这一假设可能被违反。关于重复采矿情景中的潜在轴偏移，未来的工作可以集中在数据丰富的条件下分别分析各个子盆地，以进一步完善模型并提高预测精度。