一种改进的LiDAR数据多级曲面滤波方法

王逸超，吴小忠，陈仲伟，彭康博

(国网湖南省电力有限公司经济技术研究院，长沙 410003)

0 引言

机载激光雷达(light detection and ranging，LiDAR)作为一种主动遥感测量技术，它可以不受光照等条件的影响获取高精度的密集点云，近年来得到了快速发展，已广泛应用于DTM生产、城市建模、电力线检测等各种领域[1]。系统获取的三维点云数据中除了地面点通常还包括建筑物、树木、低矮植被、桥梁等非地面点，分离地面点和非地面点的过程称为点云滤波，对于大多数应用，点云滤波是机载LiDAR点云数据处理的核心步骤，会直接影响建模和特征提取的精度。点云滤波通常通过设定阈值来进行地面点和非地面点的判断，大多都假设地形坡度变化在一定范围内[2]。滤波阈值选择的标准一般为坡度、高差距离、点到面的距离，数据组织形式可以基于单点的、格网、三角网等。

目前的滤波方法按照概念大致可分为3种，基于坡度的方法[3-4]、基于形态学的方法[5-6]和基于参数曲面的方法[7-9]。Vosselman提出了一种基于坡度的滤波方法，通过比较激光点和周围邻近点的坡度，若小于一定的坡度阈值，则将其当做地面点[3]。Susaki提出一种自适应坡度的滤波方法[4]。但基于坡度的方法随着坡度的增大滤波的效果也会变差。Zhang等提出了一种渐进形态学滤波方法，根据地形坡度在不同窗口尺寸内设定不同的高差阈值[5]，但其阈值由固定的坡度计算得到，在地形突出的地方容易引起误判。Hui等提出了基于多级克里金插值的形态学滤波方法[6]。Axelsson提出了渐进加密三角网滤波方法，并应用于TerraSolid软件[7]。Zhao等提出了一种改进的渐进加密三角网滤波算法用于植被地区的滤波[8]。加密三角网滤波方法适用范围较广，但其实现复杂，参数设置多，时间复杂度较高。Kraus和Pfeifer提出一种迭代线性最小二乘插值算法用于植被地区的点云滤波，随后扩展了该算法用于城市地区建筑物和树木的过滤[9]。此外，目前已有一些利用多回波和波形信息进行滤波的一些研究，不同地物的反射回波具有不同的回波宽度，因此可利用回波宽度作为阈值或权重辅助进行滤波[10-11]，但其在实际生产过程中的作用还待进一步验证。

可以看出，大多数算法并非一次完成滤波操作，大多通过迭代的方式进行。而相比只利用点云高差的算法，利用局部平面(曲面)的算法更加鲁棒，滤波效果也更好，但是算法的效果在不同的地形上差别很大，通常都不太能够处理复杂城市地区和地形不连续的地方[12]。另外，大多数算法都假设最低点为地面点，地面点若包含低位粗差点则会导致生成的DTM在局部下降。由于低位粗差点分布不规则，在无人工干预的情况下很难精确剔除所有的粗差点[6]。

针对上述问题，本文提出了一种改进的多级曲面滤波方法，该方法结合了形态学滤波和多级插值滤波的思想。在粗差点去除阶段采用一种简单快速有效的基于坐标直方图统计的粗差点检测方法。在种子点选择阶段，对初始种子点进行形态学运算进一步筛选得到可靠的拟合点。通过建立不同分辨率的数据集，由粗到细逐级迭代得到最终的地面模型。

1 改进的多级移动曲面滤波算法

算法的总体流程图如图1所示。包括粗差点检测、网格划分、曲面拟合，自适应阈值设置等步骤。在粗差点检测阶段本文提出了一种新的粗差点检测方法，可以在保留正常点的同时快速有效的去除低位粗差点。在曲面拟合滤波阶段，通过增加边界和空白区域的插值处理以处理边界和数据缺失地区的滤波运算，对初始选择种子点进行一次形态学开运算得到更为可靠的地面点种子点。

图1 总体流程图

1.1 粗差点检测

机载LiDAR测量得到的点云数据通常会包含粗差点，一般为高程异常点。高位粗差点通常可在滤波过程中去掉，而低位粗差点对滤波效果影响很大。本文采用一种统计高程直方图的方式来快速有效去除低位粗差点，低位粗差点通常数量较少，且高程上不连续，通过设定一定的高程间距统计所有点的高程分布，可以计算得到点云的高程直方图。直方图横坐标为当前块内点云的高程范围，纵坐标为落在各个区间上的点云个数。根据数据区域的地形起伏情况可以对数据进行分块处理，分别统计每块内点云的高程值直方图，若数据区域内有异常点，则其对应的直方图区间的点云个数数量较少，与其他正常点的区间分布存在明显差异，设定低位粗差点的判别条件如下：若第i个间隔的点云个数小于一定阈值，且位置i小于最大点云个数所在间隔位置，则该间隔内的点为低位粗差点。

图2显示了低位粗差点对高程直方图分布的影响，图2(a)显示了无低位粗差点的高程直方图分布，图2(b)显示了有粗差点的高程直方图分布。可以看出低位粗差点会使直方图分布不连续，最大值位置偏移，通过设置合适的阈值可快速检测去除低位粗差点。图3显示了使用本文的粗差点检测方法与统计平均点间距法去除粗差点的效果对比。本文方法中高程间距值为0.5m，最小点云个数为10，统计平均点间距方法中临近点搜索个数为6，标准差系数值为1.0。可以看出本文基于统计高程直方图的粗差点检测算法可以有效地去除低位粗差点，而统计平均点间距的方法在去除低位粗差点的同时也会去除部分有效的激光点。

图2 不同区域点云高程直方图分布

图3 粗差点检测方法对比

1.2 曲面拟合滤波

去除粗差点后，将点云数据按照网格大小划分网格，每个网格内存储落在网格内的点的索引。将网格内的最低点作为初始种子点，通过形态学开运算和空白区域插值得到可靠种子点。算法首先设置较大的网格大小，得到一个粗略的地面模型，然后逐渐缩小网格大小不断细化，得到最终的地面模型。曲面拟合滤波的步骤如下：

①建立网格：根据当前网格大小建立网格，计算每个点的格网号并在格网内存储该点的索引。对每个当前网格，寻找3×3的临近窗口范围内的最低点，对每个网格内的最低点进行一次形态学开运算以保证其尽可能的为地面点。在边界和空白窗口处进行插值处理保证临近窗口的完整性，以更好地拟合地形，同时避免种子点个数太少导致曲面拟合失败的情况。其中，高程值插值为周围窗口内的最低点，坐标插值为该网格中心点坐标位置。

②通过步骤1得到的9个地面点，建立二次多项式曲面方程，拟合方程如下：

z=a0+a1y+a2y2+a3x+a4xy+a5x2

(1)

二次多项式曲面有6个参数，观测量通常多余6个，可通过建立误差方程：

v=Bx-l

(2)

最小二乘可以求解得到曲面参数向量。

③计算高差阈值，地面点和非地面点的划分是由激光点到拟合曲面的距离决定的，不同窗口尺寸和不同的地形变化需要设置不同的高差阈值。通常，地形起伏较大的区域，拟合曲面和真实地形有一定的偏差，因此需要设置较大的阈值，反之，则需要较小的阈值。另外窗口的尺寸越大，地形的起伏的可能性就越大，因此本文考虑窗口尺寸和曲面拟合均方根误差的影响，窗口尺寸和均方根误差越小，拟合的曲面越接近于真实地形，高差阈值就需要设置得相对较小，反之需要设置较为宽松的高差阈值。阈值设置如下：

T=Tmin+Wratio*(Tmax-Tmin)

(3)

式中：Wratio为窗口尺寸和高程方差归一化加权之后的值；Tmin、Tmax为设置高差阈值的最小、最大值。

④更新网格大小，窗口尺寸的选择对于滤波效果的影响同样很大。窗口尺寸的选择，一般有线性增长和指数增长2种方式，线性增长方式可以很好地保持地形特征，但是会增加滤波的迭代次数。考虑到点云一般数据量较大，以指数方式增长窗口大小可以降低大范围地区点云滤波的时间。通常，点云的平均点间距可以作为窗口尺寸的最小值，而窗口尺寸的最大值，则应接近于建筑物的大小，为了保证窗口尺寸最大时，窗口内的最低点为地面点，从而使大型建筑物在滤波过程中去除。滤波过程中每一级窗口的大小可以由初始值和增长方式计算得到。

2 数据及实验分析

2.1 实验数据

ISPRS委员会提供了用于滤波实验的基准数据集，参考数据由人工滤波操作生成。数据集由Optech ALTM扫描仪在Vaihingen实验区和Stuttgart市区的8个地区采集得到。数据集共包含15个样区，其中城市地区点云平均间距为1.0～1.5 m，农村地区平均点间距为2.0～3.5 m。样区地物类型丰富，包括大型建筑物、低矮植被、陡坡等，本文选择其中的4个典型样区进行实验验证。样区12带有小型目标物(汽车等)；样区21具有狭窄的桥梁；样区31地形不连续且带有低位粗差点；样区52为带有低矮植被，不连续且坡度较大的山地。

2.2 实验结果与分析

点云滤波效果的定量评价指标通常是计算滤波的三类误差，I类误差、II类误差和总误差。其中I类误差为所有地面点中地面点被当做非地面点的比例，II类误差为所有非地面点中被当做地面点的比例，总误差为分类错误的点在总点数中的比例。总误差反映了滤波算法的总体性能。我们在VS2015平台上利用c++语言实现了本文算法，并计算I类误差、II类误差和总误差作为滤波精度的指标，与其他方法进行分析比较。

表1显示了本文方法对不同测区数据滤波结果与其他8种经典的滤波方法的对比结果。其他滤波方法结果数据来源于文献[12]及其附录。从表中可以看出，本文方法的滤波的总误差整体较小，且优于大部分方法，尤其在地形平坦的地区。不同地区的滤波精度差别较大，这跟测区地形的复杂度有关。

表1 滤波结果比较

其中测区12、21、31地形相对平坦，坡度变化不大，滤波效果较好，测区21、31的总误差在1.5%左右，在所有滤波算法中误差最低。随着坡度增大，滤波误差增大，这与现有的滤波方法的结论保持一致。测区52的总误差为9.23%，相比于其他地形平坦的测区，总误差增大了不少，但仍然好于大多数算法。可能的原因是由于在地形起伏较大的地方，拟合曲面与实际的地形表面存在差异，且受窗口大小的影响较大，相对来说，基于三角网的滤波效果要好一些，可能的原因是三角形比格网更适合表达起伏地形。

值得注意的是，没有一种滤波算法在所有地形情况下滤波效果都达到最佳，用户需要根据地形情况选择合适的算法。另外参数的设置也会影响最终的滤波精度，实验表明，高差阈值的最小值通常在0.3 m到0.5 m之间，最大值在1.5 m到2 m之间。经过两到三次迭代通常会达到比较好的滤波效果。但是通过设置参数不可能完全消除误差，I类误差和II类误差具有相关性，I类误差减少，II类误差就会增大，反之亦然。因此，用户必须在减少I类误差和II类误差之间做出选择，对于在后续的质量控制过程中，II类错误要比I类错误更容易发现和被人工修复。但是对于滤波算法来说，大多都会尽量减少II类误差，因为一些凸起或者较低的非地面点对于自动化地形重建的效果影响很大。

图4显示了参考数据和本文算法进行点云滤波的实验结果及地形重建的效果图，图4(a)为标准的滤波参考数据集，其中红色点为非地面点，蓝色点为地面点。图4(b)显示了滤波后的地面点，对比图4(a)可以看出滤波算法在去除了非地面点的同时保留了绝大部分的非地面点。图4(c)显示了滤波后的地面点通过反距离权重法利用Surf软件进行插值重建的地形效果，可以看出重建的地形效果接近于真实的地表形状。

图4 点云曲面拟合滤波

3 结束语

针对机载LiDAR点云滤波处理中存在的一些问题，提出并实现了一种改进的多级移动曲面点云滤波方法，该方法结合了形态学运算和曲面插值方法的优点，兼顾了算法运行速度和精度。实验表明，该方法能够处理各种地形的数据，在平坦和地形起伏区域均能一定程度上提高点云的滤波精度，具有良好的鲁棒性和适用性，可适用于大多数地形复杂的地区点云的滤波处理。本文的贡献在于在数据预处理阶段采用了一种基于统计直方图的方法可以快速有效地去除低位粗差点，在地面种子点选取的过程中结合形态学开运算以获取更加可靠的地面点，同时增加了对边界和空白区域的插值处理以及自适应窗口阈值的设置。如何进一步提高算法在地形坡度剧烈变化地区的滤波精度是本文下一步研究的重点，利用多源数据和额外的辅助信息进行滤波是进一步提高滤波精度的一个研究方向。