一种基于立体栅格自动分析3D点云配准误差的方法

摘要

本发明公开了一种基于立体栅格自动分析3D点云配准误差的方法，包括如下步骤：1）获得源点云；2）获得变换矩阵；3）对源点云进行下采样；4）对源点云进行滤波；5）获得目标点云；6）提取特征点；7）计算特征点的特征向量；8）寻找源点云和目标点云的对应点；9）计算对应点之间的平均距离；10）当寻找的对应点数目减少到阈值下，则进入11，否则修改寻找对应点的控制变量，返回8；11）加权平均即得到激光雷达3D点云配准的误差。本发明可以定量比较不同配准算法构建3D场景模型的准确性，帮助配准算法在实际应用中确定合适的点云规模，可用在对配准的精度要求较高的激光雷达三维测绘和结构测量领域。

说明书

技术领域

本发明涉及一种3D点云的数据处理方法，特别涉及一种基于立体栅格自动分析3D点云配准误差的方法，可用在对配准的精度要求较高的激光雷达三维测绘和结构测量领域。

背景技术

三维激光扫描技术可快速高效的获取被测对象表面的三维坐标，目前被广泛的应用于测绘工程，结构测量，文物保护等领域。但由于三维激光雷达扫描范围的局限性，有时需要多测站多角度扫描，或者多个激光雷达同时扫描。因此，需要通过3D点云数据融合和配准技术，来将多个点云配准到同一个坐标系下，获得物体表面完整的信息，实现对三维场景真实立体的再现。

针对点云配准，目前行业内的标准做法是在测量期间放置标志物。但是由于手动放置标志物该过程非常费时，并且当这些标志物不可避免地遮挡了一些对象和场景或者分布在边界不清的点云集合中，都会导致配准的失败。基于以上原因，越来越多的研究者投入激光雷达点云自动配准的研究中来。其中ICP算法无需标记，且产生了很多高效的变体，应用广泛。

但是由于激光雷达采样物体表面产生的点云常常是稀疏不均的，并且雷达采样的误差以及移动的物体，植物等反射表面产生的噪声都会造成配准精度下降。多站配准还会引起误差传递，更重要的是，针对不同的激光雷达采集的不同场景，配准算法的效果都有明显差异。在激光雷达三维测绘和结构测量领域，比如对空间位置的测量，器件面积和体积的测量，都要求配准算法具有高精度的特性。王忠武等(图像融合中配准误差的影响【J】，测绘科学，王忠武，赵忠明，2010(02))实例表明,配准误差对融合质量的影响较大,在遥感图像融合处理中,配准误差越小越好。因此，在配准算法变体繁多的情况下，针对不同的场景选择最合适的配准算法就必须研究配准算法的准确性。

除此之外，点云数目大小也会对配准算法产生影响：

其一是对算法效率的影响：点数越多，计算效率低。但是点数越少，并不一定能提高算法效率；

其二是对配准准确率的影响，点数少，则一些包含关键信息和特征的点可能被滤除，从而导致配准时出现错误对应点对增多的情况。

目前计算配准算法误差的方法大部分都只是在自身配准方法的基础上计算，因此不具有通用性和可比较性。

发明内容

在激光雷达三维测绘和结构测量领域，比如对空间位置的测量，器件面积和体积的测量，都要求配准算法具有高精度的特性，因此，为了在配准算法变体繁多的情况下，选择精度最高的算法构建3D场景，本发明提供一种基于立体栅格自动分析3D点云配准误差的方法，利用全局特征子寻找对应点，再加权平均得到最终的误差。该方法可适应稀疏密度不均的点云，并且能够自动确定邻域大小，保留足够多的特征信息，与理论误差的相关性也得到了验证；还可以帮助配准算法在实际应用中确定合适的点云规模；可以用于但不限于激光雷达三维测绘和结构测量领域。

本发明为解决上述技术问题采用以下技术方案：

本发明提供一种基于立体栅格自动分析3D点云配准误差的方法，包括如下步骤：

步骤1，通过激光雷达获得源点云；

步骤2，获得配准完成时的变换矩阵；

步骤3，对步骤1中的源点云进行下采样；

步骤4，对步骤3中下采样后的源点云进行滤波；

步骤5，将滤波后源点云中的每一点乘上变换矩阵，获得目标点云；

步骤6，计算步骤4中滤波后的源点云和步骤5中获得的目标点云中的每一点的曲率，并根据曲率的大小提取特征点，得到源点云特征点集合以及目标点云特征点集合；

步骤7，计算步骤6中提取的源点云和目标点云特征点的点特征直方图特征向量；

步骤8，寻找步骤6中提取的源点云特征点和目标点云特征点的对应点，基于随机采样一致性算法RANSAC对寻找到的所有对应点对集合进行筛选，得到对应点对集合，其中随机采样一致性算法RANSAC的拒绝阈值为激光雷达的分辨率；

步骤9，计算步骤8中对应点对集合中每对对应点两点之间的距离的平均值；

步骤10，若当前对应点对集合中的对应点对数少于第一次寻找得到的对应点对数的40％，则进入步骤11，否则将当前随机采样一致性算法RANSAC的拒绝阈值减去激光雷达分辨率的1/50作为新的拒绝阈值，返回步骤8；

步骤11，将各次迭代时步骤9得到的距离平均值进行加权平均，即得到最终激光雷达3D点云配准的误差。

作为本发明的进一步优化方案，步骤3中对源点云进行下采样，具体为：

3.1，构建一个能够包围源点云数据的最小的三维长方体空间，将该长方体空间均匀划分为若干个三维体素栅格；

3.2，对每个栅格内的每个点，利用利用kd-tree搜索其K近邻，计算邻域法向量，并求出该点和其邻域法向量之间的夹角，其中，K的取值为栅格内的点数；

3.3，对每个栅格，计算该栅格内所有点和其邻域法向量之间夹角的平均值，若平均值大于50°则该栅格的边长减小为原边长的一半，否则该栅格的边长增大为原边长的两倍；

3.4，根据3.3中的栅格边长，重新划分3.1中的长方体空间，重复3.2至3.3，直至所有栅格的边长保持不变；

3.5，选取每个栅格内的质心作为关键点来代替该栅格内的其他点，从而形成关键点集合，完成源点云的下采样。

作为本发明的进一步优化方案，步骤4中对下采样后的源点云进行滤波，具体为：

4.1，对每个关键点利用kd-tree搜索其K近邻，计算该关键点与其所有近邻点的距离的平均值，得到关键点集合对应的一个平均值集合；

4.2，计算4.1中平均值集合的标准差以及平均值，其中，该平均值集合的平均值即为全局距离平均值；

4.3，对任一关键点，若该关键点与其所有近邻点的距离的平均值超出全局距离平均值一倍标准差，则剔除该关键点，得到一个新的关键点集合，从而完成对下采样后的源点云的滤波。

作为本发明的进一步优化方案，步骤6中计算滤波后的源点云和目标点云中的每一点的曲率，并根据曲率的大小提取特征点，具体为：

6.1，分别计算滤波后的源点云中所有点曲率的平均值以及目标点云中所有点曲率的平均值；

6.2，选取滤波后的源点云中曲率大于平均值的点作为特征点，得到源点云特征点集合；

6.3，选取目标点云中曲率大于平均值的点作为特征点，得到目标点云特征点集合。

作为本发明的进一步优化方案，步骤8中寻找源点云特征点和目标点云特征点的对应点，具体为：

8.1，遍历源点云特征点集合，计算每个源点云特征点与目标点云特征点集合中每个点对应的特征向量的距离，并计算所有源点云特征点与目标点云特征点集合中每个点对应的特征向量的距离的平均值；

8.2，遍历目标点云特征点集合，计算每个目标点云特征点与源点云特征点集合中每个点对应的特征向量的距离，并计算所有目标点云特征点与源点云特征点集合中每个点对应的特征向量的距离的平均值；

8.3，若l_ij＜20％l_p且n_ji＜20％l_q，其中，l_ij表示目标点云特征点集合中第i个特征点m_i与源点云特征点集合中第j个特征点q_j对应的特征向量的距离，l_p表示所有目标点云特征点与源点云特征点集合中每个点对应的特征向量的距离的平均值，n_ji表示源点云特征点集合中第j个特征点q_j与目标点云特征点集合中第i个特征点m_i对应的特征向量的距离，l_q表示所有源点云特征点与目标点云特征点集合中每个点对应的特征向量的距离的平均值，i∈{1,2,…,MN}，MN表示目标点云特征点集合M中的特征点数，j∈{1,2,…,QN}，QN表示源点云特征点集合Q中的特征点数，则点m_i与q_j为对应点；

8.4，基于随机采样一致性算法RANSAC对8.3中得到的所有对应点对集合进行筛选，去除错误的对应关系，得到对应点对集合。

作为本发明的进一步优化方案，步骤11中将各次迭代时步骤9得到的距离平均值进行加权平均，即得到最终激光雷达3D点云配准的误差，具体为：

11.1，各次迭代时步骤9得到的距离平均值按照迭代的先后顺序进行倒序排序；

11.2，按照11.1的排序依次分配权值：1/2、1/4、1/8、……、1/2ⁿ，其中，n表示迭代次数；

11.3，根据11.2分配的权值对距离平均值进行加权平均，获得最终激光雷达3D点云配准的误差。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明提供的一种基于立体栅格自动分析3D点云配准误差的方法，采用可变大小的立体栅格进行下采样，保留了更多的特征信息；在立体栅格内进行滤波，提取特征点，计算点特征直方图特征向量实现了自动确定邻域大小，避免邻域过小不能反映局部特征，邻域过大细节特征扭曲而造成的对应点提取不准确；同时通过不断循环严格化寻找对应点的控制变量以及最后阶段加权平均都使该方法能更好的应用于实际环境中，帮助在激光雷达三维测绘和结构测量领域提高对物体建模的精度。

附图说明

图1为本发明的方法流程图；

图2为velodyne-16雷达在室内长廊采集的一帧源点云和目标点云；

图3为需要进行下采样的源点云；

图4为图3中的源点云用本发明进行下采样的结果；

图5为图4中的源点云用本发明滤除离群点的结果；

图6为稀疏点云对计算配准误差的影响(圆形点：采样点。三角形：雷达。虚线和实线表示雷达在不同位置。圆圈是雷达的扫描轨迹。由于稀疏点云点云采样点之间距离大，所以配准之后采样点本身存在距离，其中后面长方条的采样点仅表示随着物体离雷达中心距离越远，点云更加稀疏)；

图7为计算理论误差时构造的三种点云；

图8为理论误差与本发明计算出的误差对比结果；

图9为50帧点云理论误差与本发明计算出的误差的相关关系；

图10为本发明提出的方法分别计算出配准失败与成功的误差的对比图；

图11为针对velodyne-16雷达，ICPN与ICPF算法误差的比较；

图12为针对velodyne-16雷达，ICPF算法在点云数目不同时配准的误差；

图13为针对velodyne-16雷达，ICPF算法在点云数目不同时配准消耗的总时间。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

本发明一种基于立体栅格自动分析3D点云配准误差的方法，如图1所示，包括如下步骤：

步骤1，通过激光雷达获得源点云；

步骤2，获得配准完成时的变换矩阵；

步骤3，对步骤1中的源点云进行下采样；

步骤4，对步骤3中下采样后的源点云进行滤波；

步骤5，将滤波后源点云中的每一点乘上变换矩阵，获得目标点云；

步骤6，计算步骤4中滤波后的源点云和步骤5中获得的目标点云中的每一点的曲率，并根据曲率的大小提取特征点，得到源点云特征点集合以及目标点云特征点集合；

步骤7，计算步骤6中提取的源点云和目标点云特征点的点特征直方图特征向量；

步骤8，寻找步骤6中提取的源点云特征点和目标点云特征点的对应点，基于RANSAC(Random Sample Consensus,随机采样一致性算法)对寻找到的所有对应点对集合进行筛选，得到对应点对集合，其中随机采样一致性算法RANSAC的拒绝阈值为激光雷达的分辨率；

步骤9，计算步骤8中对应点对集合中每对对应点两点之间的距离的平均值；

步骤10，若当前对应点对集合中的对应点对数少于第一次寻找得到的对应点对数的40％，则进入步骤11，否则将当前随机采样一致性算法RANSAC的拒绝阈值减去激光雷达分辨率的1/50作为新的拒绝阈值，返回步骤8；

步骤11，将各次迭代时步骤9得到的距离平均值进行加权平均，即得到最终激光雷达3D点云配准的误差。

下面通过具体实施例对本发明的技术方案做进一步阐述：

本实例中，采用的是velodyne公司生产的最小型的16线3维激光雷达，雷达在不同的地点和不同的角度采集15帧点云，得到的其中一帧源点云和目标点云如图2所示。

上述步骤3中，由于激光雷达采样物体表面产生的点云常常是稀疏不均的：对距离而言，在离雷达中心越远的地方点云越稀疏，而离雷达中心较近的地方点云又较稠密。对角度而言，不同角度扫描得到的点云数据在各自重叠区域内点云密度存在差异，这都会造成算法在搜索邻域时的半径大小不具有通用性。因此，在尽可能多的保留点云特征的前提下，通过对点云进行下采样，减少密度的差异性。

如图3所示的源点云(方位角：90°、仰角：45°)进行下采样具体方法为：首先，求出能够包围输入的三维点云数据最小的三维长方体空间，将该长方体空间均匀划分为一个个三维体素栅格；然后，在每一个栅格内，利用kd-tree搜索每个点的K近邻(K的取值为栅格内的点数)，进行法向量估计。该点和其邻域法向量之间夹角的平均值可以表示该点处法向量的变化程度。变化程度大的点表示该区域起伏变化明显，所含信息更多，故应减少立体栅格的边长，这样有更多的栅格来保留较多的数据点，若变化程度小，则说明该区域较为平坦，故应增大立体栅格的边长，该区域内栅格数目减少仅保留少部分数据点。体素栅格尺寸的变化除了可以保留更多的特征信息外，对点云密度的平均也有明显效果。然后，在每一个栅格内选取质心来代替体素内其它点，在这里选取质心而不是更常用的体素中心，是因为质心可以更好的保留原始空间点的布局。最后，根据关键点所在栅格的大小更新每个关键点的K邻域大小，实现了在后续操作中自动确定邻域大小的目的，避免了邻域过小不能反映局部特征，邻域过大细节特征扭曲而造成的对应点提取不准确的现象。图4是对图3所示的源点云进行下采样的结果(方位角：90°、仰角：45°)。

激光雷达扫描得到的点云密度不均匀，雷达采样的误差以及移动的物体，植物等反射表面产生的噪声都会造成点云中存在一些稀疏的离群点。点云中这些点一般分布在雷达精度之外。例如：VLP雷达精度在3cm和100m之间。因为雷达近处不能放置物体，故实例中采集的点云会产生一些极端大值。这些离群点在计算点特征描述子时，会导致计算值异常，因此需要对点云进行离群点滤波。

上述步骤4)中对点云滤波的方法为：首先，对每一个关键点利用kd-tree搜索其K近邻；然后，计算出该关键点与其所有邻域点的距离和，再取平均得到平均距离；接着，计算出所有关键点与其邻域点的平均距离，计算这组数据的标准差和平均值，这组数据的平均值称为全局距离平均值；最后，当某一关键点与其邻域点的平均距离超出全局距离平均值1个标准差时，该关键点则被去除。图5为去除离群点的结果(方位角：0°、仰角：90°)。

上述步骤6中提取特征点的方法为：分别计算滤波后的源点云中所有点曲率的平均值以及目标点云中所有点曲率的平均值；选取滤波后的源点云中曲率大于平均值的点作为特征点，得到源点云特征点集合；选取目标点云中曲率大于平均值的点作为特征点，得到源点云特征点集合。

对于在点云集合中的每个点，首先，利用kd-tree搜索其K近邻，计算其邻域内K个点的质心，如式(1)所示，其中，p_s表示邻域内的第s个点的坐标，表示质心的坐标；然后，计算该点与K邻域的协方差矩阵和对应的特征值，如式(2)所示，其中，Cov表示该点与K邻域的协方差矩阵，表示与该协方差矩阵对应的第r个特征向量，λ_r表示与第r个特征向量对应的第r个特征值，由于该点的坐标为x、y、z三维坐标，故特征向量和特征值只有三个，即r的取值为r∈{0,1,2}；再后，根据式(3)计算该点对应的曲率；最后，求取所有点曲率的平均值，曲率大于平均值的点即选为关键点。

上述步骤8中寻找源点云特征点和目标点云特征点的对应点的方法为：已知目标点云特征点集合M和源点云特征点集合Q的点特征直方图特征向量，首先，针对集合M中的第i个特征点m_i，计算m_i与集合Q中第j个特征点q_j对应的特征向量的距离l_ij，i∈{1,2,…,MN}，MN表示目标点云特征点集合M中的特征点数，j∈{1,2,…,QN}，QN表示源点云特征点集合Q中的特征点数；然后，遍历集合M中的每一个点，得到一组距离值，计算这组距离值的平均值l_p；接着，针对集合Q中的第j个特征点q_j，计算q_j与集合M中第i个特征点m_i对应的特征向量的距离n_ji；然后，遍历集合Q中的每一个点，得到一组距离值，计算这组距离值的平均值l_q；最后，当l_ij＜20％l_p且n_ji＜20％l_q，就将该点对(m_i,q_j)记做对应点对。由于噪声的影响，不是前述所有的对应关系都是合理的，因此，采取基于RANSAC(Random>

上述步骤11得到最终激光雷达3D点云配准的误差的方法为：首先，将各次迭代时步骤9得到的距离平均值按照迭代的先后顺序进行倒序排序。因为随着拒绝阈值的严格化，寻找的对应点的准确性也不断提高。所以根据排序依次分配权值为：1/2、1/4、1/8、……、1/2ⁿ，其中，n表示迭代次数。最后，根据分配的权值对距离平均值进行加权平均，获得最终激光雷达3D点云配准的误差。

在立体栅格内进行滤波，提取特征点，计算点特征直方图特征向量实现了自动确定邻域大小，避免了邻域过小不能反映局部特征，邻域过大细节特征扭曲而造成的对应点提取不准确的现象。

以下是本实例得出的结果：

由于激光雷达采集到的两帧之间的点，本身就存在几厘米的距离，如图6所示，其中，圆形点为采样点，三角形为雷达，虚线和实线表示雷达在不同位置，圆圈是雷达的扫描轨迹。由于稀疏点云点云采样点之间距离大，所以配准之后采样点本身存在距离，其中，后面长方条的采样点仅表示随着物体离雷达中心距离越远，点云更加稀疏。同时，随着目标物体离雷达的距离越远，相差的距离也越大。已知VLP16雷达激光发射器每次旋转的角度为0.2°，在一个30m的长廊，相邻两帧采样点之间距离最大可达30*0.2/180*pi＝10.47cm(这只是最坏的假设，实际发生的概率非常小)。故实际算法的误差是比使用本发明提出的方法计算出的误差要小的，但是由于激光雷达点云稀疏是采样时不可避免的因素，故本发明提出的方法只能作为一个相对尺度，衡量算法之间的误差，不能作为一个绝对尺度衡量源点云和配准后的点云之间的距离。

为了比较配准算法之间的准确性，现将本发明的方法与其他两种现有方法进行比较。第一种配准方法为ICPF，是张建华等人提出的基于3D点云点特征直方图特征实时三维空间定位方法(申请号：CN201610659484.X)，将该方法改为用于两帧点云的配准，其中，特征描述部分用PCA方法代替ISS算法，并加入下采样和滤波操作。第二种配准方法是为ICPN，即将第一种方法的计算点特征直方图特征向量部分改为计算法向量特征。

为验证该方法的合理性，实例构造了50种理想矩阵(包含旋转和平移)和三种点云，如图7所示，将理想变换矩阵和实际变换矩阵(第一种配准方法得到的变换矩阵)相减，然后将相减得到的矩阵各个元素相加，最后结果即为理论误差。实际配准中，由于不能事先获得理想变换矩阵，故不能利用矩阵计算误差；所以用同样的数据利用本发明提出的方法测量误差，结果如图8所示。对两组数据进行相关分析，如图9所示，其相关系数达到了0.8226，属于显著相关，所以该方法可以衡量VLP16激光雷达实际应用中配准的误差。

本实施例发现该方法分析误差的结果可以分为三个等级衡量配准情况，相对尺度小于0.001，配准实际误差小于10^-5，属于配准结果理想的等级I。相对尺度小于0.25，配准实际误差小于0.01，属于配准结果一般的等级II。而配准失败时，相对尺度在0.6以上。从图10可以看出，配准失败利用该相对尺度和理想误差计算出的数值与配准成功的数据有明显的区分度。

实例采用VLP16雷达获得两帧连续的点云，两点云平移距离较小(5cm)，在ICPF和ICPN都配准成功的基础下，采用本发明提出的方法比较了两算法的准确率。在准确率方面，如图11所示，两算法基本相当。主要原因是在输入点云情况较好的情况下，算法准确率主要是由细配准算法来决定，而两配准算法的细配准方法都为ICP方法。

本发明除了可以针对激光雷达选择最优的配准算法构建3D场景模型外，还可以优化配准算法，寻找最合适的点云规模。

算法的时间复杂度是衡量算法效率的重要指标，点云的规模决定了算法中基本操作重复执行的次数。然而点云数目过少，计算特征描述子时易出错，从而导致配准失败；

针对第一种配准方法，实例将15组不同的源点云，分别采用网格化体素滤波器得到下采样后的3组点云，点数分别为10000点，7000点和5000点。比较它们与没有下采样的源点云(20000点)各自耗费的时间以及准确率。

其中与原点云的比较结果如图12所示，可以看出，点云进行下采样可以非常明显的提高配准算法的效率。然而，这种明显的影响效果在点云数目减少到7000点之后逐渐消失；原因在于点数减少时，特征描述子估算的时间上升。

点云规模除了影响算法的效率之外，对配准的准确性也有影响。实例采用本发明的方法分别计算了点数为10000点，7000点以及5000点点云配准的准确性。从图13可以看出点云数目为5000点，配准准确率最差。

因此，针对VLP16雷达实际应用中，采用7000点是综合效率和配准准确率的最佳方案。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。