改进SfM算法的多视卫星影像三维重建研究

摘要：针对现有多视卫星影像三维重建方法在精度和效率方面的不足，提出了一种改进的增量式三维重建方法，有效提升了三维重建结果的精度和完整性。通过引入透视成像模型局部近似有理函数模型的思想，在重建区域内采用透视相机模型代替有理函数模型，以简化复杂的函数关系求解，解决了传统SfM算法不能直接处理线阵卫星影像的问题。在稀疏重建阶段，设计并引入伽马校正进行色调映射，并结合SIFTGPU与RANSAC算法实现了快速特征匹配，显著提高了特征提取与匹配的速度和准确率。在稠密重建阶段，提出了基于参考面优化和视差约束的深度图重建方法，克服了传统方法在数值精度方面的局限。在MVS3DM公共数据集上的实验结果表明：本文方法在大多数情况下的重建精度和完整性均优于其他同类方法，显示出更好的重建效果。研究结果验证了所提出方法在提升多视高分辨率卫星影像三维重建性能方面的有效性，为进一步提高卫星影像三维重建的精度和效率提供了新的思路。

引用：[1] 欧远征,叶江,尹海平,等. 改进SfM算法的多视卫星影像三维重建研究[J]. 测绘科学, 2024, 49 (11): 74-87. 

引言

随着遥感技术的迅速发展，高分辨率遥感卫星已经成为全球范围内不可或缺的工具。高分辨率的遥感数据不仅大大丰富了我们对地球表面的观测，对现实世界进行真实、实时描述和管理成为可能，也为数字城市、智慧城市等概念的实现提供了重要支持，其中重建高精度城市三维模型变得至关重要。利用卫星影像进行多视影像三维重建按照采用的重建策略的不同可分为两类：基于立体匹配重建和基于运动恢复结构(structure from motion, SfM)重建 [1] 。基于立体匹配的方法是双目立体的直接扩展，通过独立处理多对影像，将多个双视图重建结果进行融合 [2-3] 。而第二类是一种计算机视觉方法，同时解决所有图像的三角测量问题，以完成三维重建 [4-6] 。立体匹配在成对多视图立体中，通过找到密集的对应关系并利用传感器相机模型对它们的位置进行三角测量来单独处理成对的影像，而运动恢复结构算法直接重建多个图像的信息。

目前，基于有理函数模型(rational function model, RFM)模型的重建算法通常先完成双视角重建，然后基于DSM融合算法完成多幅视图的重建工作。其重建流程主要包括：采用视差 [7-8] 或光流 [9] 获得图像之间的稠密对应关系；通过光束法平差优化光学卫星图像有理多项式系数(rational polynomial coefficient, RPC)模型的估计残差；基于RPC模型进行三角化恢复图像对应的三维场景结构。文献 [10]利用捆绑调整优化成对光学遥感图像附带偏差补偿的RPC模型，以满足模型之间的最大一致性，并使用RPC模型三角化恢复对应地面场景的三维结构。ASP依赖从RPC模型中的大量参数从卫星图像计算多个双视图的3D点云，最后融合点云生成DSM [11] 。S2P利用两视图之间的指向误差校正替代了复杂的非线性捆绑调整，但需要基于RPC模型进行复杂的高程迭代，以实现成对光学遥感图像的三维重建 [12] 。在此基础上，文献 [13]分析了影响两视图重建精度的主要因素，并依据这些因素对数据集中的图像对进行启发式排序，最后选择性地融合了S2P生成的独立的两视图重建结果。MicMac采用多分辨率匹配方法来限制视差搜索范围，通过生成粗分辨率视差图，作为高分辨率匹配阶段的先验信息，以提高计算速度并减少噪声 [8] 。需要指出的是，尽管上述基于RPC模型的卫星图像三维重建方法已经展现出较高的重建精度，然而这些重建方法都难以避免求解关于RPC模型的复杂三阶多项式，且估计RPC模型需要至少39组地面控制点及相关图像的对应关系，在重建过程需要耗费大量的运算时间 [14] 。

计算机视觉领域提出了许多基于光学图像的三维场景重建方法 [15-20] 。然而，这些方法通常专注于处理面阵相机拍摄的图像，无法直接处理卫星上搭载的推扫式线阵相机所获取的图像。换而言之，计算机视觉领域的三维重建方法并不适用于卫星图像的三维重建任务。为了克服这一限制，研究人员不断探索如何将视觉领域的重建算法应用于遥感图像。文献 [21]提出了一种利用地面控制点和图像点之间的真实对应关系直接拟合仿射相机模型的方法，通过这种方法可以实现光学遥感图像的仿射重建；然而，该方法对整幅光学遥感图像进行操作，因此重建精度较低。文献 [1]利用已知的RPC模型拟合弱透视投影模型，基于COLMAP实现了计算视觉重建管道的光学遥感影像三维重建。文献 [14]提出了一种基于仿射成像模型的增量式光学卫星图像重建方法，该方法不再依赖于RPC模型，只需预先获得至少4个地面控制点即可实现场景重建。

综上所述，尽管近年来文献中提出了许多基于光学卫星图像的三维重建方法，使得关于光学卫星图像的三维重建技术取得了快速的发展，然而这些方法需要反复地执行捆绑调整优化以减少近似模型误差，或者需要RPC模型中的大量参数参与重建计算，往往耗费大量的时间和CPU计算资源。而计算机视觉领域的SfM算法，能够在GPU硬件并行加速下，从不同条件下拍摄的照片中快速完成三维重建，为大量级的多视卫星影像的三维重建提供了一种新的思路。

因此，本文在传统增量式SfM算法的基础上，提出了一种改进的增量式多视高分辨率卫星影像三维重建方法。该方法在数据预处理阶段，引入透视投影模型局部范围内代替有理函数模型，避免了复杂的卫星成像模型求解，同时解决了SfM算法无法直接处理线阵卫星影像的问题。在稀疏重建阶段，设计引入伽马校正对影像进行色调映射，并结合SIFTGPU算法与RANSAC算法的特征匹配，提高了影像特征提取速度和匹配准确率。在稠密重建阶段，本文提出了基于参考面优化和视差约束的深度图重建方法，解决了传统重建方法的重建结果数值精度不高的问题，实现了多视高分辨率卫星影像的高精度三维重建。

1 研究方法

1.1 重建流程

摄影测量与计算机视觉领域的一个主要区别是相机模型的不同，卫星遥感领域通常使用RPC模型，而计算机视觉则采用更为简单的针孔相机模型。本文采用文献 [1]的方法，通过拟合弱透视原理将RPC模型局部简化为针孔相机模型，在此基础上，本文提出了一种改进的增量式SfM卫星影像的城市三维重建方法，步骤见图1。

该方法首先对多视高分辨率卫星图像进行裁剪和伽马校正的预处理，并将卫星成像RPC模型转换透视相机模型，以准备进行下一步的稀疏重建；其次在稀疏重建的特征匹配阶段使用SIFTGPU算法实现CPU和GPU混合并行加速特征提取，并引入随机采样一致(random sample consensus, RANSAC)算法对匹配结果进行筛选，以消除误匹配，进而提升匹配的精确性和效率。随后，在局部优化相机位姿和三维空间点后，通过全局光束平差法进行全局优化，并对优化后的内在矩阵进行分解以消除影像的倾斜误差。最后，在稠密重建阶段重新参数化深度以生成深度图，并对深度图进行中值滤波优化进而生成DSM,最后融合所有DSM,实现对多张高分辨率卫星图像的三维重建。

1.2 数据预处理模块

由于卫星影像图幅较大，处理整幅影像将会耗费大量时间，因此选取特定的区域进行实验，在(纬度、经度、海拔)坐标系中，感兴趣区域由2D边界框指定。为了沿着z轴限制兴趣区域，本文使用了30 m精度的SRTM数据作为兴趣区域的高程参考，将输入的卫星图像被裁剪到该感兴趣区内。将裁剪后的图像通过伽马校正将具有高动态范围(high dynamic range, HDR)的卫星影像映射为低动态范围(low dynamic range, LDR)。

伽马校正通过调整伽马值和比例因子来实现，其中伽马值用于控制图像的非线性变化，而比例因子则用于调整图像的整体亮度，而不影响像素之间的相对差异。

1.3 基于GPU优化的SIFT算法

由于传统的SIFT算法速度较为缓慢，为了加速计算，在特征提取阶段使用SIFTGPU算法加速特征提取与匹配 [22] 。该方法利用了GPU并行处理的特点，以加速SIFT算法的多个步骤，具体的实现过程如图2所示。

在特征点匹配过程中，通常会出现一些错误匹配。为了有效地识别和去除这些错误匹配，引入了RANSAC算法。通过随机选择一组数据点来拟合一个变换模型，然后利用该模型对其余的匹配点进行测试。如果一个点符合该模型，则将其视为内点，否则将其视为外点。如果内点的数量超过了预设的阈值，那么该模型被认为是良好拟合的，并且基于这些内点重新拟合一个更加精确的模型。如果内点的数量未达到阈值，那么就需要重新进行随机采样，并且重复上述的步骤，直到满足模型设定的条件或达到了最大迭代次数，本文实验过程中的最大迭代次数为2 000次。

1.4 光束法平差

为了改善稀疏场景结构以及精确求解相机矩阵，本文对稀疏三维点云和初始的相机矩阵进行光束法平差(bundle adjustment, BA) [23] 。光束法平差是一种细化一组相机参数的方法，通过最小化一组稀疏3D点的重投影误差来提高其全局一致性和准确性。为了最大限度地减少投影模糊引起的失真，向稀疏三维点的坐标中添加正则化项，保证不会偏离其精确坐标太远。

1.5 倾斜误差校正

在增量式三维重建方法中，内在矩阵中的偏斜参数通常被假定为零，因此现有的稠密重建方法通常不会对其进行建模。本文通过分解内在矩阵来实现倾斜校正。

通过分解内在矩阵，应用矩阵 A 的逆对图像进行变换，并使用矩阵 B 作为新的无偏斜内在函数，对生成的图像扭曲以原始分辨率对图像重新采样，来解决偏斜参数为零的问题。

1.6 基于参考面优化和视差约束的深度计算

在稠密重建阶段，需要计算每个视图的深度值进而重建深度图。然而，由于卫星相机与地球之间的距离较远，导致深度的平均值较大而方差相对较小。当使用普通GPU进行稠密重建时，可能会造成数值精度不高。为了解决这个问题，通过优化参考平面，并结合视差约束，对深度值进行重新计算 [25] 。

2 实验数据处理与分析

2.1 实验数据与评价指标

2.1.1 实验数据

为了评估本文方法的有效性，使用国际公用的MVS3DM数据集中的50张WorldView-3全色(PAN)影像和机载激光雷达数据，MVS3DM基准数据集中的多日期卫星图像Digital Globe提供，包括2014—2016年，时间跨度为14个月，包含在阿根廷圣费尔南多附近约100 km2区域，全色影像地面采样距离约为30 cm; 实况机载激光雷达(LIGHT LASER DETECTION AND RANGING,LiDAR)数据的地面采样距离约为20 cm, 可作为后续实验的验证数据。

由于数据集图像采集时间的跨度较长，因此数据集中的图像均可划分为多组图像，同时对其进行处理，并从数据集中挑选了部分复杂场景区域，其中包括建筑密集、存在大量阴影、大型建筑以及地面对象复杂等四种复杂场景图像，为了更加公平地评估算法的重建效率和重建性能，本文分别使用S2P、COLMAP方法和本文方法进行三维重建，并对重建结果进行对比分析。

2.1.2 评价指标

三维重建质量分析与评价是保证三维重建质量是否达到应用要求的前提与基础。本文采用了综合的精度评价指标，首先，将基准数据集中提供的LiDAR点云与重建生成的点云映射到相同的地理网格中，生成地面高程图；然后，对这两个高程图进行配准对齐；最后，统计配准后的两个高程图之间的像素值差异。评价三维重建质量主要从目视检查和构建精度评价指标两个方面进行。目视检查方法通过检查三维点云或三维点云格网化后的DSM,定性分析DSM表达三维重建的准确程度，并统计DSM中的各种误差。定量评测过程采用重建精度、完整度以及均方根误差(表1)作为评价指标。其中重建精度表示中值高度误差(median error, ME),完整度(completeness, CP)表示误差小于1 m的重建点云与整体点云的比值。

2.2 实验1

本节对数据集中的时间跨度较短的一组图像进行重建对比，图像分组方法按照影像成像时间接近和影像重叠度较高方法，本组图像选取了4张成像时间接近的遥感影像(增量式重建方法至少需要3张),选取了MVS3DM数据集中的8个场景(图3)进行实验。针对上述不同重建方法的重建结果，选取Challenge场景的重建结果(图4)进行展示，重建结果高程和重建误差对应的彩色条状图如图5所示，可推广至其他场景。

对影像时间跨度较短的一组图像进行三维重建，结果如表2所示，在各个测试场景下，S2P重建方法在各个场景的重建结果中孔洞更少，而且细节保留较好，在大多数重建场景的CP、ME、以及RMSE均优于其他重建方法，重建结果表现最好。COLMAP方法在各项指标上不及S2P方法和本文方法重建结果。本文方法在绝大多数场景上重建结果精度优于COLMAP方法，但仅在MasterSeq.1场景的完整度优于其他方法，而其他场景精度不及S2P方法重建效果，但本文方法的总体精度较其他方法相比更接近S2P方法的重建精度。本文方法属于增量式三维重建方法，在少量图像进行的重建优势还并未体现出来。

2.3 实验2

城市环境中场景复杂多变，MVS3DM数据集中的8个场景点不足以代表城市复杂场景，为进一步验证本文方法的三维重建精度，在原数据集中8个场景的基础上，共选取了20个复杂场景区域，分为4类：建筑物密集区域(建筑物的数量相对较多，密度较大)、存在阴影区域(高大楼体导致阴影)、地物复杂区域(树木、水体和建筑等多种地物)以及存在大型建筑区域(厂房、机场等大型建筑)4种类型复杂场景的卫星图像(图6)进行三维重建。

针对图6所示的不同复杂城市场景，选取了其中的一个存在阴影场景并对重建细节局部放大进行展示(图7),重建结果高程和重建误差对应的彩色条状图如图8所示。从三维重建的整体效果来看，S2P、COLMAP和本文方法的重建结果在建筑物轮廓方面，本文方法重建结果更接近真实情况，且孔洞最少，而COLMAP重建结果存在一定程度的轮廓失真。

针对4类不同复杂城市场景三维重建结果(图9),重建结果高程和重建误差对应的彩色条状图如图10所示，每类场景选取1个重建结果进行展示：在地物密集区域的第一组场景中，本文方法在重建结果中保留了微小地物细节，未出现明显损失。在存在阴影区域的第二组场景中，阴影和遮挡导致误匹配的概率增大，但本文方法成功地恢复了更多的地物细节，而S2P方法和COLMAP方法重建结果则存在更多的孔洞。在地物复杂的第三组场景中，包含大量树木和水面，增加了误匹配的风险，进而影响了三维重建的精度。而在第四组存在大型建筑场景中，存在大型建筑和茂密的树木，在本文方法重建结果成功保留了大型建筑的轮廓和细节。

图９ 部分场景三维重建结果及误差

在各个测试场景下，所有场景三维重建结果统计(表3)表明，本文方法三维重建结果在CP、ME、和RMSE的3项指标上较其他方法表现出一定的优越性，本文方法在大多数情况下取得了最低的ME与最高的CP。在建筑物密集场景类型下，建筑物密集1和建筑物密集2场景的ME、CP和RMSE指标分别为：43.6%、0.987 m、2.305 m和23.3%、1.391 m、2.345 m; 在存在阴影场景类型下，存在阴影1场景的ME、CP和RMSE指标为：66.9%、0.401 m、2.942 m; 在地物复杂场景类型下，地物复杂3场景的ME、CP和RMSE指标为：51.6%、0.783 m、2.749 m; 存在大型建筑场景类型下，大型建筑1、大型建筑3和大型建筑5这3个场景的ME、CP和RMSE指标分别为：67.5%、0.505 m、2.671 m; 67.1%、0.548 m、1.956 m和73.6%、0.847 m、1.416 m。本文所提方法在各个复杂场景中的上述7个场景的重建精度均优于其他方法。相比其他重建方法，本文方法在对大规模卫星影像进行三维重建更具有优势，重建结果精度更高。

3 结束语

针对城市区域不同复杂场景，使用S2P、COLMAP和改进的增量式多视高分辨率卫星影像三维重建方法，在MVS3DM数据集上进行三维重建实验。实验结果表明，本文方法重建结果的可视化效果优于其他方法，重建结果孔洞更少，地物边界更接近真实情况，尤其在大量级图像重建更具有优势。其中最优重建结果的CP、ME和RMSE三项指标分别为73.6%,0.847 m、1.416 m; 完整度较S2P和COLMAP方法分别提高了21.1%和36.0%;中值高度误差分别提高了9.5%和11.4%;均方根误差分别提高了9.3%和8.6%。本文方法在大多数场景重建结果中取得了最优的中值高度误差和最高的完整度，在均方根误差指标表现出与其他方法类似或更好的结果，表现出良好的三维重建性能，验证了本文方法的有效性。

本文方法与S2P相比，本文方法直接线性求解三维点云坐标，避免了复杂的高程迭代；此外，在融合步骤中只需对每个图像块生成的高程图进行融合，减少计算资源消耗的同时提升了效率。与COLMAP相比，本文方法在特征提取阶段使用RANSAC算法剔除误匹配，并在稠密重建中使用参考面优化和视差约束的深度图重建方法重新估计深度图，提高了三维重建精度。综合而言，本文方法更适用于处理大规模城市的长时间跨度卫星影像。