赵夫群

(1.咸阳师范学院 教育科学学院，陕西 咸阳 712000; 2.西北大学 信息科学与技术学院，陕西 西安710127)

0 引 言

点云配准就是为两个点集求解最优空间变换，使其能够达到最佳配准的过程。由于多个点云的配准可以通过两个点云的配准来实现，因此通常所说的点云配准都是指两个点云的配准。点云配准已经在3D重建、医学研究、目标识别以及颅面复原等领域[1～6]得到了广泛的应用。

目前，应用最为广泛的点云配准算法是由Besl P J等人[7]提出的迭代最近点 (iterative closest point,ICP) 算法，该算法简单，易于理解和实现，但对两个点云的初始位置要求较高，且要求两个点云具有包含关系。因此，国内外学者提出了很多改进的ICP算法: Zhu J等人[8]提出了一种基于双向距离的仿射ICP算法，提高了算法的配准精度和抗噪性；Han J等人[9]提出了一种基于分层搜索方案、预警机制和启发式逃生方案的加强ICP算法，能够快速准确地实现点云配准； Du S等人[10]提出了概率ICP算法，提高了算法的收敛速度和抗噪性； Du S等人[11]提出了尺度迭代最近点(scaling iterative closest point,SICP)算法，解决了点云的尺度配准问题； Zhao L等人[12]提出了皱纹—感知点云配准方法，实现了大尺度的点云配准；周文振等人[13]提出了一种基于K—means聚类的改进ICP算法，实现了精确的室内地图配准。但上述算法未解决算法收敛区间和精度之间的矛盾。

为此，本文提出一种由粗到细的点云配准算法。在粗配准阶段，采用尺度参数可变的点云配准算法，可以解决点云配准中收敛区间与配准精度之间的矛盾。在细配准阶段，通过设置旋转角约束和动态迭代系数来改进ICP算法，可以解决由旋转角变化过大引起的配准效果不佳的问题，并提高算法的迭代收敛速度。

1 粗配准

在粗配准阶段，首先采用大尺度参数测度函数作为配准测度函数，并逐渐减小尺度参数，使得算法在结束时能够获得足够精确的配准参数。

1.1 测度函数

E(Tr)=∑σ(d(Tr(mi),dj)),

i=1,…,NM,j=1,…,ND

(1)

式中mi和dj分别为点云M和D中的三维坐标点，Tr(mi)为对点mi进行空间变换，d(Tr(mi),dj)为Tr(mi)和dj之间的欧氏距离。函数σ定义为:当t=0时，σ(t)=1;否则，σ(t)=0。

对式 (1)求最大值max(E(Tr))即可实现点云配准。为了式 (1) 寻优求解，采用高斯核函数对其进行平滑，并采用平滑后的函数作为点云配准的测度函数。该高斯核函数的定义为ρh(t)=exp(-t2/s2),s为尺度参数。那么平滑后的计数函数为

i=1,…,NM，j=1,…,ND

(2)

式(2)即为点云配准的测度函数。只要尺度参数s选取合适，平滑后的计数函数Eh(Tr)即为原计数函数E(Tr)的近似。本文通过BFGS(Broyden-Fletcher-Goldfarb Shanno)拟牛顿优化算法搜索其全局极大值，从而确定点云配准的最优参数，实现点云粗配准。

1.2 变尺度配准算法

主要思想为：在算法的起始阶段，采用大尺度参数测度函数作为点云配准的测度函数，以免算法陷入局部极小；随着优化过程的推进，逐渐减小尺度参数，可在算法结束时获得精确的配准参数。具体步骤为：

1)设置初值：尺度参数s0，迭代次数k=1，初始刚体变换q0=(R0,t0)，R0为旋转矩阵，t0为平移矢量；

2)利用BFGS进行目标函数的极大值搜索；

3)通过目标函数的极大值位置获得刚体配准的参数Rk和tk；

4)若满足终止条件则令刚体变换的最终旋转矩阵为R*=Rk，平移矢量为t*=tk，算法结束;否则,转步骤(5)；

(3)

6)执行k=k+1操作，转步骤(2)。

2 细配准

细配准采用改进的ICP算法实现，在ICP算法中加入角度约束和动态迭代系数，以解决由旋转角变化过大引起的配准效果不佳的问题，并提高算法的收敛速度。

2.1 ICP算法

s.t.RTR=In,det(R)=1

(4)

式中R为旋转矩阵;t为平移矩阵。

ICP算法基本步骤如下：

1)根据第k步的已知刚体变换Rk和tk，将点云M进行Rkmi+tk变换，并建立两个点云M和D间的相关性ck+1(i)，其数学描述为

i=1,…,NM

(5)

2)计算点集M和D的刚体变换，其数学描述为

(6)

2.2 改进的ICP算法

ICP算法的改进过程分2个步骤:加入旋转角约束和加入动态迭代系数。

定义旋转矩阵R为R=RxRyRz且

(7)

式中θx,θy,θz为旋转角，令θxb,θyb,θzb为旋转角均值；Δθx,Δθy,Δθz为旋转角偏差；θxb-Δθx,θyb-Δθy,θzb-Δθz为旋转角的下界；θxb+Δθx,θyb+Δθy,θzb+Δθz为旋转角上界。

加入旋转角约束后，M和D的配准问题可描述为

(8)

在算法中加入动态迭代系数h。通常，动态迭代系数h取大于等于0的整数，当h不同时，算法的收敛速度也不同，根据不同h值对应的迭代收敛曲线可知,随着动态迭代系数h的增大，算法的收敛速度越来越快，但当h增大到一定的程度，收敛曲线会出现振荡，当h=4时，算法不再收敛。因此，要为h设定取值范围，根据实验经验，建议其取值范围在1～4之间。

在ICP算法中加入动态迭代系数h的步骤如下：

1)计算刚体变换矢量q=[R|t]T以及qk的相邻2次迭代的变化量Δqk；

2)由刚体变换矢量Δqk更新基本ICP算法中的Mk共h次，即执行Mk=Δqk(Mk)共h次。

加入旋转角约束和动态迭代系数后，该改进ICP算法的具体步骤如下：

1)给定刚体变换初值q0=[R0,t0]T，R0为初始旋转矩阵，t0为初始平移矢量，令M0=R0M+t0，迭代次数k=0，动态迭代系数h=0；

2)估计旋转角θx,θy,θz的边界，即θx∈[θxb-Δθx,θxb+Δθx],θy∈[θyb-Δθy,θyb+Δθy]和θz∈[θzb-Δθz,θzb+Δθz]；

3)利用式(5)建立点云M和D相关性ck(i)；

4)利用奇异值分解(singular value decomposition,SVD)的方法计算旋转矩阵Rk+1和平移矢量tk+1，则qk+1=[Rk+1,tk+1]T;

5)计算qk+1的相邻两次迭代的变化量Δqk+1;

6)利用式(8)计算Mk+1=Rk+1M+tk+1;

7)判断均方根(root mean square,RMS)误差，若RMSk+1-RMSk>ε，则执行h=h+1操作;否则执行h=0操作RMS定义为

(9)

8)判断动态迭代系数h，若h>0，则通过执行Mk+1=Δqk+1(Mk+1)共h次来更新点集Mk+1；

9)若满足|RMSk+1-RMSk|<ε或k>Stepmax则算法终止;否则,转到步骤(10)，ε为预设的阈值，Stepmax为最大迭代次数;

10)令k=k+1，并转步骤(2)。

3 实验结果与分析

点云配准的实验数据源于Stanford 3D Scanning Repository，如图1所示。

图1 待配准的点云

基于本文配准算法，首先采用变尺度点云配准算法实现粗配准，然后分别采用ICP 算法和提出的改进ICP 算法实现细配准、配准结果如图2所示，算法的匹配点数、匹配误差、迭代次数和耗时等运行参数如表1所示。

图2 配准结果

点云类型点云大小算法配准点数配准误差迭代次数耗时/s兔子20128,20048粗配准ICP改进ICP5566671297940.05960.02890.01353222103.22.11.1龙20920,17418粗配准ICP改进ICP5017611389070.06580.02990.01653524133.62.41.3

在采用改进ICP算法的细配准过程中，旋转角偏差取值为Δθx=Δθy=Δθz=10°，动态迭代系数取值为h=3。

从图3的配准结果来看，变尺度的点云粗配准算法可以将两个点云初步对齐；ICP算法和改进ICP算法可将两个点云进一步对齐，但改进ICP算法比ICP算法的配准误差更小，具有更高的配准精度。

从表1可见，改进ICP算法比ICP算法具有更多的配准点数、更低的配准误差、更少的迭代次数和耗时。与ICP算法相比，改进ICP算法配准精度和速度均提高了约50 %，是一种精度更高、速度更快的点云配准算法。因此，提出的由粗到细的点云配准算法是一种有效的点云配准算法。

4 结 论

本文由粗到细的点云配准算法解决了算法中收敛区间与配准精度之间的矛盾和旋转角变化过大引起的配准效果不佳的问题，是一种精度高、速度快的点云配准算法。该算法具有一定的抗噪性，但未考虑尺度各异的点云配准问题。在今后的研究中，还要综合考虑点云配准中更多的影响因素，如噪声因素、尺度因素、外点的干扰因素等，设计更加合理有效的配准算法，进一步提高点云配准的精度、速度、鲁棒性和抗噪性，扩大点云配准技术的应用领域。