摘 要：针对空中三角测量的处理时间随图像数量的增加呈指数级增长，处理大规模数据集需要大量时间的问题，该文提出了一种测区智能分区与合并算法。该算法可以自动地将无序的图像集分割成若干个有重叠的子集，将每个子集可以并行处理，根据重叠图像的连接点，将各个子集的重建结果合并在一起。实验结果表明，该算法不仅保证了时效性，同时流程也简单，处理速度很快。在一个有10处理节点的集群系统中，该算法成功地处理了大规模航空影像数据集，重建的时间效率和精度满足实际生产要求。与不分区进行比较，时间可以节省至少一半。
0 引言
随着航空影像技术的发展，人们获取图像的方式越来越多，运动恢复结构（structure from motion，SFM）问题[1]的规模也越来越多。增量式SFM方法是目前比较流行的稀疏重建方法[2-4]，它从两张或三张图像开始解算，然后持续将其他图像添加到当前已解算模型中，并对结果进行平差优化。该过程是重复进行，直到求解完所有的图像。这种方法因为每次迭代都会进行参数优化，所以比较稳健。然而，该方法优化是通过光束法平差实现，计算量很高。随着图像数量增多，会出现两个问题：一是内存空间不够，大型任务不能在一台计算机上执行；二是处理时间也呈指数式增长。为解决这些问题，本文提出一种测区智能分区与合并算法。整个场景的无序影像会自动分割为若干个有重叠的子集。然后在集群处理系统中对每个子集进行并行处理。最后通过区块合并算法将每个子集的结果合并为一个完整的重建结果。
通过查阅相关文献发现，很少有文献讨论大规模SFM划分方法。文献[5-6]是与本文主题比较相近的论文。文献[6]的作者还发表了一些其他密切相关的论文，如文献[7-8]。文献[5-6]的工作流程与本文算法类似，但区块划分算法完全不同。前者假设拍摄时间接近的照片在空间上也很接近，因此他们根据图像获取的顺序对场景进行细分。这种假设在实践中并不总正确。后者使用归一化割(Normalized Cut) [9-11]方法细分场景。但在这种情况下使用Normalized Cut有两个问题：第一个问题是Normalized Cut分区没有重叠；第二个问题是每个分区的大小不一致，需要谨慎控制。因此，他们必须采用两步法,图形分割（graph division）和图形扩展（graph expansion）来对场景进行划分：graph division步骤从整个图开始，迭代地使用Normalized Cut，将任何不满足尺寸约束的子图分割成两个平衡的子图，直到没有子图违反尺寸约束为止。然后，如果子图的重叠比率小于阈值，则进行graph expansion步骤，将图中丢弃的边和相关顶点随机地添加到其中一个连通的子图中。graph expansion之后可能会违反尺寸约束，因此需要对这两种操作进行迭代，直到两个约束都满足为止。相比之下，本文的算法更简单有效，可以很自然地处理重叠和分区大小的问题，而仅需要使用矩阵带宽缩减算法对图像进行重新排序。
本文区块合并算法使用了类似增量式SFM方法的工作流程，运行更为稳健。通过选择好的起始子块和次优子块，在每个合并步骤后进行全局光束法平差，因此可以保证整个合并过程稳健进行。如果将每个子块中的增量SFM称为细粒度SFM，则可以将块合并过程称为粗粒度SFM。实验结果表明，本文算法具有良好的鲁棒性，系统的鲁棒性甚至比ContextCapture[12]和Agisoft PhotoScan[13]的SFM模块要好，这两个软件是比较优秀的三维重建技术的代表。本文系统的稀疏重建结果与ContextCapture和PhotoScan的比较将在后面的章节中给出。
1 总体工作流程
本文稀疏重建的总体工作流程如图1所示，整个场景的图像集合是输入数据。第一个处理步骤是特征提取，可以分布式执行。特征提取的输出又是每个图像的特征。这些特征进一步用于图像检索，找到可能有重叠的图像。这些图像对是需要通过图像匹配算法进行匹配的候选匹配对。在匹配任务时也可以分布式执行。通过图像匹配，可以得到图像的邻接信息。当且仅当一个图像有足够的同名点时，该图像被视为与另一个图像相邻。图像邻接矩阵是区块分割算法的输入。区块分割算法将整个块划分为若干个较小的子块。它决定哪些图像应该分配给一个子集以及两个子集之间的重叠图像。然后用增量式SFM算法对这些子块进行分布式处理。在得到每个子块的稀疏重建结果后，采用区块合并算法从这些子块重建结果中得到完整的重建结果。
2 区块分割算法
理想的区块分割算法具备以下特点：不依赖辅助信息，如全球定位系统/惯性测量单元（GPS/IMU）信息和成像时间信息等，这些信息有助于分块，但依赖这些信息会降低算法的通用性。为避免过多的信息依赖，本研究将原始图像当作本算法工作流程的唯一输入。在区块分割算法执行中必须满足要求：具有较好连接的图像应划分在一起，以便在后续的步骤中进行鲁棒重建，而没有连接的图像不应分配到同一个子集中。而对于无序图像，使用特征匹配工作流程来获取图像之间的连接关系。
通过特征匹配，可以得到图像的邻接矩阵。假设共有n个图像，则图像邻接矩阵A是一个对称的n×n矩阵。如果图像i和图像j有连接，则矩阵的元素A[i，j]的值为1，否则该值为0。如果两个图像的同名点超过500个，则视为有连接。在实际应用中，一幅图像通常只连接到一定数量的图像，因此图像的邻接矩阵往往是稀疏矩阵，如图2（a）所示。最初图像是按字母顺序排列，这个初始排序列表用L0表示。在这种排序下，具有连接的图像不一定彼此非常接近。从邻接矩阵中亮像素（值为1的元素）的分布可以看出这一点。邻接矩阵中亮像素的分布受图像排序的影响。
本文区块分割算法的基本思想是重新排序图像，使具有强连接的图像在排序列表中彼此非常接近。然后就可以根据排序列表来分块。当具有连接的图像在排序列表中彼此非常接近时，亮像素的分布将沿着邻接矩阵的主对角线变得非常紧凑，如图2（b）所示。通过矩阵的行和列的重新排列，也即通过调整图像的顺序，可以实现这一目的，这与矩阵带宽缩减的目的非常相似。所以可以使用矩阵带宽缩减算法[14]来重新排序图像。图3显示了图像重新排序的过程示例。首先，图3（a）以无向图的形式给出了7个图像的连接关系，相应的邻接矩阵如图3（b）所示。矩阵的行数和列数等于图中的图像数。最初，图像的顺序是1、2、3、4、5、6和7。重新排序后，图像的顺序是5、3、2、4、1、6和7。最初，邻接矩阵的带宽是6。重新排序后，带宽减少到4，并且在新的排序列表中，具有连接的图像彼此更加接近。
图3 图像重排序的例子
本文设计的区块分割算法流程：初始图像排序列表用L0表示，对应的初始图像邻接矩阵用A0表示。子集的图像数量用s表示，子集之间重叠图像的数量用o表示，参数s和o可以由用户指定。然后算法执行以下步骤。
1） 利用矩阵带宽缩减算法将矩阵A0置换成块对角矩阵。矩阵置换后的新图像排序列表用L1表示。在新的排序中，连接较好的图像彼此间距离更近，而连接不好或没有连接的图像彼此间的距离会更远。
（2） 根据用户的设置，将L1的前s张图像分配给子集1。然后，使用矩阵带宽缩减算法对列表L1的[s—o，end]范围内的图像的邻接矩阵A1重新排序，以获得这些图像的新排序列表L2，这些图像包含没有被分配给任何子集的图像以及o个重叠图像。
（3） 将L2的前s张图像分配给子集2。然后，利用矩阵带宽缩减算法，对在列表L2的[s—o，end]范围内的图像邻接矩阵A2再次进行重新排序，得到新的排序列表L3。
4） 重复前面过程，直到剩下的图像不超过s个，这些剩余的图像被分配到最后一个子集。
综上所述，从初始的图像邻接矩阵开始，本文算法重复使用矩阵带宽缩减算法对图像进行重新排序，并每次将前s个图像分配给一个新的子集。基于这个特点，本文算法被命名为迭代带宽缩减分区算法。该算法的特点是：①除图像本身以外，不需要其他信息；②用户可以指定子集的大小和子集之间重叠图像的大小；③可以自动将连接较好的图像组合在一起；④区块分割算法本身很简单而且可以节省时间，可以在1~2 s内完成100 000个图像的自动分区。
3.区块合并算法
当区块分割后，每个子集的稀疏重建可以在集群上并行执行，因此,每个子区重建的坐标系是独立的。为了得到一个完整的重建结果，必须进行区块合并，重新整理连接点的观测，并对子区之间公共的参数进行统一优化。
本文中，区块合并算法的基础是图像id和连接点id。在整个区块中，每个图像和每个特征点都有一个唯一的id。基于这些id，可以找到不同子块之间的连接。研究中的区块合并算法的工作流程类似增量式SFM的工作流程（图4）。图4是两个子区进行合并的例子，根据图像id，可以找到它们之间公共的图像，然后对相机参数和图像参数进行平均。根据特征点id，可以识别出公共的连接点。对于子集1中的连接点p1和子集2中的连接点p2，可以将p1和p2合并为一个连接点p，因为p1和p2在image3和image4上共享两个相同的观测。合并之后，连接点p有4个观测，分别位于图像1、图像3、图像4和图像2上。整理好p点后，删除p1和p2两点。
在本文中，区块合并算法的工作流程如下。
1）找到最佳的两个子块。如果两个子块共享的连接点数量最多，则视为最佳。
2）进行空间相似变换。在步骤1）中找到的子块之一被视为参考块，另一个子块被视为目标块。然后将目标块的坐标系转换为参考块的坐标系。接着采用三维空间相似变换作为变换模型，根据子块间共享点的坐标计算子块之间的变换参数。
3）参考块和转换后目标块数据合并。该部分包括相机参数的合并、图像参数的合并和连接点的合并（如图4所示）。相机参数合并是将两个块中公共的相机参数进行平均，相机参数包括焦距、主点坐标和畸变参数。图像参数合并是将两个块中公共的图像参数进行平均。其中，图像参数包括位置参数和姿态参数。连接点合并不仅对点的坐标进行平均，还需要对点的观测进行重新整理。数据融合后，由两个子区的重建结果生成一个新的大些的重建结果。
4）对新生成的重建结果进行全局光束法平差。本文中，步骤3）得到的参数值是初值，为避免退化，需采用全局光束法平差对其进行优化。
5）在光束法平差步骤之后，新生成的重建结果被视为参考块。然后找到下一批最佳的k个子块。如果子块与参考块有最多的公共连接点，则该子块被视为最佳。找到的k个子块被视为目标块。k取值越大，需要进行的光束法平差次数越少，合区速度越快，但可能造成空间相似变换传递误差的增大，影响合区的效果。目前本文代码实现中的k取值为10，该值仅供读者参考，其最优性还没有进一步测试。
6） 然后对参考块和k个目标块进行空间相似变换、数据合并和全局光束法平差。之后，这些子块合并成为新的参考块。
7） 重复步骤5）和6），直至不再存有目标块。最后，可以得到一个完整的重建结果。
本文区块合并算法进行全局光束法平差的次数非常少，近似为（n/s）/k。n是图像总数，s是每个子集的大小。那么n/s大约是子块的数目。k是一次添加到参考块中的目标块数。例如，若n=30 000，s=500，k=10，则执行全局光束法平差的次数大约是6次。因此，可以发现区块合并算法所需时间较短。
4 实验结果
4.1区块分割结果
为便于可视化，将不同子集的图像用不同的颜色显示，子集之间的公共图像以这些子集其中的一种颜色显示。为了比较，本文也给出了使用Normalized Cut方法进行分割的结果。使用的Normalized Cut代码来自文献[15-16]。从实验结果可以发现，本研究的分块算法分割区块效果更优。
第一个数据集是一个只有38张图像组成的小数据集，如图5（a）所示。这些图像是环绕一个石头物体进行拍摄。在实际作业中，类似这么小的数据集不需要进行分区（这里使用它只是为了展示本文算法的特点）。该区块被划分为3个子集。第二个数据集是航空图像数据集，包含2 282张图像。[图5（b）],展示了将区块划分为6个子集的结果。第三个数据集也是航空图像数据集，包含4 427张图像[图5（c）]，图中展示了该区块被划分为10个子集的结果。从第一个数据集可以发现，Normalized Cut方法更倾向于在图像连接较弱的地方对图像进行分区。从第二个数据集可以发现，Normalized Cut方法有时会将图像分配到错误的分区中。从第三个数据集可以发现，Normalized Cut方法的有些分区形状非常不紧凑。造成这些现象的原因是由Normalized Cut方法的目标函数引起的。而本文的划分算法试图得到最小的矩阵带宽而不是最小割。从这些案例可以看出，本研究的分块算法可以将具有强连接的相邻图像分配到同一个子集中，并且子集的大小和形状都是均匀分布的。
图 5 3个数据集的区块分割结果
4.2时间效率比较
本小节比较了3种稀疏重建处理策略时间效率：第一种是增量式SFM策略，不使用区块分割算法，所有任务使用一台计算机进行处理；第二种策略是增量式SFM策略，使用区块分割算法，但要求所有这些子任务都在一台计算机上按顺序进行处理；第三种策略是增量式SFM策略，使用区块分割算法，并且在一个有10台工作节点的计算机集群处理系统上以并行方式处理这些子任务。
实验结果见表1。由于时间较长，最后8个数据集没有被第一种策略处理，因此相应的表单元格用“X”标记，表示数据集信息不可用。同样，最后4个数据集也没有被第二种策略处理。通过对比可以发现，第二种策略比第一种策略具有更高的时间效率，因为它将指数时间成本转化为近似线性成本。第三种策略比第二种策略更节省时间，因为它将工作流从顺序处理转换为并行处理。第三种策略的时间开销还受子任务数和工作节点数的影响。可以看到，区块分割对时间效率的提升非常明显，处理时间可节省一半以上。tpd这个数据只有347张影像，在3种方式中都只分了一个任务，所以3种方式处理时间相同。
表1 3种稀疏重建处理策略的时间效率比较
4.3稀疏重建稳健性比较
除了时间效率，系统稳健性也很重要。本文使用另外两个软件ContextCapture和PhotoScan进行比较，对表1中列出的20个图像数据集开展测试。所使用的ContextCapture版本是v4.3.0.506，PhotoScan的版本是1.4.1build5925。ContextCapture使用增量式SFM工作流程，PhotoScan使用层级的SFM (hierarchical SFM)工作流程。测试中发现，ContextCapture和PhotoScan经常出现稀疏重建效果不理想的现象，如图6所示。
这个数据是“njgt”数据集。它包含3 804张图片，是无人机从两个不同的飞行高度拍摄的图像。3个稀疏重建结果的比较如图6所示。在ContextCapture的稀疏重建结果中可以看到一些图像漂移，并且整个结果是弯曲的。PhotoScan的稀疏重建效果也不好，出现一些错误的相机位置，且有很多错误的点。综合对比，本文的算法结果更优。
(a) 本文算法的重建结果
(b) ContextCapture的重建结果
(c) PhotoScan的重建结果
图6 “njgt”数据集稀疏重建结果的比较
4.4 稀疏重建的精度
除了时间效率和系统稳健性，还需考虑重建的精度。为了检验本文分布式空中三角测量结果的精度，本文使用一部分地面控制点对自由网结果进行有控平差。利用其他地面控制点作为检查点，检查定向精度。
本文用gg-20k这个数据作为例子来说明定向精度情况。这个数据总共有21 957张航空影像，这些影像是由一架装有5台索尼相机的无人机拍摄的。该数据集的地面采样距离为3 cm，拍摄场景为城区。空三连接点的数量为4 393 284个，连接点的平均观测数为22.04。用13个地面控制点进行有控平差，用另外13个地面控制点作为检查点。这些检查点的精度如表2所示，满足1:500测图要求。
表2 gg-20k数据13个检查点的精度情况
5 结束语
针对大规模空中三角测量存在的问题，本文提出了一种区块分割和区块合并的智能算法，以提高大规模空中三角测量时间效率。区块分割算法主要是对图像进行重新排序，使具有强连接的图像在排序列表中彼此非常接近，然后根据排序列表来分块。实验表明，本算法分割结果良好，可以得到紧凑的子集，同时该算法具有更高的时间效率。由于本文分块算法的输入是图像的邻接矩阵，错误的图像匹配结果产生的错误图像连接关系会对算法产生负面影响。在后续研究中，将考虑一些更稳健的方法对图像匹配结果进行滤波，以获得良好的算法输入。
原标题：大规模空中三角测量的智能分区与合并算法
作者：骆奇峰1, 丁华祥1 ,鲁路平2
(1. 广东省国土资源测绘院，广州 510500；
2. 武汉大学 地球空间信息技术协同创新中心，武汉 430079）
作者简介：骆奇峰（1969—），男，高级工程师，主要研究方向为水准测量、GPS 测量、数字化地形图测量及其他基础测绘。
- END -
