基于机载LiDAR的点云数据滤波方法

黄代军

(浙江省测绘科学技术研究院 浙江杭州 310030)

机载LiDAR技术经过20多年的发展，其硬件设备已经取得了长足的进步。作为主动式数据采集手段，机载LiDAR技术已经广泛应用在场景数据采集中，如建筑物数据采集、电线巡检数据采集和中小型区域地形植被数据采集等。相比于传统的航空摄影测量，机载LiDAR不仅抗干扰能力强，而且采集的点云数据信息丰富，可以全天候作业。点云数据的优点在于可以真实反映场景中地物的细节信息，如强度信息和位置信息等，但也存在一些缺点，如采样粒度不均匀、冗余高和数据量大等。

点云数据的后处理问题一直是测绘行业研究的难题之一，相比于点云采集的硬件设备，点云数据的处理技术仍然相对滞后，其中点云数据的滤波更是点云数据处理中的关键点和难点之一[1]。经过几十年的发展，许多专家学者提出并验证了多种点云滤波算法。目前，常用的点云数据滤波算法主要包括基于不规则三角网滤波算法、移动窗口滤波算法、数学形态学滤波算法和移动曲面拟合滤波算法等。本文在传统的移动曲面拟合滤波算法的基础上，引入偏度概念，对原始算法进行了改进，并采用国际摄影测量与遥感学会提供的数据，通过试验验证改进算法的可行性，同时与原始算法进行对比，判断两种算法的效果。

1 常用的机载LiDAR点云数据滤波方法

1.1 基于不规则三角网滤波算法

Terrasolid软件使用基于不规则三角网滤波算法。基于不规则三角网滤波算法的机理：首先，对所有数据进行计算得到初始参数，再通过选取地面种子点构成不规则三角网；其次，设置条件，通过向不规则三角网中添加满足条件的地面点来加密不规则三角网；最后，当地面没有满足加入三角网的点集时，结束运算。

基于不规则三角网滤波算法仅能计算一些连续表面环境的点云，具有一定的局限性，对于城市中的一些不是连续表面中的点会因为不能满足算法设置的条件而被遗漏。为了解决这一问题，李卉等[2]使用镜像的方法将遗漏的点加入到不规则三角网中。此方法虽有一定的可行性，但是理论依据不足。为此，他们在算法中融入了一种融合区域增长思想，大大减少了漏选点。

1.2 移动窗口滤波算法

移动窗口滤波算法的机理：首先，在某个范围内取一个比较大的窗口，找到最低点；其次，设置一定的阈值，如果最低点与窗口中最高点的高差不大于此阈值，那么可将此点视为地面点。通过不断减少窗口大小和阈值，寻找到更为精确的地面点。同时，对利用不同窗口大小得到的地面点进行配权，以加权平均的方式得到DEM成果。

移动窗口滤液算法在使用时，如果设置的窗口太大，可能导致表面不连续的点被滤掉，窗口太小有可能会将建筑物等识别为地面点。设置的阈值过小，可能会将表面不连续的点识别为非地面点；设置的阈值过大，可能会将植被等点识别为地面点。所以使用移动窗口滤波算法的关键在于参数的选取是否合适。

1.3 数学形态学滤波算法

数学形态学滤波算法在数字图像处理中的应用相当广泛，其机理是通过具有一定形态的结构实现图像中对应形状的提取，达到图像识别和图像分析的目的。

数字形态学滤波算法包括两个不同的运算流程：开运算与闭运算，通过膨胀与腐蚀两个基本运算实现。先进行腐蚀操作，再进行膨胀操作的是开运算；先进行膨胀操作，再进行腐蚀操作的是闭运算。开运算可以实现对机载LiDAR点云的数据滤波。

数学形态学滤波算法受窗口尺寸选取影响较大，在实际使用前，需要对点云数据所在环境具有一定的先验知识，才能设置合适的窗口尺寸，从而实现对机载LiDAR点云数据的滤波[3]。

1.4 移动曲面拟合滤波算法

移动曲面拟合滤波算法可以实现对数据的滤波，计算速度快且简单明了。

移动曲面拟合滤波算法的机理：首先，在种子区域内找到最低的3个角点作为初始地面点，对其进行二维排序，通过求解平面方程的方式拟合初始平面，在平面方程中加入相邻角点的平面坐标，得到拟合高程值；其次，比较拟合高程与观测高程的大小，如果二者之间的差值在给定的阈值范围内，则此相邻角点为地面点，否则为非地面点；再次，将判断为地面点的相邻角点与组成初始平面的3个角点拟合为1个新的曲面；最后，按照上述方法继续判断相邻角点是否为地面点，当判断得到的地面点数量达到6个时，在增加一个新地面点的同时去掉一个最远的角点。

重复上述操作，通过简单的地形曲面拟合方式完成对整个测区的过滤[4]。

2 基于偏度平衡的移动曲面拟合滤波算法

本文在移动曲面拟合滤波算法的基础上，引入统计学中的偏度概念。偏度(skewness)，也称偏度系数，是对变量的概率分布密度函数曲线偏离对称度的统计，用以表现概率密度曲线尾部的相对长度。

假设有一随机变量X，同时拥有三阶和四阶矩阵，那么偏度Sk可表示为[5]

(1)

标准正态分布的偏度为0，当Sk<0时称为负偏态，Sk>0时称为正偏态；Sk的绝对值越大，表明其偏离程度越大。本文引入偏度概念的目的在于剔除干扰点。

试验前，假设数据中的地面点概率密度分布函数服从标准正态分布，非地面点呈偏态分布，对数据进行干扰，即通过不断地从点云中剔除非地面点的方式使得偏态分布更接近于正态分布。同时，为了使样本的分布逐渐趋于平衡，还需要不断地对偏度值进行校正。

移动曲面拟合滤波算法实现原理较为简单，但对滤波数据的要求较高，只有在点云密度足够高的前提下才能得到较好的效果。该算法通过判断相邻角点是否为地面点的方式不断拟合新的曲面，同时不断加入新点，拟合地形不断变化，误差不断累积与传递就会导致滤波精度降低。

另外，该算法的阈值不可以随着地形起伏的变化而自动调整，导致对坡度较大的地形特征点不能进行有效分离；其次，该算法对场景中的所有点都是单独进行比较，会大大降低滤波效率。本文主要是对移动曲面拟合滤波算法中的误差累积进行改进，改进后的算法流程为[6]：

1)对机载LiDAR点云数据进行预处理后，对点云区域进行虚拟网格划分，并为不同网格中的激光点建立高程序列索引。

2)确定初始窗口尺寸，在相邻3×3窗口中选取6个高程最低点，构成点集合；然后采用偏度校正的方式剔除虚拟网格中高程最低的非地面点，剩余的点组成种子点集。

3)使用最小二乘法通过种子点拟合点云数据区块。

4)遍历场景中所有点云，计算遍历的激光角点的拟合高程与实际高程的差值Δhi，设置一定的阈值Δh。如果某个点的Δhi在阈值外，那么可认定此点为非地面点，反之则为地面点，同时初步标记地面点与非地面点。设置合适的阈值对点云滤波非常重要，直接影响滤波精度。本文阈值设定的方式为

Threshold=(Hmax-Hmin)×Scale

(2)

式中：Threshold为高差阈值，Hmax和Hmin分别为区域内点云的最大高程与最小高程，Scale为窗口尺寸。

5)扩大窗口尺寸，重复步骤3)和步骤4)，得到分类更加精细的激光角点。当窗口尺寸达到阈值时，结束算法运行，实现滤波。

基于偏度平衡的移动曲面拟合滤波算法的具体流程如图1所示。

图1 基于偏度平衡的移动曲面拟合滤波算法流程Fig.1 Flow Chart of Moving Surface Fitting Filtering Algorithm Based on Skewness Balance

3 试验验证与分析

3.1 数据来源

为了验证基于偏度平衡的移动曲面拟合滤波算法的可行性，本文采用国际摄影测量与遥感学会公布的CSite2首次回波数据及利用不同滤波方法得到的参考滤波数据作为试验数据对改进后算法进行验证与分析。

国际摄影测量与遥感学会提供的数据包含激光脉冲首次、末次的反射强度和反射位置，其中，Csite2的点密度可达到0.67点/平方米。

3.2 误差指标

滤波后的数据除了需要定性分析外，还需要使用误差指标对滤波结果进行定量分析。目前常用的误差指标有两种：

1)激光点的误分率。误分率可分为总误差、一类误差和二类误差。其中，一类误差主要是在进行地面点分类时剔除部分地面点产生的误差；二类误差是在不能有效剔除地物点时产生的误差；总误差则是一类误差与二类误差通过加权求和得到的误差。

改进后算法的适用性可利用一类误差和二类误差进行反映，通常是选择有效的方法降低二类误差，以得到更加全面的地形数据[7]。本文主要使用激光点的误分率(表1)对改进算法的效果进行验证。

表1 滤波误分率Tab.1 Filtering Error Rates参考结果滤波结果地面点个数地物点个数参考数据点个数误分率/%地面点个数abe=a+bp= b / e地物点个数cdf=c+dq= c / f滤波后点数g=a+ch=b+dm=a+b+c+dr=(b+c)/m

表1中，a和d分别为在滤波过程中正确分类的地面点和地物点的个数，b为将地面点错分为地物点的个数，c为将非地面点错分为地面点的个数，e为正确地面点的个数，f为正确地物点个数，m为场景中所有的LiDAR点个数，p为一类误差，q为二类误差，r为总误差。

2)在不同区域选取离散样本点，计算滤波后地面模型的高程值与原始数据高程值的均方根误差。将滤波前后高程值的均方根误差和高程差值作为评价指标对改进算法的滤波效果进行评价，均方根误差越小，高程差值越小，表明滤波效果越好，精度越高[8]。

3.3 试验过程

为了更加直观地反映基于偏度平衡的移动曲面拟合滤波算法的滤波效果，使用CSite2中的样本数据Samp31和Samp23对改进算法进行验证。由于样本数据Samp31和Samp23中既有原始数据又有滤波参考信息，因此可同时进行定性和定量比对。

通过MATLAB软件来运行和实现改进算法，并对传统算法与改进算法得到的试验结果进行比对。在试验过程中，滤波效果受网格大小的影响较为严重，经过试验比较，采用2 m×2 m的网格得到的效果较好。

在每个网格中选取3个高程最低点组成种子点集，并通过偏度平衡过滤掉网格中高程较高的点。由于网格面积较小，可能会覆盖到建筑物表面，从而将建筑物点加入到地面点中，在进行平面拟合时由于含有非地面点，导致误差累积。剔除高程较大值后，得到所有最低点高程，对不同最低点密度进行统计得到种子点直方图，如图2所示。

图2 种子点高频直方图Fig.2 High Frequency Histogram of Seed Points

通过种子点进行曲面拟合，以Samp31数据为例得到接近于地面点的拟合曲面如图3所示。

图3 拟合曲面Fig.3 Fitting Surface

使用改进滤波算法与基于不规则三角网滤波算法对Samp31和Samp23试验数据进行滤波，结果如图4所示。

图4 不同滤波方法试验结果Fig.4 Experimental Results of Different Filtering Methods

图4(a)是采用TerraScan软件的算法对Samp31数据进行滤波处理得到的地面点集，可以看出其算法基本分离出了地面点，视觉上可以明显看到地物点和地面点的区别，但实际是地面点中仍然含有少量地物点。

图4(b)是采用基于偏度平衡的移动曲面拟合滤波算法对Samp31数据进行滤波处理得到的地面点集，可以明显看到空白区域为建筑物所在区域，过滤掉存在的小型植被区域点云。图4(c)是采用Terrasolid软件的算法对Samp23数据进行滤波处理得到的地面点集，图4(d)是采用基于偏度平衡的移动曲面拟合滤波算法对Samp23数据进行滤波处理得到的地面点集，可以看到提取效果明显好于对Samp31数据进行滤波处理的效果，但是仍存在一类误差和二类误差。其中，图4(d)中的点云更能反映出地面区域特征，更明确地表达出建筑物及植被的信息。

3.4 结果分析

为了定量分析基于偏度平衡的移动曲面拟合滤波算法的滤波效果，本文使用误分率这一评价指标对滤波效果进行评价。

由于点云数据量大，采用人工选取数据进行滤波效果对比的方式不实际，因此，本文选择与Csite样本数据中标记的参考滤波数据进行对比的方式，分别使用移动曲面拟合滤波法和基于偏度平衡的移动曲面拟合滤波算法对Samp31数据和Samp23数据进行滤波处理，并定量评价处理效果，结果如表2所示。

表2 滤波效果Tab.2 Filtering Effects 参考结果Samp31数据滤波结果地面点/个地物点/个参考数据点/个误分率/%Samp23数据滤波结果地面点/个地物点/个参考数据点/个误分率/%地面点/个15 23674513 1905.6516 34365713 1904.98地物点/个1 19912 18816 1057.441 05311 28416 1056.54滤波后点/个16 43512 93329 2956.6417 39611 94129 2955.84

根据表2可知，移动曲面拟合滤波法和基于偏度平衡的移动曲面拟合滤波算法的滤波效果都在10%以内，一类误差小于二类误差，二类误差和总误差稍偏大，但都在可接受范围内，并且基于偏度平衡的移动曲面拟合滤波算法的一类误差、二类误差与总误差均小于移动曲面拟合滤波算法。试验结果表明，基于偏度平衡的移动曲面拟合滤波算法适用于城市地区点云滤波，结果较为可靠，精度更高。

4 结 语

本文首先介绍了几种常用的点云滤波方法，在分析移动曲面拟合滤波算法的基础上，提出了一种基于偏度的移动曲面拟合滤波算法。该算法在确定种子点的过程中加入偏度概念，以降低种子点误差，达到在上一次滤波完成之后继续对非地面点和地面点进行滤波，有效分离非地面点与地面点，从而提高滤波精度的目的。同时，使用国际摄影测量与遥感学会提供的数据对改进算法进行验证，结果表明，基于偏度的移动曲面拟合滤波算法的一类误差和二类误差都在10%以内，且滤波效果优于移动曲面拟合滤波算法，可以减少一类误差、二类误差和总误差。基于偏度的移动曲面拟合滤波算法适用于城市地区点云滤波，可取得较为可靠的滤波效果和更高的滤波精度。