一种改进的形态学LiDAR点云滤波算法

周钦坤,岳建平,李子宽,杨 恒

(河海大学地球科学与工程学院, 江苏 南京 211100)

机载激光雷达(LiDAR,light detection and ranging)技术可快速获取大范围的空间三维坐标数据,目前已被广泛应用于城市三维建模、DEM快速生成、森林资源调查。经过多年其发展成熟,但数据后处理还相对滞后,针对如何从点云中区分地面点与非地面点,国内外众多学者提出了多种算法。其中Vosselman[1]提出了基于坡度的滤波算法,通过比较两临近激光脚点间高差与阈值的大小关系来确定点云属性;张小红等[2]提出了移动曲面拟合滤波算法,用一个简单的二次曲面去过滤整个点云表面;文献[3-4]中通过融合多种方法,也取得了较好的效果;Lindenberger[5]提出了形态学滤波算法,对某一固定窗口内的点云进行开运算,通过计算运算前后高差变化来确定点云属性。由于采用固定窗口,该算法仅能去除固定大小的地物,为此,Zhang等[6]在原算法基础上进行了改进,通过逐步改变窗口大小,并设置与窗口大小相适应的高差阈值以提高滤波精度。

Zhang等[6]提出的渐进形态学滤波能去除不同尺寸的地物,但是需要用户预设的参数过多,严重影响了算法的自适应性。另外,在一个标记网格中往往有多个高程不同的激光点云,如果简单的将这些激光点云归为同一类,那么点云的误判率会很高。基于此,众多学者对该算法进行了改进[7-10]。研究通过增加迭代终止条件,免除最大窗口尺寸阈值的预设工作,减少算法设置参数的个数,并进一步对网格中的点进行精分类以降低格网内点云的分类误差。

1 渐进形态学滤波

1.1 滤波理论基础

渐进形态学滤波算法基于集合运算,腐蚀与膨胀是最基本的2个算子。在进行数据处理时,常用2种算子的组合,一般有开运算和闭运算,运算可表达成如下形式:

开运算:(f∘g)=(f⊖g)⊕g,

闭运算:(f·g)=(f⊕g)⊖g,

其中:∘为开运算符号;·为闭运算符号;⊖为腐蚀符号;⊕为膨胀符号。

在用形态学对机载激光雷达点云进行滤波时,主要使用开运算,可以去除比窗口尺寸小的地物,而保留比窗口尺寸大的地物。

1.2 滤波基本原理

渐进形态学滤波的基本思想是,逐步增大滤波窗口对LiDAR数据进行开运算,以去除LiDAR数据中的不同大小的树木及建筑物等,如图1所示。

图1 渐进形态学滤波Fig.1 Progressive morphological filtering

通过变化窗口大小,地物点均被去除,但部分地面特征点也被滤除,为了在滤除非地面点的同时,保留地面点,引入了高程差阈值。对于任一点p,设dhp为前后2次滤波p点处的高程差,dhT为此次过滤的高差阈值,若dhp

2 改进的渐进形态学滤波

2.1 免预设最大窗口尺寸的迭代处理

在传统渐进形态学滤波算法中,最大窗口尺寸阈值控制着整个算法的迭代,当窗口增大到大于窗口尺寸最大阈值时,迭代停止,点云被分为地面点与非地面点2类。窗口阈值过大或过小,都会影响滤波精度,并且窗口大小一般为3、5、7、9……这样的离散数列,当设置的最大窗口阈值恰好不等于这些数时,传统渐进形态学滤波算法将以小于该阈值的窗口大小进行最后迭代,地面点就会残存大量非地面点。鉴于上述问题,依据当窗口达到最大建筑物尺寸时,能最大程度的去除非地面点,继续增大窗口,本窗口下的地面点数将近乎不变的思想,对传统算法进行改进,根据前后2次迭代地面点比率大小来判断迭代是否终止,免去最大窗口尺寸设置工作。具体步骤如下:

(1) 原始点云预处理,剔除点云噪声点与点云格网化。

(2) 进行迭代,以较小窗口过滤点云数据,统计迭代过后的地面点数量c1。

(3) 增加迭代次数,以更大的窗口过滤点云数据,统计此次迭代过后的地面点数量,记为c2。

为了验证迭代处理的效果,选取ISPRS提供的Samp12进行实验,实验结果如图2所示,在原始图像2(a)中,黑色区域为地面点集,白色点云为建筑物、植被等非地面点;图像2(b)、图像2(c)和图像2(d)中,白色点云为地面点,黑色区域为滤除的非地面点。明显看出图像2(c)框内存在未滤除的非地面点,2(d)中框内的部分地面点被错误滤除,而依据迭代终止条件停止迭代的效果较好,从而说明在此次迭代停止是恰当的,即根据迭代终止条件而非最大窗口尺寸来判断算法是否停止是有效的。

图2 迭代处理滤波结果Fig.2 Iterative processing filtering results

2.2 网格内点云精分类

如前所述,地面点网格中点的归属问题需进一步讨论。若地面点网格内仅有一点落入,该点毫无疑问为地面点,如若多个点落入同一地面点网格,除去最低点,其余点应根据判别条件进行精确分类。通过坡度阈值来判别这类网格内点的类别,以网格内最低点为参考点,计算参考点与该网格中其他点的坡度,认为大于坡度阈值的为非地面点而小于坡度阈值的为地面点。考虑到同一属性点云的空间连续性,相邻同属点间的坡度较小,坡度阈值ths可在0.8～1.2中取值,ths即为同一属性相邻点的最大坡度,不同属性相邻点的最小坡度,某点的坡度计算公式为

(1)

其中:x、y、z为某点的空间三维坐标;xmin、ymin、zmin为当前网格中最低点坐标。

利用ISPRS提供的Samp21数据验证坡度阈值设置效果,点云分类效果对比见图3。图3(a)为Samp21原始点云数据,其中有桥梁、房屋、植被等非地面点云。图3(b)框中悬浮点云应为非地面点却被误判为地面点,主要原因为该处点云密度较大,网格中往往存在多个激光点,其中最低点为地面点,其余点大部分为植被点,未设置坡度阈值的传统渐进形态学滤波将这些点均视为地面点。而设置坡度阈值的改进算法将这些点进一步分类,滤除了非地面点,滤波效果见图3(c)。

3 实验结果与分析

3.1 实验数据与评判手段

为验证改进算法对各种类型地物去除的有效性,使用ISPRS网站提供的LiDAR点云数据对该算法进行测试,该网站提供了城市与乡村点云数据各4景,并选取样本进行了人工判读,点云被精确分为地面点与非地面点。研究选取了Samp12、Samp24、Samp31、Samp42、Samp52、Samp53以代表各类地物,并采用ISPRS提出的误差评判标准对滤波结果定量分析,滤波误差定义见表1。

图3 点云分类效果对比Fig.3 Comparison of point cloud classification effect

表1 滤波误差定义

在表1中,a为地面点被正确分类的个数;b为地面点被误判为非地面点的个数;c为非地面点被误判为地面点的个数;d为非地面点被正确分类的个数。

3.2 滤波结果与定量分析

统计上述Samp12、Samp24、Samp31、Samp42、Samp52、Samp53样本的分类结果,计算各自第Ⅰ类误差、第Ⅱ类误差和总误差,具体分析见表2,其中G表示地面点;N表示非地面点;E代表误差大小。并将该算法的误差与ISPRS提供的8种经典点云滤波算法的误差进行比较,各种算法第Ⅰ类误差、第Ⅱ类误差和总误差见图4。

表2 滤波结果定量分析

由图4可见,新算法在第Ⅰ类误差和总误差方面表现较好,与Axelsson算法和Pfeifer算法接近,而第Ⅱ类误差较大,但3类误差均在10%左右,具有较好的鲁棒性。根据经验,被误分为地面点的非地面点云往往高于真实地面点云,在后续处理中,可人工干预去除这种分类错误,所以可适当增加第Ⅱ类误差以减少第Ⅰ类误差。

4 结论

研究在分析渐进形态学滤波的基础上,提出了一种不用设置最大窗口阈值的滤波方法,并进一步对网格内点云精确分类。该方法减少了参数的设置工作,不用根据不同的研究区域来设置不同的窗口阈值,具有一定的自适应性,提高了滤波算法的自动化程度,且通过去除地面点网格内的悬浮地物点,在一定程度上提高了点云分类精度,在第Ⅰ类误差和总误差方面精度较高,但第Ⅱ类误差依然较大。在后续工作中,会尝试通过改进高差阈值进一步减小第Ⅱ类分类误差。