无人机(unmanned aerial vehicle，UAV)具备高时效性、高灵活性、操作简便性等特点，已成为摄影测量与遥感领域的重要观测平台[1]。结合优视摄影或贴近摄影等几何感知路径规划方法[2-3]，无人机倾斜摄影测量系统可获取厘米级甚至毫米级空间分辨率的影像数据，保证其广泛应用[4-7]。无人机平台灵活度更大、机载POS(positioning and orientation system)系统的精度更低，导致影像定位定向难度增加。如何实现高效、高精度的无人机影像定位定向是拓展无人机平台应用范围、提升摄影测量与遥感产品可靠性的重要环节。

当前，运动恢复结构(structure from motion，SfM)已成为无人机影像定位定向的关键技术[8]，主要分为增量式SfM、全局式SfM和混合式SfM[9]。其中，增量式SfM采用迭代局部或全局光束法平差(bundle adjustment，BA)进行参数解算，具备很好的抗噪声能力，其定位定向精度更高，被广泛用于商业和开源软件系统。然而，针对大数据量、大重叠度和高分辨率的无序无人机影像，增量式SfM顺序性约束所导致的BA迭代优化误差累积和效率下降，使其很难满足大场景无人机影像的定位定向处理需求。

现有增量式SfM效率提升方法主要关注影像匹配和几何解算。对于影像匹配，有效匹配对检索是现有研究的关注点。基于数据采集先验知识或无人机平台POS数据，有序影像(航线排列规则)匹配对检索可采用时间顺序约束或空间邻近约束实现，具备非常高的检索效率[10-12]。但是，受限于特定场景或辅助数据的精度，无法适应复杂场景无序影像(航线排列不规则)，如优视摄影测量技术[2]。基于内容的影像检索(content-based image retrieval，CBIR)，特别是词袋模型(bag-of-words，BoW)，是现有SfM软件匹配对检索的标准模块，其核心思想是利用编码本构建影像特征的BoW向量，将匹配对检索转化为BoW向量的最近邻搜索问题[13-14]。然而，随着影像分辨率提升和数据量增加，词袋模型检索的效率会显著下降，主要原因是影像特征数量和编码本尺寸增加导致构建影像索引的时间消耗增加。因此，为了提升无序影像匹配效率，有效匹配对检索依然有待提升。

增量式SfM几何解算效率提升的研究包括BA优化和大场景处理[9]。对于BA优化，可利用硬件加速提升BA解算效率，比如CPU/GPU混合加速技术[15]，或挖掘相机参数内在约束或相机与三维点的稀疏投影关系，实现BA模型简化和算法优化[16]。这类方法不能从本质上解决大场景增量式SfM的弊端。近年来，基于分块-合并策略的并行化SfM技术得到了越来越多关注，其核心思想是将大场景分割为BA优化易解决、几何连接关系强的小场景，可有效降低BA迭代优化的误差累积和效率下降问题[17-19]。文献[20]提出分区优化混合式SfM，设计不稳定分区检测和可靠分区合并算法，实现大场景无人机影像增量式SfM处理。文献[21]提出顾及场景连通性的混合式SfM，其核心是引入场景关联度因子实现场景分割及分区连通性实现连接边扩展。相关研究很好地促进了大场景增量式SfM的适用性。但是，子场景稳健融合依然是分块-合并策略所需要解决的难点，主要体现在场景合并过程中：①随着子场景逐渐合并，参考模型的尺寸越来越大。如何快速、低内存消耗地实现公共三维点查找是第一个因素。②由于遮挡、重复或弱纹理影响，子模型的三维点不可避免地包含噪点，且基于公共三维点配准的子场景几何变换估计缺乏严格的几何意义。因此，如何最大限度降低噪点的影响，实现稳健子场景合并是第二个因素。

针对上述问题，本文提出一种联合全局描述子和图索引的无人机影像并行化SfM方法。针对影像特征数量大、编码本尺寸增加所导致的匹配检索效率低的问题，设计了联合VLAD(vector of locally aggregated descriptors)全局描述子和HNSW(hierarchical navigable small world)图索引的高效影像检索方法，显著降低了影像局部特征数量、提高了高维向量最近邻搜索效率，进而加速影像匹配。针对子场景合并中同名点搜索效率低、内存消耗大、几何解算精度低的问题，设计了基于按需匹配图和双向重投影误差的子场景合并方法，避免了参考模型增长导致的公共三维点查找效率低、内存消耗高、几何变化估计精度低的问题。利用多组无人机影像进行试验，验证了本文方法的可靠性。

本文方法的总体技术流程如图1所示。输入数据仅为无序无人机影像I＝{i1，i2，…，in}。首先，通过影像子集和特征子集选择，构建影像检索编码本。其次，利用编码本将局部特征聚合为VLAD全局向量[22]，并基于HNSW图索引[23]构建检索库，实现影像匹配对高效、自适应检索。然后，利用检索影像对引导特征匹配，进而构建加权拓扑连接图表示的无人机影像场景图，实现场景图分块和并行化SfM三维重建。最后，基于按需匹配图和双向重投影误差实现子场景合，构建了满足大场景无人机影像的并行化SfM技术流程。

图1   无人机影像并行化SfM技术流程

影像匹配对检索的目标是查找空间重叠的影像对，引导和加速特征匹配。本文提出了联合全局描述子和图索引的高效影像检索方法，其核心思想是一方面利用VLAD全局描述子减少局部特征点数量，加速影像索引构建；另一方面利用HNSW图索引实现高维全局描述子索引和最近邻搜索，最终加速影像匹配对选择。同时，考虑到现有方法常采用固定数量或比率确定检索影像数，存在匹配对数量过多或不足的问题，进一步提出自适应阈值匹配对选择算法。具体实现如下：

(1) 编码本构建。VLAD向量聚合依赖编码本。考虑到无人机影像高重叠度和高分辨率特性导致编码本构建计算代价高的问题，本文采用影像子集和特征子集选择策略生成训练特征，即对于影像的高冗余性，随机选择比例p的影像作为训练数据；对于特征点的高冗余性，将SIFT特征按照尺度从大到小降序排列，并选择尺度最大的h个特征点。其基本原理是影像的大尺度特征点能够覆盖到整个尺度空间，实现对影像的特征描述[24]。在获取训练特征后，采用k均值聚类[25]算法构建词汇树，利用K-means＋＋算法确定k个初始聚类中心；然后根据最近距离准则，将特征描述子划分为k个子集，并更新初始聚类中心；最后依据上述步骤，对每个子集进行迭代的聚类操作，直到达到终止条件，得到所需的k个单词数量的编码本C′＝{c1，c2，…，ck}。

(2) VLAD向量聚合。BoW统计单词在影像中的出现频次，构建k维特征向量。随着待索引影像数量增加，需要依赖更大尺寸的编码本C，导致编码本构建和影像索引的时间迅速增加。VLAD累加影像特征与最近邻单词的向量残差，组成k×d维的特征向量，可将影像局部特征聚合为高维全局向量。其中，d为特征描述子的维度。因此，VLAD可显著降低对编码本的尺寸要求。对于k个单词组成的编码本C′，VLAD的向量聚合原理：遍历距离单词ck最近的影像局部特征{di}，i＝1，2，…，N，计算两者特征描述子残差并累加求和，得到对应单词ck的VLAD向量组成部分vk＝(vk，1，vk，2，…，vk，d)。各个元素vk，j计算公式为

(1)

式中，j是特征描述子维度索引，j＝1，2，…，d；ak(di)是指示函数，当特征di包含在单词ck内时，ak(di)＝1，否则，ak(di)＝0。基于上述编码，k个单词组成的编码本C可将影像表示为k×d维的VLAD特征向量。

(3) HNSW索引与检索。针对高维向量检索的高时间代价，本文利用HNSW图索引构建VLAD高维向量索引结构，实现高效的影像匹配对检索[23]。HNSW的基本原理是将待索引向量构造成图结构G＝{V，E}。其中，V和E分别表示待索引向量组成的顶点集合及其连接关系组成的边集合。为了实现高效索引和检索，HNSW要求构图时所有节点需存在邻居(连接顶点)，且邻居数要尽可能少。同时，距离相近的顶点须互为邻居。其中，最大邻居数量M是图索引结构构建时间成本和检索精度的重要因子。图2展示了HNSW图构造与图索引原理。图构造时，规定M＝3，6个待索引向量通过5次构造实现图索引，分别插入顶点A、B、F、C、E和D。当图中顶点数大于M时，通过图检索(图2(b))查找最近邻的M个顶点，直接连接待插入顶点与M个最近邻顶点。图检索时，随机选择一个顶点作为进入点，查找其距离查询向量最近的邻居，并以此邻居作为当前进入点；迭代上述步骤，直到满足终止条件，得到查询向量的检索结果。

图2   HNSW图索引与图检索原理

(4) 自适应检索。基于HNSW构建的索引结构，现有方法常采用固定数量或比率确定检索影像数，存在匹配对数量过多或不足的问题。本文提出自适应阈值匹配对选择算法，其核心思想为，空间邻近的重叠匹配对相似性值较大，且随着重叠区域减小而显著降低；空间距离远的非重叠匹配对相似性值很小，且没有显著的变化。图3为相似性值分布，本文采用式(2)所示的幂函数对数据点进行拟合

(2)

式中，a和b是幂函数系数；x和y分别是影像编号和相似性值。可以看出，幂函数能够很好地拟合数据点。然后统计数据点相似性值的均值μ和标准差σ，并定义式(3)所示的截断函数

(3)

式中，k′>0为调节系数。联立式(2)和式(3)可解算出临界值x*，得到匹配对数量。

图3   词汇树检索影像相似性分布

影像匹配对自适应检索可获取具有空间重叠的初始影像对。基于VLAD全局描述子的影像检索不可避免地包含错误匹配对。本文进一步利用局部特征匹配优化影像匹配对，即采用SIFTGPU进行特征匹配，并结合RANSAC基础矩阵估计实现粗差剔除。初始影像对中特征匹配数量大于阈值Tm的影像对保留为正确匹配对P。为了保证影像匹配对P的正确率，本文后续试验均设定Tm＝50。

(1) 场景图构建。正确匹配对P进一步用于构建无向加权图G表示的场景图，建立影像间的拓扑连接关系。假定I＝{ii}和P＝{pij}分别表示n张影像和m个匹配对，无向加权图G＝(V，E)由顶点集合V和边集合E组成。为n张影像中的每一张ii添加一个顶点vi，组成顶点集合V＝{vi}；为m个匹配对中的每一个pij添加一条连接顶点vi和vj的边eij，组成边集合E＝{eij}。同时，为边eij赋权值wij，用于量化影像匹配对的重要性。为了保证后续场景分块和合并的可靠性，边权值wij的设定十分重要。从增量式SfM三维重建的角度看，两张影像的匹配点数量越多、像平面分布越均匀，代表其能够获取更好的位姿解算精度。为此，本文综合特征匹配数量及其分布计算边权值wij，具体见式(4)

(4)

式中，winlier和woverlap根据特征匹配点数量及其覆盖区计算；Rew控制两者权重。因此，wij是winlier和woverlap的线性加权，其计算方法见式(5)—式(6)

(5)

(6)

式中，Ninlier和Nmax_inlier分别是影像匹配对pij的匹配点数量和所有匹配对P中最大匹配点数量；CHi和Ai分别是匹配点凸多边形面积和对应影像面积，如图4所示。根据试验结果，本文的后续试验设定Rew＝0.5。

图4   边权值中重叠面积加权

(2) 分块并行化SfM。根据构建的场景图G，可以实现大场景分割与并行化SfM重建。大场景分割的关键在于子场景尺寸均衡及其强几何连接，以保证各个子场景SfM重建的时间消耗相当且有利于后续子模型合并。与文献[19]类似，本文定义了用于场景分割的尺寸约束Csize和完整性约束Ccomplete。前者限制了子场景的最大影像数量，即Csize≤Tsize；后者则保证子场景与相邻子场景的最小公共影像比例，即Ccomplete≥Tcomplete。同时，为了避免子场景公共影像数量过多导致SfM整体效率下降，本文定义了子场景的最大公共影像数量约束Cimg≤Timg。基于上述约束，场景分割步骤包括：①大场景分割。基于尺寸约束Csize，利用归一化割(normalized cut，NC)算法[26]进行大场景的初始划分，其通过剔除部分连接边实现场景分割，并满足所剔除连接边的权值和最小，进而得到若干连接紧凑的子场景{Gi}。同时，保留子场景之间剔除的连接边集合。②子场景扩展。顺序遍历，如果其对应的子场景Gi和Gj不满足Ccomplete和Cimg约束，则将连接边按照边权值降序排列，并顺序将添加到尺寸较小的子场景 G*＝min(Gi，Gj)，直到满足上述约束。通过上述大场景分割和子场景扩展步骤，可实现场景图G的分割。由于子场景SfM重建相互独立，因而可实现并行化处理。图5展示了场景分割与子场景SfM的重建结果。

图5   场景分割与SfM重建

分块并行化生成的子场景具有独立的坐标系，需通过子场景合并生成最终的全局模型。子场景合并采用SfM子模型之间的对应几何元素建立相似变换矩阵实现，如公共影像位姿和三维点坐标。假设待合并的子模型ms和mr分别表示源模型和参考模型，且存在n个对应几何元素{si}和{ri}。那么，子模型ms到mr的变换关系T＝[λR|t]，可通过最小化目标函数(式(7))得到

(7)

式中，R和t分别是旋转矩阵和平移向量；λ是缩放因子；e2是在变换矩阵T下si到ri的均方误差。式(7)所表示的最小化问题可通过SVD(singular value decomposition)算法求解[27]。考虑到子场景对应三维点数量远大于影像数量，本文将三维点作为对应几何元素。

上述子场景合并需要考虑两个关键因素。由于遮挡、重复或弱纹理影响，SfM子模型的三维点不可避免地包含噪点。如何最大限度降低噪点的影响是第一个因素。为此，本文采用RANSAC稳健估计框架进行相似变换模型估计。考虑到SfM子模型不具备绝对尺度，难以确定合适的阈值进行模型验证。本文利用双向均方重投误差计算RANSAC的几何估计残差ei

(8)

式中，m表示参考模型中影像数量；l表示源模型中影像数量；Pc(Cj，Tsi)表示三维点si在相机Cj上的投影点；xij表示三维点si在相机Cj上对应的影像点。本文设定最大残差为1.8像素，以区分变换T下的内点和外点。

另一方面，随着子场景逐渐合并，参考模型的尺寸变得越来越大。因此，如何快速查找公共三维点是需要考虑的第二个因素。现有算法一般利用子场景之间的公共影像[19]或连接点合并策略[28]建立子场景的公共三维点。但是，这类方法要么减少了候选三维点搜索范围；要么需要很高的内存消耗，无法满足大场景需求。因此，本文提出按需匹配图(on-demand match graph)，以建立特征匹配之间的映射关系。假设待合并的两个子模型分别称为源重建ms和参考重建mr。基于SfM重建的连接点，每个子模型已包含了二维特征点f和三维点之间的连接关系。因此，查找公共三维点变成了建立两个子模型的特征匹配映射关系。

按需匹配对的核心思想是源模型ms的图像数量远少于参考模型mr的图像数量，那么可只使用与源模型图像相关的特征匹配来构建映射关系。如图6所示，源模型和参考模型使用图像i1、i2、i3和i4构建匹配图，其特征匹配的映射关系采用黄线表示。基于建立的匹配图，公共三维点集合CP的查找步骤为：①对于源模型中的每个三维点，找到其特征点；②对于O中的每个特征点，利用匹配图的映射关系，从参考模型中找到其相关的特征匹配点；③对于M中的每个特征匹配点，根据参考模型的连接点关系找到其对应的三维点；④如果CP中不存在公共三维点对，则将其添加到CP中。基于建立的按需匹配图，上述步骤可高效、低内存消耗建立公共三维点。

图6   公共三维点搜索

利用3组无人机影像进行试验，试验区正射图如图7所示。数据采集均采用多旋翼无人机，详细信息见表1。第1组数据包含密集低层建筑物，采用大疆精灵4 RTK无人机采集了3743张影像，部分影像如图8(a)所示。第2组数据包含复杂建筑结构，结合优视摄影测量技术，使用DJI M300 RTK无人机采集了4030张影像，部分影像如图8(b)所示。第3组数据集包含低层建筑和茂密植被，采用5台索尼ILCE-7R相机组成的倾斜摄影设备采集了21 654张影像，部分影像如图8(c)所示。

图7  3组试验数据正射图

表1  试验数据详细信息

图8  3组试验数据

(1) VLAD聚合编码本尺寸k和HNSW最大邻居数M分析。利用数据集1分析编码本尺寸k和最大邻居数M对检索效率和精度的影响。检索效率指VLAD特征聚合、HNSW图索引构建和影像检索的时间消耗；检索精度是正确匹配对与所有匹配对的比值。每张影像检索30张影像，匹配点数超过15的影像对记为正确匹配对。编码本尺寸k取值为32、64、128、256、512和1024，试验结果如图9所示。图9结果表明：①影像检索时间消耗随着k增加而指数增长，从45.7 s增加至175.5 s，主要原因是VLAD特征聚合时间和描述符维度增加，导致图索引构建和影像检索计算量增加。②检索精度随着k增加而线性增加，从0.81增加到0.94。综合检索效率和精度，后续试验设置k为256。另外，最大邻居数M取值为6、8、10、12、16、32和64，试验结果如图10所示。可以看出：①检索效率的变化趋势分为两部分。当M取值从6增加到16时，检索效率基本保持不变；当M值从16增加到64时，检索效率急剧下降。②检索精度的变化趋势分为3个阶段。M值从6增加到8，检索精度明显提高；检索精度在8～16范围内保持不变；检索精度在16～64范围内逐渐降低。因此，后续试验参数M设置为16。

图9   编码本尺寸k对检索算法的性能统计

图10   最大邻居数M对检索算法的性能统计

(2) 基于按需匹配图子模型合并的内存和时间消耗分析。利用数据集1对比分析按需匹配图子场景合并方法的内存和时间消耗，如图11所示。方案1利用所有影像建立匹配图[28]；方案2利用两个待合并子场景的影像构建匹配图[19]。方案1子场景合并的内存消耗最高；方案2在开始阶段的内存消耗与本文方法相当，但其内存消耗不断增加，最终与方案1相当。本文方法的内存消耗在整个合并过程中最低，与方案1相比减少了1.05 GB。另外，本文方法的时间消耗也最小，在整个子场景合并过程中基本保持不变。因此，基于按需匹配图能够实现更小内存消耗、更高效的子场景合并。

图11  子场景合并的内存和时间消耗

(3) 影像匹配对检索与并行化SfM。表2统计了3个数据集的检索效率和精度。本文方法的检索精度较高，达到了90.6%、89.8%和93.5%。从检索效率上看，检索时间分别为1、2、1.8和9.6 min，每张影像的平均耗时为0.040、0.039和0.034 s。因此，本文方法可实现线性时间复杂度的影像索引和检索，满足大规模无人机影像匹配对选择需求。最终，本文方法从3组数据分别选取了59 014、65 743和353 005个匹配对，构建的影像拓扑连接图如图12所示。

表2   影像匹配对检索统计结果

图12   匹配对构建影像拓扑连接图

基于影像拓扑连接图进行大场景分块和并行化SfM重建。场景分割算法的最大尺寸阈值Tsize＝500；重叠完整度阈值Tcomplete＝0.5；重叠影像数阈值Timg＝50，结果如图13所示。图13(c)采用优视摄影测量进行采集，影像位置和倾角随地物几何特性而调整，导致分割结果从侧面看存在交叉重叠区域；图13(e)仅利用下视影像展示分块结果，河流区域下视影像缺失导致分块结果空白。可以看出，3个数据集划分为8、9和44个子场景，进而可进行并行化SfM重建，并最终合并得到全局模型，如图13(b)、图13(d)和图13(f)所示。可以看出，重建三维点可覆盖整个试验区域。表3统计了三维重建的精度和完整度。精度指BA优化重投影误差；完整度可由重建三维点数量或成功定位定向的影像数量体现。可以看到，3个数据集的精度分别为0.542、0.668和0.752像素，参与重建的影像张数分别为3724、4029和21 568张。总体而言，平均99.5%的影像实现定位定向，验证了本文方法的可靠性。因此，本文方法可实现稳健的场景图创建和分块重建。

图13   场景分块与并行重建

表3   本文SfM三维重建方法的统计结果

注：完整度指标括号外数值为重建三维点数量，括号内数值为成功定位定向的影像数量。

(4) 商业和开源软件的对比分析。本文方法与商业和开源软件进行对比分析。开源库包括ColMap(版本号：3.8)和DboW2；商业软件包括Agisoft Metashape(版本号：1.8.2)和Pix4Dmapper(版本号：4.4.12)。试验电脑配有4个2.40 GHz英特尔Xeon E5-2680 CPU，1个10 GB显存NVIDIA GeForce RTX 3080显卡，64 GB内存。

匹配对检索的统计结果如图14所示。与BoW和Dbow2相比，本文方法效率最高，3个数据集的时间消耗为1、2、1.8和9.6 min。对于数据集3，BoW和Dbow2的时间消耗分别为1 335.5和2 848.3 min。对于检索精度，BoW的性能最好，3个数据集的精度分别为90.3%、92.1%和97.6%。本文方法的精度为90.6%、89.8%和93.5%，略低于BoW方法，但均高于Dbow2。综上所述，本文方法能够实现高效、精确的匹配对选择；与BoW相比，它在3个试验数据集上的加速比达到36～108。

图14   3组数据匹配对检索结果

表4统计了三维重建结果。效率是不包含特征匹配的三维重建时间消耗；精度和完整度定义与表3相同。ColMap采用增量式SfM重建引擎；BoW和Dbow2采用本文的并行化SfM重建引擎；Metashape使用多尺度匹配和GNSS进行匹配对选择；Pix4Dmapper使用词汇树影像检索进行匹配对选择。可以看到：①仅采用增量式SfM的ColMap效率最低。采用本文的并行化SfM，BoW、Dbow2和本文算法的效率显著提升，加速比达到30倍以上，明显优于Metashape和Pix4Dmapper。因为内存不足，Metashape和Pix4Dmapper无法重建数据集3。ColMap重建数据集3的时间超过5 d，在此没有统计。②Pix4Dmapper的精度最高，其次是BoW、Dbow2和本文方法。对于重建完整度，除了Pix4Dmapper较低外，其余4个方案的完整度都相近，主要原因是Pix4Dmapper检索正确率较低。值得注意的是，尽管数据集2的影像数仅比数据集1多300张，但其重建时间消耗多4.4倍，主要原因是优视摄影的不规则数据采集方式增加了重建和合并难度。图15展示了6种方法数据1的重建结果。可以看出，本文方法生成的点云可覆盖整个试验区域。

表4   SfM三维重建对比统计结果

注：①对于数据3，ColMap的重建时间超过5 d；Metashape和Pix4Dmapper出现内存不足。②完整度指标括号外数值为重建三维点数量，括号内数值为成功定位定向的影像数量。

图15   6种方法数据1的重建结果

为了对比SfM重建的绝对精度，利用数据2中26个地面控制点进行绝对定向平差试验和精度分析。其中，3个均匀分布点作为控制点；其他点作为精度检查点。表5和图16统计了模型点坐标与地面点坐标之间的偏差。结果表明，Pix4Dmapper具有最高的精度，在X、Y和Z方向上残差分别为0.013、0.016和0.019 m。虽然BoW在垂直方向上精度排名第二，残差为0.036 m，但其水平精度低于本文方法。本文方法在X方向和Y方向上的残差为0.029和0.026 m。由于匹配对选择的精度较低，Dbow2的绝对定向精度在X和Z方向上最低，如图16(a)和图16(c)所示。综上所述，本文方法能够提供相当精度的相对定向和绝对定向三维重建模型。

表5   绝对定向残差统计结果

图16   检查点残差分布

针对无人机影像增量式SfM的匹配对检索代价大、迭代优化误差累积和效率低的问题，本文设计了联合全局描述子和图索引的并行化SfM方法。其核心思想是将大量局部特征聚合为高维VLAD全局向量，并基于HNSW高效图索引结构，实现自适应影像匹配对检索，引导特征匹配和场景图构建。联合按需匹配图的同名三维点搜索和基于双向重投影误差的稳健子模型合并，实现高效的并行化SfM三维重建。利用真实场景的无人机影像进行试验，结果验证了本文的无人机影像并行化SfM方法的有效性。后续研究将深入探索GPU并行计算技术，实现进一步的加速。