CN104240251A - 一种基于密度分析的多尺度点云噪声检测方法 - Google Patents

Abstract

本发明提供一种基于密度分析的多尺度点云噪声检测方法，所述方法主要包括两步：第一步利用多尺度的密度分析算法初步判断可能为噪声的点；第二步利用三角网约束将上一步中误检测为噪声的点重新归为非噪声点。利用本发明能够有效检测出来自激光雷达的点云和影像匹配得到的点云中包含的孤立噪声和簇状噪声,在激光雷达点云滤波、影像匹配中的误匹配点检测、光束法平差中的噪声检测等领域具有良好的应用前景。

Description

一种基于密度分析的多尺度点云噪声检测方法

技术领域

本发明属于测绘科学与技术领域，主要应用于激光雷达点云噪声的检测以及匹配点云噪声的检测，尤其是涉及一种基于密度分析的多尺度点云噪声检测方法。

背景技术

点云数据的处理一直是数字摄影测量与计算机视觉的重要研究内容。当前获取密集的三维点云主要有两种方式：1)采用LiDAR(Light Detection AndRanging)系统(及激光雷达系统)直接得到三维点云；2)利用影像匹配得到点云。但两种来源的点云数据都存在一定数量的噪声点。噪声的存在对点云数据的处理产生很多影响。如LiDAR点云滤波处理过程中，很多算法假设地面点为局部最低点，而在点云中存在局部低噪声的情况下，这种假设则是不适用的；另外噪声点对后续的大规模点云分层渲染、建筑物建模以及密集匹配所得同名点云的后续处理精度(如影像匹配、地物识别等)都会造成较大影响。因此，噪声点的检测和剔除是获取准确密集点云数据一个重要的处理环节。

从空间分布特性来说，噪声可归纳为两大特点：1)典型的高程异常的孤立点，表现为局部高噪声点和局部低噪声点；2)非典型的簇状噪声点(数量不多的噪声点的聚类)——在一定局部范围内不能明显高出邻近地面点和地物点的高程异常点，表现为低矮噪声，形态呈点簇状或块簇状。许多学者对点云噪声处理进行过研究，现有的点云去噪的算法可大致概括为以下三类：1)分析、利用高程值的分布进行去噪；2)利用数学形态学的方法进行去噪；3)根据点云密度进行去噪。方法1和方法2都是通过局部或是全局的高程值比较，并判断阈值来检测噪声，这类方法通过选取合适的阈值可以对孤立噪声有效检测但对簇状噪声的检测结果并不理想。方法3基于密度的方法是将噪声假设为位于低密度区域(相对于正常的点云密度)的单一物体或物体聚类，因此检测噪声的过程也就是检测低密度区域的过程。但现有的一些基于密度的去噪方法中，往往不能兼顾处理多种点云数据的不同类噪声。且通常密度分析的方法存在着处理簇状噪声失败的问题，主要是因为簇状噪声在局部小范围内其密度和分布表现都与非噪声点类似，因而需要找到一种既可以检测出孤立点，又可以有效检测出簇状噪声点的算法。

为此，本发明提出一种新的点云去噪算法，可适用于这两类数据中所包含的噪声点的去除。算法主要包括两步：第一步利用多尺度的密度算法去除孤立噪声和小的簇状噪声；第二步利用三角网约束将第一步中误检测为噪声的点重新归为正常点。

发明内容

本发明的目的是提供一种基于密度分析的多尺度点云噪声检测方法，能够有效检测机载激光雷达数据获取的点云与影像匹配获得的点云中包含的孤立噪声和簇状噪声，从而克服不同种类点云噪声类别差异大、表现各异难以检测的问题，该方法通过进行多尺度的密度分析点云，逐步的剔除点云数据中包含的噪声，最后通过构建三角网，并用三角网约束降低上一步误检测的噪声点。其技术方案如下所述：

一种基于密度分析的多尺度点云噪声检测方法，包括以下步骤：

(1)输入点云数据并设定尺度层级，不同尺度层级用level^l表示，其中l＝0～L，其中最小尺度层级为level0级，最大尺度层级为levelL级；

(2)采用空间六面体进行点云数据的三维分割，建立离散点云三维网格和最大范围包围盒，所述最大范围包围盒根据输入的点云三维坐标计算得到，通过局部密度分析方法标记噪声点为0并剔除，得到新的非噪声包围盒；

(3)在上一步得到的非噪声包围盒基础上，降低一级尺度层级，以1/2倍的三维网格大小逐个进行邻域的局部密度分析，进一步剔除较小的噪声点或噪声点聚类并标记为0，获得新的非噪声包围盒；迭代进行，直到完成最小尺度层级噪声点检测；

(4)完成尺度层级迭代后，检测标记为0的噪声是否为真正的噪声点，利用非噪声点构建狄罗妮三角网，检测标记为0的噪声点到对应三角面的距离，当距离小于给定的阈值时则认定为非噪声点，所述阈值设定为最精细层级的三维网格的在Z方向的尺寸。

进一步的，步骤(3)中，对特定尺度层级下，对噪声点的检测是以某一单元体为中心一定邻域范围形成空间范围做为一个单位，所述单元体代表每一个分割而成的三维网格，在一个指定的搜索空间内进行遍历，计算该邻域范围内的密度与搜索范围内的密度的关系，同时判断该邻域范围内点的数量，即领域范围内的点的数量过少的情况下也将判定为噪声点，以此两条件来判断在此三维网格内的点是否为噪声点。

进一步的，步骤(4)中，利用非噪声点构建狄罗妮三角网，检测已经标记为噪声的点到对应三角面的距离，当距离小于给定的阈值时则认定为非噪声点。

进一步的，步骤(1)中，在点云数据的数据量满足处理速度的情况下，所述L的取值为3≦L≦5；在点云数据的数据量庞大的情况下，L≧5，或者对点云数据分块进行处理。

进一步的，步骤(1)中，所述最小尺度层级level0级的三维网格是在最小的单元的基础上进行的适应性增强，所述最小的单元是定义原始点云中XYZ三个方向的点间距大小为d_x、d_y、d_z，即以d_x、d_y、d_z构成的三维网格中包含的点数应约为1，所述适应性增强是对d_x、d_y、d_z进行放大：

d_x0＝d_x·sX；d_y0＝d_y·sY；d_z0＝d_z·sZ；

2≤sX,sY,sZ≤5

其中sX,sY,sZ为放大因子，增加放大因子之后，在level0的d_x0，d_y0，d_z0三维网格中约包含sX·sY个点；

其他尺度层级levell的3个方向的间距为d_xl＝d_yl＝d_zl，

d_xl＝d_x0·2^l；d_yl＝d_y0·2^l；d_zl＝d_z0·2^l。

进一步的，步骤(2)中，对每一个尺度层级l的三维盒范围都以包围盒的左下角：P_minl(x_min,y_min,z_min)，和右上角P_maxl(x_max,y_max,z_max)确定；

同理，为在实际操作中能够处理边界的点云数据，利用公式

$\left\{\begin{array}{c}{X}_{\min }=x_{\min }-R\·d_{\mathrm{xl}};Y_{\min }=y_{\min }-R\·d_{\mathrm{yl}};Z_{\min }=z_{\min }-R\·d_{\mathrm{zl}};\\ X_{\max }=x_{\max }+R\·d_{\mathrm{xl}};Y_{\max }=y_{\max }+R\·d_{\mathrm{yl}};Z_{\max }=z_{\max }+R\·d_{\mathrm{zl}};\end{array}\right.$

将三维包络范围坐标往外拓展两个边界点的坐标变为P_minl(X_min,Y_min,Z_min)，P_maxl(X_max,Y_max,Z_max)；其中X_min,Y_min,Z_min为左下角的三维坐标,X_max,Y_max,Z_max是右上角的三维坐标，所述R为整数，代表向外拓展的程度。

进一步的，步骤(4)中，对已经判定为噪声的点，即标记为0的点进行检测，如果满足以下条件则重新标记为非噪声点：

Distance＜T_dis＝k·d_z0；1≤k≤3

其中，Distance为待判定点到其所在的三角形面的距离，T_dis为阈值参数，满足此条件的点将重新标记为非噪声点1，并将参与到下一次构建三角网，此过程将不断迭代，创建新的三角网，直到没有点再被判定为非噪声点。

由于本算法能够有效检测出激光雷达点云和匹配点云中包含的噪声，故本算法在激光雷达点云滤波，误匹配检测，光束法平差中的噪声检测等领域具有良好的应用前景。

附图说明

图1是所述基于密度分析的多尺度点云噪声检测方法的流程示意图；

图2a、2b、2c是所述基于密度分析的多尺度点云噪声检测原理的示意图；

图3a是r窗口范围与R搜索空间的概念示意图；

图3b是r窗口范围与R搜索空间的概念立体示意图。

具体实施方式

本发明提供了一种基于密度分析的多尺度点云噪声检测方法，如图1所示，该方法包含以下步骤：

(1)输入点云数据，设定最大尺度；

(2)采用空间六面体进行点云数据的三维分割，建立离散点云三维网格和最大范围包围盒；通过局部密度分析方法标记噪声点为0并剔除，得到新的非噪声包围盒；

(3)在上一步得到的非噪声包围盒基础上，降低一级尺度层级，以1/2倍的三维网格大小逐个进行邻域的局部密度分析，进一步剔除较小的噪声点或噪声点聚类并标记为0，获得新的非噪声包围盒；迭代进行，直到完成最小尺度层级噪声点检测；

(4)完成尺度层级迭代后，需要检测标记为0的噪声是否为真正的噪声点：利用非噪声点构建狄罗妮三角网，检测标记为0的噪声点到对应三角面的距离，当距离小于给定的阈值时则认定为非噪声点。

本方法涉及的基于点云数据逐层分割方法原理如图2a、图2b所示。

其步骤如下所述：

首先定义原始点云中XYZ三个方向的点间距大小为d_x，d_y，d_z，即以d_x，d_y，d_z构成的三维网格中包含的点数应约为1，即用将点云进行三维剖分，剖分获得的网格中最小的单元的大小为d_x*d_y*d_z，在这个最小的单元(体素)中应该大致包含一个点。

本发明定义不同尺度层级用levell表示，其中l＝0～L，L为间距最大的最大尺度层级，本方法将最小尺度层级level0级的三维网格大小的计算做一定适应性增强，如公式(1)：

d_x0＝d_x·sX；d_y0＝d_y·sY；d_z0＝d_z·sZ；

                                           (1)

2≤sX,sY,sZ≤5

其中，sX,sY,sZ为放大因子，为顾及地物表面的起伏情况，取值不宜过大，上式的取值为多次实验中获得的经验参数。增加放大因子之后，可以认为在level0的d_x0，d_y0，d_z0三维网格中约包含sX·sY个点。

其中，上面所述的在最小的单元(体素)中应该大致包含一个点，以及认为在level0的d_x0，d_y0，d_z0三维网格中约包含sX·sY个点都是为了计算方便，从而这样确定。

对于其他尺度层级levell的3个方向的间距为d_xl＝d_yl＝d_zl，根据公式(2)计算得到：

d_xl＝d_x0·2^l；d_yl＝d_y0·2^l；d_zl＝d_z0·2^l      (2)

各尺度层级的三维包围盒范围将从大尺度层级以较大的网格间距首先建立，逐级剔除噪声点，对剩余的点云以小尺度层级、直到level0级的最小网格间距建立。

对于每一个尺度层级l的三维盒范围都以包围盒的左下角：P_minl(x_min,y_min,z_min)，和右上角P_maxl(x_max,y_max,z_max)确定。同样，为在实际操作中能够处理边界的点云数据，用公式(3)将三维包络范围坐标往外拓展两个边界点的坐标变为P_minl(X_min,Y_min,Z_min)，P_maxl(X_max,Y_max,Z_max)；

其中X_min,Y_min,Z_min为左下角的三维坐标,X_max,Y_max,Z_max是右上角的三维坐标

$\left\{\begin{array}{c}{X}_{\min }=x_{\min }-R\·d_{\mathrm{xl}};Y_{\min }=y_{\min }-R\·d_{\mathrm{yl}};Z_{\min }=z_{\min }-R\·d_{\mathrm{zl}};\\ X_{\max }=x_{\max }+R\·d_{\mathrm{xl}};Y_{\max }=y_{\max }+R\·d_{\mathrm{yl}};Z_{\max }=z_{\max }+R\·d_{\mathrm{zl}};\end{array}\right.---\left(3\right)$

其中，R为正整数树枝上等同于公式(5)中的R，代表向外拓展的程度。本发明中我们将每一个分割而成的三维网格命名为单元体。在某一尺度层级l下的某一单元体的位置坐标可用C(u,v,w)_l表示，其中u,v,w可以通过公式(4)确定：

$u=\left(\frac{x-X_{\min }}{d_{\mathrm{xl}}}\right);v=\left(\frac{y-Y_{\min }}{d_{\mathrm{yl}}}\right);w=\left(\frac{z-Z_{\min }}{d_{\mathrm{zl}}}\right)---\left(4\right)$

本方法涉及的基于基于邻域的局部密度分析的噪声检测算法原理如图3a、3b所示。其步骤如下：

噪声点的标记：

假设输入的点云数据集为W＝{W_i}，而其中包含的噪声点集为P＝{p_i|x_i,y_i,z_i}，i＝0～n-1，n为点云总个数。噪声点的检测按从大尺度层级到小尺度层级逐级地进行。设第l级的点集为N_l＝{N_i|x_i,y_i,z_i}_l，即为点云数据集W＝{W_i}减去已经检测到的噪声点后的数据集。点集W＝{W_i}中的检测出的噪声点将被标记为0，而得到的非噪声点N_l将被标记为1。在检测初始时所有的点都标记为1。

r窗口范围和R搜索空间：

对特定尺度层级下，本方法对噪声点的检测以某一单元体为中心一定邻域范围形成空间范围做为一个单位，在一个指定的搜索空间内进行遍历，计算该邻域范围内的密度和包含的点的数量，并以此来判断在此三维网格内的点是否为噪声点。单位体的邻域范围本发明称之为r窗口范围，指定的搜索空间本发明称之为R搜索空间。

本发明定义r窗口范围为：以某个单元体C(u₀,v₀,w₀)为中心，并以r为半径的三维窗口范围(即该三维窗口的边长为2r+1)作为其邻域的空间范围。R搜索空间为以C(u₀,v₀,w₀)为中心，定义一个以更大的半径R定义的三维搜索空间(即该三维窗口的直径为2R+1)。注意：此处的R与公式(3)中的R在数值上是相等的。如图3a、3b分别以二维和三维形式展示R搜索空间与r窗口范围的定义，其中图3a中r＝1且R＝2，图3b中r＝1，R为大于2的某一正整数。对于处于levell层级的某一个单元体C(u₀,v₀,w₀)，所有满足公式(5)的单元体C(u,v,w)都在R搜索范围内。此处，R为整数，与公式(3)中的R数值相等。

$\begin{array}{c}{u}_0-R\≤u\≤u_0+R;\\ v_0-R\≤v\≤v_0+R;\\ w_0-R\≤w\≤w_0+R;\end{array}---\left(5\right)$

为表义明确将需要使用到的概念在表1中定义。

以某一尺度层级levell的某一个单元体C(u₀,v₀,w₀)为中心，需要执行四个步骤来判断该单元体内的点是否为噪声点。

1)对单元体C(u₀,v₀,w₀)，计算如果则不处理该单元体，开始计算下一个单元体；若则进入步骤2；

2)遍历获得R范围内的每一个有效单元体C(u,v,w)的NP(C_u,v,w,r)计算并以此利用公式(6)计算C(u₀,v₀,w₀)的局部密度

$\mathrm{\ϵ}(C_{u_0,v_0,w_0},r,R)=\mathrm{NP}(C_{u_0,v_0,w_0},r)-\mathrm{ANP}(C_{u_0,v_0,w_0},r,R)---\left(6\right)$

3)本发明提出的算法同时使用两个条件判定C(u₀,v₀,w₀)内的点是否为噪声点，两个条件只要有一个被满足时即接受为噪声点，并标记为0；非噪声点则标记为1。

条件1：满足阈值定义如公式(7)(8)所示：

$N_{\mathrm{threshold}}=\mathrm{MaxNP}{(C_{u_0,v_0,w_0},r)}_l\&CenterDot;\mathrm{rate};0.0<\mathrm{rate}<1.0---\left(7\right)$

$\mathrm{MaxNP}{(C_{u_0,v_0,w_0},r)}_l=\mathrm{sX}\&CenterDot;\mathrm{sY}\&CenterDot;r\&CenterDot;2^l\&CenterDot;2^l---\left(8\right)$

其中rate无法检测出很多噪声，过高则容易将非噪声判定为噪声。为levell层级下C(u₀,v₀,w₀)的r邻域范围内包含的点数的最大估算值。

条件2：满足 $\mathrm{\&epsiv;}(C_{u_0,v_0,w_0},r,R)<-\mathrm{Kt}\&CenterDot;\mathrm{\&delta;}(C_{u_0,v_0,w_0},r,R)$

Kt的取值依据正常分布的概率密度函数，如果数据分布是正常的，则位于均值Kt倍的标准差范围内的数据的比例称为误差函数(error function，ERF)，具体形式如公式(9)：

$\mathrm{erf}\left(\mathrm{Kt}\right)=\frac{2}{\sqrt{\mathrm{\&pi;}}}{\&Integral;}_0^{\frac{\mathrm{Kt}}{\sqrt{2}}}{e}^{{-x}^2}\mathrm{dx}---\left(9\right)$

99.7％的正常分布的数据位于高于或低于均值3倍标准差的范围内。故Kt＝3且满足条件2时，C(u₀,v₀,w₀)内包含的点可以被认定为噪声点。

4)重复步骤1)～3)，直到该尺度层级下的包围盒中所有的单元体都被处理。

本方法涉及的误判噪声点的检测算法原理如图2c所示，其步骤如下：

虽然以上步骤能够去除大量噪声，但是可能还是会出现一些非噪声点因分布过于稀疏而满足以上步骤3)中两个判定条件之一而被认定为噪声点。为将这些误检测为噪声点重新判定为非噪声点，本发明首先利用已经判定为非噪声点的点(标记为1的点)构建一个狄罗尼三角网，然后对已经判定为噪声的点(标记为0的点)进行检测，如果满足以下条件则重新标记为非噪声点：

Distance＜T_dis＝k·d_z0；1≤k≤3   (10)

其中，Distance为待判定点到其所在的三角形面的距离，T_dis为阈值参数。满足此条件的点将重新标记为非噪声点1，并将参与到下一次构建三角网。此过程将不断迭代，创建新的三角网，直到没有点再被判定为非噪声点。

下表是局部密度算法中的概念定义：

本发明提出了一种新的点云去噪算法，可适用于这两类数据中所包含的噪声点的去除。算法主要包括两步：第一步利用多尺度的密度算法去除孤立噪声和小的簇状噪声；第二步利用三角网约束将第一步中误检测为噪声的点重新归为正常点。由于本算法能够有效检测出激光雷达点云和匹配点云中包含的噪声，故本算法在激光雷达点云滤波、误匹配检测、光束法平差中的噪声检测等领域具有良好的应用前景。

Claims (7)

1.一种基于密度分析的多尺度点云噪声检测方法，包括以下步骤：

(1)输入点云数据并设定尺度层级，不同尺度层级用level^l表示，其中l＝0～L，其中最小尺度层级为level0级，最大尺度层级为levelL级；

(2)采用空间六面体进行点云数据的三维分割，建立离散点云三维网格和最大范围包围盒，所述最大范围包围盒根据输入的点云三维坐标计算得到，通过局部密度分析方法标记噪声点为0并剔除，得到新的非噪声包围盒；

(3)在上一步得到的非噪声包围盒基础上，降低一级尺度层级，以1/2倍的三维网格大小逐个进行邻域的局部密度分析，进一步剔除较小的噪声点或噪声点聚类并标记为0，获得新的非噪声包围盒；迭代进行，直到完成最小尺度层级噪声点检测；

(4)完成尺度层级迭代后，检测标记为0的噪声是否为真正的噪声点，利用非噪声点构建狄罗妮三角网，检测标记为0的噪声点到对应三角面的距离，当距离小于给定的阈值时则认定为非噪声点，所述阈值设定为最精细层级的三维网格的在Z方向的尺寸。

2.根据权利1所述的基于密度分析的多尺度点云噪声检测方法，其特征在于：步骤(3)中，对特定尺度层级下，对噪声点的检测是以某一单元体为中心一定邻域范围形成空间范围做为一个单位，所述单元体代表每一个分割而成的三维网格，在一个指定的搜索空间内进行遍历，计算该邻域范围内的密度与搜索范围内的密度的关系，同时判断该邻域范围内点的数量，即领域范围内的点的数量过少的情况下也将判定为噪声点，以此两条件来判断在此三维网格内的点是否为噪声点。

3.根据权利1所述的基于密度分析的多尺度点云噪声检测方法，其特征在于：步骤(4)中，利用非噪声点构建狄罗妮三角网，检测已经标记为噪声的点到对应三角面的距离，当距离小于给定的阈值时则认定为非噪声点。

4.根据权利1所述的基于密度分析的多尺度点云噪声检测方法，其特征在于：步骤(1)中，在点云数据的数据量满足处理速度的情况下，所述L的取值为3≦L≦5；在点云数据的数据量庞大的情况下，L≧5，或者对点云数据分块进行处理。

5.根据权利1所述的基于密度分析的多尺度点云噪声检测方法，其特征在于：步骤(1)中，所述最小尺度层级level0级的三维网格是在最小的单元的基础上进行的适应性增强，所述最小的单元是定义原始点云中XYZ三个方向的点间距大小为d_x、d_y、d_z，即以d_x、d_y、d_z构成的三维网格中包含的点数应约为1，所述适应性增强是对d_x、d_y、d_z进行放大：

d_x0＝d_x·sX；d_y0＝d_y·sY；d_z0＝d_z·sZ；

2≤sX,sY,sZ≤5

其中sX,sY,sZ为放大因子，增加放大因子之后，在level0的d_x0，d_y0，d_z0三维网格中约包含sX·sY个点；

其他尺度层级levell的3个方向的间距为d_xl＝d_yl＝d_zl，

d_xl＝d_x0·2^l；d_yl＝d_y0·2^l；d_zl＝d_z0·2^l。

6.根据权利1所述的基于密度分析的多尺度点云噪声检测方法，其特征在于：步骤(2)中，对每一个尺度层级l的三维盒范围都以包围盒的左下角：P_minl(x_min,y_min,z_min)，和右上角P_maxl(x_max,y_max,z_max)确定；

同理，为在实际操作中能够处理边界的点云数据，利用公式

$\left\{\begin{array}{c}{X}_{\min }=x_{\min }-R\&CenterDot;d_{\mathrm{xl}};Y_{\min }=y_{\min }-R\&CenterDot;d_{\mathrm{yl}};Z_{\min }=z_{\min }-R\&CenterDot;d_{\mathrm{zl}};\\ X_{\max }=x_{\max }+R\&CenterDot;d_{\mathrm{xl}};Y_{\max }=y_{\max }+R\&CenterDot;d_{\mathrm{yl}};Z_{\max }=z_{\max }+R\&CenterDot;d_{\mathrm{zl}};\end{array}\right.$

将三维包络范围坐标往外拓展两个边界点的坐标变为P_minl(X_min,Y_min,Z_min)，P_maxl(X_max,Y_max,Z_max)；其中X_min,Y_min,Z_min为左下角的三维坐标,X_max,Y_max,Z_max是右上角的三维坐标，所述R为整数，代表向外拓展的程度。

7.根据权利6所述的基于密度分析的多尺度点云噪声检测方法，其特征在于：步骤(4)中，对已经判定为噪声的点，即标记为0的点进行检测，如果满足以下条件则重新标记为非噪声点：

Distance＜T_dis＝k·d_z0；1≤k≤3

其中，Distance为待判定点到其所在的三角形面的距离，T_dis为阈值参数，满足此条件的点将重新标记为非噪声点1，并将参与到下一次构建三角网，此过程将不断迭代，创建新的三角网，直到没有点再被判定为非噪声点。