分区优化混合SfM方法

许 彪，董友强，张 力，孙钰珊，刘玉轩，查 冰，韩晓霞

1. 中国测绘科学研究院，北京 100830； 2. 北京建筑大学测绘与城市空间信息学院，北京 100044； 3. 北京市建筑遗产精细重构与健康监测重点实验室，北京 100044； 4. 俄亥俄州立大学工程学院土木环境和大地测量系，美国 哥伦布 43210

运动恢复结构(structure from motion，SfM)是利用一组具有重叠的无序影像恢复相机姿态的同时获得场景三维结构信息的处理过程[1]。SfM对参与排列影像的初始位姿、拍摄方式、相机参数等均没有严格的限制条件，具有良好的通用性。典型的SfM方法包括两大类：全局式SfM[2-8](global SfM,GSfM)和增量式SfM[1,9-18](incremental SfM，ISfM)。

GSfM是将所有影像整体求解，先旋转再平移，显著减少了优化的次数，但当影像姿态、尺度或者畸变较大时会出现解算不稳定，甚至解算失败[19]。不同于GSfM,ISfM以初始影像对为基础，渐进式地增加影像并扩充场景结构，稳健性较高[19-22]，当前大规模影像的SfM解决方案也多采用ISfM框架，比如COLMAP[9]、VisualSfM[22]等。但ISfM也有其缺点，主要包括：①对初始像对的依赖性高，选取不同的初始像对可能会得到差异较大的重建结果；②大场景处理效率低，增量过程中每次添加单张影像，为了消除累积误差，每次添加后都要执行光束法平差，消耗了大量的时间；③难以兼容分布式并行处理。ISfM算法属于串行模式，后续处理依赖前步结果，难以利用分布式并行计算节点提高效率。随着影像获取手段、效率的不断提高，场景影像数据量不断增加，以城市级三维建模的实际应用来说，几万张甚至几十万张影像的超大场景已变得很普遍，传统的ISfM方法越来越难以满足大规模影像三维重建的应用需求。

文献[23—24]尝试兼顾两种算法各自的优点，采用全局与增量的混合模式，提出基于分区与融合的混合式SfM方法，并将其应用于商业软件3DFZephyr。首先，将输入影像划分为多个重建单元，每个重建单元内仅包含部分影像数据；然后，利用传统的ISfM方法对每个分区单元进行重建；最后，将不同分区单位的重建结果进行融合得到完整的三维模型。通过将大影像集划分为多个计算单元，无论是分区重建还是融合过程，均可将其分配至多个计算节点实现分布式并行。这种做法不仅提高计算效率，而且避免了过多的参数的同时求解。但是，目前影像分区与融合方法研究还在起步阶段，一旦出现错误分区或不可靠分区重建结果，会严重影响最终重建精度，限制了当前混合式SfM方法在不同影像数据集上的应用。

为充分挖掘混合式SfM的优势，实现大规模无序影像的稳定快速高精度稀疏三维重建，本文提出分区优化的混合式SfM方法。具体内容包括：

(1) 提出一种基于影像关联度的分区方法，无须GPS/INS等其他辅助信息，仅利用影像间的匹配结果计算得到的关联度实现场景分区。

(2) 提出一种改进的快速增量式SfM方法(improved-ISfM)，每次重建多张影像，实现分区内影像的快速重建。在重建过程中，提出一种不可靠分区的自动检测方法，自动剔除不稳定分区，并将这些分区内的影像进行动态调整，重新划分入相邻分区。

(3) 提出一种稳健可靠的分区融合方法。利用分区间的公共地面点计算两两分区间的空间三维线性变换参数，利用变换参数将小分区合并成大分区，大分区合并成更大的分区，逐步获得完整的场景重建结果。同时，针对不同分区中可能包含同一相机拍摄影像的情况，采用了一种先融合场景后融合相机的策略，实现了相机参数、影像姿态和恢复三维信息的统一。多组具有挑战性的典型数据试验表明所提方法能够自适应于不同影像类型、不同场景和不同规模的影像数据，大大提高了重建效率的同时保证了较高的精度，尤其适用于包含几万甚至几十万张影像的大场景影像数据集。

1 基于分区优化的混合式SfM

借鉴摄影测量解析空三分区“接边”的思路，本文提出一种基于分区优化的混合式SfM方法。解析空三通常在两个“接边”区域量测公共的地面控制点并重叠少量影像，以提高接边处的加密精度，而本文方法即使在分区间不存在重叠影像时，也可完成重建过程。具体算法流程如图1所示，其中图1(a)为本文所提方法的整体技术流程，图1(b)为图1(a)中所使用的ISfM方法流程。不同于传统ISfM的是，本文的ISfM中分区动态调整的影像是从后方交会步骤开始。

图1 整体算法流程Fig.1 Overallflow chart of theproposed method

1.1 基于影像关联度的场景分区

场景分区将整个场景影像划分为多个存在关联的小块，分区过程无须飞行架次、GPS/INS等辅助信息，仅利用影像间的匹配结果。首先，利用SIFT[25]算法对影像对进行匹配，利用文献[26—27]所述方法找到待匹配影像对，并根据匹配点个数和点位分布计算两两影像间的关联度得分。其次，将场景中每张影像看作一个独立的分区，合并其邻域内关联度得分最高的同名影像，完成初次划分，此时每个分区内应包含至少两张影像。然后，取不同分区间影像关联度得分的最大值作为两个分区的关联度得分，依次合并关联度最大的两个分区，重复此过程，获得初始分区结果。最后，剔除弱连接分区，动态调整后完成整个场景分区。如图2所示，图2(a)为分区过程，场景影像经初次划分后，得到n个分区，根据关联度得分将这n个分区进行重新划分，合并关联最好的分区后得到包含m个分区的初始分区结果，经弱连接分区剔除后最终得到图2(b)所示的6个分区，图中每个方框圈定的区域代表一个分区，区域面积的差异描述为各分区包含不同的影像数，分区间无重叠影像。

图2 场景分区Fig.2 The scene partition

算法具体过程描述如下：

(1) 计算影像间关联度得分，分数的高低表征影像与其邻域的关联程度。设影像对同名像点数量hnum，按像幅划分格网(如10×10)，如果同名点落在某个格网内，则格网内数值累加1，所有同名点划分完毕后，累加值小于阈值的格网为无效格网。如图3所示，蓝色圆圈为有效格网，粉色圆圈为无效格网。统计有效格网内像点的外接矩形面积S(图3中红色框选区域)，S与像幅面积的比值即为scoredistribute。按照式(1)计算影像间的关联度得分corr

图3 点位分布得分Fig.3 Point distribution score diagram

(1)

式中，pairhleft、pairhright为左右影像构成的像对中最多同名点数；w1和w2为权重(本文选取0.7和0.3的经验值)。

(2) 影像初划分。统计每张影像关联度得分最高的一张同名影像，按分数从高到低排列，遍历排序结果，依次将未分区的影像与得分最高的同名影像划分到同一个分区。

(3) 遍历循环完成初始分区。分区间同名影像关联度得分最大值定义为两个分区的关联度得分，类似步骤(2)方法将小分区逐步合并为大分区，重复此过程，得到场景的初始分区结果。为了避免分区内影像过多或过少的情况出现，预设分区内影像数阈值范围(如15～75)，如小于阈值，将其划归至关联度得分最高的分区中，如大于阈值，回溯分区影像汇聚过程，将其拆分成两个或更多个满足阈值的分区。

(4) 剔除不可靠分区。统计两两分区间的公共地面点数量及分区与其邻域的点数最大值，排序所有分区的统计最大值，选取0.25倍的中值作为阈值，如某分区与其邻域的公共点数小于阈值，视为弱连接分区将其剔除，影像重新划分至其他分区中。至此，得到场景分区结果。

图4为一组包含372张低空摄影影像的数据。图4(a)为影像摄站点的平面位置图，粉色圈定的两张影像为ISfM选取的初始像对；图4(b)为场景分区结果。从结果可以看出，整个场景共划分为7个分区，不同分区利用不同颜色加以区分，每个分区内圈定的两张影像为选定的初始像对。

图4 实际数据的场景分区结果Fig.4 Partition result of actual data

1.2 独立分区执行改进ISfM

不同于传统的ISfM每次添加一张影像，本文通过每次添加多张影像，提出了Improved-ISfM算法，大大提高了每个分区内重建的效率。整体流程如图1(b)所示，具体如下：①选取两张影像作为初始像对，估计影像相对位姿和结构信息；②统计所有剩余影像后方交会的有效点数并排序，点数大于0.25倍中值的影像作为备选影像，然后统计备选影像进行前方交会得到的有效点数并排序，大于0.5倍中值的影像将增加至场景，完成新添加影像的位姿估计(image registration)；③前方交会得到更多的结构信息(triangulation)，增大场景的覆盖范围；④光束法优化使得重投影误差和最小；⑤过滤掉不可靠的外点(outlier filtering)，通常情况下，光束法优化和外点过滤均需循环执行多次，光束法优化后过滤外点，利用过滤结果再进行光束法优化，提高重建精度。循环②—⑤步直至无法增加更多影像至场景中为止。

考虑到影像匹配精度、地形起伏、飞行方式、相机检校特性及其他不确定性因素，为保证完整重建结果的正确性，需进一步检测相互独立的分区重建结果，剔除可能存在的不稳定项，具体包含以下几个方面：

(1) 统计误差过大的分区。选用点位平均误差MSE作为统计指标，即分区内有效点点位误差的平均值，取所有分区MSE中值的2.0倍作为阈值，删除MSE大于阈值的分区。

(2) 相机检校参数不可靠的分区。如果分区结果正确，不同分区中同源相机的检校参数(包括焦距、主点、畸变参数)应表现出大体一致的结果，剔除包含异常检校结果的分区。

(3) 成功排列影像占比低的分区。如果某分区内成功排列的影像数与总数比值小于阈值(如0.5)，认为重建结果不可靠，将其剔除。

(4) 弱关联分区。实际应用中，考虑到匹配粗差、相机检校、几何条件等因素的综合作用，不能保证重建后的每个分区间依然存在足够的公共地面点，参照2.1节中剔除弱连接分区策略，剔除可能存在的孤立分区。

将被剔除分区涉及的影像和其他未成功排列的影像追加至其他分区中，实现分区的动态调整，增加特定影像至指定的分区后，从后方交会开始执行ISfM流程完成影像排列。上述算法策略不仅降低了ISfM对初始像对的高依赖性，而且通过选取多个初始像对降低了重建失败的风险。与此同时，由于分区动态调整后额外重建所需的计算量不大，不会降低整体效率，避免了传统ISfM中因重选初始像对而必须放弃所有已完成重建结果的问题。

1.3 光束法融合场景分区

场景不同分区重建结果存在坐标系和比例尺不一致的问题，需通过融合获得统一坐标系和比例尺下的完整重建结果。本文利用重建后分区间的公共地面点，采用循环的方式，将小分区逐步融合成大分区，完成分区融合，具体过程如图1(a)右侧融合循环部分所示。

(1) 计算两两分区间的空间三维相似变换参数。通过分区间公共的地面点计算两两分区的三维线性变换7参数，统一待融合分区的坐标系和比例尺。

(2) 确定待融合的分区。任一分区可能邻接多个分区，且分区间的公共地面点数量通常不一致，需设计合理的策略，在保证算法稳健性的同时确保每个分区均参与合并。视每个分区与其邻域有效点数累加和为关联性指标，按从小到大顺序排列，简记为S1、S2、…、Sx，优先融合关联度指标小的分区。首先，选取S1作为基准，有效点数最多的Si与其融合，标记S1和Si为已处理；然后，对于尚未标记的Sj，设其有效点数最多的为Sk，若Sk已标记，则将Sj追加至Sk的融合列表中，否则Sj与Sk融合，并标记为已处理分区；最后，不断重复该过程，直至所有分区被标记，得到了融合任务列表。

(3) 确定融合分区影像位姿和结构信息坐标初值。步骤(1)中仅计算了两两分区间的变换参数，假如融合任务仅涉及两个分区，可任选其中一个为基准，将另一个变换至基准分区即可。当存在3个及以上分区需融合至一起时，选取关联性指标最大的分区为基准，如果其他分区与基准存在有效的变换参数时，直接将其变换至基准分区。当某分区与基准分区无法直接变换时，则根据空间相似变换的刚体特性及两两间的变换参数，以传递的方式通过多次变换将其转换至基准分区。

(4) 光束法优化。场景分区时，多张同相机的影像可能会被划分至不同的分区中，导致同一相机在不同分区中得到不同的检校参数。为了保证融合后参数的唯一性，遵循先融合场景再融合相机的原则。融合场景时每个分区的相机检校参数相互独立，光束法优化影像位姿和结构信息坐标。融合相机则是将同一相机的多套检校参数合并为一套，即选取或重新估计一套检校参数作为初值，再进行光束法优化，使得相机检校参数唯一。前述过程为一个融合循环，不断重复该过程，直至所有分区融合成完整场景为止。

需要说明的是，为了提高融合算法的稳健性同时保证分区融合结果的精度，本文采用了两种优化方案。①步骤(3)获取初值时，首先确定影像的位姿，然后通过多片前方交会计算结构信息的初值，用于提高初值精度。②步骤(4)光束法优化时，采用执行多次的策略，每执行一次光束法优化后，还需进行外点过滤及新插入2D、3D点的处理。一方面，降低外点对光束法优化结果的影响；另一方面，场景中新增的2D、3D点可进一步改善平差网的几何条件，有助于提高平差精度。

图5为图4所示数据的融合过程。7个分区各自完成重建处理后，利用分区间的公共地面点进行两次融合。图5(a)为影像的融合过程，第1次融合时将左侧图中的①④、②③⑥、⑤⑦融合成3个分区，再将右侧图中的①②③融合成完整的场景。图5(b)为结构的融合过程。为便于显示，已将每个分区的结果转换至比例尺一致且原点统一的坐标系下。

图5 实际数据描述的分区融合过程Fig.5 Fusion process of actual data

需要注意的是，因受多种因素的影响或干扰，融合可能会发生失败。为了保证算法的稳健性，出现融合失败时，按前文所述分区动态调整的策略进行处理。实际上本文的分区仅仅是在逻辑上指定了影像的“归属”，每张影像仍然可以灵活归入其他存在关联的分区中。

1.4 分布式并行

无论是分区ISfM还是融合分区，每一步均是相互独立的过程，因此可将每个处理过程分配至多个计算节点，实现分布式并行处理，以提高大场景的处理效率。基本的并行调度方案设计包括主机端和计算节点端。主机端负责数据组织、场景分区、确定融合任务列表、计算任务划分、分区动态调整与结果回收。计算节点端负责接受主机端指定的任务，包括单独分区重建和分区间融合两部分，同时将计算结果返还至主机端。

2 试验与结果分析

2.1 试验数据和环境

为了验证本文算法的有效性和实际性能，选取低空摄影、多视倾斜摄影和近景等共计5组试验数据。如图6所示，试验场景覆盖了环拍办公楼、狭长河道、城区和城区山地混合地形等多种地形，同时影像数量、连接点平均重叠度和地面分辨率均不相同。场景A为办公楼，利用SONY NEX-7相机环绕拍摄获取的204张近景影像，GSD为1.5 cm；场景B为长春某校园，利用大疆精灵3拍摄获取了1292张低空影像，GSD为4.8 cm；场景C为10 km长的汉北河分支河道，利用大疆精灵4拍摄13架次获取了4099张低空影像，GSD为5 cm，因存在拍摄间隔，按飞行架次定义了13组相机；场景D为城区和山地混合的四川省某市，利用SONY ILCE-5100z相机获取了29 754张5视倾斜影像，GSD为2 cm；场景E为建筑物密集的江苏省某市城区，利用SONY ILCE-5100z相机获取了121 506张5视倾斜影像，GSD为5 cm。

图6 本文试验数据影像截图Fig.6 Screenshot of test data sets

表1 试验数据概况Tab.1 Basic information of test data sets

所有试验均在如下配置的计算机上进行：操作系统为Windows10 64位专业版，处理器为Intel(R) Xeon(R) CPU E5 2.4 GHz 2处理器24核，内存为64 GB，显卡为Quadro K2200。

2.2 分区数量对算法效率的影响

为了评估本文方法中分区影像数大小对处理效率的影响，分别选取了[15,75]、[25,100]和[50,125]3种预设范围进行试验。统计结果如表2所示，每组试验场景3种分区大小的效率对比如图7所示。

图7 不同分区大小下处理效率对比Fig.7 Comparison of processing efficiency under different partition sizes

表2 不同分区大小下处理效率统计Tab.2 Comparison of processing efficiency under different partition sizes

(1) 本文方法重建精度与分区影像数预设范围无关。试验中，3种不同预设范围下，5组场景的重建精度一致，其中，场景C的1和2试验最大差异为0.01像素，场景分区大或小仅是将高关联度得分的同名影像划分多与少的问题，每个分区的内部几何条件均较为理想，并没有改变整个网形的几何条件，即便存在部分初始划分不合理的情况，分区动态调整机制的使用，也会将不合理部分修正，不会引起精度的明显差异。

(2) 分区影像数预设范围与处理效率无关。对于影像数较少的场景A、B、C来说，可认为处理效率一致。效率差异最大的为场景E的试验1和2，相差87 min。以试验1为参考，试验2效率提升5%，但是随着预设范围的增大，处理效率并没有得到进一步提升。其原因主要是因为分区与融合作为一个整体性的算法，不能简单地按两步拆分去看待，分区数的减少使得分区重建的处理时间增多。尽管融合循环及每次循环的融合次数减少了，但更多影像的融合使得光束法优化计算量也更大，还不足以明显提升效率。

2.3 对比试验分析

为验证本文方法的有效性，设计了3组对比试验。为了验证improved-ISfM方法的精度和效率，试验1对比了采用ISfM方法的著名开源软件COLMAP和本文提出的improved-ISfM方法；试验2对比了本文方法(包含整个算法流程)和improved-ISfM的精度和效率，由于本文算法内的分区重建内采用了improved-ISfM方法，因此结果的不同来源于所采用分区融合框架，该试验验证了分区融合策略的性能；试验3将本文方法与商业软件AgisoftPhotoscan(V 1.6.1)和3DF Zephyr Aerial(V 4.519)进行了对比，其中3DFZephyr Aerial采用了混合式SfM，可以与本文的混合式SfM方法进行更直接的对比。为方便起见，将以上两个软件分别简称为Photoscan和3DF。

2.3.1 试验1：COLMAP与Improved-ISfM对比

表3给出了对比试验结果。从结果可以看到，improved-ISfM精度明显优于COLMAP。在重建连接点更多的情况下，improved-ISfM对场景A、B、C的处理效率分别是COLMAP方法的6.93、4.42和7.93倍。效率提升的主要原因是improved-ISfM在每次增量过程中都加入了更多的影像和结构信息，进而改善了用于光束法优化的网形条件。这种做法在减少迭代次数加快收敛的同时，降低了外点过滤和光束法优化的循环次数。图8描述了improved-ISfM中连接点平均重叠度与每次增量增加影像数的关系。可以看到，算法效率与连接点平均重叠度和摄影方式两个因素有关。连接点平均重叠度越大，表征影像间的关联性越强，每次增量处理能够插入更多的影像，处理效率也会更高。

表3 COLMAP与Improved-ISfM重建对比Tab.3 Reconstruction comparison between COLMAP and Improved-ISfM

图8 Improved-ISfM平均每次增加影像数与连接点平均重叠度的变化关系Fig.8 The correlation between the average number of images addedeachtime and average overlap of connection points of improved-ISfM

2.3.2 试验2：本文算法与improved-ISfM对比

为了验证分区融合策略的性能，分别应用improved-ISfM方法以及本文方法分别处理整个数据集。表4给出了对比试验结果，从中可以看到：

表4 本文方法与improved-ISfM重建对比Tab.4 Reconstruction comparison between proposed methodandimproved-ISfM

(1) 本文方法能够适用于各种类型的场景。试验数据包括近景、低空摄影、倾斜摄影3种类型，包括山区、城区、河道、环拍建筑物、植被覆盖等场景，影像数从204张到12万张，尤其是试验场景C，狭长形状且飞行架次间存在弱连接区域，所提算法均能得到可靠的重建结果，说明了算法具有良好的稳健性，适用于各种类型的场景数据。

(2) 本文方法与improved-ISfM精度基本保持一致。最大的精度差别仅有0.05像素(场景A)，这种细微的差别主要是由初值差异、坐标系定义、部分临界误差像点、计算舍与误差、RANSAC算法等因素造成的。

(3) 对于包含影像数较少的场景，本文方法与improved-ISfM效率差异不大。如图9所示，对于场景A和B，两种方法的单位处理时间对比分别为1∶1.25和1∶1.31，相差较小。其主要原因是，ISfM增量处理过程中，每次新增影像和结构信息后，光束法优化所需计算开销有限，因此对整体效率的影响不显著。

图9 不同场景影像数下本文方法与improved-ISfM处理效率对比Fig.9 Comparison of processing efficiency under different number of scene images

(4) 对于大数量影像集，本文方法具有明显的效率优势。随着场景影像数的不断增加，本文方法效率显著提升。相较于improved-ISfM方法，本文方法在场景C、D、E上效率分别提升1.79、2.85和5.05倍。可以预见的是，当场景影像规模进一步增大时，本文方法的效率优势会更加明显，这是由于增加的影像仅会增加分区和融合的时间开销，不会导致处理时间的显著增长。

2.3.3 试验3：本文方法与商业软件对比试验

分别利用Photoscan和3DF处理前述试验场景，其中场景A、B和C，利用3.1节所述配置计算机进行处理。由于Photoscan和3DF在此配置下无法进行重建，故利用内存扩大至128 GB且其他配置相同的计算机进行处理。此时，Photoscan依然无法获得场景E的重建结果，而3DF无法获得场景D和E的重建结果。

试验过程中，Photoscan参数设置采用该软件推荐的默认设置，具体描述为：accuracy(high)、generic preselection、key point limit(40 000)及tie point limit(4000)。3DF利用高级设置，具体参数为：关键点密度(高)、匹配类型(精确)、匹配阶段深度(高)、照片排序(未排序)。表5给出了对比分析试验结果，每个指标最好的结果进行了加粗显示。由表5可以看出：

表5 不同方法对比试验结果Tab.5 Comparative experimental resultsof different methods

(1) 本文方法兼具精度和稳健性优势。仅有本文方法成功处理了所有试验场景，而在所有方法都成功重建的影像集上，本文方法精度与Photoscan基本一致，优于3DF软件。

(2) 本文方法效率优于Photoscan。从总体处理时间来看，本文方法在场景B和D上消耗的时间稍长；但相对于Photoscan，本文方法在场景A、B、C、D上单位时间内可处理的连接点个数分别提高2.59、2.44、1.12和1.34倍，而单位时间内可处理的像点个数分别提高1.55、2.45、1.09和1.33倍，重建效率更高。

(3) 从效率、连接点数量和精度几方面来看，本文方法都明显优于同样采用分区融合策略的3DF软件。尤其是在效率方面，本文方法在场景A、B、C上效率相对于3DF分别提高了2.09、3.26和4.58倍，随着影像数的增多，效率提升更为明显。一方面是因为improved-ISfM方法提高了分区重建的效率；另一方面，本文方法充分利用了分区间的关联性指标，而且每次融合任务中能够包含更多的分区，大大减少了融合循环次数，因此降低了时间消耗。另外，本文先融合“场景”再“相机”的策略，提高了光束法优化的收敛速度。

图10为本文分区融合算法对试验场景的重建结果，显示软件为FugroViewer(V 3.3)。

图10(a)、(e)的上部，图10(b)、(c)、(d)的左侧为影像摄站点的平面位置；图10(a)、(e)的下部，图10(b)、(c)、(d)的右部为结构信息。

图10 试验场景重建结果，包括影像摄站点平面位置和结构信息Fig.10 The reconstruction results including plane positions of cameras and reconstructed 3D points

3 结 语

本文提出了一种基于分区优化的混合式SfM方法，以应对大规模无序影像稀疏三维重建。处理过程中，无须GPS/INS等其他辅助信息，影像数量越大，本文方法的效率优势越大。多组不同规模、不同场景的无序影像试验验证了本文方法的有效性和稳健性，其中包括一组包含12万张5视倾斜影像的数据，在这组数据上仅有本文方法重建成功，商业软件AgisoftPhotoscan和同样运用了混合式SfM的商业软件3DF Zephyr Aerial均重建失败。

本文方法为大规模无序影像稀疏三维重建提供了可靠、实用的解决方案，可以促进相关应用的发展。后续工作主要包括两方面：一是改进融合算法，当场景中分区结果比较稳定时，将多个分区同时融合减少计算冗余；二是将所提方法分布式并行化，充分利用多计算节点进一步提高重建效率。