CN103679655A - 一种基于坡度与区域生长的LiDAR点云滤波方法 - Google Patents

Abstract

本发明公开了一种基于坡度与区域生长的LiDAR点云滤波方法，主要包括以下步骤：步骤1：获取LiDAR原始点云数据；步骤2：粗差剔除；步骤3：规则网格化。步骤4：地面种子点的选取；步骤5：区域生长。步骤6：将满足条件的新地面点所在的格网为中心，格网中高程最低的点作为新的种子点按照步骤4进行生长，直到没有点满足条件，停止生长。步骤7：根据步骤4和步骤5，依次对种子点集中种子点进行区域生长，直至结束。本发明有效的限制了高差阈值的大小，而且有效的降低了整个方法的计算量，能够更快的将复杂城市的LiDAR点云数据中地面点与非地面点进行分离，同时方法更加简单，方法获得结果的误差相对较小。

Description

一种基于坡度与区域生长的LiDAR点云滤波方法

技术领域

本发明属于对地面观测领域，涉及一种基于坡度与区域生长的LiDAR点云滤波方法。

背景技术

自从1998年美国前副总统戈尔提出了“数字地球”的概念，在全世界引起了强烈反响数字地球既是地球科学技术与空间技术，信息科学技术等现代科学技术交融的前沿领域，更是科学技术与社会发展以及社会科学紧密结合的结晶。作为数字地球的一部分，数字城市的建设和应用受到越来越广泛的关注，城市基础设施的三维信息在城市规划、市政管理、公共安全等方面都显现出重要的应用价值。建筑物是数字城市中最基本的信息之一，LIDAR系统具有直接描述测区表面三维信息的能力，可以快速、高效地将城市建筑物的特征信息采集到计算机当中，为快速进行建筑物信息提取带来了新的方法。

机载激光雷达(LiDAR，Light Detection and Ranging)是近些年才发展起来的新的探测技术，是集激光测距技术、惯性导航系统、动态差分GPS定位技术于一身的对地观测系统，代表了对地观测领域的一个新的发展方向，通过位置、角度、距离等观测数据直接获取目标对象表面点的三维空间坐标信息，实现地表提取和三维场景重建，被广泛的用于地面数据探测和模型的恢复、重建等应用中并显示出了巨大的前景，并逐渐成为三维城市数据模型获取的一种重要方式。机载激光雷达在获取目标地物三维空间地理坐标的同时，能够获得地表物体的垂直结构形态，同时配合上系统配备的高分辨率数码相机获取的地面相片更增强了对地面目标地物的分析和判别能力。与摄影测量技术相比，激光脉冲为主动遥感方式，对天气的依赖性较小，不受阴影和太阳角度的限制，避免了投影带来的信息损失，高程精度优势明显，提高了自动化程度，大大缩短了成图周期，可快速完成数字高程模型及数字正射影像图的大规模生产任务。

在数字城市建设中，往往需要对建筑物进行提取和重建，而建筑物特征提取的关键是数据滤波，由于地理环境的复杂性，数据滤波一直是研究的难点，滤波的最终结果是分离出地面点与非地面点。

机载激光雷达数据的滤波原理可以分为两类：第一类是根据激光脚点数据的回波强度信息进行的滤波。基于回波强度信息的滤波原理是参考了航空影像的分类原理，即不同的地物类型对激光的反射系数是有差异的，其决定了激光回波能量的大小，进而决定了该类地物表面的激光脚点的回波强度信息，将其可以转化为灰度图像，可呈现出不同的灰度值。而对于一片测区，可以将激光点云数据以其强度信息为依据转化为灰度图像，然后分析其灰度值的分布，确定出地面点的灰度范围，从而将地面点从原始点云数据中过滤出来。尽管该方法一定程度上发挥了LIDAR点云自身优势，但仅仅使用回波信息进行点云分类难以适合大部分地形。第二类是基于高程突变。基于高程突变的滤波原理是一种现在应用最为广泛的滤波原理。一般地，地形是连续的，地势起伏也是缓慢而平滑的，距离越近的两个点的高差应该越小，如果地表局部区域出现了高程突变，即邻近的两个点高程相差较大时，一般都被认为是因为此处地物类型发生了变化从而导致了局部地形的高程不连续，比如此处有房屋、树木等，该原理认为高程较高的点是发生变化的地物表面的点，高程较低的点被认为是地形表面的点。

现有的LiDAR数据滤波方法有很多，大致分为3类：形态学方法、基于内插的方法和基于曲面拟合的方法。每种滤波算法都有其适用的地形及其不足之处。2004年，GeorgeSithole等选择了8种典型的算法进行比较，分析了这些算法的优缺点。这8种算法包括：Elmqvist的主动轮廓线法、Sohn的基于区域法、AxelssOn的渐进TIN加密法、Pfeifer的迭代最小二乘滤波算法、Brovelli的样条插值法、Roggero的改进的基于坡度法、Wack的改进区域最小滤波法、Sithole的坡度滤波法。最后结论表明由于地形的复杂性，没有一种算法能够适合任何地形条件，需要对算法作更多的改进，例如使用更多的上下文信息，结合场景分解。

基于坡度的滤波算法最早由Vosselman提出，随后Sithole对其进行改进，提出了一个基于核函数的滤波方法，使得该方法适用于陡峭的地形；Roggero也对Vosselman的方法进行了改良，利用局部线性回归估算地面坡度提取初始地面。基于坡度的算法是经典滤波算法之一，该算法原理简单，实现方便，且充分考虑到地形坡度问题，适用区域较为广泛。但该算法存在一些不足之处：1)对于面积较大的地物上的点，例如大型建筑物上的点，其与最邻近的地面点之间的距离越大，高差阈值就会随距离增大而增大，即放宽了拒绝其为地物点的条件，会使得一些面积较大而又相对较低的地物上的点无法滤除。2)基于坡度的滤波算法是通过寻找点的邻接关系，逐点计算两点的高差来判断地面点，导致算法的计算量较大、速度较慢。

发明内容

发明目的：本发明的目的在于针对现有技术的不足，提供了一种有效降低了整个方法的复杂程度，计算量小，计算速度快的基于坡度与区域生长的LiDAR点云滤波方法。

技术方案：本发明提供一种基于坡度与区域生长的LiDAR点云滤波方法，包括以下步骤：

步骤1：获取LiDAR原始点云数据；

步骤2：对步骤1中获取的LiDAR原始点云进行粗差剔除；由于LiDAR原始点云数据中包含一些粗差低点，会影响到地面种子点的选取，因此为了使种子点尽量可靠，首先要粗差剔除。本发明用一个中心点的高程值与给定距离范围内每一个点的高程值比较，如果中心点明显低于其他点，这个点将作为粗差点被剔除。经过粗差剔除后，消除了一些高程异常点的影响。

步骤3：用规则格网来组织经过步骤2后的LiDAR点云数据，将所有LiDAR点云的坐标与规则网格建立起点云索引机制。规则格网的建立有助于快速搜索种子点的八邻域。其中，判断点云所在规则网格的公式如下：

X＝INT(x-x_min)／n   (1)

Y＝INT(y-y_min)／n   (2)

式中，(X，Y)为网格编号，(x，y)为点云的平面坐标，(x_min，y_min)为点云数据集中最小坐标，n为规则网格的大小。格网的大小视点密度而定，一般能保证每个规则网格内有1、2点即可。

步骤4：选取种子点集，地面种子点的选取是滤波的关键，直接影响到滤波的精度。本发明以数据区中最大建筑物的尺寸为新的规则格网大小对数据进行划分，选取每个新的规则格网中的最低点为地面点，将每个规则网格中的地面点组成地面种子点集。

步骤5：以地面种子点集中的第i个(i从1开始)地面种子点所在的步骤3中生成的规则网格为中心网格，进行8邻域搜索，将中心网格相邻的8个规则格网及中心格网中除了第i个地面种子点之外所有的点定义为待定点，分别判定每个待定点与第i个地面种子点的高差是否满足基于坡度的滤波算法的高差函数阈值，如果满足，将该待定点标记为地面点，否则，该该待定点为地物点；基于坡度的滤波算法的高差函数是一种基于两点间高差值和距离值的滤波函数。这个函数不仅表示了地形，还包括激光点的信息，可以由预估的测区地形最大坡度百分数及激光点的标准偏差确定。其中，所述基于坡度的滤波算法的高差函数为：

${\mathrm{\Δh}}_{\max }\left(d\right)=\mathrm{Sd}+1.65\sqrt{{2\mathrm{\σ}}_z},d\≤d_{\max }---\left(3\right)$

根据上述函数可以计算出高差函数阈值，其中S表示最大地形坡度百分数，S为经验值，对于较为平坦的区域，坡度值设为0.3，地形起伏较大区域的坡度值取0.5，d表示待定点与种子点之间的水平距离，d_max为待定点与地面种子点之间的最大水平距离，σ_z为原始Lidar数据的高程误差。则地面点集合就可以定义为：

$P_T=\{P_i\∈P|\∀P_k\∈P:z_i-z_k\≤{\mathrm{\Δh}}_{\max }\left(d_{\mathrm{ik}}\right)\}---\left(4\right)$

$d_{\mathrm{ik}}=\sqrt{{(x_i-x_k)}^2+{(y_i-y_k)}^2}---\left(5\right)$

P_k(x_k，y_kz_k)表示地面点，z_k为该点的高程值，P_i(x_i，y_iz_i)表示待判断点，z_i为待判断点的高程值。其中，最大距离d_max对高差函数有影响，因此需要定义邻域对激光点测试的贡献，d_max优选为10m。

步骤6：将步骤5中获得的地面点所在的步骤3中生成的规则网格为中心，规则网格中高程最低的点作为新的地面种子点按照步骤5进行生长，直到步骤5中没有新的地面点产生则停止生长。

步骤7：以种子点集中的第i+1个地面种子点所在的步骤3中生成的规则网格为中心，重复步骤4和5，遍历种子点集中的所有种子点。

工作原理：本发明结合区域生长的原理，对基于坡度的滤波算法做出了一些改进：使用规则格网代替Delaunay三角网建立点云间索引结构，不进行插值，在不损失原始精度的基础上，能够有效减小算法的复杂度；根据数据区的最大建筑物尺寸寻找出种子点集，以种子点集中的点为中心进行8邻域搜索，判断待定点与已知种子点的高差是否满足高差函数阈值，并将满足条件的地面点作为新的种子点进行生长。改进后的算法通过不断更新种子点，使得种子点与待定点之间的距离保持在设定的格网大小之内，限制了高差阈值的大小。

有益效果：与现有技术相比，本发明有效的限制了高差阈值的大小，而且有效的降低了整个方法的计算量，能够更快的将复杂城市的LiDAR点云数据中地面点与非地面点进行分离，同时方法更加简单，方法获得结果的误差相对较小。

附图说明

图1：本发明的流程图；

图2：ISPRS提供的样本数据。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

如图1所示，本发明的基于坡度与区域生长的LiDAR点云滤波方法为：

步骤1：获取LiDAR原始点云数据；为了检验算法的稳定性和可行性，本发明采用国际摄影测量与遥感学会(ISPRS)提供的专门用于滤波算法测试的样本数据进行试验。如图2所示，选取其中的Site1和Site2数据中的6个样本数据进行滤波实验。该数据点间距为1.0～1.5m，每平方米约采样0.67个点，点云的高程误差基本在0.2～0.3m之间。该实验区属于典型的城市区域，地面变化较为平缓，地物复杂，含有大型建筑物、不规则建筑物、道路桥梁和小型隧道等。实验过程中取规则格网大小为1.5×1.5m²，_σz为0.3，对于较为平坦的区域，坡度值设为0.3，地形起伏较大区域的坡度值取0.5，最大建筑物的尺寸从点云数据中直接量取。

步骤2：粗差剔除，实验数据的LiDAR点间距为1m，以每个点为中心点，搜索2m半径内的领域点，并比较中心点与领域点的高程值，如果中心点明显低于其他点，这个点将作为粗差点被剔除。

步骤3：规则网格化。对粗差剔除后的点云数据根据公式

X＝INT(x-x_min)／n   (1)

Y＝INT(y-y_min)／n   (2)

和每个点云的坐标与网格建立起索引机制。其中n取1.5×1.5m²。

步骤4：地面种子点的选取。本发明通过RealWorks软件量取数据区中最大建筑物的尺寸sm²，并利用步骤3中的公式以最大建筑物的尺寸sm²为网格大小对数据重新进行网格划分，选取每个重新划分的格网中的最低点判为地面点，并以这些点作为地面种子点集。

步骤5：区域生长。从种子点集中的第一个种子点开始，以该地面种子点在步骤3中的生成的规则网格为中心，根据步骤3中的规则网格进行8邻域搜索，将中心网格相邻的8个规则格网及中心格网中除了第i个地面种子点之外所有的点定义为待定点，分别判定每个待定点与第i个地面种子点的高差是否满足高差函数阈值，如果满足，将点标记为地面点，否则，该点为地物点，其中，基于坡度的滤波算法的高差函数为：

${\mathrm{\Δh}}_{\max }\left(d\right)=\mathrm{Sd}+1.65\sqrt{{2\mathrm{\σ}}_z},d\≤d_{\max }---\left(3\right)$

根据上述函数可以计算出高差函数阈值，其中S表示最大地形坡度百分数，S为经验值，对于较为平坦的区域，坡度值设为0.3，地形起伏较大区域的坡度值取0.5，d表示待定点与种子点之间的水平距离，d_max为待定点与地面种子点之间的最大水平距离，σ_z为原始Lidar数据的高程误差。则地面点集合就可以定义为：

$P_T=\{P_i\∈P|\∀P_k\∈P:z_i-z_k\≤{\mathrm{\Δh}}_{\max }\left(d_{\mathrm{ik}}\right)\}---\left(4\right)$

$d_{\mathrm{ik}}=\sqrt{{(x_i-x_k)}^2+{(y_i-y_k)}^2}---\left(5\right)$

P_k(x_k，y_k，z_k)表示地面点，z_k为该点的高程值，P_i(x_i，y_i，z_i)表示待判断点，z_i为待判断点的高程值。其中，最大距离d_max对高差函数有影响，因此需要定义邻域对激光点测试的贡献，d_max设为10m。

步骤6：将满足条件的新地面点所在步骤3中生成的规则网格为中心，规则网格中高程最低的点作为新的种子点按照步骤5进行生长，直到没有点满足条件，停止生长。则种子点集中的第一个种子点区域生长完成。

步骤7：根据步骤4和步骤5，依次对种子点集中种子点进行区域生长，直至结束。

本发明采用ISPRS建立的评价体系对滤波结果进行定量分析。一般地，滤波误差可分为三种：I类误差(Type I Error)、II类误差(Type II Error)和总误差(Total Error)。I类误差是将地面点分类为非地面点的误差，II类误差则是将非地面点分类为地面点的误差，对I、II类误差加权求和可得到总误差。I、II类误差反映了算法的适应性，总误差则反映了算法的可行性。ISPRS报告认为应该侧重减I类误差，因为在质量控制中，需要人工修改错误，而修改II类误差比修改I类误差容易得多。3类误差的计算公式如下：

TypeI＝b/(a+b)×100％

TypeII＝c／(c+d)×100％      (6)

Total＝(b+c)/e×100％

其中，e为总点数，a为正确判断的地面点，b为地面点错判为地物点，c为地物点错判为地面点，d为正确判断的地物点。经验表明，在人工干预时，与I类误差相比，II类误差带来的分类错误更容易被修复。在适当增加II类误差的前提下，滤波重点应放在减少I类误差上。表1为本发明算法与Sithole的坡度滤波算法、Roggero的改进基于坡度滤波算法的I类误差统计表，表2为本发明算法与Sithole的坡度滤波算法、Roggero的改进基于坡度滤波算法的II类误差统计表，表1为本发明算法与Sithole的坡度滤波算法、Roggero的改进基于坡度滤波算法的总误差误差统计表。

表1

表2

表3

从表1～表3可知，本发明算法的I类误差明显小于Roggero提出的改进的基于坡度算法和Sithole提出的坡度滤波算法，II类误差比较大，但控制在一定的范围之内，总误差基本是三种算法中最小，说明本发明的滤波方法误差比较小。

Claims (3)

1.一种基于坡度与区域生长的LiDAR点云滤波方法，其特征在于：包括以下步骤：

步骤1：获取LiDAR原始点云数据；

步骤2：对步骤1中获取的LiDAR原始点云进行粗差剔除；

步骤3：用规则格网来组织经过步骤2后的LiDAR点云数据，将所有LiDAR点云的坐标与规则网格建立起点云索引机制，其中，判断点云所在规则网格的公式如下：

X＝INT(x-x_min)／n

Y＝INT(y-y_min)／n

式中，(X，Y)为网格编号，(x，y)为点云的平面坐标，(x_min，y_min)为点云数据集中最小坐标，n为规则网格的大小。

步骤4：选取种子点集，以数据区中最大建筑物的尺寸为新的规则格网大小对数据重新进行划分，选取每个新的规则格网中的最低点为地面点，将每个新的规则网格中的地面点组成地面种子点集。

步骤5：以地面种子点集中的第i个地面种子点所在的步骤3中生成的规则网格为中心网格，进行8邻域搜索，将中心网格相邻的8个规则格网及中心格网中除了第i个地面种子点之外所有的点定义为待定点，分别判定每个待定点与第i个地面种子点的高差是否满足基于坡度的滤波算法的高差函数阈值，如果满足，将该待定点标记为地面点，否则，该待定点为地物点；其中，所述基于坡度的滤波算法的高差函数为：

$\mathrm{\Δ}{h}_{\max }\left(d\right)=\mathrm{Sd}+1.65\sqrt{{2\mathrm{\σ}}_z},d\≤d_{\max }$

其中S表示最大地形坡度百分数，d表示待定点与种子点之间的水平距离，d_max为待定点与地面种子点之间的最大水平距离，σ_z为原始Lidar数据的高程误差。

步骤6：将步骤5中获得的地面点所在的步骤3中生成的规则网格为中心，规则网格中高程最低的点作为新的地面种子点按照步骤5进行生长，直到步骤5中没有新的地面点产生则停止生长。

步骤7：以种子点集中的第i+1个地面种子点所在的步骤3中生成的规则网格为中心，重复步骤4和5，遍历种子点集中的所有种子点。

2.根据权利要求1所述的基于坡度与区域生长的LiDAR点云滤波方法，其特征在于：所述步骤2中LiDAR原始点粗差剔除的方法为：设定一个中心点，将中心点的高程值与给定距离范围内每一个点的高程值比较，如果中心点明显低于其他点，将中心点将为粗差点剔除。

3.根据权利要求1所述的基于坡度与区域生长的LiDAR点云滤波方法，其特征在于：所述步骤5中的d_max为10m。

Images