基于全局优化的实时高精度模型重建

许新傲，李艺璇，钱佳铭，冯世杰，左 超*

（1.南京理工大学 电子工程与光电技术学院 智能计算成像实验室（SCILab），江苏 南京 210094；2.南京理工大学 智能计算成像研究院（SCIRI），江苏 南京 210019；3.江苏省光谱成像与智能感知重点实验室，江苏 南京 210094）

1 引 言

目前，光学三维测量技术由于其在非接触、测量速度、测量精度等方面的突出优势，已被广泛应用于工业制造、逆向工程等许多领域［1-5］。在主动式光学三维测量方法中，颇具代表性的测量方法是条纹投影轮廓术（Fringe Projection Profilometry，FPP）。它不仅硬件搭建简单，还具有测量精度高、适用性广、实时性能好等优点，是目前最受欢迎的三维测量技术之一［6-8］。但是受相机视野限制，单次结构光三维测量只能获得被测物体一个角度的三维信息，想要获取物体全方位的三维模型，需要从不同角度获取物体的三维信息。Dana solav等人［9］通过在被测物周围增加额外的相机来获取物体各角度的三维信息。Epstein等人［10］通过使用平面镜在不同角度搭建虚拟相机与投影仪来实现高速全景三维形状测量，但这增加了系统的成本与复杂度且标定工作量与工作难度激增。

随着计算机视觉、计算机图形处理等领域的不断发展，研究人员围绕三维点云配准方向做出了许多出色的工作［11］。Rusu 等人提出了点特征直方图（Point Feature Histograms，PFH）以及快速点特征直方图（Fast Point Feature Histograms，FPFH）作为鲁棒且具有位姿不变性的多维特征，用于描述3D 点云数据集中点周围的局部几何结构［12-14］。迭代最近点算法（Iterative Closest Point，ICP）的提出［15］为点云配准工作的进一步发展打开了大门，并成为迄今为止最流行的点云配准方法。许多研究者围绕ICP 算法开展了一系列的工作，从计算速度、收敛判断、迭代方案、全局对齐等问题出发对ICP 算法进行改进与创新，衍生出了许多优秀的变体［16-22］，如Segal 等人提出的一种广义的ICP 算法（Generalized-ICP，G-ICP）［20］，Koide 等人提出的体素化广义迭代最近点算法（Voxelized Generalized-ICP，VGICP）［21］。Zhou等人提出一种面向局部重叠三维曲面的快速全局配准算法（Fast Global Registration，FGR），降低了算法对点云位姿初始化的要求与计算成本，同时也更加稳定［22］。点云配准方法的进步促进了对如何更便捷地获取物体完整三维模型问题的研究。Szymon 等人［23］提出一种基于结构光三维扫描仪与ICP 算法的实时360°三维模型获取技术。它允许用户手动旋转物体，并查看不断更新的模型，但是条纹边界编码策略获取的点云密度相当稀疏，整体三维模型成像效率较低，且仅通过ICP 算法对点云进行配准的配准方案限制了它的精度。此外，以即时定位与地图构建（Simultaneous Localization and Mapping，SLAM）为首的快速建图领域迅速成长并衍生出了诸如ORB-SLAM2、Kinect Fusion 等基于RGBD 相机进行快速建图的经典之作［24-26］。然而上述算法更加注重实时性且所使用的RGBD 相机的分辨率较低，难以满足高分辨率下的物体高精度三维重建的需求。

为了能够实时无约束地获取物体完整的高分辨三维模型，Qian 等人［27］基于FPP 系统提出了基于由粗到细配准策略的高分辨率实时360°三维模型重建方法，能够在不借助任何辅助仪器的情况下完成高精度无约束的全方位三维模型重建。该方法通过2D-3D 关键特征点匹配实现粗配准方案，虽然能够保证相邻帧点云间的高速配准，实现物体高精度三维模型的在线拼接工作，但是在应对抖动等不稳定情况下，由于特征点的错误匹配，会影响到相邻帧点云间位姿关系初始化的准确性，并且对待测物体完成全局扫描后，在重复扫描区域容易出现点云重影问题。

对此，本文提出了一种基于全局优化的实时高精度模型重建方法，针对粗配准中点云位姿初始化不稳定和全局扫描后点云重影问题进行了优化。首先，继承了Qian 所提方法中的粗配准方案，基于透视点问题求解方法（Perspective-n-Point，PnP）的非迭代变体EPnP 方法，对获取的点云数据进行实时处理得到相邻帧点云间的位姿关系。然后，引入点云法向量信息对点云初始位姿进行二次优化，能够增加粗配准点云位姿初始化的稳定性，为后续精细配准环节提供符合条件的初始位姿。在精细配准环节中，提出了一种基于图优化的全局点云位姿优化算法，在原有ICP 成对配准方案基础上，增加了点云全局位姿优化步骤，消除了连续ICP 配准后出现的累计误差，克服了物体旋转一周后出现首末两帧点云重影的问题。

2 原 理

2.1 快速三维成像

为了实现实时点云配准，必须快速获取物体各个视角的三维信息。考虑高分辨率实时3D建模对精度的要求，本文采用了FPP 中以高分辨率与精度见长的相移轮廓术（Phase-shifting Profilometry，PSP）。但是PSP 需要对时间上连续的多幅条纹图案进行分析，运动会造成每帧条纹图案间的相位误差，导致重建出的三维图像出现波纹，所以选择三步相移法，即对物体投影3 幅条纹图案，尽可能降低运动对三维成像的影响，此时相机拍摄到的条纹图像强度可表示为：

其中：(u，v)为相机中的像素点坐标，A(u，v)为背景光强，B(u，v)为调制度光强，ϕ(u，v)为当前像素点的相位，n为相移指数（n=1，2，3）。

通过3 幅相移图像，利用最小二乘法计算出物体的相位ϕ(u，v)：

由于式（2）解出的相位是截断的，即具有2π相位跳变，为了获得连续分布的相位信息，需要使用相位展开算法去除包裹相位中存在的不连续跳变，同时考虑运动敏感性，本文选择的方案是基于多视觉的立体相位展开算法（Stereo Phase Unwrapping，SPU），其原理如图1 所示。

图1 立体相位展开原理Fig.1 Principle of stereo phase unwrapping method

投影仪依次投出3 幅条纹图In(u，v)并由3 个黑白相机捕获。此时对于相机1（Camera 1，C1）中任意坐标像素，其绝对相位都有N种可能，每一个都可以与根据投影仪与C1间的坐标参数重建出一个3D 候选点。将这些3D 候选点投影到相机2（Camera 2，C2）上可以得到N个2D 候选点，通过相位相似性检验，寻找这些候选点中包裹相位值与最为接近的候选点作为C1的对应点。但在噪声和系统误差的影响下可能导致误匹配，即错误的可能更接近且投影图案的条纹频率越高越容易发生误匹配现象。为了提高相位匹配的准确度，利用C2坐标轴对N个候选点进行一次筛选并将相似度最高的2 个候选点投影到相机3（Camera 3，C3）中，并重复C1与C2的约束操作。配合自适应深度约束策略（Adaptive Depth Constraint，ADC），筛选出最佳匹配点及其对应的相位信息，通过投影仪与C1的标定关系，可以求出3D 点的坐标p(x，y，z)，经过当前帧所有点（s个）的合集可表示为当前帧点云M{p1，p2，p3，...，ps}。

2.2 实时点云配准方法的基本框架

通过上一节中提到的方法获取物体高精度的三维形貌信息，对此可以生成用于后续实时模型重建的信息流，该信息流包含物体三维点云信息、二维深度信息以及二维彩色图像。所提出的基于全局优化的实时高精度模型重建方法的整体框架如图2 所示，该方法主要分为粗配准线程与精配准线程两部分，通过两个线程之间的有序配合实现实时高精度点云配准的效果。

图2 所提出的实时高精度模型重建方法整体框架Fig.2 Overall framework of the proposed real-time high-precision model reconstruction method

首先，粗配准线程计算流程如下：

Step1.实时获取被测物体的二维深度图像与二维彩色图像，分别计算他们的二维关键特征点。对相邻帧作二维关键特征点配准，同时将二维关键特征点映射到三维坐标，利用EPnP 方法计算相邻帧点云间的变换矩阵。

Step2.对相邻帧点云间变换矩阵的范数设置阈值来筛选信息流并获取关键帧点云。计算关键帧点云法线信息，生成FPFH 特征，并以FPFH特征为依据对相邻关键帧帧点云间的变换矩阵进行优化，筛选矫正部分误差较大的变换矩阵。

Step3.以第一帧点云为基准，根据计算得到的变换矩阵计算各关键点云的位姿矩阵。最后将优化筛选后的有效关键帧点云及其位姿发送给精配准线程。

Step4.重复Step1～Step3，直至完成对被测物体的全方位扫描重建，结束粗配准线程。

在接收到粗配准线程发送的点云与位姿数据后，精配准随之进行相应工作。精配准线程计算流程如下：

Step1.将经过配准后的点云集合设为源点云，将最新接收的点云设为目标点云，利用ICP算法，在已初始化的位姿基础上进行迭代，对两者间的位姿作进一步优化，并将优化后的位姿矩阵同步到各点云数据中。

Step2.以优化后的关键帧点云的位姿为顶点，相邻关键帧点云变换关系为边，添加并参与构建超图。

Step3.重复执行Step1～Step2，直至点云粗配准线程结束。在最后一帧关键帧点云位姿添加进超图后，将最终的点云配准结果更新到系统显示界面。

2.3 粗配准位姿初始化方案

本文所提出的点云位姿初始化方法结合了PnP 算法与FGR 算法。其中，PnP 算法是求解3D 到2D 点对应关系的方法，因其突出的速度优势被广泛运用于SLAM 领域。利用此方法可以根据已知的3D 空间点及其三维坐标估计出对应相机的位姿。如果可以建立两个相机拍摄到的图像中关键点的对应关系且计算出对应的3D 坐标，便可以计算出两个相机间的位姿关系。考虑到效率与精度的平衡问题，本文使用了EPnP 方法。作为复杂度为O(n)的非迭代PnP 问题的解决方案，EPnP 方法配合尺度不变特征变换（Scale Invariant Feature Transform，SIFT）特征点匹配算法可以得到大量的特征点并能够快速计算出相邻帧点云间的初始变换矩阵。然而，在全局扫描过程中会出现由于抖动等因素导致的特征点匹配错误的问题，仅通过EPnP 方法对相邻帧点云位姿进行初始化并不能稳定输出满足ICP 要求的点云初始位姿。因此在筛选出关键帧点云后，需要使用FGR 算法对关键帧点云位姿进行二次优化，确保后续精配准中ICP 成对配准的精度并避免ICP 迭代陷入局部极小值点的问题发生。

其中：E(T，L)为邻接关键帧点云对应点间距离函数，计算出函数最小值即优化完成；T为目标变换矩阵；L为对应集合H的线性过程；lm，m'为L中的元素；Ψ(lm，m')为先验误差，受参数μ控制；参数μ控制残差对目标有显著影响的范围，随着μ的减小，匹配问题的目标函数变得更加清晰，配准更加精确。

为了求得邻接关键帧点云间的最佳变换矩阵T，需要求解目标E(T，L)的最小值。为了使E(T，L)最小化，式（3）对每个lm，m'的偏导必须为0：

将lm，m'带入式（3）中即可求解得到邻接关键帧点云间的最佳变换矩阵T。

2.4 全局点云精细配准

获得关键帧点云初始化位姿后，需要对关键帧点云位姿作精细配准以得到高精度的点云模型。原方案以最新关键帧点云作为目标点云，以一系列已配准点云整体作为源点云，使用ICP 算法进行配准，配准得到的点云整体又将作为下一次配准的源点云。这是ICP 算法用于多点云配准的常见方案，但该方法并不能对多次配准产生的累计误差做出有效处理。随着关键帧点云数量的增多，每次ICP 配准时的细小误差发生累计。当系统再次扫描到已生成点云模型的区域后，叠加后的累计误差已不可忽视，并以点云重影的形式表现出来。

本文提出了相邻关键帧点云ICP 配准与全局点云优化［28］相结合的点云精配准方案，继承原方案利用ICP 算法对点云位姿的两两配准处理原则，但是仅保留点云配准后的关键帧点云位姿而不保存配准后的整体点云，并将配准后的点云及其位姿加入超图中进行图优化处理。基于全局图优化算法的具体流程如图3 所示。将经过精细配准后的关键帧点云位姿被依次加入超图中，至所有点云成对配准结束后，通过对邻接与次邻接点云间的位姿变换关系进行迭代计算，最终得到整体收敛的全局点云位姿。

图3 图优化算法的流程示意图Fig.3 Schematic diagram of graph optimization algorithm

经过ICP配准后的点云位姿已趋于真实位姿，邻接关键帧点云间已有了优秀的匹配精度。此时，将优化后的点云位姿设为顶点（Vertex），全部点云位姿（k个）构成顶点集V={v0，v1，v2，...，vk}。结合实际配准中关键帧点云间的点云重合面积，使用ICP 算法计算邻接关键帧点云间的变换矩阵Ti，i+1=t(vi，vi+1)，vi∈V与次级邻接关键帧点云间的变换矩阵Ti，i+2=t(vi，vi+2)，vi∈V，将它们设为连接顶点间的二元边（Edge），全部变换矩阵构成边集E={T0，1，T1，2，...，Tk-1，k，将点集与二元边集组合起来便构成了二元超图G={V，E}，成功将全局点云位姿优化问题以图的形式表现出来，表达式如式（8）所示：

其中：xi为关键帧点云的真实位姿，zi≈h(xi)为关键帧点云位姿的计算值，ei=zi-h(xi)为两者间的误差，Ωi为预设的信息矩阵，代表对ei中每个元素对应的误差项相关性的一个预计。

对每条边的误差ei求目标函数的优化解，设初始点为给它一个增量Δx，于是第i条边的目标函数可以表示为：

对目标函数进行二阶展开可以得到：

其中：Ci为常数，bi为目标函数一阶展开项系数，Hi为目标函数二阶展开项系数。为了找到Δx使目标函数到达极小值，只需令目标函数求关于增量Δx的导数为0，得到一个线性方程组：

通过求解线性方程组，可以得到单次迭代的Δx最优解。对所有边进行反复迭代，即可得到全局点云的最优点云位姿。

3 实验与分析

3.1 系统环境配置

为了验证所提点云配准方案的有效性，构建了一个多视角条纹投影三维成像系统，包括1 个LightCrafter 4500Pro 投影仪（投影帧率为100 Hz，分辨率为912×1 140）、3个黑白相机（Basler acA640-750 μm，分辨率为640×480）以及1 个彩色相机（Basler acA640-750uc，分辨率为640×480）。投影仪用于投影48 周期的相移条纹图案，3 个黑白相机用于完成基于立体相位展开三维重构功能，彩色相机用于获取物体表面纹理信息。我们在Visual Studio 2017平台调用OpenCV、PCL、QT 为主的第三方库进行系统的软件编写，并在计算机（内存为32 GB，主频为英特尔i7-11700k 2.40 GHz，显卡为NVIDIA 3080 10 G）上进行测试。

3.2 粗配准阶段点云位姿初始化效果对比

对大卫模型进行扫描，沿逆时针方向水平转动模型。在实时点云配准过程中改变大卫转速并在大卫旋转接近一周时，对大卫的转动施加一个剧烈抖动，观察此时点云配准显示窗口中的实时更新的大卫三维点云数据并进行记录，记录的大卫实时点云如图4 所示。为更加清晰具体地表现本文所提方法的优势，以Qian所提方法为对照［27］，分别针对粗配准阶段点云位姿初始化效果及全局点云优化效果设计对比实验。

图4 点云位姿初始化效果对比图。（a）现有方法位姿初始化效果；（b）本文所提方法位姿初始化效果；（c）为（a）对应的局部放大图；（d）为（b）对应的局部放大图。Fig.4 Comparison of initialization effects of point cloud position and pose.（a）Pose initialization effect of existing methods；（b）Pose initialization effect of the method proposed in this paper；（c）Corresponding partial enlargement of（a）；（d）Corresponding partial enlargement of（b）.

图4 中，大卫的旋转速度与箭头颜色深度成正相关。从系统实时显示窗口更新的点云图像数据可以看出，现有方法中点云位姿初始化受大卫转动状态影响较大，当在大卫转动过程中施加一个剧烈抖动（Shake）时，该方法会出现点云配准错误。相比之下，本文提出的点云位姿初始化方案具有更好的鲁棒性。在模型转动过程中出现剧烈抖动和速度突变的极端情况下，本文所提方法能够稳定地为精配准线程提供准确的点云初始位姿以进行后续的优化操作。

3.3 点云全局优化效果对比

对大卫模型进行扫描，沿逆时针方向水平转动模型并保持较低且稳定的转速，确保在对大卫模型进行全局扫描后能够得到较多的关键帧点云数量以放大累计误差对全局点云位姿的影响。配准结束后，将大卫点云模型导入MeshLab 中，将点云进行10 倍降采样后，对使用两种方法获取的点云进行对比。现有方法得到的点云规模为2、296、930，降采样后的点云规模为175、262，本文所提方法得到的点云规模为1、655、354，降采样后的点云规模为168、594。两种算法所得到的实验结果与局部细节放大对比如图5 所示。从图5可以看出，现有方法的点云位姿初始化方案虽然能够获取大卫的完整三维模型，但对该方法重构的点云模型进行10 倍降采样后，该方法并不能有效地处理连续点云配准产生的累计误差问题。而本文所提方法能够有效地消除累计误差，实现全局一致、无重影的点云配准，即使在10 倍降采样后依然能保证重建的大卫点云模型拥有连续且平滑的表面。

图5 不同点云全局配准的结果对比。（a）现有方法导出的完整点云模型；（b）对图5（a）进行10 倍降采样的结果；（c）本文所提方法导出的完整点云模型；（d）对图5（c）进行10 倍降采样的结果；（e），（g），（i），（k）图5（b）、（d）对应的局部放大图；（f），（h），（j），（l）图5（e）、（g）、（i）、（k）对应的三角化显示。Fig.5 Comparison of point cloud global registration effect.（a）Complete point cloud model derived by the method proposed in this paper；（b）10 times downsampling result of（a）；（c）Complete point cloud model derived by existing method；（d）10 times downsampling result of Fig.5（c）；（e），（g），（i），（k）Corresponding partial enlarged pictures of Fig.5（b），（d）；（f），（h），（j），（l）Triangular display corresponding to Fig.5（e），（g），（i），（k）.

3.4 实时高精度模型重建结果

分别对一个大卫模型和一个海盗模型进行扫描，以任意速度沿逆时针方向水平转动模型，在系统显示窗口中能够实时显示模型在当前位姿的点云重构过程，最终实现对待测物体的全方位三维建模。将模型扫描一周后，结束粗配准线程，待精配准线程对模型点云全局优化结束后，模型点云的最终配准结果将更新到系统点云配准的显示窗口中。在拼接过程中，分别对相邻帧点云实时配准耗时与相邻关键帧精细配准耗时进行统计，粗配准部分数据如表1 所示，精细配准部分数据如表2 所示。

表1 相邻帧配准算法耗时Tab.1 Time-consuming of adjacent frame registration algorithm ms

表2 相邻关键帧精细配准耗时Tab.2 Time-consuming of adjacent key frame registration s

对于大卫模型，单次粗配准平均耗时约为0.6 s，单次精细配准平均耗时约为2.5 s，总模型重建耗时为73 s；对于海盗模型，单次粗配准平均耗时约为0.5 s，单次精细配准平均耗时约为2.7 s，总模型重建耗时为79 s。模型的三维点云最终配准结果如图6 所示，其中图6（a）、（c）分别为对大卫模型和海盗模型离散点云的显示结果，图6（b）、（d）为对大卫与海盗模型点云作三角剖面化处理后的显示结果。

图6 点云配准结果展示。（a）大卫模型完整点云模型；（b）图6（a）对应的三角化显示；（c）海盗模型完整点云模型；（d）图6（c）对应的三角化显示。Fig.6 Display of point cloud registration results.（a）Complete point cloud model of the David model；（b）Corresponding triangulation display of Fig.6（a）；（c）Complete point cloud model of the pirate model；（d）Corresponding triangulation display of Fig.6（c）.

3.5 完整三维模型重建精度实验

为了验证本文所提方法的精度，我们对经过江苏省计量局校准过的标准陶瓷球进行实验分析，陶瓷球的半径分别为25.396 7 mm 和25.398 9 mm，球心距为100.116 6 mm。实验结果如图7 所示。由于精度球表面没有足够的特征点，因此将一个具有复杂面型结构的海盗石膏像放置于陶瓷球后，两者位置关系如图7（a）（正视图）和图7（c）（俯视图）所示。图7（b）和图7（d）为待测物体旋转一周后的三维重建结果。在Geomagic 中对陶瓷球作球面拟合，得到的误差分布如图7（e）和图7（f）所示。完整陶瓷球的重构半径分别为25.366 3 mm 和25.417 7 mm，误差分别为31 μm和18.8 μm；测得的球心距为100.200 0 mm，误差为83.4 μm。根据实验结果可知，与现有方法相比，本文所提方法重建的完整三维模型的整体精度提升约40 μm。

图7 精度球与海盗模型的重建结果分析。（a）、（c）精度球与海盗模型间的位置关系；（b）、（d）旋转一周后精度球与海盗模型的重建结果；（e）、（f）球面拟合后的误差分布。Fig.7 Analysis of reconstruction results of precision ball and pirate model.（a），（c）Positional relationship between the precision ball and the pirate model；（b），（d）Reconstruction results of precision ball and pirate model after one rotation；（e），（f）Error distribution after spherical fitting.

4 结 论

本文提出了一种基于全局优化的实时高精度模型重建方法，利用基于条纹投影技术的立体相位匹配方法实时获取高精度物体三维点云信息，改进原有实时360°点云配准方案，利用二维关键特征点与法向量特征对三维点云进行快速、准确的位姿初始化。此外，引入图优化技术对整体点云作全局配准处理，有效克服了点云重影问题，高效稳定地获取被测物体的三维点云模型。该方法能够在系统中高效且稳定地获取被测物体的三维点云模型，运行稳定，拼接精度高，实时全方位模型重构误差约为84 μm，且多次实验测量与理论相符。本工作进一步探索了非接触式测量在工业中的应用，实现了兼顾鲁棒性和高精度的三维模型重建，未来有望应用于缺陷检测、逆向工程等工业领域。