机载LiDAR点云数据的多级滤波方法*

时培强，江 虹

（西南科技大学 信息工程学院，四川 绵阳 621010）

0 引 言

激光雷达（Light Detection and Ranging，LiDAR）是一种应用全球定位系统（GPS）﹑惯性测量单元（IMU）﹑激光测距仪和计算机技术于一体的系统[1]。与传统航测相比，机载LiDAR具有精度高﹑速度快﹑受外界因素影响小﹑自动化程度高等技术优点，且能提供高精度﹑高密度的三维信息[2]。LiDAR在三维数字城市建模﹑城市规划﹑海岸线检测﹑电力和公路选线巡线﹑森林资源的管理和评估等领域，具有广阔的应用前景[3]。

激光雷达数据由地面点和地物点组合而成。滤波是指从离散的点云数据中将地面点和非地面点进行分离[4]。由于点云数据呈现出的无规律性和不确定性，滤波一直是激光雷达数据处理的首要问题，且其精度直接影响后续工作的准确性。

典型算法主要分为：三角网滤波算法[5-6]﹑曲面拟合滤波算法[7-10]﹑数学形态学滤波算法[11]﹑移动窗口法和基于坡度的滤波算法[12]等。ISPRS第Ⅲ委员会对不同的滤波算法进行了全面分析，结果表明：多数滤波算法对特定地形可以得到理想的效果，但对于多种复杂地形滤波效果不佳[13-14]。三角网滤波算法整体适应性较好，但计算量和复杂度相对较大，且会丢弃部分地面点；曲面拟合算法对地形特征保护较好，算法相对简单，但种子点的不当选取会影响滤波精度；数学形态学算法简单易实现，计算效率高，需要根据实际数据确定窗口大小；移动窗口算法的滤波适用性相对较差；基于坡度的滤波算法主要依赖初始坡度的设定。实际研究中发现，可以将均值限差法﹑曲面拟合算法与三角形角度限制法结合，通过偏度和峰度的变化对阈值进行相应调整，从而提高滤波效果。该方法以均值限差法获取可靠性较高的地面点为基础进行曲面拟合，通过对不断更新的区域进行曲面拟合，完成“粗分类”“细分类”滤波，然后将剩余点与邻近已确定的地面点通过三角形角度限制法判断，利用各方法的优点进行逐级滤波，最终使滤波结果更加合理。

1 滤波算法的基本思想

本文的滤波算法是一个逐级滤波处理过程，思路为：区域广且点云密度大时，如果对整体进行处理会影响滤波的效率和准确性，所以本文采用数据格网化，并进行区域分块。首先对点云数据进行格网化分，其次对整个区域进行分块，分块后的各区域应包含地面点，最后将分块后的区域再分成均等的若干小区域，在小区域内选择最低点，并利用均值限差法判断当前最低点是否为地面点。如果当前区域选择的地面点数大于二次曲面拟合方程求解系数个数，则对所选的地面点进行曲面拟合，并利用拟合曲面对此区域的点云数据实施Ⅱ类误差最小化的“粗分类”过程。通过对加入点前后两次偏度和峰度的计算，根据偏度和峰度的变化对该区域中的阈值进行动态调整。当满足地面阈值条件时，将点加入到地面点集合；反之，当满足地物阈值条件时，将点加入地物点集合。

经过滤波后，大部分建筑物﹑植被和高程极值点被滤除。然后，分别对小区域以“粗分类”的地面点结果为基础选择下一级地面点。在小区域内进行曲面拟合，利用拟合曲面方程对小区域进行“细分类”。根据偏度和峰度的变化对阈值进行相应调整。当满足地面阈值条件时，将点加入地面点集合；反之，当满足地物阈值条件时，将点加入到地物点集合。遍历所有分块，直到整个区域被覆盖，从而实现点云数据的初次分类。最后，对部分剩余的点通过与其邻近且已确定的地面点利用三角形角度限制判别法再次进行分类。角度阈值的设定需要根据峰度和偏度的变化作相应调整，将判断点之间形成的角度与设定的角度阈值进行比较。当满足地面阈值条件时，将其加入地面点集合；反之，加入到地物点集合，从而实现点云数据的整体分类。滤波流程如图1所示。

图1 滤波算法流程

2 算法实现

2.1 偏度和峰度

统计学中，从样本的概率密度曲线直观看来，偏度（Skewness）就是其尾部的相对长度[15]，峰度（Kurtosis）就是其尾部的厚度。通常，偏度和峰度分别用sk和ku表示。若某一随机变量X的三阶矩和四阶矩同时存在，则计算公式为：

其中N表示样本总数，xi为随机样本点，σ和μa分别表示样本的标准差和算术平均值，定义为：

偏度和峰度均是为无量纲的量。若sk＞0，则称该分布为正偏态或者右偏态；若sk＜0，则称该分布为负偏态或者左偏态。|sk|越大，表示其偏离程度越大。类似地，若ku＞0，则该分布比较陡峭；若ku＜0，则该分布比较平坦。偏度和峰度均是相对于正态分布来比较的。正态分布的偏度和峰度均为0。偏度和正态分布如图2所示。本文利用偏度﹑峰度的变化对阈值进行相应调整，通过比较新增加点云数据对偏度和峰度的影响，定位到新增点云的类型。如果前后偏度差小于0，表示新增点接近均值，更偏向于是地面点。如果偏度差大于0，表示新增点偏离均值较大，更偏向于地物点。如果前后峰度差大于0，表示新增点云与当前点云数据偏差不大，偏向于是地面点。如果峰度差小于0，表示新增点云与当前数据偏差较大，偏向于地物点。通过偏度﹑峰度的变化量对阈值进行调整，将高程差值与阈值进行比较，决定新增点云是地面点还是地物点，以提高滤波的精度。

图2 不对称分布及正态分布

2.2 点云数据格网化

数据格网化是将所测区域的点云数据用规则的格网全部覆盖，如图3所示。

图3 格网化分

首先统计所测区域点云总数N，并得到X方向的最小值XMin﹑最大值YMax以及Y方向的最小值YMin﹑最大值YMax。所需矩形区域的面积大小为Area=(XMax-XMin)(YMax-YMin)。由点云总数N和区域面积Area计算点云密度，根据点云密度确定格网间隔，使点云数据分布均匀。假设根据以上信息得到规则格网的大小为m×n，则某点的坐标信息(Xk,Yk,Zk)所对应的行列号(i, j)为：

如图4所示，根据某点（图中标注点）的坐标信息可以得到所在格网位置。每个格网需要记录该格网内包含的点数和当前点所对应的序列号，通过序列号可以快速索引到该点包含的信息。

图4 格网索引

2.3 均值限差法

均值限差法是以选中点为中心，并设定窗口大小，计算窗口内所有点的均值与中心点之间的差值，通过比较差值与设定的阈值判断中心点的取舍。本文根据非地形点的高程大于地形点的高程且非地形点的高程一般偏离均值较大的特性，计算此区域内的点云高程平均值，然后计算选中点与平均值之差，通过与设定的阈值进行比较，与平均值高差大于设定阈值的点被视为非地面点，如图5所示。由于点云数据包含一些高程极高点或极低点，而这些极点并不是真实数据，会影响滤波精度，需剔除这些误差点。本文采用此方法通过剔除极低点来选择初始地面种子点。首先选择局部区域的最低点为中心点，以此点为中心计算其周围若干格内的点云高程均值，然后得到中心点高程与均值的差值。通过比较差值与设定阈值确定此点是否为地面点，从而排除极低点的干扰。此方法可以有效得到最初地面种子点，从而为曲面拟合提供准确﹑可靠的地面点。

图5 均值限差法

2.4 曲面拟合算法

由于地形表面比较复杂，可以看作是一个有多个独立的曲面组合而成的空间曲面[16]。二次曲面拟合对地形特征保护较好，通过部分地面点可得到局部地形的总体特征，利用地形特征能够有效辨别邻域内的点是否属于地面点。每一个独立的曲面可以用一个二次曲面进行拟合，其拟合方程为：

式（8）中，(xn, yn, zn)为点云n 的3维坐标值，式（9）中X为二次多项式的系数，即待求参数。式（7）拟合方程用矩阵表示为Ax=b，其一般包含m个等式和n个未知数，且m＞n。该方程组一般无解，为了得到尽可能合适的X且使等式尽量成立，引入残差平方和函数：

式（11）中s(x)取最小值：

通过对式（11）中的s(x)进行微分求最值，得到：

当方程个数大于方程系数个数时，由此求出方程系数，将相应区域的点云数据代入此方程等到相应的高程值。同时，更新所在区域的偏度和峰度，通过偏度和峰度的变化量相应调整阈值。原始点云高程和拟合曲面高程的差值小于地面限差值的点被视为地面点；反之，差值大于地物限差值的点被视为非地面点。

当高程差值不符合限差条件时，当前点需要进行下一步的细分类，如图6所示。当所测区域的窗口较大时，此时偏度和峰度的变化量相对较大，对阈值的设定权重需要降低，设置的限差值应较小；反之，当所测区域的窗口较小时，偏度和峰度的变化量相对较小，对阈值的设定权重较大，设定的限差值应较大。阈值的动态设置可以提高算法的适应性。

图6 曲面拟合

2.5 三角形角度限制法

当水平距离较近时，垂直高度差越小，形成的空间立体角越小；反之，水平距离较近，垂直高度差越大，形成的空间立体角越大，如图7所示。

图7 角度限制法

文中通过偏度和峰度的变化量对角度阈值进行相应调整。当前判断点与邻近已确定地面点构成三角形。当二者构成的角度大于设定的角度阈值时，则当前判断点被视为非地面点；反之，当角度小于设定的角度阈值时，则当前点被视为地面点。空间立体角的计算公式为：

3 滤波实验

3.1 实验数据和评价标准

采用ISPRS第Ⅲ委员会在2003年提供的数据对本文算法的适应性﹑可行性进行验证。本文选取包含高大建筑物﹑地物﹑低矮植被等分布的地区（Sample12），包含激光点云总数52 119个；另一区域包含高大建筑物﹑植被﹑建筑物与植被混合区（Sample31），包含激光点云总数28 862个。

滤波评价标准包含Ⅰ类误差﹑Ⅱ类误差以及总误差三种。Ⅰ类误差是指将地面点误判为非地面点的比率，Ⅱ类误差是指将非地面点误判为地面点的比率，总误差是由Ⅰ类误差与Ⅱ类误差加权求和得到的[17]。其中，Ⅰ﹑Ⅱ类误差表现算法的适应性，总误差反映算法的可行性。

3.2 滤波结果和分析

利用本文所提的算法分别对Sample12区域和Sample31区域进行滤波处理。可以发现，在城市区域滤波效果较好，在地形起伏较大的山区或地物点与地面点高程相近处出现部分分类误差。如图8﹑图9所示，通过对滤波效果进行分析，可以得出以下结论：

（1）对于平坦的道路，本文算法能够比较准确地从点云数据中提取出道路附近的地面点。由于道路起伏较小，曲面拟合可以获取较为准确的原始地形信息。当道路旁低矮的植被与地面点接近时，会将部分低矮植被点误分为地面点。在保证Ⅱ类误差最小化的前提下，通过偏度和峰度的变化对阈值进行相应调整，能够比较准确地分离出地面点。

（2）对于高大建筑物和植被，本文算法能够从点云数据中将绝大部分建筑物和植被点提取出来，并可对建筑物和植被点进行准确分类。对于建筑物的起始边缘或低矮植被与地面接近的情况，此时会影响部分建筑物边缘点和低矮植被的误分类，如图10所示。随着相差距离的增大，分类结果会越来越准确。

（3）对于建筑物与植被混合区，本文算法能够将大部分点云正确分类，但仍存在部分误分类情况，如图11所示。在建筑物与植被之间会存在少量的地面点或低矮植被。对于此部分点云，在经过粗分类和细分类后不能准确分类。根据此区域的偏度和峰度的变化对角度阈值进行调整，对此部分点云使用角度限制法可以有效提高分类效果。和峰度的变化不明显，会将部分低矮植被点误分为地面点，从而影响了Ⅱ类误差。为了降低Ⅱ类误差，本文将地面点阈值区间设置较小。当地形出现明显突变时，不可避免会将此部分地面点误分为地物点，从而增加了Ⅰ类误差率。

图8 Sample12区域滤波结果对比

图9 Sample31区域滤波结果对比

图10 Sample12区域局部误差

表3是本文算法的误差与部分经典算法的误差之间的对比。可以发现，本文算法在测试区域总误差处于中等以上水平，可知所提算法具有一定的可行性，并能够得到可靠的滤波效果。

图11 Sample31区域局部误差

从表1和表2可以看出，通过对滤波误差进行分析，可以得出以下结论：

（1）不同的测试区域虽然三类误差有所区别，但三类误差中，Ⅱ类误差在不同测试区域都是最低的。这验证了本文算法是在降低Ⅱ类误差的前提下，进一步提升了整体的滤波精度。

（2）Ⅱ类误差出现的主要原因是地面上有低矮植被，当低矮植被与地面点高程差很小时，偏度

表1 Sample12区域滤波结果统计

表2 Sample31区域滤波结果统计

表3 滤波结果比较

4 结 语

本文提出一种多级相结合的LiDAR点云数据滤波算法。主要工作和结论如下：

（1）以均值限差法获取可靠性较高的初始地面种子点为基础，提高曲面拟合的准确性。以曲面拟合算法对数据进行滤波，能有效剔除地物点，同时保留原始地貌，且算法简单。

（2）以三角形角度限制法对剩余点进行判断，可以很好地利用坡度信息和已经确定为地面点的数据进行处理，进一步提高滤波的精度。

（3）通过实验验证，该滤波算法通过偏度和峰度的变化对所在区域的阈值进行相应调整，对含有高大建筑物和植被的城市区域有较好的适应性，且取得了理想的滤波效果。

尽管本文算法的可行性和实用性较高，偏度和峰度的变化在城市区域有着较高的适应性，能够对阈值进行相应调整，但对于地形起伏较大的山区效果不显著，滤波精度较低。此外，区域分块需要根据建筑物的大小划分，自适应性也需要进一步完善。以后的研究需要进一步对算法进行改进，并对已确定的地物点进行插值处理，以便更好地保持地形原貌，也可以利用回波强度信息进一步提高滤波精度。

[1] 赖旭东.机载激光雷达基础原理与应用[M].北京:电子工业出版社,2010.

LAI Xu-dong.Basic Principles and Applications of Airborne Lidar[M].Beijing:Publishing House of Electronics Industry,2010.

[2] 张小红.机载激光雷达测量技术理论与方法[M].武汉:武汉大学出版社,2007.

ZHANG Xiao-hong.Airborne Lidar Measurement Technology Theory and Method[M].Wuhan:Wuhan University Press,2007.

[3] 张玉方,程新文,欧阳平等.机载LIDAR数据处理及其应用综述[J].工程地球物理学报,2008,5(01):119-124.

ZHANG Yu-fang,CHENG Xin-wen,OUYANG Ping,et al.The Data Processing Technology and Application of Airborne LIDAR[J].Chinese Journal of Engineering Geophysics,2008,5(01):119-124.

[4] 刘经南,许晓东,张小红等.机载激光扫描测高数据分层迭代选权滤波方法及其质量评价[J].武汉大学学报:信息科学版,2008,33(06):551-555.

LIU Jing-nan,XU Xiao-dong,ZHANG Xiao-hong,et al.Adaptive Hierarchical and Weighted Iterative Filtering of Airborne LIDAR Data and Its Quality Assessment[J].Geomatics and Information Science of Wuhan University,2008,33(06):551-555.

[5] Axelsson P.DEM Generation from Laser Scanner Data Using Adaptive TIN Models[J].International Archives of Photogrammetry and Remote Sensing,2000,33(B4/1;PART 4):111-118.

[6] 曾妮红,岳迎春,魏占营等.车载LiDAR点云滤波的改进不规则三角网加密方法[J].测绘科学,2016,41(09):136-139.

ZENG Ni-hong,YUE Ying-chun,WEI Zhan-ying,et al.An Improved Irregular Triangular Network Encryption Method of Vehicle-borne LiDAR Point Clouds[J].Science of Surveying and Mapping,2016,41(09):136-139.

[7] 张小红.机载激光扫描测高数据滤波及地物提取[D].武汉:武汉大学,2002.

ZHANG Xiao-hong.Airborne Laser Scanning Altimetry Data Filtering and Features Extraction[D].Wuhan:Wuhan University,2002.

[8] 熊俊华,方源敏,邓德标.最小二乘曲面拟合的LiDAR数据滤波方法[J].测绘科学,2013,38(04):74-76.

XIONG Jun-hua,FANG Yuan-min,DENG De-biao,et al.Surface Fitting Filtering based on Least Square Method from LiDAR Data[J].Science of Surveying and Mapping,2013,38(04):74-76.

[9] 孙崇利,苏伟,武红敢等.改进的多级移动曲面拟合激光雷达数据滤波方法[J].红外与激光工程,2013,42(02):349-354.

SUN Chong-li,SU Wei,WU Hong-gan,et al.Improved Hierarchical Moving Curved Filtering Method of LiDAR Data[J].Infrared and Laser Engineering,2013,42(02):349-354

[10] 张皓,张永生,刘军等.一种基于平面拟合的LIDAR点云滤波方法[J].测绘科学,2009,34(04):141-143.

ZHANG Hao,ZHANG Yong-sheng,LIU Jun,et al.A Method for Filtering LIDAR Points Cloud based on Planar Fitting[J].Science of Surveying and Mapping,2009,34(04):141-143.

[11] 董保根,秦志远,朱传新等.关于机载LiDAR点云数据形态学滤波的几点思考[J].测绘科学,2013,38(04):23-25.

DONG Bao-gen,QIN Zhi-yuan,ZHU Chuan-xin,et al.Perspectives of Morphological Filtering for Airborne LiDAR Point Clouds Data[J].Science of Surveying and Mapping,2013,38(04):23-25.

[12] Vosselman G.Slope based Filtering of Laser Altimetry Data[J].International Archives of Photogrammetry and Remote Sensing,2000,33(B3/2;PART 3):935-942.

[13] Sithole G,Vosselman G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds[J].Isprs Journal of Photogrammetry& Remote Sensing,2004,59(01-02):85-101.

[14] Sithole G,Vosselman G.Filtering of Airborne Laser Scanner Data based on Segmented Point Clouds[C].Laser Scanning,2005.

[15] 董保根,秦志远,陈静等.无需阈值支持的机载LiDAR点云数据滤波方法[J].计算机工程与应用,2013,49(15):219-223.

DONG Bao-gen,QIN Zhi-yuan,CHEN Jing,et al.Threshold-free Method for Airborne LiDAR Point Clouds Data Filtering[J].Computer Engineering and Appl ications,2013,49(15):219-223.

[16] 刘志青,李鹏程,郭海涛等.融合强阈值三角网与总体最小二乘曲面拟合滤波[J].红外与激光工程,2016,45(04):128-135.

LIU Zhi-qing,LI Peng-cheng,GUO Hai-tao,et al.Integrating Strict Threshold Triangular Irregular Networks and Curved Fitting based on Total Least Squares for Filtering Method[J].Infrared and Laser Engin eering,2016,45(04):128-135.

[17] 黄先锋,李卉,王潇等.机载LiDAR数据滤波方法评述[J].测绘学报,2009,38(05):466-469.

HUANG Xian-feng,LI Hui,WANG Xiao,et al.Filter Algorithms of Airborne LiDAR Data:Review and Prospects[J].Acta Geodaetica et Cartographica Sinica,2009,38(05):466-469.