CN105069840A - 一种基于曲率特征的三维正态分布变换点云配准方法 - Google Patents

Abstract

一种基于曲率特征的三维正态分布变换点云配准方法。该方法包括：根据设定的曲率阀值，分别提取模型点云和目标点云的曲率特征点；将模型点云曲率特征点的表示形式转换为正态分布组合形式；通过优化目标函数求取最优变换矩阵，使得目标点云中点的坐标值准确地变换到模型点云坐标系下。相比于已有的三维正态分布变换点云配准方法，该方法排除了大量冗余点的干扰，降低了参与配准过程的点集的数量，有效地缩短了配准算法的运行时间。

Description

一种基于曲率特征的三维正态分布变换点云配准方法

技术领域

本发明属于三维重建技术领域，特别是涉及两片带重叠区域的点云数据集之间实现配准的方法。

背景技术

由于待重建场景的复杂性、激光传感器扫描范围的局限性等原因，往往需要将一个大的场景分为多个区域分别扫描。而每次扫描的相对坐标系都是以激光扫描仪的发射点(视点)为原点的新建坐标系。为了使多次扫描的点云数据集统一到同一个坐标系下，需要求出每两次扫描的点云数据集之间的变换矩阵，求变换矩阵的过程就是点云配准的过程。三维重建系统中常用的配准算法主要分为三类：基于特征的配准算法^[1–4]，迭代最近点算法(ICP)^[5–7]和正态分布变换算法(NDT)^[8–11]。

基于特征的配准算法就是利用两片点云对应的特征点对或者对应的特征向量估计变换矩阵，使得一个点云数据集准确地配准到另一个点云数据集。由于此类算法只存在一次变换过程，对应特征估计的准确性要求非常高，因此在实际应用中，基于特征的配准算法均作为一种粗配准或者初始配准，如文献[3]中提取SIFT特征完成初始配准以及文献[4]中根据k近邻搜索点设计的一种不变特征。完成初始配准得到比较理想的初始位置之后，采用精确配准算法提高配准精度。常用的精确配准算法是Besl等人在1992年提出的迭代最近点算法(ICP)^[6]。ICP算法的明显特征是在计算过程中的迭代思想，本质上是一种计算两个点集间的距离并将其最小化的方法。ICP算法对初始位置要求严格、易于陷入局部最优，而且算法时间复杂度高^[8]。

正态分布变换配准算法于2003年由Biber和Strasser提出[9]，正态分布变换配准的主要思想是改变模型的表示形式，由原来离散点的表示形式转换为一种正态分布的组合形式，新的模型描述了某一个点分别属于各个表面的概率。

发明内容

本发明目的是解决配准过程耗时长的问题，提供一种基于曲率特征的三维正态分布变换点云配准方法。

本发明采用三维正态分布变换算法在很大程度上降低了配准算法的运行时间，而且降低了实验误差，提高了实验精度。但应用到一个场景范围大、点云密集程度高的场景中，配准过程耗时长的问题依然不容忽视。通过提取曲率特征，排除了大量冗余点的干扰，降低了参与配准过程的点的数量，有效地缩短了配准算法的运行时间。

本发明提供的一种基于曲率特征的三维正态分布变换点云配准方法，包括如下步骤：

第1、曲率特征提取，提取曲率特征的原则是：曲率种子点邻域半径范围内的点为曲率特征点，邻域半径计算过程为在提取曲率特征点之前，首先设定曲率阈值ε，曲率值大于阈值的点为曲率种子点，并确定邻域半径的计算公式：

R＝kH(1)其中H为该点曲率值，k为比例参数。k的取值根据实际情况进行调节，场景特征丰富时，同比例减小比例值；反之，增大比例值；

曲率特征提取过程为：

第1.1、分别计算模型点云中点的曲率值，计算公式如下：

$K=\frac{\left|y^{\′\′}\right|}{{(1+y^{\′2})}^{3/2}}---\left(2\right)$

其中，曲面函数为y＝f(x)，y′,y″分别为函数y关于x的一阶和二阶导数。

第1.2、遍历模型点云中所有点，若曲率值大于阈值ε，计算该点的邻域半径R，提取邻域半径R内的点为曲率特征点。

第1.3、参照第1.1步和第1.2步方法，提取目标点云的曲率特征点。

第2、坐标变换矩阵估计

采用精简过程为：首先将点云模型划分为体素栅格结构，然后分别计算每个体素栅格内包含数据点的重心，用重心点代替同一栅格内其他点，依次处理完所有体素栅格。然后应用三维正态分布变换算法求取坐标变换矩阵，完成配准，具体步骤如下：

第2.1、首先，将模型点云的曲率特征点均分为一系列小立方体，每个小立方体为一个单元，边长长度一般设定为场景宽度的1/20，并且保证小立方体内点数不少于10个。

第2.2、然后对于每一个单元，用公式(3)计算均值向量和公式(4)计算协方差矩阵，每个单元的点数要大于一个设定的阈值；

$q=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}{x}_k---\left(3\right)$

$C=\frac{1}{n-1}\underset{k=1}{\overset{n}{\Σ}}(x_k-q){(x_k-q)}^T---\left(4\right)$

其中q为均值向量，C为协方差矩阵，n为此单元内点的个数，x_k＝1…n是此单元里包含的点，(x_k-q)^T中的上标T表示对该向量进行转置运算。k表示点的序号，为累加符号，表示从k＝1的项，累加到k＝n的项。位置x处的点在单元b内的概率用正态分布N(q,c)表示，概率密度函数为：

$p\left(x\right)=\frac{1}{c}\exp (-\frac{{(x-q)}^T{C}^{-1}(x-q)}{2})---\left(5\right)$

其中，q是此单元的均值向量，C是此单元的协方差矩阵，C的上标-1代表对C进行矩阵求逆运算，c是归一化常数，exp表示自然指数函数。

第2.3、需要优化的参数是旋转变量和平移变量，把所有参数用一个向量p表示。定义一个变换函数T(p,x)，即将点x通过p代表的旋转和平移变换得到新的位置坐标，在二维空间中，T(p,x)为：

$T(p,x)=\left[\begin{array}{cc}cos\mathrm{\φ}& -sin\mathrm{\φ}\\ sin\mathrm{\φ}& \cos \mathrm{\φ}\end{array}\right]x+\left[\begin{array}{c}{t}_x\\ t_y\end{array}\right]---\left(6\right)$

其中t_x,t_y是平移变量，φ为旋转变量，因此p＝[t_x,t_y,φ]。

第2.4、假设目标点云的曲率特征点集为χ＝{x₁,...,x_n}，坐标变换参数向量为p。定义目标函数s(p)，求最优变换p使得χ变换之后在模型点云上的概率最大，具体公式为：

$s\left(p\right)=-\underset{k=1}{\overset{n}{\Σ}}p\left(T\right(p,x_k\left)\right)---\left(7\right)$

等价于求目标函数s(p)的最小值；应用迭代牛顿算法优化此函数，直到最后一次迭代的参数向量模长|p|小于0.0001。

第3、将目标点云集变换到模型点云坐标系

根据第2.3步求得的坐标变换函数T(p,x)，将目标点云数据集χ＝{x₁,...,x_n}变换到模型点云坐标系，完成两片点云的配准。

本发明的优点和有益效果：

本发明设计的改进方法与原三维正态分布变换方法在设定的初始旋转偏差范围内运行时间随着偏差的增大而上升，但本发明方法的运行时间相比于原方法表现出更加明显的优势。

附图说明：

图1为配准实验点云模型；

图2为初始平移偏差对本发明方法和文献[9]中方法的配准时间及最终误差的影响；

图3为初始旋转误差对本发明方法和文献[9]中方法的配准时间及最终误差的影响；

图4表示全景点云模型；

图5表示模型点云提取曲率特征结果图；

图6表示目标点云提取曲率特征结果图；

图7表示配准之前两片点云的位置关系图；

图8表示配准之后两片点云的位置关系图。

具体实施方式：

实例例1：

为了进一步分析比较原三维正态分布变换方法与本发明方法，设计一组评估实验：待配准的模型点云与目标点云来自同一个点云数据集(点云模型如附图1所示)，人为设置一定的初始平移偏差与初始旋转偏差，在此初始偏差的基础上分别应用三维正态分布变换方法与本发明方法进行配准实验，每组参数均进行10次实验并求取均值。

第一组实验：

1.1a、设定初始偏差

为了验证初始平移偏差对配准效果的影响，保证初始旋转为0rad，设定初始平移偏差为100～500mm，并将目标点云按照相应的变换矩阵完成变换。

1.2a、提取模型点云与目标点云的曲率特征值

设置曲率阈值0.5和邻域半径比例值500，根据公式(2)计算各个点的曲率值，根据公式(1)计算邻域半径的值，保证特征点数量大概为原点云数量的1/10，分别提取模型点云和目标点云的曲率特征点。

1.3a、应用三维正态分布变换方法完成两组曲率特征点之间的配准

将模型点云的曲率特征点在空间范围内划分为一系列小立方体，小立方体的边长为50cm，根据公式(3)、(4)计算小立方体内三维点的平均向量和协方差矩阵，根据公式(5)得到三维点在该立方体空间内的概率密度函数。根据公式(6)、(7)求变换矩阵使得目标点云变换到模型点云空间坐标系内概率最大。

1.4a、比较两种方法的运行时间，以及配准完成后最终平均偏差值

附图2展示了附图1中点云数据集实验结果的曲线图，当初始平移偏差由100mm增大到500mm时，原方法的运行时间由279.32ms增大到546.54ms。而在此过程中，本发明方法的运行时间一直低于原方法运行时间的一半，误差值基本保持为一个很小的定值。分析上述实验结果可得，本发明方法与原方法有一个共同特征，即初始偏差越大，方法运行时间越长，但是本发明方法的运行时间相比于原方法有明显地缩短。而且在设定的初始偏差范围内，两种方法最后均能收敛到非常小的误差值。

第二组实验：

1.1b、设定初始偏差

为了验证初始旋转偏差对配准效果的影响，保证初始平移为0mm，设定初始旋转偏差为0.1～0.5rad，并将目标点云按照相应的变换矩阵完成变换。

1.2b、提取模型点云与目标点云的曲率特征值

设置曲率阈值0.5和邻域半径比例值500，根据公式(2)计算各个点的曲率值，根据公式(1)计算邻域半径的值，保证特征点数量大概为原点云数量的1/10，分别提取模型点云和目标点云的曲率特征点。

1.3b、应用三维正态分布变换方法完成两组曲率特征点之间的配准

将模型点云的曲率特征点在空间范围内划分为一系列小立方体，，小立方体的边长为50cm，根据公式(3)、(4)计算小立方体内三维点的平均向量和协方差矩阵，根据公式(5)得到三维点在该立方体空间内的概率密度函数。根据公式(6)、(7)求变换矩阵使得目标点云变换到模型点云空间坐标系内概率最大。

1.4b、比较两种方法的运行时间，以及配准完成后最终平均偏差值

附图3展示了附图1中点云数据集实验结果的曲线图，当初始旋转偏差由0.1rad增大到0.5rad时，原方法的运行时间由363.84ms增大到896.07ms。而在此过程中，当初始旋转偏差比较小时，本发明方法的运行时间几乎是原方法的三分之一，本发明方法的优势更加明显，而且两种方法最终均收敛到非常小的误差值。

分析上述实验结果，本发明方法与原方法，在设定的初始旋转偏差范围内，运行时间随着偏差的增大而上升，但是本发明方法的运行时间相比于原方法，表现出更加明显的优势。

实施例2:

将本发明方法应用到大范围室内场景三维数据获取实验中。实验选用的激光测距传感器为HOKUYO公司的UTM-30LX2D激光扫描仪，硬件环境为CPU：Intelcorei5-23002.8GHZ；RAM：4.00GB；IDE：visualstudio2010。

2.1、确定扫描策略

在激光扫描仪扫描范围有限的条件下，为了覆盖整个室内场景，需要制定一个扫描策略，达到在有限次扫描之后得到整个场景的点云模型。

首先为了获取三维空间的数据，需要步进电机带动激光扫描仪在第三个维度上旋转同时完成在另外两个维度的扫描。根据室内场景大小以及每次扫描装置可以覆盖的范围，将室内数据的获取分5次进行。每次选取一个视点，当次扫描完成后，再选取另外一个视点进行扫描，直到扫描过程结束。

2.2、每两片点云之间配准拼接

得到5片点云模型之后，确定其中一片点云所在坐标系为参考坐标系，每两片点云之间应用本发明方法进行配准拼接，最终统一到同一个坐标系，实验结果如附图4所示。

附图5至附图8以两个点云模型为例具体展示了配准实验的实验结果。附图5、附图6分别为模型点云和目标点云提取曲率特征点的实验结果，图中白色点为提取的曲率特征点。曲率阈值设为0.2094，求邻域半径的比例因子大小为500。附图7中标记区域表明配准之前偏差较大，附图8中模型点云和目标点云准确地配准到一起。

参考文献

[1]张政,张彩明.一种基于法向特征的点云数据配准算法,中国图学新进展2007–第一届中国图学大会暨第十届华东六省一市工程图学学术年会论文集,2007.

[2]王蕊,李俊山,刘玲霞等.基于几何特征的点云配准算法,华东理工大学学报:自然科学版,35(5):768–773,2009.

[3]王程冬,程筱胜,崔海华等.SIFT算法在点云配准中的应用,传感器与微系统,31(2):149–152,2012.

[4]邢正全,邓喀中,薛继群.基于K-近邻搜索的点云初始配准,测绘科学,38(2):93–95,2013.

[5]G.C.Sharp,S.W.Lee,ICPregistrationusinginvariantfeatures,IEEETrans.onPatternAnalysisandMachineIntelligence,24(1):90–102,2002.

[6]P.J.BeslandN.D.McKay,Amethodforregistrationof3-dshapes,IEEETrans.onPatternAnalysisandMachineIntelligence,14(2):239–256,1992.

[7]S.Rusinkiewicz,M.Levoy,EfficientvariantsoftheICPalgorithm,inProceedingsofthe3rdInternationalConferenceon3-DDigitalImagingandModeling,2001:145–152.

[8]T.Zhiqiang,Z.Ke,Y.Chenglong,etal,ImprovedICPregistrationin3-Dmodelreconstruction,TransactionsoftheChinaWeldingInstitution,1:024,2013.

[9]M.Magnusson,A.Lilienthaland,T.Duckett,Scanregistrationforautonomousminingvehiclesusing3D-NDT,JournalofFieldRobotics,24(10):803–827,2007.

[10]M.Magnusson,A.Nuchter,C.Lorken,etal,Evaluationof3Dregistrationreliabilityandspeed:AcomparisonofICPandNDT,inIEEEInternationalConferenceonRoboticsandAutomation,2009:3907–3912.

[11]A.Das,s.L.Waslander,ScanregistrationusingsegmentedregiongrowingNDT,TheInternationalJournalofRoboticsResearch,33(13):1645–1663,2004.

[12]代星,崔汉国,胡怀宇.基于曲率特征的点云快速简化算法,计算机应用,29(11):3030–3032,2009.

Claims (1)

1.一种基于曲率特征的三维正态分布变换点云配准方法，其特征在于该方法包括：

第1、曲率特征提取，提取曲率特征的原则是：曲率种子点邻域半径范围内的点为曲率特征点，邻域半径计算过程为在提取曲率特征点之前，首先设定曲率阈值ε，曲率值大于阈值的点为曲率种子点，并确定邻域半径的计算公式：

R＝kH(1)

其中H为该点曲率值，k为比例参数，k的取值根据实际情况进行调节，场景特征丰富时，同比例减小比例值；反之，增大比例值；

曲率特征提取过程为：

第1.1、分别计算模型点云中点的曲率值，计算公式如下：

$K=\frac{\left|y^{\′\′}\right|}{{(1+y^{\′2})}^{3/2}}---\left(2\right)$

其中，曲面函数为y＝f(x)，y′,y″分别为函数y关于x的一阶和二阶导数；

第1.2、遍历模型点云中所有点，若曲率值大于阈值ε，计算该点的邻域半径R，提取邻域半径R内的点为曲率特征点；

第1.3、参照第1.1步和第1.2步方法，提取目标点云的曲率特征点；

第2、坐标变换矩阵估计

采用精简过程为：首先将点云模型划分为体素栅格结构，然后分别计算每个体素栅格内包含数据点的重心，用重心点代替同一栅格内其他点，依次处理完所有体素栅格；然后应用三维正态分布变换算法求取坐标变换矩阵，完成配准，具体步骤如下：

第2.1、首先，将模型点云的曲率特征点均分为一系列小立方体，每个小立方体为一个单元，边长长度一般设定为场景宽度的1/20，并且保证小立方体内点数不少于10个；

第2.2、然后对于每一个单元，用公式(3)计算均值向量和公式(4)计算协方差矩阵，每个单元的点数要大于一个设定的阈值；

$q=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}{x}_k---\left(3\right)$

$C=\frac{1}{n-1}\underset{k=1}{\overset{n}{\Σ}}(x_k-q){(x_k-q)}^T---\left(4\right)$

其中q为均值向量，C为协方差矩阵，n为此单元内点的个数，x_k＝1…n是此单元里包含的点，(x_k-q)^T中的上标T表示对该向量进行转置运算，k表示点的序号，为累加符号，表示从k＝1的项，累加到k＝n的项；位置x处的点在单元b内的概率用正态分布N(q,c)表示，概率密度函数为：

$p\left(x\right)=\frac{1}{c}\exp (-\frac{{(x-q)}^T{C}^{-1}(x-q)}{2})---\left(5\right)$

其中，q是此单元的均值向量，C是此单元的协方差矩阵，C的上标-1代表对C进行矩阵求逆运算，c是归一化常数，exp表示自然指数函数；

第2.3、需要优化的参数是旋转变量和平移变量，把所有参数用一个向量p表示；定义一个变换函数T(p,x)，即将点x通过p代表的旋转和平移变换得到新的位置坐标，在二维空间中，T(p,x)为：

$T(p,x)=\left[\begin{array}{cc}cos\mathrm{\φ}& -sin\mathrm{\φ}\\ \sin \mathrm{\φ}& \cos \mathrm{\φ}\end{array}\right]x+\left[\begin{array}{c}{t}_x\\ t_y\end{array}\right]---\left(6\right)$

其中t_x,t_y是平移变量，φ为旋转变量，因此p＝[t_x,t_y,φ]；

第2.4、假设目标点云的曲率特征点集为χ＝{x₁,...,x_n}，坐标变换参数向量为p；定义目标函数s(p)，求最优变换p使得χ变换之后在模型点云上的概率最大，具体公式为：

$s\left(p\right)=-\underset{k=1}{\overset{n}{\Σ}}p\left(T\right(p,x_k\left)\right)---\left(7\right)$

等价于求目标函数s(p)的最小值；应用迭代牛顿算法优化此函数，直到最后一次迭代的参数向量模长|p|小于0.0001；

第3、将目标点云集变换到模型点云坐标系

根据第2.3步求得的坐标变换函数T(p,x)，将目标点云数据集χ＝{x₁,...,x_n}变换到模型点云坐标系，完成两片点云的配准。