复杂零件点云的高精度鲁棒配准方法

王素琴， 姚奕伸， 石 敏， 朱登明

(1.华北电力大学 控制与计算机工程学院 北京 102206；2.中国科学院 计算技术研究所 前瞻研究实验室 北京 100190；3.太仓中科信息技术研究院 江苏 太仓 215400)

0 引言

机械制造业的不断发展，对复杂零件的测量技术提出了更高的要求。传统的测量方法多采用人工或三坐标测量仪等精密检测设备，前者需要耗费大量的人力，后者成本高且测量效率低。同时，零件某些部位空间狭小,使得检验探头无法触及内部，从而难以准确测量。基于机器视觉的三维重建技术为复杂零件的非接触式测量提供了新的思路[1]。

基于机器视觉的三维重建技术首先利用三维扫描设备重建零件多视角下的点云，而后通过点云配准、融合得到完整的零件点云，再利用表面重建技术得到零件的网格模型，最后将网格模型与标准CAD模型进行误差分析以实现零件的尺寸测量。在这个过程中，点云的精确配准是点云融合的关键，是获取零件完整点云的先决条件。点云配准首先要找到两个点云中相匹配的若干点对，依据这些点对来求解点云之间的变换矩阵，然后将不同坐标系下的两个点云变换到同一坐标系下实现点云的融合。其中点对匹配的准确性和变换矩阵计算的鲁棒性对点云配准的精度起着关键的作用。

点云配准分为粗配准和精配准两种方法，粗配准方法根据预先确定的匹配点对集一次完成变换矩阵的计算，其配准结果依赖于匹配点对的准确性。以ICP算法及其变体[2-5]为代表的精配准算法以最近邻点作为匹配点对，在多次迭代中逐步完成变换矩阵的计算。但易受初始位置和噪声点的影响，往往会陷入局部最优解而无法达到最优配准结果。在高精度配准中，往往将两种方法联合使用，由粗配准为精配准提供良好的初值。

面向零件测量的点云配准精度要求高，因此采用粗配准和精配准结合的方式来完成零件点云的配准。其中，粗配准中的匹配点对的准确选取是实现高精度配准的先决条件，它决定了粗配准的结果，进而影响精配准的精度和鲁棒性。现有的粗配准算法[6-9]一般依据点云特征描述子来确定匹配点对，由于零件结构的特点，其表面往往存在很多特征相似的点，在匹配点对时易产生误匹配，导致最终配准结果精度较低，无法满足工业测量对点云配准的精度要求。

针对以上问题，本文提出一种面向复杂零件精密测量要求的，基于旋转平移不变性的精确几何结构的点云粗配准方法。主要创新点如下。

1)利用两点及其法线形成的夹角及两点间距离构成四元组来描述点之间的几何结构关系，该关系不因点云的旋转平移而改变。

2)为解决法线之间夹角无法描述法线相对位置的问题，提出一种基于指示向量的奇异性消除方法。

3)为解决因四元组中变量单位不同而难以度量相似性的问题，提出一种基于归一化的相似性度量方法来表征四元组之间的相似性。

1 相关工作

1.1 点云粗配准方法

Aiger等提出了4PCS算法[10]，该算法通过在目标点云和源点云之间寻找具有相同拓扑结构的四点对来计算两个点云之间的变换矩阵。该方法的时间复杂度比较高，不适用于大规模点云。Mellado等提出了Super 4PCS算法[11]，采用智能索引技术来解决4PCS算法效率低的问题。张琮毅等提出了尺度可变的快速全局点云配准方法[12]，该算法通过引入尺度因子，构造并优化能量函数来实现不同尺度点云之间的配准。Yang等使用截断最小二乘法描述配准问题,并使用一种图论框架来对尺度、旋转和平移进行级联求解[13]。Bhattacharya等利用特殊欧式群和迭代加权最小二乘法的鲁棒特性，对匹配点对集进行求解以完成点云配准[14]。陆军等提出了一种基于多尺度加权法向投影均值差提取关键点的点云配准算法[15]。Yang等提出了一种全局优化的Go-ICP方法[16]，将分支定界的思想与ICP算法相结合，在一定程度上降低了局部最小值对ICP算法的影响。张顺利等提出了一种基于自适应邻域匹配的点云配准方法[17]，利用自适应邻域信息形成的三角网格，定义残差角进行局部曲面描述，从而确定匹配点对，计算点云之间的变换矩阵。盛敏等提出一种基于曲率和颜色特征度的点云配准方法[18]，首先描述点云的曲率和颜色特征度，之后采用颜色或曲率相似性策略确定匹配点对，进而求解变换矩阵。Zhou等提出了一种全局配准算法FGR[7]，该方法依据点云特征描述子确定匹配点对集，通过构造能量函数，利用高斯牛顿法求解最佳的变换矩阵。

1.2 点云特征描述子

点云特征描述子是确定匹配点对的依据。PFH[19]特征描述子通过描述点和邻域点之间任意两点的法线分布关系来表征点的几何特征，其维度为125。SHOT[20]特征描述子将点的邻域空间进行分块，分别统计每块内的法线夹角余弦值的分布情况，其维度为352。RoPS[21]特征描述子将点的邻域点依次投影到局部坐标系的三个坐标平面，对坐标平面上的投影点形成的分布矩阵提取中心矩和香农熵来计算该点的特征，其维度为135。快速点特征直方图(fast point feature histogram,FPFH)[6]是对PFH的改进，它仅考虑描述点与其邻域点之间的法线分布关系。假设点云中含有n个点，考虑每个点周围的k个点，则PFH的时间复杂度为O(nk2)，FPFH为O(nk)。而且FPFH维度低，仅为33，点对匹配速度快，有效缓解了由于点的数量多而导致的匹配效率低的问题。因此，本文选用FPFH特征描述子来计算点云的特征向量。

1.3 点云精配准方法

在点云精配准方法中，最为著名的是由Besl等提出的ICP算法[2]，后续研究者又相继提出了多种ICP方法的变体。Granger等将归一化高斯权重加权与ICP方法结合[22]。Bouaziz等提出了Sparse-ICP方法[23]，通过使用稀疏诱导范数来解决ICP算法对缺失值和异常值敏感的问题。Pavlov等利用安德森加速思想来提高ICP算法的迭代速度[24]。Rusinkiewicz等为对称点云设计了新的目标函数[25]，简化了ICP算法，提高了计算效率。Zhang等提出了FRICP算法[26]，该算法使用Welsch函数来抑制噪声点对ICP算法的干扰，同时引入安德森加速器来加速ICP算法，目前该算法在相同点云数据集上相比其他精配准算法精度更高。尽管各种ICP变体在不同层面上对原始ICP算法进行改进，但仍未完全解决ICP算法对初始值依赖、容易陷入局部最优解的问题，它必须在有良好初值的情况下才能够实现高精度的配准。

2 算法概述

点云配准算法可以描述为：给定目标点云P和源点云Q，寻找源点云Q到目标点云P的最佳变换矩阵(R,t)，使得P和Q的对应点对之间的均方根误差最小，如公式(1)所示，

(1)

其中：R为旋转矩阵；t为平移矩阵；pi和qi是一组匹配点对。

求解变化矩阵首先需要确定点云之间的匹配点对，为了确定该匹配点对，本文首先利用点云特征描述子计算点云的特征向量，根据特征相似性原则确定两个点云中的初始匹配点对集。由于零件表面存在很多特征相似的点，因此根据特征相似性原则确定的初始点对集中包含较多的误匹配点对。这些误匹配点对会干扰后续变换矩阵的求解，严重影响点云配准的精度。为了解决零件点云误匹配点对过多的问题，本文受文献[27]启发，利用具有旋转平移不变性的精确几何结构特征对初始匹配点对集进行筛选，剔除误匹配点对。最后利用列文伯格-马夸尔特(L-M)算法对筛选后的匹配点对集进行非线性优化以计算两个点云之间的变换矩阵，完成点云粗配准。最后利用精配准算法FRICP[26]实现零件点云的精确配准。算法流程图如图1所示。

图1 算法流程图

3 初始匹配点对集的确定

首先计算点云的FPFH特征描述子，流程如下。

3)按照步骤2)的方法计算点Ps的所有邻域点的SPFH值。

4)根据步骤2)、3)，计算点Ps的FPFH值，如公式(2)所示，

(2)

其中：wk为权重，一般以点Ps到其邻域点Pt的距离作为权重。

5)为点云P及点云Q中的所有点计算FPFH值，将点云P和Q的FPFH特征记为F(P)={F(p)|p∈P}，F(Q)={F(q)|q∈Q}。

依据FPFH特征描述子并基于特征相似性原则、利用k-d树来确定初始匹配点对集。k-d树是一种用于高维空间向量搜索的数据结构。利用k-d树分别对目标点云和源点云的特征向量建立索引。对于点云P中的每一点p，其特征向量为F(p)，在点云Q特征向量的k-d树中查找F(p)的最近邻特征向量F(q)；对于点云Q中的每一点q，其特征向量为F(q)，在点云P的特征向量的k-d树中查找距离F(q)最近的特征向量F(p)。当F(p)和F(q)互为最近邻特征向量时，认为p和q是一组匹配点对，并将其加入初始匹配点对集C1。

4 精确匹配点对集的确定

为解决零件点云误匹配点对较多的问题，本文利用具有旋转平移不变性的精确几何结构特征对匹配点对集进行筛选。从初始点对集C1中随机挑选两组匹配点对(p1,q1)、(p2,q2)。由于零件为刚性物体，其点云在空间变换过程中几何结构保持不变。因此如果这两组点对均为正确的匹配点对，那么点p1、p2在点云P中所形成的几何结构特征与点q1、q2在点云Q中形成的几何结构特征应该相同。

图3 四元组示意图

然而夹角θn具有奇异性，即仅通过θn不能唯一地表示n1和n2的空间位置关系。如图4所示，n1和n2所处的空间位置不同，但其夹角均为θn。

图4 θn的奇异性问题

(3)

将点p1、p2的四元组表示为(θp1,θp2,θpn,dp)，点q1、q2的四元组表示为(θq1,θq2,θqn,dq)。在四元组的相似性计算时，由于角度和距离的量纲不同，不能直接使用基于欧氏距离的向量相似性度量方法。为了解决这个问题，本文提出一种归一化相似性度量方法，归一化是将有量纲公式转化为无量纲公式的一种方法，如公式(4)所示，

(4)

其中：s表示两个四元组的相似度，越接近0，表示点p1、p2和点q1、q2的几何结构越相似。当s满足式(5)时，即s低于阈值τ时，认为点p1、p2和点q1、q2的几何结构特征相同，即匹配点对(p1,q1)和(p2,q2)为正确的匹配点对，将其加入匹配点对集C2，

s≤τ，

(5)

其中：τ为给定的阈值，τ越接近0，约束条件越严格，筛选出的匹配点对越准确。但是由于数据采集视角的原因，目标点云和源点云的点分布情况不完全相同，因此过低的τ可能导致没有可用的点对。

5 变换矩阵求解

基于匹配点对集C2，采用非线性优化方法求解点云的变换矩阵T。常用的非线性优化方法有最速下降法、牛顿法、高斯牛顿法等。最速下降法又称一阶梯度下降法，它沿着梯度的反方向前进，从而使得目标函数收敛，但由于梯度下降法的贪心策略，经常会导致收敛速度过慢。牛顿法又称为二阶梯度法，它需要计算迭代点的海塞矩阵H，不适用于数据量大的点云。高斯牛顿法使用一阶导数矩阵与其转置的乘积JTJ来近似牛顿法中的H矩阵，解决了H矩阵计算复杂的问题。高斯牛顿法的增量方程可以表示为JTJΔx=-JTf(x)。由于零件点云中复杂的点分布情况，增量方程的参数JTJ经常会出现奇异矩阵的情况，导致增量方程无法求解，算法无法收敛。L-M算法通过引入一个非奇异矩阵来解决高斯牛顿法奇异矩阵的问题，从而可以提供更稳定、更准确的增量。

本文使用L-M算法求解点云变换矩阵，首先基于匹配点对集建立目标函数及残差函数，

(6)

其中：T是点云P和Q之间的欧氏变换矩阵；p和q是一组匹配点对。

与高斯牛顿法相似，L-M算法也需要求解增量方程，

(JTJ+μI)ξ=-JTf(T)，

(7)

其中：J为残差函数f(T)关于T的导数矩阵；I是为解决JTJ可能为奇异矩阵而引入的单位矩阵；μ为阻尼系数；ξ=(α,β,γ,a,b,c)，(α,β,γ)为L-M算法每次迭代的欧拉角，(a,b,c)为每次迭代的平移向量。由于L-M算法的每次迭代的旋转角度较小，所以可以近似认为sinα～α，sinβ～β，sinγ～γ。因此可以根据公式(8)实现变换矩阵T的计算，

(8)

T′是上一次迭代中计算出的变换矩阵。

在公式(7)中，当μ较大时，μI起主导作用，此时L-M算法接近最速下降法；当μ较小时，JTJ起主导作用，此时L-M算法接近高斯牛顿法。因此μ的取值可以在迭代过程中根据实际情况动态调整。为了确定阻尼系数μ，本文设置一个衡量指标ρ来指引μ的取值，ρ计算为

(9)

当ρ接近1时，表明函数实际的下降值与模型预估的下降值相差不大；当ρ较大时，表明函数实际的下降值比预估的大，在迭代过程中可以增大μ的值；当ρ较小时，说明函数实际的下降值比预估的小，在迭代过程中可以减小μ的值。

6 数据集构建与实验结果分析

6.1 数据集构建

数据采集系统由工业相机、投影仪和旋转台等设备组成。零件置于旋转台上，旋转台单次旋转60°，共采集零件6个视角下的点云数据。

本文的算法采用C++实现，在搭配Intel i7-9700F处理器、16 GB内存的计算机上进行实验。

6.2 点云配准评价指标

为了评估点云配准结果的优劣，本文采用均方根误差(RMSE)作为评价指标。RMSE反映了点云P和Q配准后点对之间的欧氏平均距离，RMSE值越小表示配准精度越高，其计算方法为

(10)

其中：pi和qi是一组对应点对。

6.3 实验结果及分析

6.3.1点云配准实验 为了验证本文所提出的误匹配点对筛选方法的有效性，将本文算法与现有的配准算法FGR、FRICP进行对比实验，实验结果如图5所示。其中红色点云为目标点云P，绿色点云为源点云Q，零件的初始位置如图5(a)所示。图5(b)～(d)分别展示了FGR算法、FRICP算法和本文算法的配准结果，可以看出对于零件1和2，FGR算法配准结果均不够理想，其RMSE值仅为1.122 0 mm和1.394 4 mm，主要原因是没有进行误匹配点对的剔除。精配准算法FRICP在零件1上配准结果精度较高，在零件2上配准结果不够理想。说明FRICP算法鲁棒性不高，配准效果对初始值依赖较大，配准精度不稳定。这是由于在初始位置不好的情况下，FRICP算法不能保证配准结果收敛到最优结果。本文算法由于实现了对匹配点对集的筛选，因此对零件1和2均可以实现良好的配准效果。在粗配准方法的对比实验中，相比FGR算法，本文方法对零件1配准的RMSE值降低了88%，对零件2配准的RMSE值降低了81%。详细数据如表1所示。

图5 零件点云粗配准结果对比

表1 不同方法配准结果的RMSE值

本文所提出的配准算法为粗配准算法，与精配准算法结合可以实现更高的配准精度。本文在FGR算法和本文算法的基础上进行FRICP精配准实验。实验结果如图6所示，图6(a)是点云的初始位置，图6(b)、(c)分别展示了FGR+FRICP和本文算法+FRICP的配准结果。可以看到FGR+ICP算法对于零件1配准精度较高，RMSE值达到了0.023 1 mm。对于零件2的配准效果较差，RMSE值为0.052 4 mm。而本文算法与FRICP算法结合，对零件1和2均可以实现更好的配准精度，其RMSE值分别达到了0.023 1 mm和0.007 0 mm。在搭配FRICP的精配准方法的对比实验中，与FGR+FRICP方法相比，本文方法+FRICP对零件1配准的RMSE值相同，对零件2配准的RMSE值降低了86%。详细数据如表1所示。实验结果验证了本文提出的粗配准算法对于复杂零件点云配准的有效性，相比于FGR算法可以提供更好的粗配准结果，搭配精配准算法的配准精度能够满足工业测量的精度要求。

图6 零件点云精配准结果对比

6.3.2四元组有效性实验 本文使用四元组(θ1,θ2,θn,d)来描述两点之间的几何结构特征。为了验证四元组中每一变量对配准结果的有效性，本文进行了四组对比实验，分别是：1)剔除距离d的三元组的对比实验；2)剔除角度θ1的三元组的对比实验；3)剔除角度θ2的三元组的对比实验；4)剔除角度θn的三元组的对比实验，如图7所示。其中图7(a)为使用四元组的点云配准结果，图7(b)～(e)依次为上述四组对比实验结果。

从图7(b)可以看出，在四元组中剔除变量d对配准的结果影响较大。而剔除变量θ1、θ2、θn后配准的结果相差无几，但均无法达到图7(a)所示结果的精度。

图7 四元组有效性对比实验

图8分别展示了初始匹配点对集和精确匹配点对集的效果图，其中红色的线表示错误的匹配点对，绿色的线表示正确的匹配点对集。可以看到仅依靠点云特征和特征相似性原则会产生较多的误匹配点对，而使用具有旋转平移不变性的精确几何结构方法对匹配点对集进行筛选后的精确匹配点对集中包含较少的误匹配点对。

图8 四元组筛选点对效果图

7 结语

本文针对零件表面有很多相似特征的点、容易产生误匹配点对的问题，提出了一种基于点云旋转平移不变性的点对筛选方法。针对两法线夹角的奇异性问题，提出了一种基于指示向量的奇异性消除方法。针对四元组提出了一种基于归一化的相似性度量方法。实验结果表明，本文方法与其他方法相比具有更高的配准精度，在为精配准方法提供初始值方面表现出较高的鲁棒性。同时结合点云精配准算法可以实现更高的精度，能够满足零件测量的要求。考虑工业领域零件测量通常要求实时性，因此，本文方法还需要在算法效率上采用基于openMP、CUDA等并行化方法进行优化。