基于高程统计的机载LiDAR点云三角网渐进滤波方法

陈琳，范湘涛，杜小平

（1.中国科学院遥感与数字地球研究所 数字地球重点实验室，北京 100094；2.中国科学院大学，北京 100049）

1 引 言

激光扫描测距技术（Light Detection And Ranging，LiDAR）是一种快速直接获取地表模型的技术［1］。机载LiDAR点云数据包括地面点之外，还有树木、建筑物、桥梁、车辆等信息，获取真实的数字高程模型需要剔除这些非地面点，即点云数据滤波［2］。在机载LiDAR数据处理过程中，滤波是最关键的部分之一，是进行后续高层次数据处理的基础［3］。准确的滤波结果才能保证各种应用的有效性。

经过十几年的发展，滤波算法已经取得一定的研究成果。现有滤波方法根据技术路线差异，主要分为3类：形态学方法、基于内插的方法和基于曲面约束的方法［4］。基于形态学的点云滤波方法使用与被测地物相似的结构元，进行腐蚀、膨胀等形态学运算过滤地物点数据。众多学者开展了此方面的研究［5－8］，这类方法比较适用于城市区域，结果依赖移动窗口尺寸，对非平坦地区效果并不理想。基于内插的点云滤波方法首先建立粗糙的初始DEM，然后逐渐内插加密DEM，去除地物点。其中Axelsson提出的不规则三角网加密算法是比较有代表性的方法［9－11］，该方法在城市和森林地区都有较好的滤波效果，并且已经应用到商业软件TerraScan中。基于曲面约束的方法假设复杂地形在局部地区可以用简单的曲面表达，进而排除地物点。曲面的种类有样条曲面［12］、正交多项式曲面［13］和移动二次曲面［14－15］等。此类方法计算相对复杂，拟合过程会造成误差累积，较适合平缓地形。由于地形地物的复杂性，没有一种算法可以完全适用。

2 方 法

2.1 基于高程统计的阈值分割

不同地物要素会导致局部地区点云的高程分布差异，可通过高程统计得到高程阈值进行数据分割。大部分数理统计量与数据的空间位置和数据格式无关［17］，基于统计的处理方法普适性强。

高程直方图反映点云数据的高程分布。对比高程直方图，迭代计算点云数据的高程阈值，得到的最佳阈值将点云数据分割为两个高程差异明显的数据集。当研究区地形较为平缓、地物与地面差异明显时，单次分割即可将地物与地面点有效分离；反之，需要对已经分割得到的数据集进行多次分割，逐级分离地物与地面点。其具体处理流程为：

①绘制高程直方图。点云数据的高程直方图以高程为横坐标，频数为纵坐标，反映各高程与其出现频数之间的关系。

②参数阈值计算。计算原理根据文献［18］：查找待分类点云的最大高程值gmax和最小高程值gmin，取两者平均值为初始阈值g0；g0把点云分成大于g0和小于g0的两部分，g1为两部分期望值的平均值；如果g0与g1之差的绝对值大于自定义的极小值，则将g1赋值给g0，重复上一步的阈值计算，反之，则g1是分割的最佳阈值，退出计算。其中，极小值是一尽可能小的实数，以保证最终得到的最佳阈值趋于稳定（在本实验中自定义极小值取为0.01）。

③多阈值分割。如果待分类点云中有多个不同高度的目标，可对已分割结果进行多次分割，逐步将地物点剔除。

2.2 三角网渐进滤波

三角网渐进滤波算法默认局部地区的最低点为地面点，由这些最低点创建初始不规则三角网（Irregular Triangle Network，TIN）地形，通过计算点到初始地形的距离和角度判定点的类别，其具体流程为：

①参数设置。参数包括：研究区最大建筑物尺寸、点到三角面的最大高程阈值和点到三角面节点的最大角度阈值。

1)小范围内精细化导航，避免游客迷路。AR导航更直观、准确，能让步行的用户知晓自己周边位置场景的更多细节。

②选择种子点。以最大建筑物尺寸划分格网后，查找每个格网中的最小高程点，由这些点创建初始TIN。

③加密三角网。查找激光点，如果点到它落入的三角面片的距离和角度满足阈值参数，则判定该点为地面点，并加密到三角网中，反复迭代直到没有新的点被判定为地面点为止。

2.3 基于高程统计的三角网渐进滤波

基于高程统计的阈值分割方法简单易行，但由于缺少空间约束条件，滤波得到的地面点集会混杂较多的低矮地物点。而三角网渐进滤波方法缺少先验知识，在选择种子点时格网尺寸受到建筑物尺寸限制。本文方法引入高程统计分割结果作为先验知识，在此基础上划分格网时，格网尺寸不再受建筑物尺寸约束，可以定义小尺寸的格网以提取更多初始地面点，同时保证初始值的准确度，得到更加细致的初始地形。

由于先验知识的准确度会影响后期三角网滤波的精度，对点云数据进行高程统计分析时有三个注意点。第一，研究区范围较大或地形剧烈变化时，低地处具有较高高程的地物点与高地处的地面点可能有相同的高程值（“同高异类”），需要对研究区按照高程和平面位置进行分区统计。具体操作时，将原始数据按照横向和纵向分别二等分，得到4个次级尺寸的区域，如果小尺寸区域内仍存在“同高异类”现象，继续对小尺寸区域按照横向和纵向进行二等分，直到分区内不存在“同高异类”为止；同时，如果相邻的小区域内点云高程分布相近，地形变化相对平缓，也可将相邻小区域进行合并。第二，为了避免地物边缘点混入分割结果中造成初始点选择误差，在进行分割时，可以适当调低分割阈值，舍弃部分地面点，以保证统计分割结果的准确率。第三，在微小的地形变化中，单纯依靠高程统计分割方法很难剔除低矮植被点，但由于植被是略高于局部地表的，按照选择极低点的假设条件，这些点不会影响种子点的精度。

由于高程统计分割方法在建筑物密集的城区具有较高的准确性，而在地形起伏较大的林区和山区准确性相对较低，所以基于高程统计的三角网渐进滤波方法更适合于城区点云数据。

该方法滤波主要流程：

①剔除系统粗差点。依次查找每个点一定范围内（以5m～8m为半径的圆周内）的其他激光点，如果该点比任意其他激光点都低1m以上，则将该点判定为粗差点，予以剔除。此处的1m为经验值，可以根据地形起伏情况，手动调整。

②绘制点云数据的高程直方图，统计点云高程得到阈值参数，对照直方图进行高程分割。当研究区范围过大或地形变化较剧烈时，按照地形特征将点云数据分成范围较小的多个区域，对每个子区域的点云数据执行高程统计和阈值分割后，合并地面点分割结果。

③对分割得到的备选地面点划分格网，提取每个格网内的最低点作为种子点，创建初始三角网地形。

④计算粗差处理得到的所有激光点到初始地形的距离和角度，满足阈值参数的点被加密到三角网中进行下一次迭代过程，当没有新的地面点加入到三角网时，迭代计算结束，得到滤波结果。滤波流程如图1所示。

图1 基于高程统计的三角网渐进滤波流程图

3 数据与实验

为了检验此算法，使用ISPRS公布的滤波实验数据，其中Site2数据覆盖范围广，用以解释分区原理，其他3组数据有半自动半手工分类结果作为真值，与原始三角网渐进滤波方法相比进行精度验证。

3.1 实验数据

所有数据都为城区数据，点密度为0.67pts／m2，点间距为1.0m～1.5m。Samp21中高程为288.48m～320.28m，研究区内有桥梁、公路、小隧道和植被等。samp23中高程为262.27m～348.29m，有较明显的地形起伏，研究区内有大型建筑物、不规则形状的建筑物以及汽车等人工地物。samp31中高程为226.94m～343.95m，数据中存在部分极低点，研究区内分布有高低起伏的建筑物、零散汽车以及分布于建筑物间的各种植被。Site2为大范围数据集，覆盖区域包含了samp21和samp23，高程范围为248.1m～482.63m，地形变化较大，地物复杂。

3.2 实验结果

Site2数据集包含了408921个激光点，覆盖范围广，研究区内地形变化较大，在进行高程统计前，需对该数据集进行分区以保证子区域内没有“同高异类”现象，然后对分区内数据进行独立的统计分割，再将分割结果合并，提取地面点，建立初始地形。图2分别为原始数据集、数据分区结构以及基于本文算法的滤波结果。图2中的（a）图是原始数据集按照高程大小进行渲染的显示结果。从图中可以看出，数据集的左上部分地形较平坦，建筑物与植被高差较明显，而数据集的右下部分地形变化剧烈，植被以及建筑物分布较复杂。因此，在图2（b）的数据分区结构中，右下侧被划分得更为详细，以确保高程统计方法可以准确地识别地物与地形间的高程差异。图2（c）为基于高程分割数据集进行三角网渐进滤波后得到的最终结果。该方法有效地剔除了研究区内的建筑物和植被，同时，在高程渲染图的左上角和左下角显示不明显的桥梁也被有效地剔除了。在剔除地物点的同时，此方法有效保留了研究区最右侧的高地地形。

图2 大数据集分区

图3是3个数据集的ISPRS公布的半自动半手工分类结果、经原始三角网方法和本文方法滤波得到的对比结果图。

图3中原始算法与本文算法的滤波效果都比较理想，在3组数据中桥梁、大部分建筑物和零碎植被点都被较好地剔除了。不同之处在于，本文算法在3组数据的边界建筑物滤波效果优于原始算法。

在samp21中，对比右下角的建筑物，原始三角网算法错误地保留了部分建筑物边界，并且滤波结果中建筑物的边缘也相对粗糙，本文方法较好地剔除了这部分地物点；同时，对比左下角处的地物部分，原始算法相比较本文算法，过度地剔除了一部分地面点。在samp23中，本文方法较原始三角网算法，更好地剔除了右边界处的建筑物，建筑物边沿与真值数据也很匹配，但是左边界下半部分的少量地物点被错误地保留了。在samp31中，本文方法正确地剔除了下边界左部的植被、下边界中部的建筑物以及中部的植被建筑物交错地物点。

图3 3组数据滤波结果对比图

3.3 误差分析

根据ISPRS公布的误差评定方法，对3组数据集进行误差统计。表1和表2分别是原始三角网渐进滤波方法与本文方法得到的3个数据集的统计结果。其中，以ISPRS公布的半自动半手动分类结果作为真值判断依据，a为被正确分类的地面点，b为被错分类为地物点的地面点，c为被错误分类为地面点的地物点，d为被正确分类的地物点；一类误差为＊100%，二类误差为＊100%。

在原始算法精度较高的情况下，3组数据经过本文方法滤波后，提取了更多地面点的同时，有效剔除了地物点的影响，精度均得到一定的提高。在samp21和samp31数据集中，研究区地形起伏较小，进行高程统计分割后，先验知识的准确度也相对较高，两个数据集的滤波精度提升幅度，相比较地形起伏更加剧烈的samp23数据集更大，所以研究区的地形条件是影响误差的一个重要因素。所有数据集在先验知识的基础上提取了更多的初始地面点，待分类激光点相对减少，同时建立了更加细致的初始地形，所以在三角网渐进滤波时，待分类激光点到地形的参数值更加准确，滤波精度得到提高。

表1 原始三角网渐进算法精度评价

表2 本文算法滤波结果

原始三角网渐进滤波算法中，格网尺寸的划分会受到研究区内最大建筑物尺寸影响，以免划分的格网中出现没有地面点的情况。在本文方法中，经过高程统计分割后，在很大程度上解决了格网尺寸的约束，格网可以按照用户自定义的尺寸进行划分。为了测试格网尺寸对滤波结果的影响，本文从5m～40m划分了8个不同的格网尺寸，然后统计了3组数据在不同尺寸下得到的滤波结果的两类误差，图4为3组数据的误差变化情况。研究发现，本方法得到的滤波结果，二类误差相对平稳，一类误差随着格网尺寸变化相对剧烈，并且随着格网尺寸变大，一类误差总体呈上升趋势。当格网尺寸设定在10m～20m之间时，一类误差和二类误差较小且相对平衡，可以得到比较理想的滤波结果。

图4 本文算法中两类误差随格网尺寸变化

4 结束语

本文在三角网渐进滤波方法及其改进方案［19－21］基础上，提出基于高程统计的三角网渐进滤波方法，该方法先从点云数据中提取先验知识，基于先验知识再进行三角网滤波，经过3组数据验证，均得到较好的滤波效果，比较适合于城区点云数据处理。在进行高程统计分割时，当研究区面积较大或者地形起伏较剧烈时，需要对研究区进行分区统计，以保证统计分割的精度。由于有高程统计分割得到的先验知识，在进行格网尺寸划分时，该方法可以不受原始三角网渐进滤波算法中最大建筑物尺寸的约束，自由选择尺寸，经过实验验证，此方法以10m～20m的格网尺寸得到的滤波结果更优。与原始三角网渐进滤波方法相比，此方法提高了滤波精度，但对数据进行分区时，本文方法未引入有效的评估参数，自动化水平有待进一步提高。

［1］刘春，陈华云，吴杭斌.激光三维遥感的数据处理与特征提取［M］.北京：科学出版社，2010.

［2］周晓明，马秋禾，许晓亮，等.LIDAR点云滤波算法分析——以ISPRS测试实验为参考［J］.测绘工程，2011，20（5）：36－39.

［3］刘经南，许晓东，张晓红，等.机载激光扫描测高数据分层迭代选权滤波方法及其质量评价［J］.武汉大学学报：信息科学版，2008，30（6）：551－555.

［4］黄先锋，李卉，王潇，等.机载LiDAR数据滤波方法评述［J］.测绘学报，2009，38（5）：466－469.

［5］KILIAN J，HAALA N，ENGLICH M.Capture and evaluation of airborne laser scanner data［J］.ISPRS，1996，31（B3）：383－388.

［6］CHANG K J，TABOADA A，CHAN Y C.Geological and morphological study of the Jiufengershan landslide triggered by the Chi－Chi Taiwan earthquake［J］.Geomorphology，2005，71（3－4）：293－309.

［7］梁欣廉，张继贤，李海涛.一种应用于城市区域的自适应形态学滤波方法［J］.遥感学报，2007，11（2）：276－281.

［8］隋立春，张熠斌，柳艳，等.基于改进的数学形态学算法的LiDAR点云数据滤波［J］.测绘学报，2010，39（4）：390－396.

［9］AXELSSON P.Processing of laser scanner data－algorithms and applications［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1999，54（2－3）：138－147.

［10］AXELSSON P.DEM generation from laser scanner data using adaptive TIN models［J］.International Archives of Photogrammetry and Remote Sensing，2000，33（Part4－B4／1）：110－117.

［11］AXELSSON P.Ground estimation of laser data using adaptive TIN－models［C］.Proceedings of OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models，Stockholm，Sweden，March，2001：185－208.

［12］ELMQVIST M.Terrain modeling ＆analysis using laser scanner data［J］.ISPRS，2001，34（3／W4）：219－224.

［13］BROVELLI M A，CANNATA M.Digital terrain model reconstruction in urban areas from airborne laser scanning data：the method and the example of the town of Pavia（Northern Italy）［J］.ISPRS，2002，34（2）：43－48.

［14］苏伟，孙中平，赵冬玲，等.多级移动曲面拟合LIDAR数据滤波算法［J］.遥感学报，2009，13（5）：833－839.

［15］尚大帅，马东洋，赵羲，等.一种基于移动曲面拟合的机载LIDAR点云数据滤波方法［J］.测绘技术装备，2012，14（2）：23－26.

［16］SITHOLE G，VOSSELMAN G.Experimental comparison of filter algorithms for bare－Earth extraction from airborne laser scanning point clouds［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2004，59（1－2）：85－101.

［17］杨晓云，梁鑫.基于参数统计的LiDAR数据分割算法［J］.广西民族师范学院学报，2012，29（3）：33－35.

［18］庄军，李弼程，陈刚.一种有效的文本图像二值化方法［J］.微计算机信息，2005，21（06X）：56－57.

［19］李卉，李德仁，黄先锋，等.一种渐进加密三角网LIDAR点云滤波的改进算法［J］.测绘科学，2009，34（3）：39－41.

［20］隋立春，张熠斌，张硕，等.基于渐进三角网的机载LiDAR点云数据滤波［J］.武汉大学学报：信息科学版，2011，36（10）：1159－1163.

［21］左志权，张祖勋，张剑清.知识引导下的城区LiDAR点云高精度三角网渐进滤波方法［J］.测绘学报，2012，41（2）：246－251.