植被茂密地区点云双重滤波方法研究

王云云，唐菲菲*，王章朋，肖 敏，唐天俊，王铜川

(1.重庆交通大学 土木工程学院，重庆 400074；2.上海宝冶集团有限公司，上海 201908；3.杭州钱塘测绘有限公司，杭州 310000)

引 言

机载激光雷达(light detection and ranging,LiDAR)系统是一种集激光测距技术、全球导航卫星定位系统动态定位技术、惯性导航系统姿态测量技术和计算机技术4种技术于一体的空间测量系统[1]，可获得高分辨率高精度的空间信息。目前，机载LiDAR技术已广泛应用于数字地面模型获取、道路提取、电力线提取、林业调查、3维城市模型建立、点云分类等地球空间信息学科的众多领域[2]。

尽管机载LiDAR系统硬件已发展较为成熟，然而，机载LiDAR点云数据的滤波算法仍相对比较落后，许多学者对点云滤波算法进行了研究，并提出了多种滤波算法，大概可概括为两种：第1种是利用点云数据的几何特性(3维坐标)进行的滤波；第2种是利用点云数据的回波特性(回波强度、多回波信息)进行的滤波[3-6]。经典的滤波算法都是基于几何特性提出的，如：基于坡度的方法、基于数学形态学的方法、基于曲面拟合的方法以及基于不规则三角网(triangulated irregular network，TIN)算法。近年来，硬件设备快速发展，点云的回波特性得到了学者们的重视。LIU等人利用不同介质的反射率，提出了一种激光强度信息与点云高程信息结合对LiDAR数据进行分类的算法[7]。ZHANG等人提出利用首末次回波的高程差实现了机载LiDAR数据的分类[8]。TANG等人基于点云数据回波信息，将点云以不同分辨率的体素为单位实现点云滤波，在森林地区达到了较好的滤波效果[9-10]。LIN等人利用激光脉冲的多回波特性，预先剔除一部分非地面点的激光脚点，增强了滤波的效果[11]。SUN等人提出了一种基于多回波及Fisher判别的滤波算法，验证了该方法在陡坡点云滤波中的适应性与有效性[12]。LI等人提出一种基于回波强度的k均值聚类算法，该方法可提高研究区点云滤波精度并减少地面脚点的数据量[13]。

在前人的基础上，作者提出一种结合几何特性与回波特性实现点云数据滤波的方法。针对回波特性，知道多次回波会产生不同反射阶段，导致激光的能量均有所减弱，就会使单次回波和末次回波之间强度存在差异；且植被密集区域通过多回波分离，发现单次回波的地物类型多样(道路、建筑、植被)，而末次回波的地物较单一(植被)。因此，本文中对单次回波和末次回波点云数据分别结合不同特征利用不同方法实现自适应粗滤波。另外，不规则三角网渐进加密滤波算法对各种不同地形能取得较好的滤波效果，但其种子点的选取至关重要，所以本文中引入分块思想，利用局域点云数据密度来选取种子点，使种子点选取更加准确，进而提高不规则三角网渐进加密精滤波算法的结果精度。

1 点云双重滤波

利用多回波信息进行滤波之前，需根据原始点云的回波次数和第几次回波数据序列分离出各回波信息，即回波分离。通过回波分离结果分析得出：单次回波包含地面点、大部分植被冠层点以及建筑物等多种地物的回波信号；首次回波则大概率都是植被冠层、树干等非地面点；末次回波为地表点和植被点的集合；而中间次回波数据量极其少，因此不做分析处理。基于此，本文中只选取单次回波和末次回波的激光脚点参与滤波算法，且针对单次回波和末次回波所包含地物类别的差异性，利用不同的方法分别对单次回波和末次回波进行粗滤波，提出了一种基于偏度平衡-最大类间方差的联合粗滤波方法；然后，将粗滤波后的单次回波和末次回波数据叠加，参与后续不规则三角网渐进加密精滤波处理，实现点云自适应双重滤波方法。操作流程如图1所示。

Fig.1 Flow chart of filtering algorithm

1.1 点云偏度平衡-Otsu联合粗滤波

1.1.1 基于偏度平衡的单次回波粗滤波 由于原始LiDAR点云中存在噪声点，首先利用TERRASOLID软件通过目视检查剔除部分相对明显的高位和低位噪声点，然后结合强度信息进行频数统计分析，剔除强度异常值。剔除噪声点后，对单次回波的实验数据结合强度信息利用偏度平衡理论自动计算强度阈值。

2006年,BARTELS等人最先提出了一种基于偏度平衡的滤波算法[14]，通过引入偏度统计量与峰度统计量的平衡思想，该方法凭借其非监督、无阈值的优点被广泛使用。根据统计学的中心极限定理，BARTELS提出了两个假设：(1)自然状态下，无序离散分布的点云数据中的地面点数据的概率密度函数总是呈正态分布;(2)由于植被、建筑物等非地面点云数据的存在，影响了原始点云数据中的地面点数据的正态分布，导致LiDAR点云数据的概率密度函数从原本近似正态分布呈现出偏态分布[15]。因此，将原始LiDAR点云数据的概率密度函数呈现的偏态分布曲线校正为近似正态分布的过程，就是把点云数据中的非地面点滤除来实现点云数据滤波的过程。本文中基于偏度平衡理论剔除部分非地面点，实现单次回波点云数据的粗滤波。

基于偏度平衡的点云粗滤波方法的具体实现步骤如下：(1)将去噪后的单次回波的点云数据(包括3维坐标xi,yi,zi和强度信息Ii(i=1,…,N)，N为点云总数)单独标记出来;(2)对点云强度信息进行统计，获取其最大、最小强度值Imax,Imin，并将最小强度值Imin设为初始强度阈值I;(3)计算单次回波的点云强度偏度值Sk，然后判断Sk值的正负。计算公式如下：

式中，Sk为偏度，μ3为3阶中心矩，σ为标准差;(4)若Sk<0，则激光强度阈值I自动加1，即I=I+1；若Sk>0，输出结果Ith，即点云数据的最佳强度阈值;(5)将所有反射强度值大于最佳强度阈值Ith的点云数据保留,并参与后续精滤波。

1.1.2 基于最大类间方差的多次回波粗滤波 最大类间方差法是由日本学者OTSU于1979年提出的[16]，是一种图像分割中确定阈值的最佳算法,简称Otsu法。本文中利用该方法获取最佳阈值的优势，计算首末次回波的高差阈值，进而实现末次回波点云的粗滤波。

记T为高差阈值，大于高差阈值的点所占比例为ω0,平均高差为μ0；小于高差阈值的点比例为ω1，平均高差为μ1，首末次回波高差的总平均高差为μ，计算类间方差δ，则有[17]：

μ=ω0×μ0+ω1×μ1

(2)

δ=ω0×(μ0-μ)2+ω1×(μ0-μ)2

(3)

将两式联立可得：

δ=ω0×ω1×(μ0-μ1)2

(4)

当方差最大时，可以认为获取的高差为最佳阈值，所以利用最大类间方差法可自适应计算阈值的特性，实现首末次高差的自适应阈值。其具体的算法步骤如下：(1)首先将同一束激光的首次回波和末次回波一一对应，进行标记;(2)计算同束激光对应首末次回波的高程之差，统计其最大、最小值分别为Hmax和Hmin;(3)假设首末次高差进行分类的最佳阈值记为T，则有Hmin

1.2 点云精滤波

经典不规则三角网渐进加密滤波算法是AXELSSON于2000年提出的[18]，该方法对不同场景的适应能力较强、滤波质量较好[19]。通过选取初始地面种子点构建初始不规则三角网，由于初始三角网的精度会影响后续加密过程的判断，则地面种子点的选取一定程度上会影响滤波的效果[20]。在植被茂密地区，激光点云数据分布很不均匀，因此，本文中引入分块点云密度，基于局部点云密度设置自适应格网大小，从而获取可靠地面种子点。

具体步骤如下：

(1)统计粗滤波后获取的点云数量N和点云数据所占面积S，计算参与精滤波的激光点云分布的整体平均密度ρ=N/S。

(2)点云数据按最大建筑物边长L进行分块，记每块面积为Sj(j=1,2,3,…)，计算分块后激光点云分布的局部平均密度ρj=Nj/Sj(j=1,2,3,…),其中,Nj为每块实验区的点云数量。

(3)点云分布的局部平均密度与整体平均密度ρ进行比较，若ρj>ρ,则计算规则格网大小为gj=int[M÷ρ]-B(j=1,2,3,…)；若ρj<ρ,则计算规则格网大小为gj=[M÷ρ]+B,其中,M表示每个格网中必须包含的最少点数量，B为常数。

(4)点云数据格网分配并建立索引，根据计算得到的虚拟格网大小gj(j=1,2,…)、分块边长L以及局部点云最大高差值Hj，可以将局部区域分别确定一个mj×nj×1(j=1,2,…)的虚拟格网，其中,mj和nj的计算方式为：

每块虚拟格网的编号(Rj,Cj)的计算公式为：

式中，⎣」为向下取整；xj,k为第j块第k个点的x坐标；xj,min为第j块的坐标x的最小值；yj,k为第j块第k个点的y坐标；yj,min为第j块的坐标y的最小值。

(5)遍历格网。通过索引对每个格网中点云数据的高程值z从小到大进行排序，则第1个数据值即为格网中的最小高程值对应激光点，所有格网获取的第一个数据构成地面种子点集。

(6)将最后获取的地面种子点构建初始TIN,然后不断迭代加密实现点云数据精滤波，得到地面点集。

2 实验数据与结果分析

2.1 实验数据

本文中选取了3组具有不同地形特征的数据集，如图2所示。图2a中的数据集1为平坦区域，包含建筑物、道路、少许植被等地物；图2b中的数据集2为较平坦区域，但包含大量植被；图2c中的数据集3为地形起伏区域，包含植被、少许建筑物、道路等地物。3组实验数据集的单次回波和末次回波点云总数分别为：3882510，3716640，4570200。其中对应单次回波激光脚点数分别为：3721050，3327300，3950280，则末次回波激光脚点数分别为：161460，389340，619920，且每条数据均记录了3维位置坐标x,y,z,回波号和回波强度信息。

Fig.2 Three experimental data sets

2.2 点云粗滤波结果

2.2.1 基于偏度平衡的单次回波滤波结果 实验区的单次回波点云数据中主要包含建筑物、道路以及大量植被点等地物，由于不同介质有不同反射率，如：混凝土、沥青等物质的反射率在20%左右；植被不同种类各不相同，大概在30%～60%范围；而黏土的反射率则接近于75%，则部分植被与土壤的反射率较为接近[7]，则依据不同地物的反射率存在差异的特性，可利用强度信息滤除一大部分地物点。而在植被覆盖地区，多回波中首次回波点云数据一般为树冠、树干等，且考虑到首次回波没有激光能量消耗，将其作为单次回波的一个比较对象。图3为单次回波和首次回波点云强度示意图。其中条形图表示单次回波，折线表示首次回波，横坐标轴为无量纲量强度值，纵坐标为对应强度值点云个数。

Fig.3 Schematic diagram of single and first echo intensity

从单次回波强度分布可以发现，其强度在0～20范围的数据量极少，而255数据量极大，造成强度变化幅度过大，影响整体回波强度分布规律，则将单次回波强度值为0～20和255的数据作为噪声点或异常点剔除。另外，图中单次回波有两个峰值，与首次回波对应分析可得出，峰值在200左右大概率为地面点，峰值在40左右大概率为地物点，而首次回波在100左右出现数据量急速减少趋势，可认为强度值在100之前几乎均为非地面点，对单次回波强度值为100左右的点云进行可视化分析，通过强度数据和可视化判别的综合分析，最后选取强度值在80～254范围的单次回波点云数据参与粗滤波。

通过MATLAB平台编写该实验的程序，实现基于偏度平衡的粗滤波方法。将实验区的单次回波数据序列包含3维坐标x,y,z和强度值(80～254)加载到该滤波方法程序中，作为输入数据，最后输出相应的强度阈值分别为160,145,130，结果如图4所示。从图中可以看出,利用强度信息能较好地滤除建筑物、道路及大量植被点等地物，且较好保留了地面形态细节，增加了后续种子点选取的准确性。

Fig.4 Experimental data single echo result graph

2.2.2 基于最大类间方差的多次回波滤波结果 实验数据的首末次回波主要分布在植被覆盖区域，计算同一激光束对应的首末次回波对应高差值，采用最大类间方差法对实验数据的首末次回波的高差通过MATLAB编程实现阈值自适应，得到高差阈值分别为6.5m,7.0m,6.0m。将实验数据中高差大于阈值所对应的末次回波激光脚点提取出来，参与后续滤波处理，其结果如图5所示。

Fig.5 Figure of final echo results of experimental data

基于偏度平衡的单次回波粗滤波和基于Otsu法的多次回波粗滤波分别完成后，得到对应单次回波和末次回波的粗滤波结果并进行融合，实现了偏度平衡-Otsu联合粗滤波，如图6所示。从图中能明显看到道路和建筑物被剔除，地形结构得到了较好的保留，为后续地面种子点的选取提高了可靠性。

Fig.6 Figure of coarse filtering result of point cloud

2.3 点云精滤波结果

将粗滤波后的点云数据参与改进的不规则三角网渐进加密精滤波，获取地面点集。为获取较高可靠性的地面种子点，首先将实验区域按照20m×20m进行分块划分，然后分别计算分块后的局部点云密度，将局部点云密度ρi(i=1,2,3,…)与总平均密度ρ进行比较来确定格网大小，其中每个格网中必须包含的最少点数量M和常数B值的设置与平均密度相关，分别分别设置为150,1;100,1;100,2。最后将每个格网中数据按z值进行排序，选取最低点作为地面种子点，进行不规则三角网渐进加密滤波算法实现。图7为3组数据精滤波得到的地面点俯视图。从图中可以看出,该方法较好保留了地形结构，但在植被茂密区域，获取地面点较小甚至会出现一些空洞，这是以后需要进行改进的一个方向。

Fig.7 Figure of experimental results of point cloud precision filtering

对滤波方法进行结果评价时，仅通过定性分析具有一定片面性，可结合定量分析综合评定滤波结果。3类误差是评定点云滤波效果的一个重要指标，将误差分为第Ⅰ类误差、第Ⅱ类误差和总误差，其中Ⅰ类误差指地面点误分为非地面点，Ⅱ类误差指非地面点误分为地面点，总误差为两类误差的综合。表1为精度评定表。表中,e表示参考数据的地面点，即a+b；f表示参考数据的非地面点，即c+d；p表示本文中滤波结果的地面点，即a+c；q表示本文中滤波结果的非地面点，即b+d；h为总数据[21]；且a和d分别为正确滤波的地面点和非地面点，b和c则为错误分类的地面点和非地面点。

Table 1 Precision evaluation

3类误差的计算公式如下：

利用该方法分别对3组实验数据经过精度进行评定，计算的3类误差分别如表2所示。从表中可以看出，该滤波方法得到的滤波结果综合效果较好，数据集1的地形相对平坦且植被覆盖较低，其3类误差都相对较小；数据集2则植被覆盖率极高，相比其它两组数据集的3类误差相对较大，但第Ⅱ类误差即误分为地面点的误差相对来说也是较低的，说明获取的地面点准确度相对比较高；数据集3的地形有一定坡度且植被覆盖率较高，其结果精度相对较高，也较好保留了地形细节。因此，本文中提出的双重滤波方法不仅在平坦且植被覆盖较少的地区取得较好结果，而且在地形起伏且茂密植被覆盖区域仍能较好保留地形细节，但植被覆盖率过高时，可考虑与其它多源数据进行融合分析。

Table 2 Error evaluation table of experimental results

3 结 论

本文中的方法对单次回波和末次回波的激光脚点进行粗滤波，实现阈值自动化，且降低了后续种子点选取无法自动剔除建筑物和植被等非地面点的可能性，大大提高了精滤波运算效率和滤波精度，但该方法只适用于植被高度分布相对均匀的区域，高差阈值的确定方法需要进一步研究。另外，针对植被茂密地区滤波后出现的稀疏点云或者空洞区域如何处理仍需深入研究。