基于特征提取的点云自动配准优化研究

杨高朝

YANG Gaozhao

国家测绘地理信息局 海南基础地理信息中心，海口 570203

Hainan Geomatics Center,National Administration of Surveying,Mapping and Geoinformation,Haikou 570203,China

1 引言

地面三维激光扫描技术在测绘领域的应用越来越广泛[1-4]，打破了传统测量数据获取与处理模式，推动了点云数据处理理论的发展。由于受到测量设备和环境的限制，物体表面完整数据的获得往往需要通过多站扫描才能完成。因此，为了得到物体完整的点云数据，需要对这些局部点云数据进行坐标转换，最终拼接到同一个坐标系里。商业中应用最广泛最成熟的方法就是通过在重建物体上贴标志点[5]来实现点集的配准，但有些文物古迹的恢复和逆向工程存在无法贴标签甚至无法接触的问题，因此研究点云自动配准算法不但具备重要的实用价值，还具有重要的商业价值。

目前，许多国内外学者提出了多种点云配准算法。到20世纪90年代中后期比较流行的基于几何特征的匹配算法，如几何Hush[6]算法、Spine image等，其中最常见的就是基于迭代最近点（Iterative Closest Point，ICP）的方法[7]。但ICP算法可能存在局部最优缺陷，如果两站点云没有较高的重叠度或没有很好的初姿，很可能陷入局部最优，最后得到错误的匹配结果。Blais等对点云进行随机采样以提高速度，算法精度受到影响。Okatani提出了通过考虑平移和旋转错位的特性来获得最佳矩阵的方法，该方法易陷入局部收敛。朱延娟等[8]提出了基于点云曲率和邻域曲率相似度的粗配准。张广鹏等[9]提出了主元素贴合法的全局配准。刘艳丰[10]提出了基于四点算法的全局配准。

本文在总结分析国内外学者专家理论研究的基础上，针对上述点云数据配准中所存在的问题，从研究初始匹配入手，结合改进的ICP算法，设计出一种利用点云法向量夹角变化度及特征曲率相结合的自动配准方法。

2 自动配准的主要流程

点云自动配准一般分两个过程：粗配准、精配准[11-14]。粗配准是通过特征提取[15]、信号滤波等手段，来获得不同视角数据集间相对应的控制点，从而估计出配准变换参数。精确配准采用迭代方式逐渐逼近最佳结果，由于迭代方式不能保证得到最佳收敛结果，精确配准需要粗配准为其提供一个迭代的起始位置，来增大收敛到最佳结果的概率。具体步骤如图1。

图1 点云配准流程图

首先P、Q两站点云有一定的重叠度[16]（一般特征比较明显的点云超过60%即可得到正确的配准结果），根据每个点云的法矢量大小重采样点云，为每个点云定义一个N维向量，向量值为点云与邻近点云的法向量夹角。然后根据这些法向量夹角进行对应点搜索，进行粗配准，再进行精匹配，求取转换参数。具体详细步骤见下文。

2.1 粗配准

为了缩小点云之间的旋转和平移错位，使得精确配准不致得到错误的结果，需要进行点云粗配准。在研究上述算法的基础上，提出了法向量邻域夹角搜索的粗配准算法。首先对每个点构造一个不变特征向量，即每个点的法向量与其K-近邻点法向量[17]的夹角组成的向量，以此向量作为输入量构建k[θ1,θ2,…,θk]T，其中θk为弧度值，在两点集中搜索最近点。法向量夹角如图2、图3。

图2 模型点云法向量夹角

图3 目标点云法向量夹角

在计算法向量夹角的过程中，法向量的估计是比较重要的一步。由于采集到的原始数据往往很庞大，为了加快速度，先对原始数据进行点云抽希[18]，然后采用基于局部表面拟合的方法进行法向量估计[19]（为了使配准速度更快，本文选择平面拟合方案，但曲面拟合可能效果更优，但费时较多）。为此，对于点云中的每个扫描点pi，搜索与其最相近的K个相邻点，然后计算这些点最小二乘意义上的局部平面P，此平面可以表述如下：

式中，n为平面P的法向量；d为P到坐标原点的距离。

计算该点与最近点的法向量夹角余弦值，如图2和图3所示。当法向量夹角值越大时，余弦值就越小，说明该点所在区域起伏变化较大，特征较为明显，反之则较为平坦。因为平坦部分的点的几何特征较为接近，不适合作为特征点集进行配准，所以选取适当的阈值ε=0.99（经过实验求得），保留法向量夹角余弦值小于ε的点。

设待配准的两个点云数据分别为P和Q，利用上述特征点的提取方法，分别对两个点云进行特征点提取，得到P的特征点集为Q的特征点集为，其中m和n分别为P和Q的特征点的个数。

2.2 获取匹配点对

2.2.1 对应点对的初步获取

几何特征信息种类越多，计算就越复杂，但是特征信息量少，特征识别度低，就很难收敛。经过多种特征点云的实验，最终决定选取两种基本几何特征搜索匹配点[20]。下面详细介绍点云几何特征的建立过程。

（1）对于点集Pt中的每个点pti，以该点邻域法向量夹角加权值k[θ1,θ2,…,θk]T的特征距离作为该点的特征量：

（2）对于点集Pt中的每个点pti，利用邻域点集中k个最近点，求出邻域的重心O(pti)，以该点到其邻域重心的距离值作为该点的特征量：

对于点集Pt中的每个点pti，为其在点集Qt中搜索对应点qtj，若qtj是pti的对应点或近似对应点，则对两点来说，以上提出的两种几何特征应该相同或近似相同。因此对于点集Pt中的一点pti和点集Qt中的一点qtj，只要满足以下4个条件则认为它们具有对应关系：

阈值ε1、ε2取值过大，会导致产生很多错误的匹配点对，取值过小，无法获取足够的匹配点对，根据其他学者的经验以及本人反复试验，式中取ε1=ε2=0.01，使用以上关系进行过滤并初步建立点对应关系。由于点云中本身就存在很多特征相似区域，使用多个特征关系进行对应点的建立可以很好地避免出现大量一对多的情形，再通过合理阈值的选择可以为点集Pt中的每个点pti初步找到其在Qt上的对应点。最后建立初始匹配点对数组并记为：

式中，N为匹配点对的数量。

2.2.2 错误对应点对的剔除

采用上述方法确立的对应点集中会存在错误对应点对（或称为不可靠、不兼容点对），或者会产生多对互相矛盾的点对。为此，必须引入评价标准或约束准则去除错误对应点对，这是许多ICP改进算法的关键，也是国内外诸多学者研究的重点。区别于深度图像可利用图像自身的灰度、颜色、纹理等辅助信息[21-23]，对于点云，能够利用的信息主要是点云自身的曲面几何特征。

如果以上获取的初始匹配点对都是正确匹配点对，依据几何刚性变换原理可知，对于任意两对匹配点对和应满足点对间距离的偏差在极小的允许范围之内，即可表示选取阈值ε=0.02，对每个点对计算P中除其以外与其符合刚性距离约束条件的点的数目Qi，若P中一个点对满足：

2.2.3 匹配点对搜索精度检验

为了检验匹配点搜索精度，引入度量ε5：

2.2.4 计算初始匹配参数

在求得目标点集P和Q匹配点对后，需要计算两点云之间的平移变换矩阵，比较典型的方法有最小二乘法、四元数法[24]、矩阵的奇异值分解法（SVD）[25]等。经过对比发现，最小二乘法求解收敛速度更快。根据旋转矩阵和平移向量对目标点云进行初步转换。

2.3 精配准

2.3.1 改进后的ICP算法

进行初始配准后，两站点云的初始姿态基本上接近，可以作为ICP算法的初值进行精匹配。

根据经典的ICP算法，本文在以下方面对ICP算法进行改进，以提高算法的效率：

（1）使用基于法向夹角的点云重采样[26]；

（2）应用法向投影选取对应点，该搜索过程中使用kd-tree或八叉树可以提高速度；

（3）权重设置选取距离的平方加权平均值；

（4）误匹配点剔除选取距离和法向夹角相兼容的原则；

（5）目标函数选取点到面的误差最小。

法向投影-邻域搜索将源数据投影到目标点云上，然后在目标点云上进行最近点搜索（点周围的几何特征信息也作为判断的标准如曲率等）。

误差准则函数选择：

（1）如图4，对应点对直接距离的平方和最小[27]，这一形式是最常用的误差准则。

其中，pi是目标点云的坐标向量；qki是模型点云的坐标向量；R是变换矩阵；t是平移矩阵；ε是对应点对直接距离的平方和。当ε达到极值最小时，此时的R和t就是要求的值。

图4 最近点搜索算法

（2）如图5，每个源数据点到以其对应点法向量建立的局部平面的距离平方和最小。

其中，ni是pi的法向量。其他变量的含义同上。表示R⋅qki+t投影到pi点的切平面的距离的平方。采用此目标函数，只要两个点在与切平面平行的方向上滑动，距离都不变，因此在优化过程中允许滑动，这是此目标函数获得较快收敛速度的直观原因。

图5 法向量垂直搜索算法

根据图6实验结果，可知第二种方案无论是收敛速度上还是精度上都优于第一种方案，本文选取第二种方案作为误差准则函数。

图6 精配准误差图（比较点到点与点到面的差别）

2.3.2 精匹配实现步骤

运用改进的ICP算法进行精细配准步骤如下：

步骤1构造参与算法的两个初始点集，分别记作P和X；设定阈值Δ＞0，把点云的特征加入到配准策略上来，比如点云的曲率、法向量夹角等。

步骤2初始化p0=p,k=0。

步骤4计算Rk和Tk,使

最小，ni是目标点云的法向量。

步骤5对p0进行三维变换，得到pk+1=Rkp0+Tk。

步骤6k=k+1,如果k=1,转到步骤3；如果k≥2,且则终止迭代并转到步骤7，否则转到步骤3。

步骤7最后所得到的Rk、Tk即为所要求的精细配准参数，利用所得的参数，对第一片点云进行转换，配准结束。

2.3.3 误差评定

为了评价点云自动匹配的精度，由文献[28]可得：

其中，mp为配准中误差；pi、p′i为公共点对转换前后向量值。

3 实验及结果分析

在Windows平台下利用MFC和OpenGL软件实现了算法和配准界面，如图7。采用MFC单文档结构，利用分割视图法，管理3个视图。第二个视图显示model（模型）点集，第三个是data（目标）点集，然后点击菜单打开配准窗口，可在第一个视图同时显示模型点集、目标点集和配准后的点集。提供手动配准，自动配准，Mesh（或点云）旋转平移、点云删除、重选、缩放以及点云数据简化等功能。

图7 基于MFC和OpenGL实现的配准软件

3.1 实验结果

为了验证算法的有效性，本实验利用网上下载的不同视角的两片点云，点云个数29876个。首先对点云进行压缩，保留一些特征比较明显的点，压缩后的点云为7652个，如图8所示。接着对两站点云进行粗配准，结果如图9（a），最后再对点云进行精匹配，结果如图9（b），配准结束后还可以查看本次实现的配准精度以及配准所花费的时间。为了更能显示本文方案的可用性，选取一个现今比较成熟的方案，应用同样的数据同时实验，以便和本文方案进行比较。经过对比，选取了文献[20]配准算法对扫描点云进行迭代配准，其配准效果如图10所示。为更形象地查看迭代过程及结果，根据每次迭代的残差，结合MATLAB生成线状图，如图11。

图8 两组点云初始位置

图9 本文配准效果

图10 文献[20]配准效果

图11 文献[20]和本文算法迭代配准误差图

为了更直观地分析问题，对本次实验进行统计，如表1。

表1 本文算法和文献[20]算法的比较

3.2 结果分析

为了定量比较上述两组实验结果，通过计算多视点云重叠区域的最邻近点对距离平方和的平均值的平方根，作为整体配准后的距离中误差，用于评估多视点云间的配准精度，同样的数据运用式（12）求得本文算法配准后的中误差为30.4 μm,文献[20]配准后的中误差为38.7 μm。已知原始点云的密度为0.15 mm，实验表明，本文算法要优于文献[20]的算法。另外从实验结果还可以看出，本文算法配准中误差近似为0.3倍点云密度，因此本文配准算法的精度比较理想。

在算法迭代收敛方面，使用本文算法进行配准共迭代9次满足平差终止条件。其中3个平移参数的改正数收敛较快，3次迭代时已收敛；3个旋转角参数的改正数迭代后期收敛相对较慢，第5次迭代时3个旋转角参数的改正数绝对值都小于6×10-5，而本文中设定的限差为4×10-5，直到第9次满足限差要求，迭代结束。使用文献[20]算法同样对两站点云进行配准，一直迭代13次才满足迭代终止条件，从收敛效率上可以看出，本文算法要优于文献[20]的算法。

在算法运行速度方面，本文方案粗配准和精配准一共用了6.04 s，文献[20]一共用了8.09 s。由于本文方案尽管使用的是基于点云特征的配准方法，但并没有直接计算拟合曲面曲率，直接通过法向量夹角作为三维点云刚体不变的特征来实现自动配准。因此和传统利用点云拟合曲率实现点云自动配准的方案相比，本文方案提高了运行速度。

综上所述，本文算法无论是从配准精度上还是从迭代效率上都优于文献[20]。实验配准效果图和实验获得的对比充分说明了本文算法的有效性。

4 结论与展望

本文将利用邻域法向夹角搜索进行的初始配准作为初始值，并利用改进的ICP算法对提供的初始值进行精配准。实验结果表明，改进后的ICP算法在对应点的搜索时间和配准精度等方面都优于经典的ICP算法，既解决了传统的ICP算法进行配准，算法陷入局部最小，从而得到错误的配准结果的问题，又提高了点云配准的效率。点云中本身就存在很多特征相似区域，单一的特征匹配识别度低，在实验中发现，很多特征较少的点云也会陷入局部收敛，但采用多个几何特征来搜索对应点所花费的时间也会相应增加。因此，如何在两者间找到一个平衡点并在实际生产中广泛应用是接下来需要解决的问题。