基于多尺度虚拟格网的点云自适应滤波算法

Adaptive Filtering Algorithm of Point Cloud Based on Multi-Scale Virtual Grid

摘要: 点云滤波是LiDAR点云数据处理中的重要步骤，针对传统曲面拟合滤波算法参数求解不合理，滤波阈值单一且自适应较差的问题，提出一种改进自适应阈值滤波算法。首先对预处理之后的点云数据引入虚拟格网进行分割并依据邻域格网选取种子点，然后采用混合最小二乘法对曲面拟合参数进行解算，计算真实高程与拟合高程的差值并结合k-means聚类与正态分布确定自适应阈值，最后采用多级滤波的策略逐级改变虚拟网格大小实现点云的高精度滤波。将该算法与其它经典滤波算法进行对比分析，试验结果表明：本文算法在不同场景下有良好的滤波精度，且稳定性较好。

Abstract: Point cloud filtering is an important step in LiDAR point cloud data processing. Aiming at the prob-lems of unreasonable parameter solving, single filtering threshold and poor adaptability of tradi-tional surface fitting filtering algorithm, an improved adaptive threshold filtering algorithm was proposed. Firstly, the pre-processed point cloud data is introduced into the virtual grid for segmen-tation, and the seed points are selected according to the neighborhood grid. Then, the mixed least squares method is used to solve the fitting parameters of the surface, and the difference between the real elevation and the fitted elevation is calculated, and the adaptive threshold is determined by combining K-means clustering and normal distribution. Finally, the multi-stage filtering strategy is used to change the size of virtual grid step by step to achieve high precision filtering of point clouds. The proposed algorithm is compared with other classical filtering algorithms, the experi-mental results show that the proposed algorithm has good filtering accuracy and good stability in different scenes.

1. 引言

激光雷达(LiDAR)作为一种新兴的测量手段，由于其精度高、可靠性强、可实时快速获取大范围空间三维信息与全天候作业等特点，现已被农林业与三维城市建模等广泛运用 [1]。点云滤波是点云数据处理中的重要步骤，国内外学者也对此开展了大量的研究。常见的滤波算法有基于坡度的滤波算法 [2]，基于不规则三角网的滤波算法 [3]，基于数学形态学的算法 [4]，基于曲面拟合的滤波算法 [5] 与布料模拟滤波算法 [6] 等。算法各式各样，但其算法思想都是围绕一定区域内地物点高程高于地面点这一认知而展开。

传统的坡度滤波方法原理简单、效率高但阈值的选取与自适应性有待提高。基于不规则三角网的滤波方法效果好，但对于不连续地形与地形复杂区域适用性差且对异常值敏感。基于数学形态学的滤波算法依赖于移动窗口尺寸的选择，窗口尺寸过大易过滤掉地形凸包，过小则会忽略大型建筑物。布料模拟滤波算法将点云倒置，通过多次迭代确定布料节点最终位置，最后计算点云与布料点的距离并设定阈值进行滤波。该方法效果好，但需人为设定的参数过多。基于曲面拟合的滤波算法能较大程度上还原地形特征，但滤波精度依赖于参数解算方法与种子点的选取，滤波阈值的自适应性也有待提高。针对曲面拟合算法的以上问题，文献 [7] 将曲面拟合算法与移动曲面相结合，提升了算法适应性，文献 [8] 在传统算法基础上加入窗口邻域，完善种子点的选取并采用自适应阈值进行滤波，提升了滤波精度与普适性。文献 [9] 将混合最小二乘法用于曲面拟合参数解算并利用点云曲率值计算自适应阈值，进一步减少了算法的人工参与。

在综合分析上述曲面拟合算法的优势后，本文提出一种多尺度格网曲面拟合的自适应阈值点云滤波算法。该算法使用金字塔策略构建多级虚拟格网，根据坡度因子大小设置邻域窗口选取种子点，采用混合最小二乘(least squares-total least squares, LS-TLS)进行曲面参数解算，借鉴k-means聚类思想确定自适应滤波阈值。与传统曲面拟合算法相比进一步提升了算法的自适应性。

2. 算法原理与流程

传统的曲面拟合算法仅采用6个种子点进行拟合 [10]，但该方法忽略了地形特征的连续性，拟合效果过度依赖中心窗口的种子点，而且对于不同的地形仍采用相同的窗口邻域与固定的窗口尺寸并不合理。针对以上问题，本文算法选择在三个方面进行改进：1) 依据坡度因子选择不同格网邻域，并采用金字塔格网策略，格网尺寸逐级递减，逐渐贴近真实地形。2) 采用LS-TLS将二次曲面拟合系数矩阵中的常数列与非常数列分开，只对非常数列进行修正。3) 采用k-means聚类与正态分布确定高差阈值，提升算法普适性。本文算法具体流程如图1所示。

Figure 1. Algorithm flow of this paper

图1. 本文算法流程图

2.1. 数据预处理

由于多路径效应等影响，点云数据在采集时不可避免的会存在一些噪声点 [11]，其中的低位噪声点容易被选取为种子点，影响曲面拟合效果与滤波精度。本文选择统计滤波的方法将噪声点剔除。对每一个点的邻域进行分析并计算其相邻点的距离，将平均距离在范围之外的点判定为离群点去除。

2.2. 多尺度虚拟格网的构建

虚拟格网的构建即点云格网化，忽略点云数据中的高程值并将其投影在X-Y二维平面上。格网的构建步骤如下：首先根据点云数据的X-Y最大最小值确定格网边界。其次设置初始格网大小来进行点云数据的格网ID号计算。其计算公式如下：

$\{\begin{array}{l}Col=int\left(\frac{x_n-x_{\min }}{l}\right)+1\\ Row=int\left(\frac{y_n-y_{\min }}{l}\right)+1\end{array}$ (1)

式中：Col代表行号；Row代表列号； $x_n$ 、 $y_n$ 代表点云的坐标值；l代表格网尺寸； $x_{\min }$ 、 $y_{\min }$ 代表点云坐标的最小值。

大多数移动格网将格网邻域设置为2*3或3*3，其格网尺寸也固定不变。然而考虑到地形特征的连续性与完整性，固定的格网尺寸与邻域在地形起伏较大时容易忽略地形特征点，导致拟合后的平面与实际地面相差过大。因此本文的格网尺寸与邻域根据地形坡度来进行设置。首先将邻域格网设置为3*3，计算中心格网最低点与周围8邻域格网最低点的坡度，再依据距离越远，权重越低的原则将这些坡度进行反距离加权。最后根据加权坡度值来设置格网的尺寸与邻域大小。加权坡度值计算公式如下：

$\{\begin{array}{l}slop{e}_{ave}={\displaystyle \underset{i=1}{\overset{8}{\sum }}{w}_i{s}_i}\\ w_i=\frac{\frac{1}{D_i}}{{\displaystyle \underset{i=1}{\overset{8}{\sum }}\frac{1}{D_i}}}\\ S_i=\frac{Z_{\min }-Z_{\min }\left(i\right)}{D_i}\end{array}$ (2)

式中， $w_i$ 为权重因子； $D_i$ 为中心格网最低点到第i个邻域格网最低点的距离； $s_i$ 为中心格网最低点与第i个邻域格网最低点计算得到的坡度值。领域格网如下图2所示。

Figure 2. Eight neighborhood grid

图2. 8邻域格网

根据计算出来的加权坡度值将格网尺寸与邻域大小进行调整。加权坡度值大则说明地势起伏大，格网内地面特征点多，此时保持格网邻域不变，通过缩小中心格网尺寸得到更多地面种子点；加权坡度值小则说明地势起伏较小，此时保持中心格网尺寸不变，通过扩大格网邻域搜索更多种子点。对于格网尺寸来说：如果加权坡度值未超过20˚，则选择初始格网尺寸，不做调整；若超过20˚，说明地势有起伏，则将格网尺寸缩小一半。对于邻域大小来说：如果加权坡度值超过10˚则领域格网设置为3*3；否则重新计算每个格网相对于5*5邻域格网的反距离加权平均坡度值。若计算后的加权坡度值超过5˚，则设置邻域格网为5*5，否则设置为7*7。最后选择各个格网内的最低点作为种子点进行二次曲面拟合计算。多尺度格网划分及种子点选取如下图3所示。

Figure 3. Grid segmentation and ground seed point selection

图3. 格网划分及地面种子点选取

2.3. TLS-LS曲面拟合

在选取足够的种子点之后，用这些种子点来拟合近似地形，二次曲面拟合方程如下：

$Z_i=a_0+a_1{X}_i+a_2{Y}_i+a_3{X}_i^2+a_4{Y}_i^2+a_5{X}_i{Y}_i$ (3)

式中： $a_j\left(j=0,1,\cdots ,5\right)$ 表示曲面方程的系数， $X_i$ 、 $Y_i$ 、 $Z_i$ 表示点云数据的三维坐标值。

将该式矩阵化又可表示为：

$\left[\begin{array}{c}{Z}_1\\ Z_2\\ ⋮\\ Z_N\end{array}\right]=\left[\begin{array}{cccccc}1& X_1& Y_1& X_1^2& Y_1^2& X_1{Y}_1\\ 1& X_2& Y_2& X_2^2& Y_2^2& X_2{Y}_2\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& X_n& Y_n& X_n^2& Y_n^2& X_n{Y}_n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\end{array}\right]$ (4)

式(4)可简化为

$L=AX$ (5)

式中A为系数矩阵，L为观测向量。考虑到点云在数据采集时存在不可避免的误差，即系数矩阵A与观测向量L存在误差，式(5)应该写为：

$L\approx AX$ (6)

又因为A中只是部分向量存在误差，其常数列不存在误差，因此将A分解为

$A=\left[\begin{array}{cc}{A}_1& A_2\end{array}\right]$ (7)

将对应的X分解为

$X={\left[\begin{array}{cc}{X}_1^{\mathrm{T}}& X_2^{\mathrm{T}}\end{array}\right]}^T$ (8)

分解后的混合最小二乘变量含误差(Error In Variables, EIV)模型如下：

$L+{\Delta }_L=A_1{X}_1+\left(A_2+{\Delta }_{A_2}\right)X_2$ (9)

式中，A₁为n × 1的常数列，A₂为n × 5的非常数矩阵， ${\Delta }_L$ 与 ${\Delta }_{A_2}$ 表示观测向量与系数矩阵A中非常数元素的误差；n为观测数。其参数求解过程如下：

首先构建增广矩阵 $C=\left[AL\right]$ ；然后对其进行QR分解。

$C=QR$ (10)

分解得到的Q为n × n的酉矩阵，R为7 × 7的上三角矩阵。则

$R=Q^{\mathrm{T}}C=Q^{\mathrm{T}}\left[\begin{array}{ccc}{A}_1& A_2& L\end{array}\right]=\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{1L}\\ 0& R_{22}& R_{2L}\end{array}\right]$ (11)

将此方程以分块的形式表示可得：

$R_{11}{X}_1+R_{12}{X}_2=R_{1L}$ (12)

$R_{22}{X}_2=R_{2L}$ (13)

对式(13)构造增广矩阵并进行奇异值分解可得到X₂的整体最小二乘估计值，最后将X₂回代至式(12)对X₁进行经典最小二乘可得其解。

混合最小二乘以不同的约束准则分别对X₁与X₂进行参数估计。该方法不仅考虑到了观测向量与系数矩阵的测量误差，还将系数矩阵一分为二进行参数解算。普通最小二乘法只考虑了因变量的误差，忽略了自变量的误差；整体最小二乘法考虑到了二者同时存在误差，却忽略了系数矩阵中的常数项。混合最小二乘法综合了经典最小二乘与整体最小二乘法，兼顾了自变量与因变量误差的同时将系数矩阵分开处理，采用该方法进行曲面方程参数解算是较为合适的。

2.4. 自适应高差阈值选取

大多数滤波算法在面对不同场景的点云数据时需手动设置不同的滤波阈值参数。为了减少人工对算法的干预，提升算法的普适性，需设置自适应阈值。文献 [12] 指出，在排除非地面点的影响后，大量的地面点呈现出正态分布。因此可认为，在经过了二次曲面拟合后，地面点的拟合高程值与真实高程值之差值较小，且它的高差数据应呈现正态分布；而非地面点的差值较大，且扰乱了地面点高差数据的正态分布。因此根据地面点与地物点的分层现象，先采用k-means聚类算法将两者高差数据分开，再结合正态分布思想来确定曲面拟合的高差阈值。

考虑到部分区域存在地物点与地面点较为接近的情况，例如低矮植被，其分层现象不明显。这种情况下不使用k-means聚类方法进行分类，直接计算高差值数据的均值 ${\mu }_2$ 和均方差 ${\sigma }_2$，将高差数据大于 ${\mu }_2+n{\sigma }_2$ 的点作为非地面点剔除。

在使用k-means聚类后可得到两种数据，聚类中心较小的高差数据与较大的高差数据，将中心较小的聚类结果视为地面点结果并计算出其均值 ${\mu }_1$ 和均方差 ${\sigma }_1$，将 ${\mu }_1+n{\sigma }_1$ 当作高差阈值进行滤波分类。利用正态分布确定阈值时，n可取2或3。小于该阈值的分为地面点，否则分为非地面点。然而在实际情况中可能存在格网内全是地面点的特殊情况，此时的格网拟合高差数据全部偏小。因此需设置一个最低高差阈值，若滤波格网内的最大拟合高差值小于该阈值，则滤波格网内的点全部分为地面点。通常情况下最低地物点应高于地面0.5 m，故本文将最低高差阈值设置为0.5 m。

3. 点云滤波试验及分析

3.1. 试验数据

为验证算法的有效性，选择国际摄影测量与遥感学会(International Society for Photogrammetry and Remote Sensing, ISPRS)发布的公开点云数据集进行实验，共15个样本，包括了城市及森林两种场景，其点云数据密度分别为0.67点/m²和0.18点/m²。数据集已经手动分类并贴上标签。在点云滤波之前先利用统计滤波进行数据预处理。数据预处理之后可根据点云样本数据的最大建筑物尺寸设置初始格网大小，其后的多级滤波格网尺寸则逐级缩小，例如二级格网大小为初始格网的一半，三级格网则缩小为初始格网的1/3，以此类推。所选的数据集中的15个样本数据的初始格网与数据特征如下表1所示。

Table 1. Terrain features and initial grid size

表1. 样本数据特征与初始格网大小

3.2. 试验结果及分析

为更好说明本文算法的滤波过程，选取数据集中的样本12与样本54为例进行展示。样本12为城市区域，存在大量建筑物，地势较为平坦，但含有大量零散地物。样本54为村庄区域，地势崎岖不平，建筑物旁有大量植被。共进行3次滤波，图中的白色点云为错误分类的地物点。由图中可看出，在第一次滤波后，样本中的大尺寸建筑物基本被剔除，但仍然存在部分错分情况，主要是建筑物旁的低矮植被与小物体有残留。而在经过2级与3级滤波之后，错分点明显减少，滤波处理相比第一次而言更加精细，地面细节得到了保留，滤波结果更加贴近真实地面。由此可知，多级滤波的策略逐步剔除地物点，还原真实地表，相比于单次滤波效果更佳。滤波结果如图4与图5所示。

(a) Samp12真实地面点 (b) 一级滤波后 (c) 二级滤波后 (d) 三级滤波后

Figure 4. Sample12 filtered ground point results

图4. Sample12滤波后地面点结果

(a) Samp54真实地面点 (b) 一级滤波后 (c) 二级滤波后 (d) 三级滤波后

Figure 5. Sample54 filtered ground point results

图5. Sample54滤波后地面点结果

为进一步对试验结果进行分析，本次试验采用ISPRS于2003年提出的误差评价标准对滤波结果进行定量分析。该标准将误差分为三类，一类误差为地面点错分为地物点的误差，二类误差为地物点错分为地面点的误差，一类误差与二类误差之和与总点云的比值称为总误差。表中a代表滤波后分类正确的地面点，b代表分类错误的地面点即地面点错分为地物点，c代表分类错误的地物点即地物点错分为地面点，d代表分类正确的地物点。e和f代表参考数据中地面点和地物点的个数，g和h代表滤波结果中地面点和地物点的个数，n代表总点云个数。评价标准如下表2所示。

Table 2. Filter error evaluation criteria

表2. 滤波误差评价标准

根据误差评价标准进行定量分析，可以看出sample21，sample31，sample42，sample51，sample61，sample71等样本滤波效果较好，总误差均达到了6%以下，而sample11与sample23效果较差，误差达到10%以上。其中的sample31与sample61均为不连续地形且地形有起伏变化，说明本文算法在面对不连续地形时依然有一定适用性。sample11地形坡度较大低矮植被较多且同一高程面上地面点与地物点错杂在一起，导致曲面拟合效果不佳，滤波分类存在难度。部分样本结果如下表3所示。

Table 3. Statistical table of filtering result error (Partial sample)/%

表3. 滤波结果误差统计表(部分样本)/%

3.3. 算法对比

将本文算法数据集的滤波结果与ISPRS发布的8种经典滤波算法进行算法精度对比，进一步对本文算法进行分析。由表4中可知本文算法的总误差均值为7.04，均方差为3.61，高于其中的7种经典滤波算法。其中Sample61精度最高，Axelsson [13] 算法精度最高，但算法时间长，计算量大，需迭代构建TIN计算角度与距离等进行滤波。Elmqvist [14]、Wack [15] 与Brovelli [16] 算法在森林和农村场景滤波精度不高，说明在植被较多与陡坡场景时适用性不足。Roggero [17] 算法需要提前知道样本数据信息设置阈值，不同的地形地貌需要不同的阈值，算法自适应性较差。Sithole [18] 算法在不连续地区与坡度地区滤波精度欠缺。Soho [19] 算法与Pfeifer [20] 算法精度整体较为稳定。综合来看，本文算法不论是在不连续地区、陡坡地区、多建筑物地区与低矮植被地区都能取得不错的滤波效果，滤波后能较好的保留地形细节，算法稳定性好。对比结果见表4所示。

Table 4. Comparison of filtering algorithms

表4. 各滤波算法对比

4. 结论

针对传统曲面拟合算法的不足，本文根据地形坡度构建多尺度虚拟格网与邻域选取种子点，探讨了TLS-LS方法相较于经典最小二乘法解算曲面方程的优点，借鉴正态分布思想实现高差阈值的自适应选取，减少了人工对算法的干预程度，研究了多级滤波策略逐级分离地面点与非地面点的可行性。与其它经典滤波算法的对比试验表明，本文算法在不同的场景下有着良好普适性且精度较高。但本文算法也有不足之处，在坡度变化较大、地面点与非地物点分层现象不明显且交杂在一起时，滤波效果有待提高。

基金项目

国家自然科学基金(41901225)，教育部规划基金项目(22YJAZH083)，贵州省科技计划项目(黔科合支撑[2022]一般204)。

NOTES

^*第一作者。

^#通讯作者。