倾斜摄影测量在三维重建、工程测量和灾害应急响应等领域应用前景广泛^{[1, 2]}。作为其关键技术，多视密集匹配(multi-view stereo，MVS)也成为诸多领域的研究热点。

MVS是一种在已经恢复相机方位的前提下，重建场景表面稠密的三维点云的精确3D重建技术。它可以利用多视影像冗余信息，最大限度地减少信息盲区，处理影像中的遮挡问题，提高匹配的可靠性。近年来，国内外很多学者相继提出了大量实用的算法：Goesele等^[3]提出以稀疏匹配点作为种子，进行区域增长的方法；Furukawa等^[4]提出基于物方面元实现由稀疏面元扩散到密集面元的多视匹配的算法；Campbell等^[5]提出基于离散的MRF模型来选择最终的深度估计值的方法；王竞雪等^[6]提出对像方特征点和物方平面元进行集成的密集匹配算法。密集匹配算法总体上可以分为体素法、表面演化法、特征点生长法和深度图合并法4类^[7]。其中，基于深度图的方法更适合大尺度场景，本文采用的基于PatchMatch的多视密集匹配方法，即基于深度图进行场景重建。最初，PatchMatch作为在图像块间快速寻找近似最邻近匹配的方法而被提出。Bleyer等^[8]把PatchMatch引入到双目立体视觉，可以在校正图像对的视差空间图像(disparity space image, DSI)找到每个像素的近似最优支持平面。基于PatchMatch算法，可以采用不同的匹配代价函数，确定PatchMatch算法的传播和随机搜索过程中深度信息的正确性。NCC(normalized cross-correlation)是多视匹配中一种常用的匹配代价函数，它对光照和材料变化不敏感，且匹配精度较高^[9]，但NCC在弱纹理和重复纹理区域可能会失败，在物体边界处会出现边缘膨胀，影响匹配的质量。

本文提出一种具有边缘特性的ANCC(the adaptive normalized cross correlation)作为匹配代价，实现倾斜影像的密集匹配过程，有效地改善了生成点云的质量。

1 算法主要流程

本文将每张图像作为一次参考图像，为N张图像I={l_i}产生深度图。核心算法采用了PatchMatch，通过传播和随机搜索得到深度信息。图像密集匹配流程包括：①由SFM得到稀疏点云影像; ②选择图像; ③基于PatchMatch算法计算深度图(深度和法向量赋初值、PatchMatch对深度和法向量迭代优化和得到深度图); ④融合深度图; ⑤形成三维点云。

1.1 图像选择

为了合理地选择图像，采用类似文献[10]的两层视图选择算法。在图像层，全局视图选择，为每张参考图像选择合适的邻域视图集合；在像素层，局部视图选择从全局视图选择得到的邻域视图集合中针对每个像素分别确定一张最佳匹配图像，最佳图像往往跟当前待匹配像素相关。

1) 全局视图选择。对每个参考图像R，全局视图选择可以确认R的邻域视图集合N。选择过程中主要考虑场景内容和尺度以及图像间重叠度信息，根据文献[10]中的评分函数来衡量每个候选邻域图像的质量，从而确定最佳的候选图像。将SFM(structure from motion)阶段产生的特征点作为指示器，图像共享特征点越多，图像重叠度就更高。共享特征点最多的图像也许不能提供足够的基线进行精确重建，另外，SIFT特征检测匹配了不同分辨率图片的特征点，这样的分辨率差异使得密集匹配存在问题。

对每个参考图像R，其邻域视图集合N中的每个图像V的评分函数定义参见文献[10]。邻域视图集合N中所有图像对的积的权重函数使相邻图像间有更好的视差范围，衡量图像R和V在特征点f处分辨率的相似度的权重函数使图像间有更相近的分辨率。这两种权重函数的计算方法参见文献[10]。

2) 局部视图选择。全局视图选择为参考图像确定了邻域视图集合N，N集合中候选图像的数量可以根据实际情况进行设置。在参考图像每个像素位置处，从集合N中的图像中针对当前像素选择一个最合适图像A。在参考图像逐像素匹配过程中，根据当前像素的局部视图选择标准来不断更新A图像。

给定当前像素的深度和法向，可以计算当前像素在相机坐标系下对应的空间三维点，分别得到各候选搜索图像I_j和参考图像I_i在当前的可视点下相机中心光线之间的夹角θ_ij。对每个候选图像，计算其和参考图像的相机间光轴中心的距离d_ij。根据夹角和距离的乘积θ_ij×d_ij，将邻域视图集合中的候选图像按升序排列，取乘积最小且空间三维点在图像上可见的候选图像作为当前像素最佳搜索图像，参与当前像素的匹配。

1.2 深度图计算

在图像深度图的计算过程中，本文提出一种具有边缘特性的ANCC作为匹配代价函数，采用PatchMatch算法寻找最优深度值，计算每张图像深度图。对于图像上每个像素，尝试寻找其对应的空间中三维点的局部切平面，使得参考图像和其搜索图像在该点匹配代价最小。局部切平面本质上是场景表面的切平面，通过三维点X的坐标和它的法向量n来表示。

对于参考图像I_i和搜索图像I_j, 给图像I_i中每个像素p分配一个三维平面，p对应的三维点X_i一定位于p所在的直线上，设d为深度值，x为像素p的坐标，则X_i在C_i坐标系中坐标为：

$ {X_i} = {\rm{d}}{\mathit{\boldsymbol{K}}_i}^{ - 1}x $

(1)

在本文中，参考图像和搜索图像的相机参数分别是{K_i, R_i, C_i}，{K_j, R_j, C_j}, 其中K是内参数，R是旋转矩阵，C是相机中心。像素对应的平面在相机C_i坐标系中表示为f_p={X_i, n_i}，X_i为三维点，n_i为法向量。由文献[11]可得到图像间对应点之间的单应矩阵H_ij。

与传统方法使用NCC作为匹配代价不同，本文采用ANCC作为相似性测度。文中的匹配代价ANCC，即将像素色度空间和几何空间邻近性引入权重计算, 采用双边权重方法计算对应像素块间的匹配代价。ANCC定义为：

$ C = \frac{{\sum\limits_{i = 1}^M {{w_R}{w_S}} \left( {{I_R} - \overline I {_R}} \right)\left( {{I_S} - \overline I {_S}} \right)}}{{\sqrt {{{\sum\limits_{i = 1}^M {\left| {{w_R}\left( {{I_R} - \overline I {_R}} \right)} \right|} }^2}} \sqrt {{{\sum\limits_{i = 1}^M {\left| {{w_S}\left( {{I_S} - {{\overline I }_S}} \right)} \right|} }^2}} }} $

(2)

式中, C值的范围为[-1, 1]，为了使参考图像在当前像素p的匹配代价值非负，这里的匹配代价设为m=1－C。由单应矩阵H_ij可以确定像素p在搜索图像中的像素p′，得到搜索图像中对应的像素块。其中，I_R和w_R是参考图像上以像素p为中心的像素块中像素的灰度值和权重值。I_S和w_S是搜索图像上以像素p′为中心的像素块中像素的灰度值和权重值，M为像素块中的像素数量，权重值w_R、w_S可由双边权重函数w得到，有：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;w\left( {p, q} \right) = \\ {\rm{exp}}\left(- {\frac{{ {{\left\| {p - q} \right\|}^2}}}{{2{\sigma _d}^2}} - \frac{{{{\left\| {I\left( p \right) - I\left( q \right)} \right\|}^2}}}{{2{\sigma _r}^2}}} \right) \end{array} $

(3)

式中，‖p－q‖表示像素p和q的空间距离；‖I(p)－I(q)‖表示色彩空间中灰度值的欧氏距离；σ_d和σ_r分别表示控制空间和颜色差异的参数。

${\overline I _R}$的定义为：

$ \overline I {_R} = \sum\limits_{}^i {{w_R}{I_R}/\sum\limits_{}^i {{w_R}} } $

(4)

式中，${\overline I _R}$为参考图像像素块的灰度加权平均值；${\overline I _S}$可以根据搜索图像上对应像素块中像素的权重值和灰度值采用类似方法计算得到。

使用PatchMatch算法来计算图像的深度图，分为初始化、传播和随机搜索3个过程。首先，由SFM得到稀疏三维点，通过构网和插值来初始化图像上其他像素点对应的深度值和法向值。随后，进行像素间的深度值和法向值的传播及随机搜索，通过多次迭代进一步优化像素所关联的切平面的深度值d和法向量n。

传播过程可以传播空间邻域像素的平面信息到当前的像素。传播的前提是相邻像素可能有相似切平面，即相似的深度和法向。对于高分辨率图像，像素个数足够多，则在初始化的深度和法向值中很可能至少有一个像素有接近真值的深度值和法向值。

设在图像(x, y)坐标处像素p当前的深度和法向分别为d₀、n₀，其中，x、y表示像素的坐标值。将(x－1, y)和(x, y－1)处像素的深度d和法向n传递给(x, y)处像素，检查在坐标(x, y)的像素对应的匹配代价m(d, n)和m(d₀, n₀)。如果m(d, n) < m(d₀, n₀), 则采用d和n作为当前(x, y)处像素新的深度和法向，反之则深度法向不变。在偶数次迭代中，从左上像素开始到右下像素，传播左上邻域, 即(x－1, y)和(x, y－1)，遍历所有像素；在奇数次迭代中采用相反的顺序，从右下到左上，传播右下邻域，即(x+1, y)和(x, y+1)。在传播过程后进行随机搜索过程，进一步精化平面。随机搜索可以通过测试一些随机的深度和法向值，进一步减少匹配代价m，同时避免收敛于局部最小值。

给定一个搜索范围[δd, δn]，①在随机搜索范围内随机选择参数[d′, n′]，使得d′∈ [d－δd, d+δd], n′∈ [n －δn, n +δn]; ②如果匹配代价m(d′, n′) < m(d₀, n₀), 则接受新的参数[d′, n′]，令当前的深度和法向分别为d′和n′，反之则深度法向不变；③返回第①步。上面的过程重复k次。

经过一定次数的传播和随机搜索，即可得到每张图像对应的深度信息。

1.3 深度图融合

在上述过程后，得到了图像的一系列的深度图，再将所有深度图合并表示一个场景。采用基于可见性的多视深度图融合的点云生成方法，对不同角度深度图上的同名点对的深度进行加权平均，计算融合后的深度值，对像素深度中偏离超过一定阈值的，则认为是异常值。由深度信息和相机模型信息即可得到场景的三维点云。

2 实验及分析

实验数据共分两组。第1组采用的是某地区拍摄的倾斜航空影像，采用的相机为Sony QX1，图像大小为2 728像素×1 816像素，主距为20.492 mm，像素大小为4 um。利用VisualSFM软件进行图像处理，利用SIFT匹配特征点，稀疏重建，得到稀疏点云。图 1为本文试验的航空影像。

本文实验中对于ANCC算子采用的像素块尺寸为15像素×15像素，权重中的空间和颜色差异参数σ_d=11，σ_r=11, 迭代次数k=3，全局视图选择中邻域视图集合的图像个数为4，α_max=10，对于NCC作为匹配代价函数时图像像素块尺寸为7像素×7像素，迭代次数k=3，全局视图选择参数和ANCC相同。

将第1组实验选用的全部影像融合，比较得到的点云，如图 2、图 3所示。图 2为局部三维点云；图 3为整体三维点云。其中，图 2(a)、图 2(c)、图 2(e)和图 3(a)为NCC为匹配代价; 图 2(b)、图 2(d)、图 2(f)和图 3(b)为ANCC作为匹配代价。

从结果来看，采用NCC生成的点云数量为3 663 956，点云整体存在很多漏洞，尤其弱纹理的区域；采用ANCC的点云数量为4 611 592，点云数量增加，同时漏洞减少，重建效果得到改善。在屋顶、路面等区域点云密度更高，没有出现明显的漏洞区域。出现这种情况的可能原因是图像中存在的噪声，在NCC计算中，采用中心像素块的像素信息计算相应的平均值和标准差，对像素块内的所有像素给予相同的权重，而ANCC的双边加权可以一定程度减小噪声影响，匹配结果更好。

第2组实验采用包含11张图片的Fountain-P11数据集^[12]，分辨率是3 072像素×2 048像素，分别采用和第1组实验相同的实验参数，由ANCC和NCC代价函数分别得到融合后喷泉处的点云。如图 4所示。

将密集匹配得到的点云分别和激光点云配准，采用ICP方法通过CloudCompare软件进行精度的比较，采用NCC得到的点云数量为8 930 571，RMS值为0.029 7 m；采用ANCC得到的点云数量为9 262 260，RMS值为0.028 4 m。可以看到，两组点云的整体效果都很好，精度相当，ANCC精度略高。

3 结束语

本文采用基于PatchMatch的方法进行多视密集匹配，通过一种两层视图选择方法确定搜索影像，并利用具有边缘特性的相似性测度，得到密集匹配的点云。通过试验验证了本文方法的可靠性，试验结果表明，本文的算法增加了点云的数量和点云的完整性，提高了密集匹配的质量。