1 引言

随着三维激光扫描技术的快速发展,地面三维建模的效率和精度越来越高。实际工程中,点云数据量十分庞大,在数据的采集过程中不可避免地存在着很多噪声点,点云噪声的性质和种类随点云采集环境的变化而变化。噪声会极大地影响地面三维建模的精度,因此,需要采取合适的算法对不同场景下的噪声进行剔除。

点云目标提取作为近年来的一大研究热点,国内外专家提出了较多的算法和理论^[1-14],采用提取算法的前提是剔除非研究对象点云,即滤波。这些滤波算法大致可分为4类:基于形态学,基于聚类密度,基于曲面拟合,基于统计特征。王竞雪等^[15]提出结合形态学滤波和不规则三角网(TIN)的滤波算法,解决了传统TIN的滤波精度对初始种子点选取的依赖问题;邢承滨等^[16]提出基于形态学滤波的等值线生成法,可直观地确定滤波算法中的关键参数阈值;赵明波等^[17]提出一种渐进多尺度数学形态学滤波算法,通过改进形态开运算处理大面积空白区域。文献[ 15-17]方法归类于第一类,此类方法较难解决阈值的自适应问题。牛晓静等^[18]提出自适应聚类和滤波融合的方法对点云进行去噪平滑处理;苏本跃等^[19]基于K-means聚类对RGBD点云进行外部噪声去噪和精简,可保留点云的原始特征;张巧英等^[20]将密度聚类算法应用于连续分布点云去噪中,极大限度地保留了原始点云的特征;田青华等^[21]运用改进的欧氏聚类的散乱工件点云分割方法解决了场景分割困难的问题。文献[ 18-21]方法归类于第二类,但该方法邻域搜索普遍耗时较长。冯东恒等^[22]提出用共面法向量特征去除离群面的RANSAC算法,用于去除近地表噪声和散点噪声;林万誉等^[23]提出基于最小截取二乘法的分段拟合的滤波算法,该算法适用于多曲面三维点云;赵凯等^[24]提出基于体素栅格的快速密度聚类算法用于实时处理点云。文献[ 22-24]方法都归类于第三类,更适用于小范围、表面连续的数据。李明磊等^[25]提出的孤点的加权平均距离去除法降低了对远距离离散点的误判;李仁忠等^[26]提出一种基于方法库的点云模型去噪与精简算法。文献[ 25-26]方法都归类于第四类。

三维激光扫描仪是以固定角度分辨率进行等角扫描的,周围地物、地貌的空间分布不同导致距离扫描仪越远的点云越稀疏、密度越不均匀^[27]。因此,不能简单依靠一种方法进行目标提取。本文采用大型矿山点云数据,其具有扫描范围大、点云数据多、噪声明显、地形地貌不规则等特性。为有效剔除该场景下的非地面点数据并提高运算效率,提出一种结合统计滤波与密度聚类的矿山地面点云提取算法。该算法从与聚类算法相似的特征空间入手,通过分析不同近邻点特征空间,找出最优近邻点数,最后通过对特征空间进行提取并求交得到地面点集,在保证地面点提取精度的前提下加快了运算速度。

2 原理与方法

本文基于三维激光扫描仪RIEGL VZ-1000收集的德兴铜矿点云数据,研究矿山地面点云的提取算法,该数据的角度分辨率为0.0005°,100 m的扫描精度为5 mm,图1为2018年12月该矿区的部分点云数据,共510519个点,横纵方向扫描角设置为0.02°,受扫描角度的限制,及地形地貌和随机误差的影响,数据中存在大量的非地面点和少许未扫描到的空洞。由于矿区点云的范围巨大、点数众多,并掺杂无规律的噪声点,因此提取地面点云算法的精度与剔除随机误差点的能力密切相关。通常矿山场景中点云的噪声主要是来自于大尺度远离主体的孤立噪声点和小尺度靠近主体的非地面的随机噪声点等。因此,首先需要对该地区的噪声空间分布特性进行分析及分类;其次,建立点云索引,搜索一定邻域内的点集,并运用统计滤波的特征空间表达点集;然后,基于特征空间的特征表达原则再结合二维特征密度的分布特性进行聚类求交;最终可消除随机误差的影响,并提取地面点集。

图 1. 原始的点云

Fig. 1. Original point cloud

下载图片 查看所有图片

2.1 噪声空间分布特征

本文算法主要研究的是矿山地面点的提取,首先根据非研究对象的空间分布特性对矿山点云模型中的该对象作尺度定义,本文将矿山点云模型中的非研究对象分为大尺度孤立点、中尺度径向分布点及小尺度障碍物噪声点。这三种对象在矿山点云中表现尤为明显,严重影响地表点云提取精度与后期精确地表建模,因此,剔除这三类对象的精度及效率对地面点云提取的研究具有重大作用及意义。

图2(a)为平面点云图。该三种对象的简化二维模型如图2(b)所示,横坐标为点的个数,纵坐标为距离,以O、O₁、O₂这3点为例,圆点O代表小尺度障碍物噪声点,菱形点O₁代表大尺度孤立点,三角点O₂代表中尺度径向分布点,搜索其20个近邻点的位置,结果显示:大尺度孤立点的极差、方差及均值等特性远大于中尺度径向分布点和小尺度障碍物噪声点,而中尺度径向分布点的特征值也略微大于小尺度障碍物噪声点。对这20个近邻点到原点的欧氏距离进行线性回归后可明显看出:大尺度孤立点的近2个邻域点值很小,其统计也较为接近一般点,但随着近邻点的增加,其拟合的斜率又出现间断跳跃增长;中尺度径向分布点随着近邻点数量的增加对其拟合的斜率的影响并不是很大,但其斜率大于非噪声点的;小尺度障碍物噪声点的邻域点较为特殊,会出现明显的层次感,前10个点的斜率一直很小,但延续到11~16个点又会出现波动,到17~20个点时影响更会大,其斜率的波动随着邻域点的个数增加而增加,增加量变化有一定规律。

图 2. 噪声的空间分布。(a)平面点云图; (b)邻点的距离图

Fig. 2. Spatial distribution of noise. (a) Plane point cloud map; (b) distance map of neighbor points

下载图片 查看所有图片

总的来说,大尺度孤立点的统计特征表现为极差、方差等统计特征远大于其他噪声;中尺度径向分布点的统计特征表现为邻域点距离呈现规律增长,及斜率的变化基本稳定;小尺度障碍物噪声点的统计特征表现为斜率具有间断变化,且其变化率相对稳定。

2.2 KD-tree索引与统计特征

本算法基于统计滤波算法的统计描述和密度聚类算法的空间描述。统计滤波的统计描述是基于KD-tree的k近邻点搜索,即给定查询点及邻域点个数k,就能从数据集中找到距离该点最近的k个数据,再由近邻点的统计特征变化确定截断阈值,找出噪声点并剔除。

本文借鉴这种高效率的搜索算法对传统算法进行改进,统计滤波只是单一地使用一种或两种特征值进行特征空间分析,从滤波效果来看,其有一定的噪声剔除能力,但是缺乏效率与准确性,因此本算法需要在减少运算次数的同时,提高非研究对象的剔除精度,并分析和挖掘统计数据的规律。

为了有效识别非研究对象,并表明近邻点的空间特征,根据文献[ 28-30]中的相关理论知识,建立统计特征量均值(mean)、方差(variance)、斜率(slope)、极值(range)、能量(energy)。k个近邻点距离的集合为{d_ij,i≤N₀,i∈N⁺,j≤k,j∈N⁺},其中N₀为总的搜索点个数,d_ij为第i个搜索点到第j个近邻点的距离;同时, 第i个搜索点的均值为dmean=1k∑j=1kdij,第i个搜索点的方差为dstd=1k-1∑j=1k(dij-dmean)2,第i个搜索点的斜率为dslope=∑j=1kwijdijxj/∑j=1kwijxj2,其中w_ij表示距离的权重,x_j表示近邻点标号,第i个搜索点的极差为d_range=max(d_ij)-min(d_ij),第i个搜索点的能量为 den=∑j=1kdij2。由于利用单一的特征值无法对非研究对象进行精确区分,例如:中尺度径向分布点和小尺度障碍物噪声点的均值有可能相似,随着近邻点的增加这两种噪声的斜率也难以区分,因此,需要综合考虑各种统计特征的变化情况。此外,为了应对特征值的间断分层变化,对特征值进行分段加权处理,以便更好地区分,如:在斜率的计算中,对不同近邻点设置不同的权重。

对以上特征取不同的近邻点个数时,所反映的特征信息是不同的,近邻点个数太少或太多,都会增加错分率。为了能够高效地分析出不同场景下最少错分率的近邻点个数,可先建立一个较多近邻点的KD-tree索引,避免因重复建树而导致效率降低,例如:对本文的510519个点,当近邻点个数k=100时,先储存510519×100的距离矩阵,当近邻点个数选择k=40时,直接调用即可,不需要重复建树。

2.3 二维特征密度聚类算法

以上对统计特征的分析是对目前传统统计方法的优化与改进,但也只是基础分析与应用。基于相同属性地物具有相同或相似的特征空间分布的原理,首要分析非研究对象的特征空间,并以此确定划分阈值。同时,为了更加细致地分析近邻点的统计特征变化情况,并以此剔除非研究对象点,本文结合统计滤波和二维特征密度聚类提出一种矿山地面点的提取算法。所提算法流程图如图3所示。

图 3. 算法流程图

Fig. 3. Flow chart of proposed algorithm

下载图片 查看所有图片

本文算法在KD-Tree的创建中,仅进行一次建树,计算量为n²;在统计特征值描述上,因为特征统计量增加为5个,计算量为5n;在空间聚类上,需遍历特征向量组的密度空间进行噪声剔除,再结合交并差提取地面点,这一步均只使用索引ID进行处理,计算量为3n;因此,可定义本文算法的算法复杂度为o(n²),从理论上讲本文算法所需的运行时间更短。

在统计描述中,从KD-tree索引的建立、距离矩阵的查询以及统计特征的计算上,已经对算法的效率进行了优化,节约了运算时间;在空间描述和特征密度计算上,本文对5个基础特征:均值、方差、斜率、极值、能量,分别建立10个一一对应的二维密度特征空间,特征向量组为均值-斜率、均值-方差、均值-极值、均值-能量、斜率-方差、斜率-极值、斜率-能量、方差-极值、方差-能量、极值-能量,依次标号为(1), (2), (3),…,(10)。

根据密度空间的离散情况,可保留不同倍数中误差δ的有效数据作为所需的研究对象,其中2δ、3δ的有效数据比例约为95.44%和99.74%,由此可计算出各密度空间的截断阈值ω_i,并对这10种二维特征密度聚类后的结果取交集Q_j,Q_j即为有效研究对象,可表示为

Dij=Ti,{Ti≤ωi,i∈N+}Qj=D1j⋂D2j⋂…⋂D10j,{j∈N+},(1)

式中:N⁺为正整数;i为特征向量组的标号;j为第j个近邻点;T_i为各特征向量组的密度;D_ij为T_i中小于阈值ω_i的有效点云;Q_j为各特征向量组剔除非研究对象后的有效地面点的并集。 最后,计算非研究对象点的剔除效率是否最高,是即可输出,否则重新选择近邻点个数k。

该算法中,均使用索引ID号进行交并差处理,因此,在计算中仅进行两层遍历,且只针对ID号进行索引,最后根据ID号返回确定真实坐标值。

3 实验结果与分析

3.1 特征密度空间分析

点云去噪的首要目标是对10个二维特征密度空间进行分析,分析各组特征向量之间的非研究对象的剔除能力,分析结果如图4、图5、图6、图7所示,其中直线表示各组的点的拟合趋势线。

图 4. 邻域10个点的特征密度分布

Fig. 4. Characteristic density distributions of 10 points in neighborhood

下载图片 查看所有图片

图 5. 邻域20个点的特征密度分布

Fig. 5. Characteristic density distributions of 20 points in neighborhood

下载图片 查看所有图片

如图4所示,近邻点为10时,噪声点在这些特征空间中都具有较大的离散情况,且噪声点的特征空间密度都相对较低。如:图4(1)的特征空间中,均值越大,斜率也越大,剔除非地面点后的密度空间表示出较高的密集性;图4(2)、(3)的特征空间中,虽然正常点云有一定的离散性,但非研究对象的离散性更大,例如:当某些点云极差较大时,有较高概率被判定为大尺度孤立点,同时,它们的均值和方差也较大时,这类点云就一定是需要被剔除的噪声点。图4(4)、(7)的特征空间中,噪声点的能量都很大,且与均值和斜率的变化相关性较大,正常点云的能量、均值、斜率都较低;其他的特征空间中,都有效地去除了一些点,这些特征空间密度聚类后的交集 Qj包含的是符合所有特征空间判定的有效地面点云,而除此之外的点即为非地面点,这其中包含了三大类非研究对象;对于在某一特征空间中具有特有属性的非研究对象来说,其判断结果具有较大的不确定性,对于这些无法准确区分的点云而言,可以通过调节特征值在各近邻点上所占层次权重进行聚类,以此筛选判断出有效地面点。

根据图4,10个近邻点的二维密度特征空间能够反映出绝大部分非研究对象,再结合密度阈值,即可获得有效点云对象,结果表明该算法可通过剔除非地面点云,实现目标点云的提取。

如图5所示,近邻点为20时,相比于10个近邻点的二维特征密度空间分布,20个近邻点的特征空间有一定的变化。首先,由图5(2)、(3)、(5)、(8)、(9)、(10)的特征空间分布可知,其噪声点的离散程度大于10个近邻点的,同时,密度空间的范围也有所扩大,分辨出非地面点云的能力也在增强。在图5(1)、(4)、(6)、(7) 的特征空间中,有效点的密集度良好,且其特征空间的范围呈有规律的扩增。

如图6所示,近邻点为30时,相比于图5的密度分布,30个近邻点的特征空间的变化较大,这有可能是受到间断分层式非地面点云特性的影响。首先,由含方差的特征空间可以看出,方差对噪声点离散程度的响应大于10和20个近邻点时的情况,同时,二维特征密度空间的范围呈非等差增长,其特征空间大小的增率放缓,即分辨出非地面点云的能力的增率大于邻域点个数的增率。

图 6. 邻域30个点的特征密度分布

Fig. 6. Characteristic density distributions of 30 points in neighborhood

下载图片 查看所有图片

如图7所示,近邻点为40时,相比于图4、图5、图6的密度分布,40个近邻点的特征空间中的噪声点离散性最强,且剔除非地面点云后的特征空间明显优于10和20个近邻点时的效果。综上,40个近邻点使特征空间中的噪声变得格外突出,即有效点的密度集聚度好、离散度大。但其特征空间范围的增率较大,这导致有可能出现过剔除的情况,即将部分有效地面点也剔除了。

算法实验结果如图8所示,随着近邻点的增加,首先,大部分的小尺度障碍物点、大尺度孤立点以及中尺度径向分布点被剔除了;其次,少许边缘凹陷处和空洞周围的稀疏点云也被剔除。剔除非地面点云后的三维点云效果图均能有效表征提取的矿山地面三维形态。图8中,在近邻点为40个时,部分有效地面点也被剔除,而在近邻点为10个时,该部分点周围仍有未剔除的非地面点。

综上所述,本文算法在地面点的提取上是有效的,邻域点个数的增加对提取的效果有一定的影响。为了能定量定性分析该算法的优劣性,本文将从以下4点进行阐明: 1)不同邻域点及向量组的拟合斜率差变化情况;2)不同邻域点及向量组的拟合精度;3)不同邻域点及向量组联合去噪后的地面点个数;4)均方差等误差曲线的走势。

图9(a)为不同邻域点及向量组的拟合斜率差的变化,在选不同个数的邻域点时,图9(a)中(1)、(2)、(3)、(4)、(8)、(10)的特征空间对应的拟合的斜率差变化不明显,即提取地面点的整体拟合的斜率变化不大;(5)、(6)的特征空间的拟合斜率随着邻域点的增加逐渐减小,即其特征空间中的有效地面点更偏重于横坐标上的特征值——斜率; (7)、(9)的特征空间的拟合斜率,随着邻域点的增加逐渐增加,即其特征空间中的有效地面点更偏重于纵坐标上的特征值——能量。这表明能量、斜率与地面点的相关性很大。

图 7. 邻域40个点的特征密度分布

Fig. 7. Characteristic density distributions of 40 points in neighborhood

下载图片 查看所有图片

图 8. 提取效果。(a)原始云点;(b) 10个邻域点;(c) 20个邻域点;(d) 30个邻域点;(e) 40个邻域点

Fig. 8. Extracted results. (a) Original point cloud; (b) 10 neighbor points; (c) 20 neighbor points; (d) 30 neighbor points; (e) 40 neighbor points

下载图片 查看所有图片

图9(b)为不同邻域点及向量组的拟合精度差的变化情况,在邻域点个数不同时,其中(1)、(4)、(7)的特征空间的拟合精度差变化不大,从侧面印证了剔除的非地面点云大部分为随机噪声;(2)、(3)、(5)、(6)、(7)、(8)、(9)的特征空间的拟合精度差随着邻域点的增加逐渐增加,即其提取后的拟合精度优于提取前的拟合精度;尤其是(2)、(5)、(9)的特征空间,其提取后的拟合精度远高于提取前,说明提取后的特征密度集聚度高,剔除非地面点的能力强。由此可见,剔除随机噪声能使特征空间的拟合精度增加,聚集性更好。

图 9. 实验结果及精度。(a)斜率差的变化;(b)精度差的变化;(c)地面点云个数;(d)误差曲线

Fig. 9. Experimental results and accuracy. (a) Trend of slope difference; (b) trend of accuracy difference; (c) number of ground point clouds; (d) error curves

下载图片 查看所有图片

图9(c)为各特征向量组在取不同邻域点个数时提取的地面点云个数,在取不同的近邻点个数时,(1)、(4)、(7)的特征空间的提取效果较差,只能去掉2000~4000个随机噪声点;(2)、(3)、(5)、(6)、(7)、(8) 、(9)的特征空间剔除随机噪声点的个数随着邻域点的增加而增加,即其剔除非地面点的能力也在增强;图9(c)中最后1列为本文算法的最终提取地面点云个数,并且随着近邻点个数的增加,其剔除噪声点的个数也在增加,同时增率放缓。

为了验证本文算法的有效性,采用国际摄影测量与遥感(ISPRS)的总错分误差标准,将有效地面点错误分为噪声点的数量和将噪声点错误分为有效地面点的数量之和占总和点云的比例为

Rtotal⁃error=c+da+b+c+d×100%,(2)

式中:a为有效点数量;b为噪声点数量;c为将有效地面点错误分为噪声点的数量;d为将噪声点错误分为有效地面点的数量,本实验已知有效点个数为493049个。图9(d)横坐标是近邻点的个数,纵坐标为总错分误差和最近邻点欧氏距离的均方差。根据图9(d)及表1所示的误差曲线及提取结果,均方差随着近邻点的增加逐渐降低,说明随着近邻点的增加,提取的效果越来越好;由总误差曲线可知,随着近邻点的增加总误差减小,并且在36个近邻点时达到极小值,总误差最优为0.007707,MSE为0.019633 m;由表1可知,本算法在对510519个点进行矿山地面点云提取时的运算时间较短。

3.2 提取算法验证

为验证本文算法的提取精度与效率,使用基于方法库的点云模型。由于文献[ 26]中未曾提及各方法的具体参数,本文先用控制变量法研究半径滤波、统计滤波和体素栅格滤波法在不同参数下的去噪效果。原始数据扫描角为0.02°,实验结果如图10和表2所示。

图 10. 半径滤波、体素滤波及统计滤波提取效果

Fig. 10. Extraction results of radius filter, voxel filter, and statistical filter

下载图片 查看所有图片

由图10和表2可知,当半径滤波的滤波半径与点数不同时,其噪声的剔除能力不同,均方差和分类误差均大于本文算法;当选择不同的体素栅格时,其去噪效果与数据的抽稀类似,无法准确剔除噪声且均方差和分类误差极大,无法实现有效的地面点提取;当统计滤波选择不同近邻点数时,剔除噪声能力较好,但容易将稀疏地面点误判为噪声,有效地面点保留较少,且耗时较长,在选取相同近邻点的情况下,用时约为本算法的7倍。

由于双边滤波对低频数据的过滤效果较好,但是耗时较长。经比较,确定方法库中各方法对应的参数:统计滤波中选择近邻点个数为36;双边滤波中滤波器边长为0.5 m,标准偏差为[1.0,0.1];体素滤波中滤波器边长为0.1 m;半径滤波中滤波半径为0.8 m,近邻点数为4。基于该方法库对矿山点云数据进行处理,提取得到地面点云493522个,方差为0.0196,极差为3.9511,整体误差为0.0062,耗时374.5536 s。实验结果如图11所示。

图 11. 基于方法库的提取效果。(a)去噪正面效果;(b)去噪侧面效果

Fig. 11. Extraction results based on Method-Library. (a) After denoising in front view; (b) after denoising in side view

下载图片 查看所有图片

综上,本文算法的时间效率明显优于半径滤波、体素滤波、统计滤波及方法库去噪模型,且在地面点的提取上效率较高,相比其他算法而言,对非目标对象的误判率更低。为验证本文算法的普适性,选用一期扫描参数与原实验数据相似的更大面积的点云数据,并沿用上文中遍历得到的最优近邻点个数36进行地面点提取,得到提取效果如图12所示。原始点云为5866314个,去噪后为5830264个,耗时277.7259 s。

图 12. 对大面积点云的提取效果。(a)原图正面;(b)去噪后正面;(c)原图侧面;(d)去噪后侧面

Fig. 12. Extraction results of large-area point cloud. (a) Before denoising in front view; (b) after denoising in front view; (c) before denoising in side view; (d) after denoising in side view

下载图片 查看所有图片

4 结论

本文提出了一种结合统计滤波与密度聚类的矿山地面点云提取算法,利用KD-tree搜索近邻点并计算其统计特征,同时结合二维特征密度空间的分布特性对密度空间进行聚类并提取地面点云。以VZ-1000地面三维激光扫描仪获取的大型矿山数据为研究对象,实现对大型矿山地面点云的提取。对KD-tree在不同近邻点个数下的地面点提取精度与效率进行研究,分析表明当近邻点个数为36时,算法总误差最小。此外选用另一期矿山点云数据对算法的普适性进行验证。综合分析考虑10组特征向量组的提取结果,同时进行求交集处理以提高有效点云的精度,通过对传统算法的研究可知,统计滤波的提取效果较好,但精度仍低于本文算法,且耗时约为本文算法的7倍。本文算法原理简单,易于实现,且能快速、精确地提取有效地面点云,具有一定的实用价值。但所提算法仍存在不足之处,主要表现为凹陷处边缘、大型障碍物的阴影边缘等位置的点云较为稀疏,容易造成错误划分,还需进一步深入研究,希望通过自适应调整统计特征值的权重解决该问题。