点云滤波是机载激光雷达 (Light Detection and Ranging, LiDAR)点云数据处理的重要步骤。主动式成像的LiDAR点云数据具有数据量大、数据分布离散等特点，因此数据的处理方法与其他遥感数据类型有所不同。三维点云数据滤波是后续点云分类以及植被、道路、建筑物等特定物提取等处理步骤的基础。如何实现快速、高效的点云滤波，提升滤波算法在多种不同地形的适应性，具有重要意义。

点云滤波的算法自20世纪90年代以来，一直是热门的研究方向。典型的算法包括基于数学形态学的算法、基于坡度的算法、基于三角网的算法、基于曲面拟合的算法和基于聚类分割的算法等。根据ISPRS的结论，多数传统算法针对特定的地形，能够达到令人满意的效果，但是在复杂和混合地形中表现不佳^[1]。

数学形态学算法最早被用于剖面分析的点云滤波，后续研究者对形态学方法进行改进以提升形态学滤波方法的准确率。Sun等^[2]提出一种基于迭代多尺度形态学开重建的方法，对含有各种复杂地形地物特征的城市场景均有较强的适应性。此类方法主要问题是在地形起伏较大区域难以兼顾一些较小且复杂的地形，渐进窗口难以适应局部地形坡度的阈值。

基于坡度的算法，确定两个相互临近的点，依据两点的坡度是否大于一定阈值，对于较高点进行判定。算法需要频繁寻找临近空间点，在面临大量的定位、查找和修改工作时，运算量较大，导致点云数量较大时算法效率偏低。

基于曲面的方法在一定范围内采用二次或者其他曲面形式来逼近当前窗口内的地形表面，依照点到曲面的高差来判定是否为地面点。Zhang等^[3]提出的布料模拟滤波实质上也是使用质点弹簧模型，去找到最优的地形表面。

基于三角网的方法通过构建三角网，对点云进行定位和判断。在ISPRS的第III工作组的8种滤波算法对比^[1]当中，三角网渐进加密 (Progressive TIN Densification, PTD)算法获得最小的最高总误差以及最高的Kappa系数。后续的研究改进方向主要在：优化构网和查找速度^[4]、优化种子点选取策略 (避免地物点和低点噪声被选入种子点)、优化在山脊和地形断裂区域的表现^[5]以及减少外部参数等方面。

基于分割的方法，将点云分割成块之后，基于块对象的关系来进行滤波。Wang等^[6]采用连通点集的方式基于局部表面凹凸性进行点云的分割，对散乱点云具有良好的效果。此类方法对分割精度要求较高。针对不同的地形采取的分割策略和参数设置等都对滤波结果优劣产生影响。

基于分类的方法，训练分类器将原始点云分为地面点和非地面点两类。传统的机器学习的方法使用信息向量机^[7] (Information Vector Machine, IVM)、相关向量机^[8](Relevance Vector Machine, RVM)、随机森林 (Random Forest, RF)等采集样本特定的特征训练分类器对点云进行分类。Rizaldy等^[9]使用全卷积网络在ISPRS滤波数据集和荷兰ANH-3数据集上均得到较好的滤波效果。Janssens等^[10]评估了PointNet在点云滤波的应用效果，基于分类效果提出了改进方向。Zhao等^[11]将迁移学习引入机载激光雷达点云分类，有效提升分类精度，较大程度减少了训练时间。

基于扫描线对点云进行分割，避免了针对点直接进行对象分割的巨大的计算量，同时能够在一定程度内保留点云之间的拓扑信息。Ye等^[12]在进行多方向扫描线分割的基础上定义了地面显著性的概念，提出基于地面显著性进行近似全局优化的半全局滤波 (Semi-Global Filtering, SGF)方法。

混合的方法是指结合上述多种经典或改进算法以提升算法在不同区域的综合表现，补偿某一类算法在特定区域的错误情况。Liu等^[13]结合强阈值三角网渐进加密滤波和总体最小二乘曲面拟合滤波方法，提升了算法的可靠性和自动化程度。

文中在Ye等^[12]的地面显著性计算的基础上对异形建筑物区域进行了改良以提升显著性计算准确性，并通过对呈现不同地面显著性的网格采取不同的分类面并进行滤波，整合形成最终滤波结果，能够较好地处理大型建筑物和各种复杂地物，同时对有起伏区域使用拟合方法计算分类面，更好表达真实的局部地形，对多种地形环境的点云能够达到较为稳定的滤波效果。

文中在以虚拟格网组织激光雷达三维点云的基础上，使用多方向扫描分段的方式进行分割并计算优化后的地面显著性值，以此初步判断地形类别，针对不同的地形采取包括曲面拟合、显著性插值等手段获取更准确的滤波分类面，以完成点云滤波过程。

虚拟格网是使用相同规模的矩形网格单元将点云在二维上进行划分，如图1所示。多数情况下网格单元在XOY平面上的长和宽可设置为相等。对于激光脚点，根据其在二维平面的x、y坐标，将其划归到所属网格当中，根据点坐标可以计算出所属网格行列号，通过网格ID也能索引到所包含的所有激光脚点序号。点云数据以网格方式组织后，为了节省计算量，一般采用单个网格内的最低点高程代表网格单元高程。格网化公式如下：

式中：i和j为格网的行列号；x和y为点云的水平坐标；$\min x$和$\min y$为点云的水平边界值；d为网格间距。在一些点云处理算法中，会保存单个网格内部的所有三维点序号以便于索引。

图  1  虚拟格网

Figure 1.  Virtual grid

地面显著性γ求解的目标是为了将高程区分度明显的地物网格和地面网格以不同的地面显著性值分离，此算法使用基于多方向扫描的分割算法来计算网格单元的地面显著性。

显著性求解的第一步是对诸多格网的分段过程，即阶跃探测 (如图2所示)，旨在分离在高程上有显著差异的网格单元，将相邻且高程接近的网格划分至同一个分段。

图  2  多方向高程阶跃探测

Figure 2.  Elevation mutations detection in multiple directions

扫描方向以点云XOY平面坐标为基准，自y轴正方向起算，顺时针旋转，每45°为一个扫描方向，共8个扫描方向，各方向扫描过程相互独立。在8个方向内对于每个方向的扫描线依照高程进行分段，此处所使用的高程为单个网格内所有点的高程最低值，后续的显著性计算中所使用的高程与此相同。使用一定的高差阈值ΔH，在单个扫描方向的某一条扫描线上，从起算的网格单元开始沿扫描方向逐个向前遍历，每当下一处格网的网格高程与当前网格高程的高差超过高差阈值ΔH，或下一处格网为空时，则将起算网格到当前网格的所有网格划分为一个分段，然后从下一个非空网格重新开始起算。分段的主要目的是区分典型的建筑物，因此高差阈值ΔH一般设置为低于测区内最低矮建筑物的高度即可，多数情况下使用1.0~1.5 m的高差阈值即可有效进行网格高程的区分。

地面显著性求解是一个多方向累计递减的阶跃探测过程。基于各方向的高程变化，对分割段的显著性进行运算。

单个方向的显著性判断是基于分段的，即在某一方向的某次判断中，满足高程条件的分段内部所有网格单元同时修改显著性；多方向叠加过程是则基于网格单元的，对网格单元在8个方向的显著性计算叠加得到最终的显著性值，显著性叠加过程与其所处分段无关，与同一分段的其他网格单元也无关，网格单元之间的运算过程独立。

求解显著性之前，需要将所有非空的网格的地面显著性的值初始化为1。在上一步网格分段的基础上，对每个方向再次以相同方式进行扫描。在某个方向的单条扫描线的扫描中，依据网格分段过程中得到的分段信息，定位位于不同分段邻接处的网格单元。在邻接单元处，若当前分段末尾网格相对于下一分割段的起算网格的高程差值大于一定阈值，则将当前分段的所有网格的地面显著性值减去1/8。这一过程在8个方向上进行叠加，当某一片网格高程在8个方向均显著高于区域的网格高程，则地面显著性将最终计算为0，若某一格网单元在8个方向的邻接分段的高程值均未显著低于当前分段，则该区域的最终地面显著性计算值为1.0。值得注意的是，分段是一维的概念，指某扫描方向上某一扫描线上的若干连续网格单元的集合，同一单元各方向的分段并不一致，有别于对象分割中的分割块。

扫描线方向上分割段发生高程阶跃的情况有图3所示3种：第一种是由高到低，此时Segment 2的所有网格单元的显著性减去1/8；第二种是由低到高，因此Segment 4内的网格单元的显著性不做处理；第三种同样是由高到低，但由于Segment 5只有一个网格单元，于是将Segment 6的末尾单元与Segment 7起算单元做高差计算，若高差较小，则对Segment 4和Segment 6的网格显著性不做处理。若Segment 6 的高度显著高于Segment 7，则将Segment 4和Segment 6中的网格单元的显著性值均减去1/8。如图4所示为点云显著性渲染结果，高程分层设色图从蓝色到红色代表由低到高，显著性灰度图纯黑色代表显著性为0，白色代表显著性为1，颜色深代表地面显著性低。

图  3  地面显著性计算

Figure 3.  Ground saliency calculating

图  4  点云高程分层图（a）及地面显著性灰度渲染图（b）

Figure 4.  Layered map of point cloud elevation(a) and gray rendering map of ground saliency(b)

正常情况下高程显著高于周围地区的建筑物应当具有较低 (接近于0)的显著性值，然而实地的建筑物在形态和分布上较为复杂，地面显著性的计算值往往远高于0。如中心凹陷状建筑，当屋顶的高差阈值大于显著性分割时的高差阈值时，在扫描时凹陷处有较多的方向处于较低水平，中心处显著性会远高于0甚至达到1。凹陷处过高的地面显著性数值会对网格类别划分的过程产生误导，需要进行消除。增大基本网格尺寸可在一定程度上减少此类现象的发生，但大格网会导致细小地形起伏被削平 (类似渐进形态学滤波效果)，不宜作为主要优化手段。为保证建筑物分离，此处使用一种简单的邻域窗口搜索策略，如图5所示。

图  5  特殊形状建筑物屋顶

Figure 5.  Roofs of special shape building

对于一处地面显著性为1或接近1的网格，若在其一定大小的邻域窗口内，地面显著性显著小于1的非空网格占绝大多数，则将其地面显著性设置为窗口内低显著性网格的显著性的均值。建筑物中心的凹陷处导致的异常显著性值，使用此种方法基本能够有效去除。

这种建筑物凹陷在绝大多数情况下不会呈现广泛的分布，因此在优化的时候考虑使用邻域信息。此处采用的优化方式类似数字图像处理领域的数学形态学腐蚀方法。使用一定大小的窗口对网格显著性值进行类似形态学的腐蚀操作。如图6所示，将地面显著性高于阈值γ_thre的网格单元看作二值图像的前景像素，将小于阈值的网格看作单元背景像素。使用形态学腐蚀操作之后，孤立的高显著性单元将会被视作噪声滤除。这种过滤操作针对网格单元，即将孤立的高显著性单元的显著性值置为周围低显著性邻域网格单元的显著性均值。

图  6  中心凹陷建筑物区域显著性分布

Figure 6.  Ground saliency distribution of center sunken buildings area

图6中的红色虚线框内区域即屋顶中部有凹陷的建筑物的网格显著性。通过使用不同窗口大小的窗口进行类似形态学腐蚀的操作，将孤立的高显著性区域的显著性设置为邻域的均值，从而使建筑物区域的显著性保持一致性。

一些类似四合院的建筑，中心处实际上为地面，由于上述显著性优化过程导致其地面显著性最终计算值较低，但对后续的滤波结果影响不大，在后续的分类面插值求解过程中弥补这一过程。

在点云求取到地面显著性后，根据不同的地面显著性数值，来表达当前网格内最低点为地面点的可能性：数值接近1，网格最低点为地面点的可能性较高；数值为0，则为非地面点的可能性较高。对于地面显著性较高的网格点，其本身的高程可作为滤波分类面的参考值，而对于地面显著性较低的点，由于最低点为非地面点的可能性较大，因此要找到一个尽量贴近真实地面的分类面，其分类面的高程值需要根据临近的地面点寻找参考。

在对不同的网格点处以不同的分割手段前，需要先确定一个合理的阈值将二者分开。在正常情况下，由于二维网格化过程使用网格最低点作为网格高程，因此建筑物屋顶处的网格高程应接近屋顶面高程，少数情况下屋顶下有噪声点干扰，应在预处理过程中剔除，而房屋附近的地面点网格高程则接近地面高程，因此建筑物屋顶分段在大多数方向可以找到高程显著低于其的邻近地面进行分段，故建筑物屋顶分段内的网格地面显著性计算值一般接近或等于0；斜坡区域则有少数方向临近分段高程显著低于自身高程，故地面显著性在0.5左右；平坦地面点在无低噪声点和地面凹陷的情况下各个方向的临近分割段极少有低于自身高程的情况，地面显著性则接近或等于1。

包含地面点的网格点，自身的高程是确定当前网格单元分类面的主要参考数据；而建筑物网格点，需要寻找附近的地面点作为参考。进行显著性数值的分割时，首先应确保将建筑物点分离开来，分离的核心在于地面显著性划分阈值的确定。具体的显著性阈值设置方法，详见实验部分。

在使用合适的分割阈值将格网分开后，不同的网格点的后续处理方式如下所述。

对于平坦地面点，使用网格最低高程作为分类面高程。事实上，此类区域同样可以基于网格邻域点使用曲面拟合方案得到当前区域的近似分类面，但是对于较平坦的地面区域，拟合得到分类面与直接使用网格最低高程效果相近，对于结果影响极小，因此亦可直接使用网格最低高程作为水平分类面高程以减少不必要的计算量。或者判定网格最高点与最低点高程差较小时使用最低点高程作为水平分类面高程。

对于有起伏的区域或是斜坡区域的分类面获取，比较合适的是采用拟合的方式。此处的分类面拟合方法，基本思路类似于移动窗口最小二乘曲面拟合滤波方法。从${\rm{3}} \times {\rm{3}}$窗口的网格单元的最低点中各选一个点，共计选取小于等于9个点，纳入当前网格邻域点集合中，根据点数不同，采取不同的处理策略。在确定当前网格非空的前提下，具体分以下几种情况：点数小于3时，无法使用拟合方案，则仍使用建筑物网格区域的内插方法获取分类面高程，具体过程参考下述建筑物网格点插值方法。当点数大于3且小于6时，此时点数足够拟合一个平面，但不足以拟合一个二次曲面，若低噪声点已被去除，此时取其中最低的3个点作为拟合平面的参考点。点数多于3个点，可以使用剩下的点作为检查点。当点数大于等于6时，即拟合二次曲面作为分类面。从中选取高程最低的6个点，作为参考点拟合二次多项式曲面。若点数多于6个点，可以从剩下的点中选取检查点，以保证拟合的分类面的准确性。

对于建筑物区域的网格点，分类面则根据临近地面点信息来生成。在这一过程中，地面显著性和距离因素均可以作为辅助判断依据，文中采用了基于地面显著性的网格点插值手段。

对于地面显著性分割后确定为建筑物的疑似点，在获取分类面时采用如下方法：

(1) 从当前网格向周围8个方向延伸搜索，直至每一个方向均找到地面显著性大于阈值且高程低于当前网格高程的网格单元。若某一方向搜索距离大于距离阈值仍未找到符合条件的网格，则当前的方向判定为失效方向，不作为分类面插值的参考。对于符合条件的网格，记录该方向找到的网格高程和距离当前格网的距离 (格子数)。

(2) 以当前网格单元为中心，在8个方向搜索到的显著性值较高的网格中，找到距离最近的网格。若同一距离下有多个点，则使用显著性值最大的网格高程值作为分类面高程。

(3) 对网格水平分类面的高程，设定一定的高差阈值，对当前网格内所有点进行类别划分。高程显著高于分类面一定高度的，判定为非地面点，在此之下的判定为地面点。

文中研究以Visual Studio 2015为开发平台，系统环境为Windows 7 64位。实验数据为广州某地区机载激光雷达点云数据，共9个子测区。单个点云数据的点数为500万左右，点云密度为每平方米10~20点。数据源位于城市区域，地物种类复杂多样，测区内包含高层建筑、密集房屋、植被、桥梁、道路等多种地物。网格数据组织过程中，需保证网格单元尺寸大于平均点距，实验中网格间距为1.0 m。各子测区地形分布见表1。

用于验证参数设置和结果可用性的标准是国际摄影测量与遥感协会(ISPRS)的评价体系(交叉表法)：两类误差和总体误差指标（如表2所示）。一类错误率指拒真错误，即地面点被误分为非地面点的比例，二类错误指纳伪错误，即非地面点被误分为地面点的比例。总体错误率即所有的类别错分点占总点数比例。错误率越低代表整体滤波效果越好。同时使用了Kappa系数作为定量评价的标准之一，Kappa系数越大，表明滤波效果越好。

在进行地面显著性求解之后对网格点进行类别划分过程中，针对阈值的设置，设置了从0.375~0.875等一系列的数值，比较分割的结果和最终滤波的结果，以探索分割阈值的最优值。当γ_thre2为1.0时代表无点使用局部最低高程，此时的划分代表二类划分，在可靠地面点区域全程使用拟合方式求取分界面，计算量较大，故未列出γ_thre2为1.0的方案。

使用不同的显著性划分阈值，对多组数据分别进行滤波得到结果后，在不同的显著性阈值下，统计一类错误率、二类错误率和总体错误率。错误率的统计结果为在多组数据测试下的平均值。具体结果如表3所示，E1、E2、Et、K分别表示一类错误率、二类错误率、总体错误率和Kappa值。

表3中，Site1和Site2分别代表平坦区域和地形起伏区域的统计值。根据统计结果，错误率最小时所取γ阈值在平坦区域和起伏地形区域表现不尽相同。平坦和建筑物较密集区域，使用较大γ_thre1阈值时二类错误率和总体错误率更低；在地形有起伏区域，显然γ_thre1采用0.5或更低的数值时具有更好的滤波效果。综合各种测区，γ_thre1在采取0.5附近的数值时，对多数情况下分割的结果比较合适，对不同的地形都具有不错的适应性，可以作为算法的默认阈值。当整个测区的地形普遍偏向于平坦或者起伏时，仍可以通过调整γ_thre1获得更好的滤波表现。

γ的设置对最终的错误率具有明显影响，在其他参数采用默认的相同值的情况下，使用不同γ阈值下的错误率差距可能超过1%。对于斜坡区域和平地区域的γ阈值分界γ_thre2，在斜坡区域和非地面区域的γ_thre1相同的情况下，不论是在城区还是山区，γ_thre2的设置对于最终错误率的影响并不明显。在山区γ_thre2的数值对滤波精度的影响相较于城市更大，但总体而言差距仍不算显著，故大多数情况下γ_thre2无需进行调整，可以作为固定的参数存在，这里主要讨论γ_thre1阈值的设置对结果的影响。在平坦区域，增大γ_thre1，错误率明显降低。在此类区域，建筑物分布密集且广泛，此时较为稳健的做法应当是使用插值的手段而非拟合，因此增大γ_thre1的值，更多的过渡区域将会使用显著性插值手段寻找分类面。在起伏地形区域，减小γ_thre1，一类和总体错误率则更低。此时使用插值的手段从附近的地面点插值的高程可能会与当前区域的地形不一致，曲面拟合获取地形分类面方式则在这些区域的表现更加稳健，有利于表达起伏地形的真实情况。从统计结果来看，错误率表现与算法设计的思路基本一致。

实验使用改进渐进形态学算法 (Progressive Morphological Filter, PM)、比较成熟的商业化软件TerraSolid 中三角网渐进加密算法(PTD)和布料模拟滤波算法 (Cloth Simulation Filter, CSF)作为对照。错误率和Kappa值如表4所示。

多数情况下，PTD算法和文中表现出了更低的一类错误率，其中PTD表现最为优秀，即地面点保留较为完整。文中算法一类错误率略低于PTD算法，但差距不大。CSF算法在大部分情况下，能够保持较低一类的错误率，部分情况下错误率会剧增，尤其是在地形较复杂的情况下，很难找到较为合理的参数适应整个测区的情况。PM算法在保持地面点上表现不如其他3种算法。

在二类错误率的表现上，各算法的表现因测区而异，总体而言各算法差距不大，具体情况仍与测区地形密切相关。在地势平坦区域，文中算法和CSF算法表现出最好的二类错误率，其中，当地形较为复杂，例如有高架、桥梁和连续建筑等地物时，PTD算法表现更加稳定，文中算法与PTD算法精度相当。总体而言，相较于其余几种算法，文中算法在二类错误率的表现上最为优秀。

在总体错误率上，文中算法在9个测区中有5个测区总体错误率优于PTD算法。CSF算法在多数情况下错误率低于PTD算法，但差距不大。在地形较为复杂或特定地形下，CSF算法错误率显著上升。PM算法在部分平坦区域表现也比较优秀，但实际应用时参数多，多数需调整以达到最佳效果，算法自适应性差。Kappa系数的统计表现与总体错误率接近，4种算法在绝大多数区域的Kappa系数均在80%以上；文中方法和PTD算法仍在不同地形区域获得了最好的Kappa数值表现。

算法效率上，实验使用的数据测区点数在400~700万点之间，单个测区点数500万左右。4种算法中，PM算法效率最高，滤波时间在1 s内(不包含文件I/O花费时间)。文中算法在效率上低于PM算法，滤波时间在2~6 s之间。而PTD算法和CSF算法花费时间较长，尤其是当迭代条件较为苛刻时，滤波时间显著增加，多数测区点云滤波时间超过10 s。

图7为滤波结果和带有真实类别的点云对比，测区1地形简单，地势平坦，地物和地面界限分明，文中算法在此测区表现出良好的滤波效果，植被、建筑等地物均能够较好地和地面分离。测区9地形复杂，地物和地面界限不分明，且有高架、斜坡、密集植被(植被下地面点稀少)和低矮人工地物。除了极为低矮的人工地物和植被(图7(d)红色区域附近)由于高程参数无法适应，导致被算法归类为地面点之外，其余地物均有效分离。同时，位于高架桥下的少数低矮植被由于无法找到足够的有效参考地面点，未能完全和地面分离。

图  7  测区1、9滤波结果对比。(a) 测区1滤波结果；(b) 原始点云；(c) 测区9滤波结果；(d) 原始点云

Figure 7.  Comparison of filtering result of site 1 and 9. (a) Filtering result of site 1; (b) original point clouds; (c) filtering result of site 9; (d) original point clouds

综合各个测区来看，算法在不同地形下均能够达到较好的滤波效果，较少出现集中和成规模的错误滤波情况，大多数滤波错误可以通过调整阈值参数宽容度来降低或消除。同时，算法在效率上在同类算法中也处于较高水平，能够适应大量点云数据的生产应用需求。

对于多数滤波算法在单一地形表现较好，对混合和复杂地形表现不够，需要针对某些特定地形做改进和补偿的现状下，文中采用扫描线分割的方式替代传统方法“渐进加密”和“全局优化”的思路，大大减少了计算量。基于地面显著性对不同地形进行粗分类并对显著性计算结果进行优化，针对不同地形有选择地使用曲面拟合或基于显著性插值的方法计算分类面并滤波，有效避免了单一算法在不同地形频繁进行针对性修改和参数调整的过程。在多种地形点云的滤波测试中，均得到稳定可靠的滤波结果，且需要调整的外参数极少；相对于经典算法，在滤波准确率和自适应程度上都具有明显的优势。下一步改善的目标是针对地形粗分类机制的改进，使用自适应程度更高乃至完全自适应阈值的方式区分地形，并对目前算法未能解决的地形断裂区域进行优化。