一种机载LiDAR点云自适应坡度阈值去植被算法

An Adaptive Slope Threshold Airborne LiDAR Point Cloud Devegetation Algorithm

摘要: 针对利用机载LiDAR点云数据制作高精度DEM时易受地形及茂密植被干扰问题，提出一种改进的基于坡度和地形的滤波算法：新算法按一定尺度将点云格网化后，引入地形判断因子，判断格网内地形并自适应调整格网的坡度阈值，从而避免过渡平滑或忽略较小地形特征导致的地形失真，以及在陡坡断崖处滤波效果不佳的问题。对FARO SINGAPORE PTE. LTD公司提供的机载LiDAR点云数据处理结果表明，相对于经典坡度滤波本文算法优势明显，在地形复杂且植被密度大的林区也能很好地保留地面点并滤除植被等非地面点，得到真实的地形信息。在陡坡断崖地形区域，本文算法总误差较布料模拟算法降低了2.68%，滤波效果较好。

Abstract: Aiming at the problem that the use of airborne LiDAR point cloud data to produce high-precision DEM is susceptible to terrain and dense vegetation interference, an improved filtering algorithm based on slope and terrain is proposed. After the point cloud is grid based on a certain scale, the terrain judgment factor is introduced to judge the terrain in the grid and adaptively adjust the gradient threshold of the grid, so as to avoid the problem of terrain distortion caused by smooth transition or ignoring small terrain features, and poor filtering effect on steep slopes and cliffs. The processing results of the airborne LiDAR point cloud data provided by FARO SINGAPORE PTE. LTD show that the algorithm in this paper has obvious advantages compared with the classical slope fil-tering. In the forest area with complex terrain and high vegetation density, the ground points can be well preserved and filtered out. Non-ground points such as vegetation can get real terrain infor-mation. In the area of steep slope and cliff terrain, the total error of the algorithm in this paper is reduced by 2.68% compared with the cloth simulation algorithm, and the filtering effect is better.

1. 引言

机载LiDAR因具备主动性、穿透性和高效性等特点，使其在林区数字高程模型(Digital Elevation Model, DEM)生产中显现出巨大技术优势。然而，林区复杂的地形和高植被覆盖度导致的DEM不精确问题制约其进一步发展，因此，如何利用高效去植被滤波算法构建林区高精度DEM逐渐成为近年来机载LiDAR技术应用中的研究热点问题之一 [1]。目前国内外学者提出了多种植被滤波的方法，其中包括基于形态学的滤波方法 [2] [3]、基于坡度的滤波方法 [4] [5] [6]、基于插值的滤波方法 [7] [8] [9] 以及基于分割的滤波方法 [10]。文献 [11] 提出一种以窗口化和地形坡度为基础的植被茂密区域点云滤波算法，该算法将大量点云数据格网化，将非坡度因素引起高程差异的相邻两点中高程值较大的点归类为非地面点，该方法在中矮植被覆盖区域有较强适用性，但对于地形复杂林区的适用性不强。文献 [12] 提出一种密集低矮植被区LiDAR点云滤波算法，此算法在低矮植被区域效果良好。文献 [13] 针对机载LiDAR数据采用单一渐进式形态学滤波难以生成高质量DEM的问题，提出一种基于空间向量投影的后处理组合滤波方法，可有效滤除对DEM构建精度影响较大的近地面点，但在地形破碎、具有陡峭不连续地形区域滤波效果较差。文献 [14] 充分考虑点云数据组织形式以及地表形态，借鉴基于坡度滤波的基本原理，实现坡度值的自适应更新。该算法在地势平坦城区效果良好。文献 [15] 为了抑制突出地形特征引起的遗漏误差，提出了一种基于多级克里格插值的改进形态学算法，该算法本质上是渐进形态学滤波算法和多级插值滤波算法的结合，在平坦城市地区以及农村地区效果良好。文献 [16] 针对布料模拟算法易受到崎岖地形影响的问题，提出一种基于地面点归一化的适应地形的滤波方法，该方法先高程归一化后使用布料模拟滤波从而消除地形影响。综上所述，国内外越来越多的学者开始重视在地形崎岖、植被覆盖度高的林区的机载LiDAR点云滤波算法研究，并进行了卓有成效的探索。但算法依旧受到地形崎岖的影响，本文在以往基于坡度的、渐进式数学形态学滤波算法基础上引入地形因子，提出一种基于格网化及坡度阈值的自适应去植被滤波算法，经验证，在具有崎岖地形的高植被覆盖度林区，可以获得良好的滤波效果。

2. 自适应坡度阈值去植被滤波算法

现有坡度滤波方法，在处理由非地形坡度引起高程差异的相邻两点时，将相对较高的一点视为植被点给予滤除，较低的一点视为地面点加以保留，经过不断比较和搜索，最终达到滤除植被的目的。这种算法对于地形单一和植被稀疏地区，滤波效果良好，但在具有复杂地形和植被茂密的地区，这种采用相对统一的坡度阈值的算法难以获得高质量去植被滤波结果。本文对现有坡度滤波方法进行改进，即在将点云格网化的同时，引入地形判断因子，针对不同地形起伏情况，通过地形判断因子调节坡度阈值决定一点是植被点或地面点(非植被点)，从而实现在复杂地形条件下的林区机载LiDAR点云自适应去植被滤波的目的。

本文算法是基于格网化及坡度变化的思想进行的：首先将点云格网化，在X-Y平面规划一定尺寸的格网，该格网总体上覆盖所有的点云在X-Y二维平面的投影，通过一定的计算将点云分配至每一个格网之中：

$\begin{array}{c}H=\mathrm{int}\left(\frac{x_i-x_{\min }}{L}\right)+1\\ V=\mathrm{int}\left(\frac{y_i-y_{\min }}{K}\right)+1\end{array}\}$ (1)

在式(1)中，H表示点处于的格网行序号，V表示点处于的格网列序号，x_i表示i点的X坐标，y_i表示i点的Y坐标，x_min和y_min分别表示整个点云中坐标最小值，L和K分别表示每一个格网的长度与宽度即格网的尺寸大小。格网化完成后，找到每个格网中最低点并视为地面点，计算每个格网中的点与相邻所有格网中最低点构成的向量与X-Y平面夹角 ${\theta }_i$，计算公式为：

$\cos {\theta }_i=\frac{\stackrel{\to }{u}\cdot \stackrel{\to }{v}}{\left|\stackrel{\to }{u}\right|\cdot \left|\stackrel{\to }{v}\right|}$ (2)

在式(2)中， $\stackrel{\to }{u}$ 为格网点与相邻格网中最低点构成的方向向量， $\stackrel{\to }{v}$ 为相邻格网中最低点所在的X-Y平面法向量。本文算法具体实现步骤如下：

1) 对点云数据预处理后，将其初步格网化，找到每个格网中最低点，遍历每一个点云格网(如图1所示)，将格网中的最低点与8邻域格网中最低点计算夹角 ${\beta }_i$ (见图2)，并根据格网中点云数据分布情况设置为 $h=2\times {\beta }_{i\max }$。其中 ${\beta }_{i\max }$ 为最大夹角。

Figure 1. Traversal of direction vector at the lowest point of grid

图1. 格网最低点方向向量遍历

Figure 2. β_i Schematic diagram of X-Z section

图2. β_i在X-Z剖面示意图

2) 判断 ${\beta }_i$ 并引入地形影响因子 $d_i$，当所有夹角都为正即认为此格网包含或处于高程较高的坡顶区域,此时 $d_i=1+L÷10$。当所有夹角都为负即认为此格网包含或处于高程较低的谷底区域，此时 $d_i=1-L÷20$。当所有夹角有正有负时即认为该格网所在区域地形变化不大较为平缓，则 $d_i=1$。对于每个格网设置坡度阈值为 $h_i=d_i\times h$。

3) 遍历并计算格网中的每一个点与相邻所有格网中最低点构成的方向向量和X-Y平面的夹角绝对值 $\left|{\theta }_i\right|$，与(3)中坡度阈值 $h_i$ 作比较，当所有夹角都小于坡度阈值时则认为该点是地面点，否则认为该点为非地面点。

4) 将点云数据再次以更小的尺寸格网化，当格网内点云数量少于5个时，将此格网视为噪声或非地面点；执行(2)、(3)、(4)步，直到滤波结果满意时结束算法。

3. 实验及结果分析

3.1. 实验结果

本文实验所用数据为FARO SINGAPORE PTE. LTD公司提供的机载LiDAR点云学习数据。如图3所示，该实验区域地形复杂，整体植被覆盖度高，区域A (坐标范围X：5150~5250 m，Y：1050~1150 m)地势平缓、植被点较少但地形数据缺失较多；区域B (大致坐标范围X：5250~5300 m，Y：1050~1200 m)地形绵延但植被点数较多；区域C (坐标范围X：5150~5250 m，Y：1150~1200 m)地形陡峭包含有陡坡和断崖且植被覆盖度高，区域内点云数量极多。由原始点云图像(图3)可以清晰看出，在植被覆盖度高的低谷处即图中右部区域，地面点云的数量极少，需要剔除大量的植被点而尽量保留真实地面点。

Figure 3. Original point cloud image

图3. 原始点云图像

本文实验根据测区情况，先后设置格网尺寸为5 m*5 m，2 m*2 m，进行两次迭代之后得到实验结果。第一次结果如图4(a)所示。第二次缩小格网尺寸后重复迭代计算得到滤波后结果如图4(b)所示。

(a) (b)

Figure 4. Overall filtering results of study area with terrain factor. (a) Results after the first filtering; (b) Results after the second filtering

图4. 引入地形因子的研究区整体滤波结果。(a) 第一次滤波后结果；(b) 第二次滤波后结果

根据图5结果，相对于原始点云数据，第一次滤波之后大量的树冠等高层植被点被剔除，而地面点与低矮的植被点被保留，在地形低洼区域，以及点云数量较少的区域剩余植被点较多、较明显。在第二次滤波之后，低矮植被点被大量剔除，地面点保存完好，滤波效果明显。

(a) (b)

Figure 5. Filtering results of two schemes. (a) Filtering effect of scheme I algorithm; (b) Filtering effect of scheme II algorithm

图5. 两种方案滤波结果。(a) 方案一算法滤波效果；(b) 方案二算法滤波效果

3.2. 结果分析

为对滤波效果以及精度进行分析，采用了三种种实验方案分别滤波后，以可视化、生成DEM以及统计试验区内代表区域分类情况等方法进行对比分析，方案一：经典基于坡度的滤波方法，方案二：布料模拟算法(参数设置为布料格网尺寸：0.2，阈值：0.8)，方案三：本文算法。

由以上结果分析，可以明显地看出：1) 方案一的滤波算法保留了较多低矮植被，特别是在点云数量稀少且植被点数量相对较多的区域保留大量的植被点，而方案二的算法不仅在整体区域剔除大量低矮植被点，且在边缘保留了真实地形。2) 利用方案一的滤波算法结果构建的数字高程模型效果较差(图6(a))，特别是在总体点云数量稀少、植被较高且覆盖度大、地面点较少的区域有大量的尖锐小山峰，没有较好地反映地形特征。利用方案二滤波后的点云构建的数字高程模型(图6(b))效果良好，能够极好地反映真实的地形特征，且模型较平滑，大量的植被点被剔除，仅在原来总体点云数量稀少，林木覆盖度极大且地面点较少的区域少量错误剔除了地面点。可见滤波效果好、质量高。

(a) (b)

Figure 6. Digital elevation model filtered by two schemes. (a) DEM model generated by scheme I; (b) DEM model generated by scheme II

图6. 两种方案滤波后的数字高程模型。(a) 方案一生成的DEM模型；(b) 方案二生成DEM模型

因本文实验数据每个点没有属性信息，故采用抽样法选取三个区域内30 m*30 m范围的点云数据，对这些点进行手工分类后根据2003年国际摄影测量和遥感学会(ISPRS)发布的误差评价标准对本次实验滤波结果进行精度评价。该评价方法将滤波后点云分为四类：a) 被正确分类的地面点个数；b) 被错分类为非地面点的地面点个数；c) 被错分类为地面点的非地面点个数；d) 被正确分类的非地面点个数。该方法共计算三类误差，分别为：I类误差、II类误差、III类误差。其中I类误差表示地面点被错误分类为非地面点的比例、II类误差表示非地面点被错误分类为地面点的比例、III类误差是总误差，表示总体点云被错分的比例。三类误差计算方法如下：

$\text{I}类误差=\frac{\text{b}}{\text{a}+\text{b}}$ (3)

$\text{II}类误差=\frac{\text{c}}{\text{c}+\text{d}}$ (4)

$\text{III}类误差=\frac{\text{b}+\text{c}}{\text{a}+\text{b}+\text{c}+\text{d}}$ (5)

为评价滤波的精度，分别选取三种方案滤波后的三个30 m*30 m的区域手动分类并统计正确以及错误分类的情况。其中区域A地势平缓区域且地形数据缺失较多、区域B地形绵延区域且植被覆盖度较高、区域C地形崎岖且植被覆盖程度高，三个区域内点云分类情况如表1所示。

Table 1. Manual classification of point clouds in three regions

表1. 三个区域点云手动分类情况

根据该三个区域内点云分类的统计情况，计算三个区域内两种滤波算法的三类误差(保留两位小数)并统计如表2所示。

Table 2. Statistics of three types of errors in three regions

表2. 三个区域内三类误差统计情况

对本次实验数据中选取的三个不同地形代表性区域进行统计计算三类误差并对比分析方案一和方案二可以得出：1) 在区域A，三种方案均剔除了较多地面点，一类误差都较高。总体上，方案二精度最高、方案三次之、方案一最低；在区域B ,三种方案都保留了较多非地面点，二类误差都较高；在区域C，三种方案的二类误差都较大，总体上方案三精度最高。在该区域方案三总误差较方案二降低了2.68%。2) 经过误差对比、滤波效果分析，方案二即布料模拟算法在地形缺失区域以及陡坡断崖处效果不佳，而本文算法在该类区域效果及精度优于布料模拟算法。3) 方案三滤波后点云数据经过与滤波前原始数据比对分析，被错误分类为非地面的地面点多处于地形凸起及点云缺失区域，如整体数据区域边缘的凸起及水体等点云数据缺失区域的边缘等。而被错误分类的非地面点多处于地形低洼谷底区域。通过三个代表区域的分类情况分析可得出：方案三算法在不同地形区域均有较好的效果，特别在植被覆盖度高且地形起伏大的区域滤波效果有显著提升。

4. 结论

本文针对以窗口化和单一地形坡度为基础的点云滤波算法难以获取定性复杂植被茂密区域DEM问题，提出引入地形判断影响因子并判断格网内地形起伏的方法以达到自适应设置每个格网中最大坡度阈值的目的。本文实验结果表明，引入地形判断因子对传统基于坡度的机载LinDAR滤波算法具有明显的改善作用，在陡坡断崖处的滤波效果优于布料模拟算法，新算法能很好地滤除植被等非地面点，生成的DEM模型能很好地反映真实地形，且滤波后三类误差均较小，尤其在地形起伏大且植被密集区域，本文算法优势明显。但是在地形判断因子计算以得到更为精确的自适应阈值方面还需研究，提高算法的适应性。

基金项目

贵州省科技厅科技支撑计划：高分四号卫星和CYGNSS协同的贵州多云山区土壤湿度智能反演方法研究(黔科合支撑[2022]一般204)；贵州大学研究生创新基地建设项目：测绘科学与技术研究生创新实践基地建设项目(贵大研CXJD[2014]002)。

NOTES

^*通讯作者。