1、引言

机载激光雷达（LightDetectionandRanging,LiDAR）数据滤波是目前点云数据处理领域中的关键技术，其精度直接影响数字高程模型（DigitalElevationModel,DEM）等产品的质量。因此，滤波算法研究具有重要意义。

机载LiDAR系统可以快速获取地表高精度、高密度的3D点云数据及其强度信息。但是，目前的滤波算法往往仅利用地面的几何信息（如梯度[1,2]、连续[3,4,5,6]等），这种数据利用方式未考虑点云强度信息在滤波中的辅助作用，未能充分发掘点云数据提供的所有信息，势必影响滤波的准确性。另外，已有滤波算法设计针对地形起伏不大、地形连续、地物较为单一的简单场景，由此导致算法无法适用于地形起伏较大、地形不连续、地面附着物形态各异、规模不一、数量众多的复杂场景[7]。

如果能综合利用点云数据的几何及强度信息并针对复杂地形特点设计算法，必将提高滤波结果的准确性并促进机载LiDAR在复杂场景下的滤波应用。综合利用所有信息的关键首先在于对机载LiDAR数据的表达方式。但传统的数据结构，如栅格格网[8,9]、不规则三角网（TriangulatedIrregularNetwork,TIN)[9]、点云[11]和八叉树[12]应用于机载LiDAR数据表达时均存在局限性：栅格格网和TIN为2.5D数据结构，用其表达3D点云将导致信息损失并进一步影响滤波结果的完整性；点云并未明晰表达激光点间的空间结构及拓扑信息，由此导致滤波算法设计困难、效率低下；八叉树各节点的尺寸不一，各节点间的邻接关系难以建立，这同样增加了算法设计的难度。为了克服上述限制，采用一种更为简单的真3D数据结构表达机载LiDAR数据显得尤为重要。同时，针对已有的滤波算法存在的无法适应复杂场景的问题，有必要开发一个适应范围更广泛的算法以更好地处理复杂多变场景。

因此，本文提出了一种基于灰度体素结构分割模型下的机载LiDAR3D滤波算法。该算法首先将机载LiDAR数据规则化为灰度体素结构，其中灰度为各体素内激光点的平均强度值的离散化表示；然后基于体素间的空间连通性和灰度相似性准则，将灰度体素结构分割并标记为若干个3D连通区域；最后依据地面的高程特性检测地面对应的3D连通区域。该算法的优势在于：(1)灰度体素结构的真3D特性使其具备完备表达3DLiDAR点云的能力且该结构同时融合了机载LiDAR数据的高程（各体素间的邻接关系为高程接近的隐含表达）和强度（各体素依据其内的激光点强度赋值）信息，有助于综合利用LiDAR数据的所有信息从而提高滤波精度；(2)3D滤波算法。利用地面目标在灰度体素结构中的“空间上具有连通性、灰度上表现出一致性”特性构建滤波模型，基于灰度体素结构分割模型、针对复杂场景中地面的几何及辐射特征设计，明显优于仅利用地面几何信息的滤波算法[1,2,3,4,5,6]的设计，可用于更复杂场景的滤波应用。(3)滤波结果为3D地面体素形式，可直接用作3D地面模型。

2、基于体素分割的3D滤波

算法包括机载LiDAR数据规则化、灰度体素结构分割及地面目标检测3个步骤。

2.1机载LiDAR数据规则化

规则化就是将包含场景目标的空间范围依据体素分辨率划分成3D格网并依据格网单元（即体素）内激光点的强度属性对各个体素赋值为不同灰度等级的过程。其中，空间范围由包围机载LiDAR数据的轴向平行包围盒确定；体素分辨率依据激光点间的平均点间距由式（1）确定；体素赋值方式如式（2）所示，将LiDAR数据中的各激光点映射到3D格网中，并对含有激光点的体素赋值激光点强度均值、不含有激光点的体素赋值0，进而把上述体素值离散化到{0,1,∙∙∙,255}得到体素灰度值。

式中：（Δx,Δy,Δz）为体素沿x、y和z方向的分辨率；Axy(Axz,Ayz）是LiDAR数据的XY(XZ,YZ）平面投影所得2D点集沿轴向最小外接矩形的面积；n为点数，式中取最小值是为了使构建的体素结构与LiDAR数据间的精度损失更少。

式中：i是激光点的索引；（xi,yi,zi）是第i个激光点的坐标；（ri,ci,li）是第i个体素的坐标（行、列和层号）。

另外，机载LiDAR数据中通常包含异常数据，如图1所示。该数据源于低空飞行物、车辆等目标对激光束的直接反射或树干等对激光束的多次反射，为非真实目标信息，其存在会影响体素结构构建的效率及准确性，须在体素结构构建前予以剔除。剔除方案：统计LiDAR数据中各激光点高程的频次，并通过直方图可视化显示；确定与真实目标对应的最高（低）高程阈值Th(Tl）；对各激光点，若其高程高于Th或低于Tl，则判为异常点，予以剔除，否则保留。

图1Samp31异常数据分布侧视图

2.2灰度体素结构分割

分割的目标就是将3D连通且灰度特征相似的体素合并到一个3D连通区域。假设2.1节所得灰度体素结构中共包含l个3D连通区域，分割的任务就是将l个标签分配给体素结构中的各个体素，以便属于同一个3D连通区域的体素具有相同的标签，而属于不同的3D连通区域的体素具有不同的标签。详细方案：依次扫描灰度体素结构中的非0值体素，直到扫描至第k个（k=1,2,…）未被标记的体素。假设该体素的体素值为u、标签L1,…,Ld-1（其中，d是3D连通区域的索引，1≤d≤l）已被使用，选择一个新的标签Ld。调用程序LABEL(k,u,Ld）（图2）标记与第k个体素3D连通（6、18、26或其它邻域尺度的连通）且体素值相近（相似性准则将在后文确定）的体素的标签为Ld。继续扫描灰度体素结构中的未被标记的非0值体素，直到所有的体素都被标记，得到若干个3D连通区域。

图23D连通区域标记LABEL(k,u,Ld）程序流程

在上述分割过程中应用不同的邻域尺度及体素值相似性准则会得到不同的分割结果并由此影响后继的滤波精度。为此，需面向应用确定最优值以保证地面的提取质量。最佳邻域尺度的确定见实验3.2节。最佳的体素值相似准则方案如下（以实验数据Area3为例）：统计灰度体素结构中的非0值体素的灰度值的频率，并可视化显示，如图3所示。

图3Area3的非0值灰度值频率直方图

由图3可知：地物目标的的灰度分布呈现多峰性（4个峰，分别位于1、28、80和133附近）且为多峰正态混合分布。则可借助高斯混合模型[13]对该分布进行拟合，估计所得4个高斯分量的参数记做（μ2，σ2），其中，j为高斯分布的索引，j=1、2、3、4。为了保证属于同一正态分布的体素被分割到一个3D连通区域，可利用图1中的峰值和谷值（位于10、61和103附近，分别记做v1,v2和v3，可利用局部极小值法确定）确定各个高斯分布的灰度范围。以第二个高斯分布为例，若令μ2-m2l×σ2=v1，μ2+m2r×σ2=v2,则可确定m2l和m2r，由此可确定第二个高斯分布的灰度范围[μ2-m2σ2，μ2+m2σ2]，其中，乘数m2=average{m2l,m2r}。同理可确定其它高斯分布的灰度范围，并记做[rjl,rjr]。但是，相邻的高斯分布间可能存在重叠。为了避免分布重叠且有利于后继的滤波处理，各高斯分布的范围分别设为（0,r2l),[r2l,r2r],(r2r,r3r],(r3r,r4r]。原因在于：第2个高斯分布对应地面目标（该先验知识可通过目视灰度体素模型获取），第1和第3个高斯分布的灰度范围设置均以第2个为基准设置以便准确地获取地面目标。而第3和第4个高斯分布的分布重叠并未处理的原因在于算法目的在于分离地面目标，此处即使植被和建筑物被混分也不影响算法的滤波精度。由此，前文所述的体素值相似性准则确定转换为了求取各个地物目标的灰度范围，图2中的与u接近条件等价于“体素值位于第k个体素所对应的目标的统计灰度范围内”。

2.3地面目标的3D连通区域检测

利用地面低于周边地物高度的特性，从分割所得的3D连通区域中分离地面目标对应的3D连通区域。详细方案如下：

对任一3D连通区域，若其轮廓线的高程低于周围地物一定的高度（高差阈值Te,Te=2m），则作为候选，否则判定为非地面连通区域。其中，轮廓线的高程指的是位于轮廓线上的体素的平均高程。周围地物的高程可用下述方案获取：利用结构元素[111;111;111]对该3D连通区域做膨胀，记膨胀处理所得外轮廓线的体素集合为Cs={va(ra,ca,la),a=1,2,…,q}，其中，s代表第s个3D连通区域，a是第s个3D连通区域的外轮廓线中各个体素的索引，q是第s个3D连通区域的外轮廓线包含的体素个数。va(ra,ca,la）∈Cs，搜寻与其平面坐标相同的非0值体素，上述体素的平均高程即是第s个3D连通区域的周围地物的平均高程。

计算各个候选3D连通区域的重心，若某一连通区域的重心高程值与其一定水平邻域（如5×5）内非0值体素的平均高程的高差小于等于某一阈值，则判定该3D连通区域对应地面，否则对应非地面。

3、实验数据、结果及讨论

3.1实验数据

采用不同密度的两组实测LiDAR点云数据作为实验数据检验提出的算法的有效性和可行性。

数据集1：国际摄影测量与遥感协会（InternationalSocietyforPhotogrammetryandRemoteSensing,ISPRS)III/3工作组提供的滤波测试数据（https://www.itc.nl/isprs/wgIII-3/filtertest/downloadsites/）[14]:Samp11、12、21、22、23、24、31、41和42，图4(a）以Samp31为例进行展示。这些区域包括了不同场景中滤波可能遇到的主要困难，如粗差点的影响、复杂不规则形状的大型建筑物、地物与地面相连、过街天桥、低矮植被等。点云数据记录了首、末次回波及其强度信息，点云密度为0.67点/m2。另外，该数据已被手工分为地面和非地面点两类，可用作标准数据定量评价提出的算法的精度。

图4实验（点云）数据示例

数据集2:ISPRSIII/4工作组提供的目标分类算法测试数据[15]（数据需申请，http://www2.isprs.org/commissions/comm3/wg4/data-request-form2.html):Area2和Area3（网站未提供与Area1对应的地面真实数据，无法进行定量精度评价，因此该数据未被采用），图4(b）以Area3为例进行展示。其中，Area2为被树木围绕的高层城市住宅；Area3为包含有独立房屋和许多周围树木的住宅区。点云数据记录了多次回波及其强度信息，点云密度为4点/m2。另外，该数据已被准确的分类为沥青地面、自然地面、建筑物、植被等类别，可用其中的地面真实数据作为参考数据定量评价提出的算法的精度。

3.2滤波过程及结果

实验以Samp31和Area3为例展示滤波过程及其结果，并对“邻域尺度”参数进行敏感性测试及分析。首先，Samp31和Area3经异常数据剔除、空间范围确定、体素分辨率计算、3D格网划分、点云映射及体素赋值等步骤处理得到灰度体素结构，如图5所示。其中，各步骤的处理结果见表1。可目视图5获取各地物目标的灰度分布情况，并作为先验知识用于随后的体素值相似性准则确定。

然后，对灰度体素结构进行分割。如2.2节所述，分割结果与各地物目标的灰度范围和邻域尺度有关。其中，前者可由2.2节所述方案确定：由图3可知第二个高斯分布对应地面目标，因此需首先确定其分布参数（μ2=33.34,σ2=12.79），然后结合谷值10和61确定乘数m2l=1.82,m2r=2.16，并由此确定乘数m2=(1.82+2.16）/2=1.99，进而确定第二个高斯分布的灰度范围为[33.34-1.99×12.79,33.34+1.99×12.79]，即[8,59]。同理可确定其它高斯分布的灰度分布范围。避免灰度分布重叠后的各个高斯分布的灰度范围分别为（0,8),[8,59],(59,103],(103,255]。后者的最优值实验确定方案如下：在相同条件下，不同的邻域尺度分别被应用到各个测试数据，对应的滤波结果的误差见表2。此处提出的算法所得地面数据是用体素表达的，而参考数据中则是离散的激光点。为和参考数据做对比进而计算提出的算法的误差，此处首先统计算法所得地面体素内的激光点，然后和参考数据对比，进而用I类误差（将地面点错分为非地面点比例，Ie）、II类误差（将非地面点错分为地面点比例，IIe）及总误差（错分的地面点的比例，Te）等指标定量描述算法误差。

图5构建所得灰度体素模型

表1灰度体素模型构建过程及结果

由表2可知：(1)数据集1中，采用6、18、26、51及56邻域时，提出的算法的平均总误差分别为0.0907、0.0765、0.0723、0.0685及0.0800。也即，从总误差指标来看，点云密度为0.67点/m2时，51邻域是最佳邻域尺度。(2)数据集2中，采用6、18、26、51及56邻域时，提出的算法的平均总误差分别为0.1029、0.1025、0.1019、0.0838及0.0962。也即，从总误差指标来看，点云密度为4点/m2时，51邻域是最佳邻域尺度。(3)邻域尺度的增加并不意味着总误差的必然降低。提出的算法的思想是：地面目标信息可以通过在灰度体素模型中定义的连通性和灰度相似度来传播。以6邻域为例，地面目标信息的传播只能由中心体素向上、下或4个基本方向传播。由此导致采用6邻域只能将平坦地区的地面体素合并到一个3D连通区域并继而被正确分离。如果能增加邻域尺度，如18、26邻域，其传播方向的增加可能会纳入更多的地面体素，从而降低滤波的I类误差（如表2中采用18、26邻域时的Ie明显低于采用6邻域算法的Ie）。但是，如果邻域尺度过大，则可能将一些非地面体素错分为地面体素，从而导致II类误差增加从而降低滤波的准确性（如表2中56邻域算法的IIe明显高于51邻域算法的IIe）。

51邻域尺度下的体素模型分割结果见图6。由图7可知：地面目标被分割为多个独立的3D连通区域。

最后，基于各个3D连通区域的特性分离地面目标。Samp11、12、21、22、23、24、31、41、42、Area2和Area3所得地面体素数分别为18517、18728、7095、15854、11253、3893、12366、3096、9905、39729、39154,Samp31和Area3所得地面体素见图7。

表2提出的算法在不同邻域尺度下的的误差统计

图6灰度体素模型分割结果（各连通区域用不同颜色表示）

图7提出的滤波算法的滤波结果

3.3滤波算法定量精度评价

为了定量评价提出的算法精度，以基于像元的评价方法[16]进一步统计了滤波结果的完整率（Rcom）、正确率（Rcor）、质量（Rq）及Kappa系数，如表3所示。

由表3可知，数据集1的平均完整率、正确率、质量及Kappa系数分别为0.9611、0.9248、0.8934及0.8484；数据集2的平均完整率、正确率、质量及Kappa系数分别为0.8490、0.8531、0.7404及0.7922。为了探究影响算法完整率和正确率的因素，考察了滤波结果的误差分布情况，如图8所示。

结合图8分析可知，影响提出的算法正确率的主要因素包括：某些和地面目标相邻且强度接近的非地面目标被错分为地面目标；算法中与地面对应的连通区域分离仅利用了地面低于周围地物的特性，可能导致错分。影响提出的算法完整率的因素包括：某些高程突变的地面被错分。

表3基于像元的评价方式下的提出算法的算法精度

图8误差分布顶视图

3.4滤波算法精度对比

针对数据集1，提出的算法同已有经典滤波算法的总误差对比见表4。其中，Axelsson[10]为2003年前发表的精度最高的滤波算法，其他则为2003年后发表的在多个测试数据中误差低于Axelsson[10]的算法。

由表4可知，提出的算法对比Axelsson[10]改进了9个样本中的3个样本精度；提出的算法对比已有的精度最高的3D滤波算法[3]精度相当，改进了9个样本中的2个样本精度。

鉴于目前利用数据集二的均为分类算法（分为建筑物、树、低矮植被、天然地面（naturalground）及沥青地面（asphaltground）等类），为了将提出的算法同上述算法（参见链接，Area2和Area3的链接分别为http://www2.isprs.org/commissions/comm3/wg4/results/a2＿detect.html和http://www2.isprs.org/commissions/comm3/wg4/results/a3＿detect.html）做对比，本文在提出的算法基础上进一步利用天然地面的激光反射强度值大于沥青地面的激光反射强度值这一特性对二者加以区分，进而统计上述检测的完整率、正确率及质量并和已有算法做对比，如表5所示。其中，天然和沥青地面的区分细则为：若某一已判做地面的3D连通区域的灰度（3D连通区域内各个体素的强度均值）大于某阈值，则为天然地面，否则为沥青地面。阈值的选取规则为：统计已判做地面的各3D连通区域的灰度频率，设置频率直方图中的谷值作为灰度阈值Ti。

由表5可知，提出的算法对比其他已有分类算法，其质量指标均为最高值，这也说明了算法的有效性。

综合上述算法精度对比实验可知，提出的算法在低密度点云数据滤波中与已有的3D滤波算法精度相当，但随着点云密度的增加，其滤波精度明显高于其它滤波算法，优势明显。

表4提出的算法和其他算法的总误差对比

表5提出的算法与其他算法的自然地面分类精度对比

3.5讨论

提出的算法的滤波结果由输入参数如阈值（Te、Ti）、目标的统计灰度范围和邻域尺度等参数决定。其中，Ti和目标的统计灰度范围是根据实际数据源设置的，使用本文给出的设置方案，可以很容易地确定“数据源”型阈值，并不能限制提出的算法的普适性。Te的设定是经验型的，因为地面通常比其周围的建筑物至少低3m，考虑到地面周围还有低矮植被，所以本文设置高度差阈值为2m。点云密度约为0.67和4点/m2时的邻域尺度可以直接使用51邻域，因为本文已证明51邻域为滤波的最佳邻域尺度。因此，本文提出的算法与地面的形态、坡度及连续性无关，可适用于不同复杂度场景的地面目标提取，具有普遍性。

4、结论

本文面向机载LiDAR数据提出了一种3D滤波算法。算法首先将机载LiDAR数据规则化为灰度体素结构，然后利用体素间的连通性和灰度相似性将其分割成多个3D连通区域，最后根据地面和其他目标间存在高差这一特性提取与地面对应的3D连通区域。该算法基于灰度体素结构设计，以3D连通区域分割理论为基础，综合利用了LiDAR点云数据包含的高程及强度信息。对比仅利用几何信息的滤波算法，其将强度信息用于滤波辅助决策，从而为地面和非地面目标的区分提供了更有效的信息，因而可适用于更复杂的场景。实验采用不同密度的ISPRS基准测试数据测试了邻域尺度参数的敏感性及提出的滤波算法的有效性，并和其他经典的滤波算法的精度进行对比。定量评价的结果表明，51邻域为最佳空间邻域尺度；点云密度为0.67(4）点/m2的数据集1(2）的滤波的平均完整率、正确率及质量分别为0.9611(0.8490）、0.9248(0.8531）及0.8934(0.7404）；提出的算法相对其他经典滤波算法在高密度点云数据滤波时表现更佳。总体而言，提出的算法有助于综合利用几何及辐射信息以提高复杂场景下地面提取结果的准确性。但是，算法仅支持高程或强度不同的地面目标的有效分离。为支持更复杂场景的地面目标分类，后续研究可尝试构建点云数据的多值体素结构模型（即体素赋值体素内激光点的色彩信息，可通过LiDAR和光学影像融合数据或者多光谱LiDAR数据等方式获取激光点的色彩信息）并在此基础上提出相应的地面目标分离算法。

参考文献：

[5]刘凯斯,王彦兵,宫辉力,等.机载LiDAR点云数据的二面角滤波算法[J].地球信息科学学报,2018,20(4):414-421.

[6]王丽英,王圣,徐艳,等.结合体元数据结构的机载LIDAR建筑物检测[J].中国图象图形学报,2017,22(10):1436-1446.

[7]张杰.机载LiDAR点云滤波及特征提取技术研究[D].徐州:中国矿业大学,2017.

[8]孙美玲.机载LiDAR数据滤波及城区汽车目标检测方法研究[D].成都:西南交通大学,2014.

[9]孙蒙,顾和和.基于微分形态学断面的机载LiDAR数据滤波新方法[J].大地测量与地球动力学,2016,36(7):591-594,599.

[12]周晓明.机载激光雷达点云数据滤波算法的研究与应用[D].郑州:解放军信息工程大学,2011.

[13]赵泉华,石雪,王玉,等.可变类空间约束高斯混合模型遥感图像分割[J].通信学报,2017,38(2):34-43.

王丽英,赵元丁.一种灰度体素结构分割模型下的机载LiDAR3D滤波算法[J].地球信息科学学报,2020,22(11):2118-2127.

基金:辽宁省教育厅科学技术研究项目(LJ2019JL015).