利用倾斜影像进行多视图立体视觉三维重建包括倾斜影像空中三角摄影测量、影像密集匹配、基于影像信息的密集点云三维构网与优化以及纹理映射等关键技术^[1-2]。与利用三维扫描仪获取的三维点云不同，利用影像获取的点云数据存在疏密程度不均匀、不够致密和精确等问题, 常常导致三维模型与真实场景在几何一致性上存在差异。因此，对重建得到的三角网格进行优化，构造完整的表面模型是基于影像进行三维重建的重要一环。

基于影像进行三角网格优化是利用投影几何原理，计算原始影像与投影影像构成的影像对间的零均值归一化相关系数(zero mean normalized cross correlation, ZNCC)，对物方顶点求导得到其梯度值，使用梯度下降法对三角网格顶点进行调整，使局部的能量函数代价最小，即通过相关系数求出的网格顶点梯度为零，从而达到优化的目的^[3]。国内外学者对三角网优化进行过许多相关研究。文献[4]同时使用最小内角最大化方法和最小权法进行优化，但是基于影像的三维重建模型优化效果不佳；文献[5]提出采用窄带图割的可见外壳和多目立体匹配融合的方法，结果证明该方法相比直接采用图割算法在整个三角网上进行优化效率更高，但未考虑如何适应大场景、大数据量模型优化的问题。文献[6]给出了利用数字线划图(digital line graphic，DLG)重建地形三维模型过程中基于三角形拓扑关系的构网优化算法；文献[7]证明了采用最新顶点平分法在自适应细化有限元网格上对离散椭圆问题的一致收敛算法。然而它们均属于二维构网优化算法。文献[8-9]分别介绍了早期的网格优化方法是通过优化边缘信息的梯度实现的，但是存在场景限制的缺点。

利用三角网格表示曲面，对三角形的顶点不断地迭代调整，以达到多视图影像相关的一致性，从而实现网格优化的目的。但是在三角网格顶点梯度的计算过程中，由于影像数量多、影像分辨率大等原因，可能会出现大量的冗余，且三角网顶点密度越大，与网格有关的计算越慢^[10-12]。事实上，在三角网格优化过程中，并不是影像中所有区域对应的网格顶点梯度都需要进行计算，如影像的平面区域中，颜色相对单一、纹理不丰富区域的顶点由于计算出的梯度过小，因此无需重复计算。同时，在常规的网格优化预处理过程中，模型表面的平面区域经过网格细分会产生大量密集的小三角形，给优化过程带来额外的计算，降低网格优化效率。因此，效率过低是优化过程中出现的最主要问题之一^[13-19]。

针对以上问题，本文提出一种自适应快速优化方法。在第一次迭代过程中，计算得出顶点梯度初值后，对网格中的每一个三角形进行评价并自适应标记为活跃区与怠惰区两个不同区域。活跃区是迭代计算过程中顶点移动量较大的区域，此区域对优化效果影响较大；怠惰区是迭代过程中顶点移动量较小的区域，对三角网格的最终优化效果影响较小，在后续迭代过程中无需重复计算怠惰区三角形顶点梯度。在将三角网进行初始标记后，进一步通过图割算法获得三角网的最优标记, 选择性地优化活跃区三角网，从而达到快速优化的目的。

1 利用影像信息进行三角网优化

三角网优化是指利用影像信息对重建后的原始模型网格进行优化处理，从而获得最优的网格模型。利用影像信息迭代计算网格顶点梯度，将三角网的顶点沿着梯度方向进行移动，使得网格顶点对应同名像点间的相关性达到最大，当计算出的网格顶点梯度为0时，迭代结束，获得最优的网格模型。

1.1   影像对的相似性测度

如图 1所示，对于网格S中的任意一个三角面t，都可在影像对上找到公共重叠的区域。对于网格中某个三角形t中存在的P点，在同名影像对I_i、I_j上有其对应的同名像元x_i、x_j，对于影像I_j上的像元x_j，利用影像I_j的内外方位元素将x_j反投到物方空间得到三角形t内的一点P，将P以影像I_i的内外方位元素进行投影得到在影像I_i视角下的像元点$I_i^{S, j} $，即通过投影可将影像I_j经物方的三角形t的诱导投影到影像I_i的视角之下得到影像$I_i^{S, j} $。

图  1  基于影像像元的三角网重投影

Figure  1.  Relation of Image Re-project

下载: 全尺寸图片 幻灯片

以ZNCC作为影像相似性的度量，定义两个像元点之间的灰度值方差C(i，j，x)为：

式中，N表示以当前像元为中心的窗口内的像元总数；I_i(x_k)是在影像I_i上像元x_k的灰度；$ I_i^{S, j} (x_k)$是在投影影像上像元x_k的灰度；$ {\bar I_i}\left( {{x_k}} \right)$、$\bar I_i^{S, j}\left( {{x_k}} \right) $是该像元在对应影像上的灰度的均值。则ZNCC $ M({I_i}({x_k}), I_i^{S, j}({x_k}))$的计算公式为^[20]：

1.2   顶点梯度计算

对于网格平面上的一点P，其对应的可视影像I_i上的像元x_i存在映射关系$ {x_i} = \prod\nolimits_i P $，在影像I_j上存在映射关系$ {x_j} = \prod\nolimits_i P $，则P = ∏_i^-1x_i。点P在网格平面S对应微小面元δS，则P点的梯度G_S (P)为：

式中，N是物方点P所属三角面的法向量; d_i是影像I_i的相机中心指向点P的方向向量。

表示影像对相似性度量在像元x_k处相对于影像坐标的偏导数，可由该点处的相关系数ZNCC对像点坐标求偏导得到。

表示影像坐标关于物方点坐标的偏导函数。

将P点的移动转换为对应面元的移动：

对于每一个顶点来说，其所参与构建的多个三角形投影在影像上后都包含多个像元点，因此在投影之后对于每一个像元点q来说，都可由其所属的三角形的3个顶点v_k来表示其位置$q = \sum\limits_k {\varphi _k}{v_k}\left( {k = 0, 1, 2} \right)$ ，且$\sum\limits_k {\varphi _k} = 1$ ，φ是三角形的某一顶点对三角形内部一点的权重。

物方上q点的梯度G_S(q)也可用三角形3个顶点的梯度G_S(v)来表示：G_S(q)= $\sum\limits_k $φ_kG_S(v_k)。因此可通过稀疏矩阵求得三角网中每一个顶点的梯度：

式中，m表示像元点的数量；n为顶点的个数；A中每一行的元素表示为：（1）每个影像点对应的三角形顶点所在的元素位置由所属三角形顶点的权重决定，（2）无关三角形顶点所对应的元素位置由0填充。以顶点的梯度为未知数，矩阵A与对应影像点的梯度为已知值，通过构建稀疏矩阵计算出所有顶点的梯度值，再对网格顶点进行统一的移动操作，即可得到优化后的网格模型。

2 网格三角形自适应标记策略

三维模型表面网格优化的实质是利用影像灰度信息，通过对灰度求其关于顶点坐标的梯度，将顶点沿三角面的法向量方向进行移动从而得到最优顶点。在计算过程中发现，在纹理丰富区域计算出的顶点梯度值较大，而在纹理不丰富的平坦区域所计算出的顶点梯度值则明显偏小。利用网格中三角形所有顶点的梯度，计算网格中每一个三角形的位移量以及优化该三角形的时间消耗，将网格划分为纹理丰富优化效果明显的活跃区和位于平坦区域优化效果不明显的怠惰区两部分。图 2所示为本文优化流程。

图  2  本文优化流程

Figure  2.  The Process of Our Method

下载: 全尺寸图片 幻灯片

2.1   构建精度-效率函数

将三角网格中的每一个三角形按照所属的区域分别标记，将所属活跃区的三角形标记为true，怠惰区的三角形标记为false。三角网格表面S上的每一个三角形t都有其标记f，即f（t）=true或f（t）=false。

构建三角形精度-效率函数u(f)：

式中，u_accuracy (f)表示网格优化的精度函数，主要由活跃区三角形优化的位移量决定，因此可由三角形顶点梯度的模来表达；；u_{time_efficient_loss}(f)表示网格优化的时间减少量函数，主要通过减少优化怠惰区三角形来达到减少优化时间的目的，因此可由怠惰区三角形的初次优化时间来表达。为了保证加速优化后的网格能够在尽可能保持优化效率的前提下，仍能够与全局优化后的网格保持最大程度上的吻合，应保证此精度-效率函数取得最大值。

2.2   网格三角形自适应标记

由于在优化的过程中每个顶点的初始移动量不同，因此每一个顶点的优化程度也不相同。以顶点梯度的模作为顶点初始移动量，对于网格表面中的某一三角形t，顶点v_i是三角形的某一顶点，i = 1，2，3。顶点v_i的初始移动量g_i可表示为：

式中，grad(v_i)是此顶点的梯度。

对于三角形t的初始移动量g_i可表示为：

每一个三角形的优化都由其3个顶点决定。三角形顶点的梯度则是三角形内部点梯度的矢量和，因此三角形的优化时间可由三角形的面积来进行量化，进而对某一三角形t的时间代价函数可表示为：

式中，v₀、v₁、v₂是三角形t的3个顶点。

对于三角网格中的每一个三角形都定义一个效率代价c_e，以每一个三角形所对应的初始移动量g_i与优化时间消耗t_c的比值作为三角形的效率代价。三角形的效率代价越高, 表示被标记为活跃三角形的可能性越大。

2.3   三角形的标记

2.3.1   获取三角形初始标记

将三角网格的每一个三角形的效率代价按升序进行排序，可以得到一条如图 3（以优化实验中教堂的c_e曲线为例）所示的曲线。图 3中，横坐标为排序后三角形在序列中的序列号，纵坐标为对应的c_e值。

图  3  效率代价曲线图

Figure  3.  c_e Curve

下载: 全尺寸图片 幻灯片

将c_e值的分子与分母分别进行归一化，以三角形的面积作为时间统计量，那么通过标记可实现的时间减少率为：

即怠惰区三角形优化时间占总优化时间的比例。标记后的精度损失率由怠惰区三角形的总位移量与所有三角形的总位移量的比值表示：

式中，lazy表示怠惰区三角形；all表示网格中的所有三角形。即怠惰区三角形的移动量相对于全局优化三角网移动量的比值。l、r的关系可用图 4表示。如图 4所示，当曲线在原点（0, 0）时，精度损失最小，时间消耗最大，此时为全局优化；随着优化的推进，所需时间逐渐减少，优化后网格的精度损失逐渐变大，需获得曲线上的最佳平衡点。根据此结果对网格三角形进行标记，保证优化效率的提升与精度损失达到最佳的平衡。

图  4  精度损失-时间减少曲线

Figure  4.  Accuracy Loss-Time Reduction Curve

下载: 全尺寸图片 幻灯片

将精度-效率公式（8）改写为：

以l、r为自变量，将式（9）改写为：

式中，w_l、w_r分别表示l、r在函数中的权重。必存在一点e(r₀, l₀)，使得该函数取得最大值。通过对函数求导，最优解点e(r₀, l₀)在图 4中函数曲线上对应的斜率为$k = \frac{{{w_r}}}{{{w_l}}} $，这个点就是将三角形分为两个属性的最优点。在本文实验中，视l、r等权，即取$k = \frac{{{w_r}}}{{{w_l}}}=1 $，通过微分法求得曲线斜率为1时对应的c_e值，获取网格三角形的初始标记f^initial，小于该值的三角形属于怠惰三角形，而大于该值的三角形属于活跃三角形。

2.3.2   图割获取最优标记

由于初始标记后的网格中某些区域存在一个或数个单独标记的三角形，为了达到最优的标记效果，使用图割算法获取每个三角形的最优标记。

对于三角网格中的每个三角形t都有其初始的标记f（t）={true, false}，以三角网中的每个三角形为图割过程中图的节点，以邻接三角形的邻接边作为所构建的节点间的边构建图，则图的能量函数可表示为：

最优的标记依赖于三角形的初始标记f^initial，E_data(f)是对于所有的三角形可能存在的不同标记的权值和，即：

而对于每一个三角形与所构建的终端节点的连接边的权值为：

本文中使用Potts模型进行平滑项的设置：

其中p、q是具有邻接关系的两个相邻三角形，则：

通过以上的图割过程得到每一个三角形的最优的标记f^optimal，由图 5~6可以看出，经过图割，纹理丰富区域的活跃三角形密度更大，而平坦区域的怠惰三角形密度更大。图 5~6中红色区域为标记的活跃区三角形，白色为怠惰区三角形。

图  5  图割前标记

Figure  5.  Before Graph-Cut

下载: 全尺寸图片 幻灯片

图  6  图割后标记

Figure  6.  After Graph-Cut

下载: 全尺寸图片 幻灯片

3 三维模型快速优化实验与分析

3.1   实验数据及结果

分别选取标准影像，如图 7~8所示的喷泉与教堂影像，以及建筑物丰富、特征明显的某测区152张无人机倾斜影像数据（见图 9）。在构建初始三角网后，采用本文所给的三角网快速优化方法进行优化重建并与全局优化方法^[21]进行对比（测试用机CPU E3-1231 v33.40 GHz，内存8 GB）。图 10~11分别为利用全局优化方法与本文方法对原始影像优化后的效果。图 12~14为针对无人机影像不同优化方法局部区域优化前后的比较。表 1~2分别为标准影像与无人机影像网格模型优化前后的时间对比。为比较本文方法与全局法三维重建两个点集的匹配（相似）程度，分别对标准影像及无人机影像数据，以Hausdorff距离为指标，统计了本文方法加速优化网格与全局优化网格间采样区的Hausdorff距离与平均距离（见表 3）。图 15为无人机影像整体优化后的效果。

图  7  喷泉影像

Figure  7.  Fountain Image

下载: 全尺寸图片 幻灯片

图  8  教堂影像

Figure  8.  Church Image

下载: 全尺寸图片 幻灯片

图  9  无人机倾斜影像

Figure  9.  Unmanned Aerial Vehicle(UAV) Image

下载: 全尺寸图片 幻灯片

图  10  全局优化效果

Figure  10.  Global Refining Result

下载: 全尺寸图片 幻灯片

图  11  本文方法的优化效果

Figure  11.  Refining Result of the Proposed Method

下载: 全尺寸图片 幻灯片

图  12  优化前的局部细节图

Figure  12.  Local Region Before Refining

下载: 全尺寸图片 幻灯片

图  13  本文方法优化后的局部网格图

Figure  13.  Local Region After Refining with Our Method

下载: 全尺寸图片 幻灯片

图  14  全局方法优化后的局部网格图

Figure  14.  Local Region After Global Refining

下载: 全尺寸图片 幻灯片

表  1  网格模型优化时间对比

Table  1.  Comparison of Mesh Refining Time

影像名称	影像分辨率/像素	网格顶点数/面数	活跃三角形个数	本文方法优化时间/s	全局方法优化时间/s	效率提升/%
喷泉	3 072×2 048	190 306/379 922	80 318	289.587	1 232.293	76.5
教堂	3 072×2 048	305 952/610 812	95 961	1 428.181	4 660.299	69.4

下载: 导出CSV 

| 显示表格

表  2  无人机倾斜影像数据模型优化时间对比

Table  2.  Comparison of UAV Image Reconstruction Model Refining Time

影像数	影像分辨率/像素	网格顶点数/面数	活跃三角形个数	本文方法优化时间/s	全局方法优化时间/s	效率提升/%
152	1 988×1 326	1 299 699/2 579 931	1 310 652	2 698.60	8 994.09	70.1
	3 976×2 652	1 552 672/3 083 567	1 738 133	6 556.32	20 076.40	67.3

下载: 导出CSV 

| 显示表格

表  3  本文方法优化结果与全局方法优化结果匹配程度/m

Table  3.  Matching Degree Between Our Refining Results and Global Refining Results/m

数据名称	网格间Haudorff距离	网格间平均距离
喷泉模型	0.005 163	0.000 168
城堡模型	0.006 370	0.000 327
无人机数据	0.002 137	0.000 175

下载: 导出CSV 

| 显示表格

图  15  整体优化后的结果

Figure  15.  Overall Refining Results

下载: 全尺寸图片 幻灯片

3.2   实验结果分析

由表 1及表 2可以看出, 本文方法相较于全局优化效率提升明显，表 3反映了本文加速优化后网格与全局优化后网格的匹配程度。由表 3可以看出，两种方法优化后网格极为接近，匹配度较高。结合表 1~2结果可以得出：本文加速优化方法相较于全局优化方法，在精度不降的前提下，优化效率提升明显。通过以上实验分析发现，在基本保持原有精度的前提下，利用本文提出的三维重建网格自适应快速优化方法，在优化的过程中既能实现对初始三角网的处理，保持优化达到预期效果外，还能有效地减少网格顶点的迭代计算时间，提高了优化效率。

4 结语

针对基于无人机倾斜影像大场景三维建模网格优化过程中效率低的问题，本文在基于影像的网格优化方法中对网格顶点梯度初始计算结果进行评价，将优化网格中的三角形标记为活跃区三角形与怠惰区三角形。对活跃区三角形顶点的梯度进行迭代计算，忽略怠惰区三角形顶点梯度的计算，减少冗余计算，达到在基本保持原有精度的前提下提升优化效率的目的。标准影像及对地观测真实无人机影像的实验结果均证明了本文方法的有效性。

虽然本文方法在实现优化效率快速提升的同时，模型的精度基本与全局优化后的精度一致，但在网格优化的过程中未能同时实现精度的提升。如何在保持网格优化效率的同时，提高网格优化精度，实现效率与精度的共同提升将是后续的研究内容。