近年来，无人机（unmanned aerial vehicle，UAV）已成为一种广泛使用的遥感数据获取工具^[1]。通过搭载非量测相机，无人机获取的高空间分辨率的光学影像具有重要的应用价值^[2]。同时，倾斜摄影测量技术可获取目标的顶面和立面信息，简化目标识别和解译的难度。由此可见，无人机倾斜摄影测量系统可以有效地发挥两者的优势。但由于载荷能力或设备成本限制，无人机平台无法搭载高精度的直接定位系统，影像定向成为无人机倾斜摄影测量成功应用的关键环节。

运动恢复结构（structure from motion，SfM）技术可以同时恢复影像姿态和场景结构信息^[3]。但SfM技术一般用于小规模场景或者低分辨率影像的定向，影像间连接点匹配问题限制了其在无人机影像定向方面的应用。一方面，无人机倾斜影像分辨率高、数据量大、重叠度高导致特征匹配的计算代价高。另一方面，不同视角倾斜影像的初始匹配结果经常包含大量的错误匹配^[4]。因此，基于严格几何模型估计的粗差检测算法，比如随机采样一致性（random sample consensus，RANSAC）算法，往往不能提供稳健、高效的验证结果^[5]。

全球导航卫星系统（Global Navigation Satellite System，GNSS）/惯性测量单元（inertial measurement unit，IMU）信息可用于计算影像地面覆盖，以尽可能减少特征匹配的时间消耗。文献[6]建立了基于重叠面积和相机方向的影像连接图，提高了特征匹配效率，但仅仅依靠重叠准则可能造成影像匹配对的大量冗余。对此，文献[7]利用重叠面积进行边加权，建立影像的拓扑连接图，将提取的骨络相机网络引入SfM，实现了无人机影像的高效稀疏重建，但倾斜影像特征匹配还受到成像角度的影响。因此，在构建影像拓扑连接图时，仅仅利用重叠面积进行边加权是不合理的。另外，拓扑连接图简化可能剔除过多的旁向影像连接。此外，尽管抗外点率能够高达70%，但是RANSAC算法的计算复杂度随外点率呈几何级数增长^[8]。文献[9]利用霍夫变换进行尺度和旋转参数估计，用于剔除初始匹配的外点，从而提高几何验证效率。文献[10]将霍夫变换引入RANSAC算法，用于搜索最可能的变换参数。因此，结合霍夫投票和RANSAC算法可以发挥两者的优势，实现高效、高精度的几何验证^[5]。然而，基于匹配点对穷举比较的点集旋转、尺度参数估计等仍会造成高时间消耗^[11]。

针对上述问题，本文提出考虑影像拓扑连接的匹配对选择算法和基于物方投影的高效几何验证算法，以解决影像匹配对选择和高效几何验证问题，最终建立无人机倾斜影像的高效定向方案。

1 影像匹配对选择

首先利用平台的POS（positioning and orientation system）数据和相机相对于平台的相对姿态数据计算影像的位置和姿态^[12]；然后通过判断影像在平均高程面的地面覆盖是否重叠得到初始影像匹配对。

由于倾斜摄影的影响，上述初始匹配对包含大量重叠区域很小或狭窄的匹配对，这些匹配对特征匹配数量少或正确匹配率低，因此，本文将剔除不满足式(1)的影像匹配对：

$ \left\{ {\begin{array}{*{20}{c}} {{W_o} \ge {W_t} \times {R_o}}\\ {{H_o} \ge {H_t} \times {R_o}} \end{array}} \right. $

(1)

式中，W_t和H_t分别为影像地面覆盖外接矩形的宽度和高度；W_o和H_o分别为重叠区域外接矩形的宽度和高度；R_o为重叠比值。

但是，仅仅采用重叠准则会产生大量冗余的影像匹配对^[7]。为进一步简化，本文设计了最大生成树影像配对选择算法（maximum spanning tree expansion, MST-Expansion）。首先利用初始匹配对构建影像拓扑连接图(reduced topological connection network, Reduced-TCN)；然后提取Reduced-TCN的最大生成树（maximum spanning tree，MST）^[13]，最简化影像连接网络；最后利用局部拓扑连接分析，逐步提高网络连接的稳定性。

1.1 基于重叠面积和交会角度加权的拓扑连接图构建

假定n张影像和m个匹配对分别表示为I={i_i}和P={p_ij}；无向加权图G=(V, E)表示为Reduced-TCN，其中V和E表示顶点和边集合。构建Reduced-TCN：为每一张影像i_i添加一个顶点v_i，构成顶点集合$V = \left\{ {{v_i}:i = 1 \cdots n} \right\}$；为每一个匹配对p_ij增加一条连接顶点v_i和v_j的连接边e_ij，构成边集合$E = \left\{ {{e_{ij}}:i,j = 1 \cdots n} \right\}$。同时，边e_ij权值为w_ij，表示匹配对的重要程度。对于影像匹配，匹配数量可用于描述影像对的重要性。对于倾斜摄影测量，成像角度对特征匹配的影响较大。因此，本文同时利用重叠面积和交会角度计算匹配对的重要性，即按式(2)计算边权值：

$ {w_{ij}} = {R_w} \times {w_{{\rm{area}}}} + \left( {1 - {R_w}} \right) \times {w_{{\rm{angle}}}} $

(2)

式中，w_area和w_angle分别表示重叠面积和交会角度权值；R_w是取值范围为0到1.0的加权比值。加权项按式(3)、式(4)计算：

$ {w_{{\rm{area}}}} = {\rm{area}}/{\rm{are}}{{\rm{a}}_{\max }} $

(3)

$ {w_{{\rm{angle}}}} = \left\{ \begin{array}{l} \cos \beta ,\beta \le {90^ \circ }\\ 0,其他 \end{array} \right. $

(4)

式中，area是重叠区域对应多边形的面积；${\rm{are}}{{\rm{a}}_{{\rm{max}}}} = \mathop {\max }\limits_{{p_{ij}} \in P} \left( {{\rm{are}}{{\rm{a}}_{{p_{ij}}}}} \right)$, 是所有匹配对的最大重叠面积；β表示影像匹配对的交会角。重叠面积和交会角度定义参见图 1(a)。

1.2 基于拓扑连接分析的TCN简化与扩展

Reduced-TCN的最简化连通图可用MST表示^[13]，用于实现最大程度的影像匹配效率提升。当利用影像匹配对重要性进行边加权时，MST表示连接边最少且结构最稳定的连接序列。因此，本文首先利用Kruskal算法^[14]从Reduced-TCN中提取MST。

拓扑连接图MST中绝大多数顶点只保留着生成方向的连接边。为了提高网络的稳定性，本文提出的MST-Expansion算法将沿着扩展方向（与生成方向垂直的方向）添加关键连接边。假设图g_t表示原始Reduced-TCN，首先利用MST的连接关系创建拓扑连接图g_m。则图g_m可看作图g_t的子图，用实线表示图g_m的边集合，然后对图g_m的每个顶点执行扩展操作。

1）根据协方差矩阵的特征向量确定扩展方向。首先从图g_m中查找当前顶点的连接顶点，如图 1(b)所示，红色和绿色圆点分别为当前顶点和连接顶点。然后，计算当前顶点和连接顶点物方坐标的协方差矩阵，且其最小特征值所对应的特征向量表示的方向为扩展方向。如果最大特征值和最小特征值之间的比值满足式(5)，则取消对当前顶点的扩展操作。

$ {e_{\max }} \le {R_e} \times {e_{\min }} $

(5)

式中，e_max和e_min分别表示最大和最小特征值；R_e为特征值比值阈值。

2）扩展区域的已有连接边搜索。首先将扩展方向向量分别沿着顺时针和逆时针方向旋转角度α，得到扩展区域边界，且规定包含扩展方向向量的区域为扩展区域。图 1(b)中蓝色虚线表示扩展区域边界。然后从图g_m中查找扩展区域中已有连接边。当扩展区域的已有连接边数大于阈值T_e时，表示当前顶点的连接稳定，则跳过当前顶点的扩展操作。

3）搜索候选扩展顶点。候选扩展顶点表示与当前顶点在图g_t中存在连接边，但在图g_m中不存在连接边的顶点。对于给定的当前顶点，首先从图g_t中查找所有连接顶点；然后逐一检查这些顶点与当前顶点是否在图g_m中存在连接边，如果不存在，则记为候选顶点。在图 1(b)中，蓝色圆点表示候选顶点。每个扩展区域存在一个对应的候选顶点集合，并按照连接边的权值进行降序排列。

4）连接边扩展。在每个扩展区域内，按照排列顺序，依次将连接边添加到图g_m中，直到新增加的连接边数量达到阈值T_e。按照上述步骤对每个顶点执行扩展操作，可以得到一个连接更加稳固的影像拓扑网络。与原始的拓扑连接图MST相比较，MST-Expansion算法能够增强影像连接网络的稳定性。

2 高效几何验证策略

几何验证通常从初始匹配中估计一个预定义的变换模型。如果待验证匹配点对与模型兼容，归为正确匹配点（内点）, 否则归为错误匹配点（外点）。当外点率较大时，经典的RANSAC算法往往不能提供高效可靠的几何验证结果。对此，本文首先使用影像的粗略POS数据简化同名点的几何变换模型；然后分步进行粗差剔除：①剔除明显的错误匹配，提高初始匹配的内点率, ②利用RANSAC算法的基本矩阵估计进行严格验证模型估计，优化最终匹配结果。

本文设计了一种分层运动一致性约束（hierarchical motion consistency constraint, HMCC）算法，用于提高初始匹配的内点率，设计思想是将同名点在像方的复杂变换模型转换为物方的平移变换。具体做法是：利用影像的粗略POS数据，将同名点投影到物方平均高程面，得到一个由起始端点和终止端点组成的方向向量，包含长度L和相对于物方坐标系X轴的方向θ两个属性（图 2）。该向量是HMCC算法的基本元素，称为对应匹配点的motion。

假定n个初始匹配点组成的motion集合为$M = \left\{ {{m_i}:i = 1,2 \cdots n} \right\}$；每个motion m_i由一个起始点${s_i} = \left( {{x_s},{y_s}} \right)$和终止点${t_i} = \left( {{x_t},{y_t}} \right)$组成；motion的方向θ和长度L可以按照式(6)、式(7)计算得到：

$ \theta = {\rm{a}}\tan 2\left( {{y_t} - {y_s},{x_t} - {x_s}} \right) $

(6)

$ L = \sqrt {{{\left( {{x_t} - {x_s}} \right)}^2} + {{\left( {{y_t} - {y_s}} \right)}^2}} $

(7)

式中，方向θ的取值范围调整为0°~360°；长度L为实数。基于以上定义，HMCC算法的具体实现如下。

1）全局方向一致性约束。其基本思想为：正确匹配motion方向变化范围较小，错误匹配motion方向随机。因此，基于motion的方向统计可以剔除明显的错误匹配。对于motion m_i，首先根据式(6)计算方向θ_i；然后得到n个motion的方向列表${d_{{\rm{list}}}} = \left\{ {{\theta _i}:i = 1,2 \cdots n} \right\}$，如图 2(a)所示，其中，方向列表显示在底部。本文基于方向列表d_list的霍夫投票策略进行错误匹配剔除，以适应初始匹配中外点率高的情况。在获取极值bin后，检查其邻域距离5以内的所有bin值，如果票数超过极值bin的20%，则该邻域bin作为候选bin。最后，将候选bin关联的匹配点标记为内点。

2）局部方向变化一致性约束。局部方向变化一致性表示motion方向变化的平滑性约束（方向变化范围设为0°~30°）。对于motion m_i，首先从集合M中搜索k最近邻子集${M_{nn}} = \left\{ {{m_j}:j = 1,2 \cdots k,j \ne i} \right\}$；然后，计算当前motion m_i与任一邻近motion m_j的方向变化，即方向θ_i和θ_j之差的绝对值，如式(8)所示：

$ {\rm{d}}{\theta _{ij}} = \left| {{\theta _i} - {\theta _j}} \right| $

(8)

式中，|·|表示取绝对值。因此，根据式（8）可以得到方向变化列表$d_{{\rm{lis}}{{\rm{t}}_i}}' = \left\{ {{\rm{d}}{\theta _{ij}}} \right\}$。在图 2(b)中，当前motion用红色向量表示；邻域motion用绿色和橙色向量表示；当前motion的方向变换列表按升序排列。为了降低错误匹配的影响，列表$d_{{\rm{lis}}{{\rm{t}}_i}}'$的中值代表当前motion m_i的方向变化。同样，利用所有motion方向变化进行霍夫投票。在获取极值bin后，检查其邻域距离3以内的所有bin值，如果票数超过极值bin的40%，则该邻域bin作为候选bin。最后，将候选bin关联的匹配点标记为内点。

3）全局长度一致性约束。其基本思想为：正确匹配motion的长度变化范围比较小，错误匹配motion的长度变化随机。因此，长度变化较大的motion可以标为错误匹配。对于motion m_i，首先根据式(7)计算长度l_i；然后得到n个motion的长度列表${l_{{\rm{list}}}} = \left\{ {{l_i}:i = 1,2 \cdots n} \right\}$；最后，根据式(9)进行z-score测试检测错误匹配：

$ z{\rm{ - scor}}{{\rm{e}}_i} = \frac{{{l_i} - \bar l}}{\sigma } $

(9)

式中，l和σ分别是长度列表l_list的均值和标准差；l_i是当前motion m_i的长度。当z-score > 3时，对应匹配点被标记为外点。

3 SfM重建总体方案

基于SfM的无人机倾斜影像稀疏重建流程如图 3所示，主要包括以下几个部分。

1）影像匹配对选择。首先计算影像地面覆盖；然后基于重叠准则得到初始影像匹配对，并构建拓扑连接图；最后利用MST-Expansion算法简化影像连接图。

2）特征提取和匹配。利用SiftGPU^[15]进行单张影像的特征点提取。为了避免大尺寸影像导致内存溢出问题，采用分块策略进行特征提取和特征匹配^[16]。

3）初始匹配几何验证。首先，利用HMCC算法剔除明显错误匹配，提高匹配结果的内点率；然后，基于RANSAC的基本矩阵估计进行严格几何验证，优化最终匹配结果；最后将同一物方点在不同影像上投影所组成的同名点集合构成连接点。

4）增量式SfM重建。如式(10)所示，相机姿态和场景几何的重建问题被描述为非线性目标函数的联合最优化问题，即最小化像方投影点与对应像点观测值的误差和。

$ \mathop {\min }\limits_{{C_j},{X_i}} \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{\rho _{ij}}{{\left\| {Q\left( {{C_j},{X_i}} \right) - {x_{ij}}} \right\|}^2}} } $

(10)

式中，X_i和C_j分别表示三维点坐标和相机参数；$Q\left( {{C_j},{X_i}} \right)$是三维点X_i在相机C_j上的投影点；x_ij是对应的影像点；$\left\| \cdot \right\|$表示L₂范数；ρ_ij是指示函数, 当三维点X_i在相机C_j中可见时，ρ_ij=1, 否则，ρ_ij=0。目标函数的最小化利用开源库Ceres Solver^[17]实现。因为未知数的较好初始值是获取全局最优解的关键，本文采用增量式SfM框架进行相机参数和场景几何的求解^[2]。

4 实验与分析 4.1 实验数据

利用4组无人机倾斜摄影数据验证稀疏重建方案的有效性，详细信息见表 1。其中，数据2包含43个地面控制点，可以用于重建模型的绝对定向精度分析。实验部分的所有处理都在一台Windows 10操作系统的台式机上完成，处理器3.4 GHz Intel Core i7-4770，独立显卡GeForce GTX 770M，显存2.0 GB。

4.2 结果与分析

MST-Expansion算法的参数设置与影像匹配对提取的统计结果参见表 2。MST-Expansion参数设置为：重叠度比值R_o=0.5；加权比值R_w=0.6；特征值比值R_e=3；扩展角度α=45°；扩展阈值T_e=1。其中，利用定向影像数量和定向残差作为指标进行重建实验，得到重叠度比值R_o和加权比值R_w的最佳取值。对于后续实验，本文约定Full-TCN为初始匹配对；Reduced-TCN表示过滤掉重叠面积很小或很窄影像对后保留的结果；MST表示最大生成树简化后的匹配对；MST-Expansion表示本研究算法的提取结果。

表 2结果表明，MST-Expansion算法分别为4组数据保留了640、664、1 516和309个匹配对；Reduced-TCN保留的匹配对数量是MST-Expansion的8倍多；与Full-TCN匹配对相比，MST-Expansion算法可以实现19~35倍比率的匹配对简化，显著地降低了后续特征匹配的时间消耗。

为了分析MST-Expansion算法的简化过程，选择数据1和数据3两个典型的倾斜影像数据集进行分析，并采用权矩阵表示拓扑连接图。其中，权矩阵的非零元素表示影像匹配对；非零元素的取值表示对应匹配对的边权值。图 4(a)~4(c)表示数据1的简化过程；图 4(d)~4(f)表示数据3的简化过程。实验结果表明，Full-TCN包含了很多小权值的匹配对（绿色和蓝色），影像平均匹配数分别为36和69。尽管Reduced-TCN将影像平均匹配数减少到16和24，保留结果中仍然存在很多权值很小的匹配对。图 4(c)和图 4(f)分别表示两组数据基于MST-Expansion简化后的拓扑连接图。可以看出，几乎所有沿着航向的连接边都被提取出来。由于具有较大的重叠面积和较小的交会角度，这些连接边都赋予了较高的权值。对于数据3，大部分被保留的连接边由一张垂直影像和一张倾斜影像组成。这些连接边对于保证影像连接网络的稳定性具有重要作用。最终，影像平均匹配数减小为2，在保证三度重叠的基础上实现影像连接网络的最简化。

图 5显示了HMCC用于数据4中一个匹配对的明显错误匹配点剔除过程，初始匹配的内点率为32.8%。首先，图 5(a)利用全局方向一致性进行投票，剔除具有随机方向的21个错误匹配，图 5(b)利用局部方向变化一致性进行投票，剔除方向变化较大的8个错误匹配；然后，利用长度约束过滤掉剩余错误匹配点，如图 5(b)中的标记为z-score测试的motion。最终，绿色motion对应的匹配为HMCC算法保留的匹配结果，其内点率为42.3%。可见，HMCC算法能够过滤明显的错误匹配，提高初始匹配的内点率。

为验证算法的效率，利用数据3设计了两组实验：①利用Reduced-TCN提取的所有匹配对，不采用分块策略进行特征提取和匹配（实验1）；②利用MST-Expansion选择的匹配对，并在特征提取和匹配阶段采用分块策略（实验2）。两组实验中RANSAC的最大迭代次数设置为100 000，并且将基本矩阵估计作为严格几何约束，其最大距离残差设置为1个像素。图 6显示了LO-RANSAC（locally optimized-RANSAC）^[18]与HMCC-RANSAC的对比结果。对于实验1，部分匹配对外点率高，导致LO-RANSAC的效率很低（时间消耗大于5 s）；相反地，HMCC-RANSAC可以高效地过滤明显错误的匹配，提高效率；LO-RANSAC和HMCC-RANSAC时间消耗分别为43.4 min和7.3 min。对于实验2，一方面，匹配对选择了具有较大重叠区域和交会角度小的匹配对；另一方面，分块策略可以进一步降低错误匹配的发生。因此，两个算法的效率相当，分别为2.1 min和1.3 min。总之，当初始匹配外点率较高时，HMCC-RANSAC算法的几何验证更高效。

在获取实验区域稳健、可靠的连接点后，基于增量式SfM可以实现场景重建。表 3统计了重建结果。结果表明，本研究的技术方案能够实现所有实验数据的完整重建。其中，数据1和2的总体时间消耗大约为35 min；数据4的时间消耗不到15 min；数据3由于具有较大的影像数量和较高的重叠度，因此时间消耗最多（大约74 min）。最终，4组数据相对定向的精度均优于0.5个像素。另外，分析匹配时间消耗与影像数量关系可以发现，本文的重建方案可以实现具有线性复杂度的影像匹配，即影像匹配的时间消耗与影像数量线性相关。

为了评价稀疏重建的绝对定向精度，本研究测量和准备了数据2中的43个地面控制点，如图 7所示。其中，编号为7、9、35和37的4个地面控制点参与光束法平差，实现重建模型的绝对定向。其余的地面控制点作为检查点，用于定向精度验证。检查点在X、Y和Z方向的单点残差如图 8所示。结果表明，在X、Y和Z方向上，检查点最大残差不超过0.1 m；结合残差的随机分布，实验结果表明定向结果不存在系统性偏差。

为了验证影像匹配对选择、分块特征提取和匹配策略对稀疏重建精度的影响，本研究设计了4种重建方案：①保留所有的重叠影像对，并不采用分块策略（Full-TCN）；②保留重叠度比值大于0.5的影像对，并不采用分块策略（Reduced-TCN）；③利用本文的算法提取影像匹配对，并不采用分块策略（MST-Expansion）；④在方案③的基础上，引入分块特征提取和匹配策略（MST-Expansion2）。采用上述控制点方案进行绝对定向，并统计4种方案的重建精度，如表 4所示。实验结果表明，过滤掉重叠面积很小或者非常狭窄的匹配对可以提高重建精度，特别是Z方向的中误差从0.074 m降到0.065 m；虽然实现了影像匹配对19~35倍的简化，MST-Expansion算法使得重建精度在X、Y和Z方向损失了0.11、0.08和0.26倍地面采样间隔（数据2的地面采样间隔为3.67 cm）；分块策略能够增加匹配点数量和改善点位分布，使得模型精度在X、Y和Z方向分别达到0.044 m、0.037 m和0.057 m。

5 结语

针对无人机倾斜影像数据量大、重叠度高、初始匹配外点率高等问题，本文设计了基于拓扑连接分析的影像匹配对选择算法和基于分层运动一致性约束的几何验证算法，并建立了无人机倾斜影像高效SfM重建的总体方案。利用4组不同倾斜设备采集的无人机影像数据，本研究从稀疏重建的效率、完整性和精度等方面进行了深入的对比分析。结果表明，MST-Expansion算法能够显著减少影像匹配对数量，减少特征匹配的时间消耗；HMCC-RANSAC算法能够简化匹配点的复杂空间关系，剔除明显的错误匹配，提高初始匹配的内点率。相对定向和绝度定向结果表明，结合分块特征提取和匹配策略，本文提出的方案能够在保证模型精度的前0提下，实现无人机倾斜影像的高效、稳健的稀疏重建。