基于自适应坡度的无人机载LiDAR点云滤波方法

摘要

本发明提供的基于自适应坡度的无人机载LiDAR点云滤波方法，包括步骤1：获取测区内的点云数据，对点云数据进行滤波；步骤2：对点云数据建立索引；步骤3：计算并保留格网内的最优平面；步骤4：选择预设窗口内的高程最低点作为地面点的初始种子点；步骤5：计算最优平面和地面点之间的距离，以及该地面点到最优平面中心的坡度，根据坡度阈值获取地面点集；步骤6：将初始地面点集中的点作为基点，采用反距离加权插值法生成数字地面模型；步骤7：通过新生的数字地面模型，获取坡度阈值的更新值，重复迭代步骤4~步骤6，获取最终的地面点集。基于自适应坡度的滤波提高了算法运行效率。

说明书

技术领域

本发明属于遥感与摄影测量中新型遥感数据应用领域，特别是涉及到无人机系统中搭载的激光扫描仪在获取地形、地貌等信息时的快速提取过程中，海量点云数据的自适应滤波处理以及滤波质量的评定方法。

背景技术

自上世纪80年代开始，激光扫描技术(Light Detection And Ranging，LIDAR)，作为一种新型的主动式遥感技术，在三维空间信息获取方面取得了重大突破^[1]。随着商用机载LiDAR系统的普及和推广，机载LiDAR已逐渐应用于三维城市建模^[2]，数字地面模型的获取^[3]，真正射影像的制作^[4,5]等多个领域。近年来，LiDAR系统在硬件和系统集成方面的发展已经十分成熟，而点云数据的后处理及数据挖掘则相对滞后，主要体现在如何高效组织与管理海量不规则离散点云，识别与提取高精度的地物类别。而几乎所有点云数据的后处理算法研究以及应用的关键步骤便是滤波处理，因此，点云滤波算法的研究已成为大多数专家学者关注的热点之一。LiDAR点云数据后处理研究的重要任务是从海量激光点集中提取有效信息，其中，相当一部分是关于数据的滤波处理，即提取出真实地形和地物特征，然后为进一步将机载LiDAR点云分成有意义的点类(道路、裸地、草地、建筑物、树木等)提供良好的条件。目前，大多数滤波算法都是根据判断某点云与周围邻近点的高程差异或几何关系作了假设，然后形成了四种滤波方法:基于数学形态学的滤波方法,基于分割的滤波算法,基于渐进加密的滤波方法和基于地表的滤波方法。

发明内容

本发明要解决的技术问题在于提供一种使算法自动化切尽量保留点云原始信息的角度设计滤波算法，以解决现有滤波算法应用时存在的问题。

本发明采用以下技术方案：

基于自适应坡度的无人机载LiDAR点云滤波方法，包括：

步骤1：通过无人机载LiDAR获取测区内的点云数据，对点云数据进行滤波，剔除高程异常点；

步骤2：使用虚拟规则格网对点云数据建立索引；

步骤3：对步骤2中形成的所有格网，利用每个格网中的所有点，使用均方根误差最小原则分别进行平面拟合，同时保留拟合的平面中保留最小均方误差不大于设定阈值的最优平面；

步骤4：预设第一次选择的窗口尺寸，选择预设尺寸的窗口内的高程最低点作为地面点的初始种子点；

步骤5：计算步骤3中的最优平面和与其最近的地面点之间的距离，以及该地面点到最优平面中心的坡度，若距离小于设定的距离阈值，且坡度小于设定的坡度阈值，则认为该平面为地面平面，将该地面平面延伸至另一格网，同时判断另一格网内的点到该地面平面的距离和平面中心的坡度，若距离小于设定的距离阈值，且坡度小于设定的坡度阈值，则判断为地面点，依次类推到其它格网，获取地面点集；

步骤6：将初始地面点集中的点作为基点，采用反距离加权插值法生成数字地面模型，该数字地面模型的格网尺寸大于步骤2中建立索引时的格网尺寸；

步骤7：计算新生的数字地面模型的均值，并从该均值数字地面模型中提取出局部坡度最大值，作为坡度阈值的更新值，同时，使用新尺寸的窗口内的高程最低点作为地面点的初始种子点，重复迭代步骤4～步骤6，直到达到设定的迭代结束条件，获取最终的地面点集。

所述步骤1中，通过基于局部临近点拟合的噪声检测方法剔除高程异常点。

所述步骤7中，迭代结束的条件为：迭代次数达到设定的阈值，或者生成的数字地面模型中的格网窗口的尺寸包含测区内最大建筑物的尺寸。

所述步骤7中，新尺寸的窗口大小为上次迭代过程中窗口大小的5倍。

本发明的有益效果：

(1)采用虚拟格网索引，实现点云搜索效率的提高，而且最大程度地保留激光点云的原始信息。

(2)基于自适应坡度的滤波是一个迭代更新参数的过程，实现了在坡度参数更新的同时，数字地形模型也进行了迭代更新，且无需人工设置阈值参数，提高算法运行效率。

(3)针对不同类型的城区点云数据，该滤波算法能够很好的顾及大尺度数据整体的地形起伏和局部细节，为激光点云后续处理以及应用提供保障。

附图说明

图1为基于虚拟规则格网的索引建立示意图。

图2为判别地面点与平面之间的坡度和距离的示意图。

图3为本发明的工作流程图。

图4为测区1的滤波前效果图，其中(a)畸变纠正后的影像数据，(b)滤波前点云灰度示意图。

图5为测区2的影像与点云数据，(a)原始影像数据，(b)对应的原始点云灰度示意图。

图6为测区1的滤波前后效果图，(a)滤波后的地面点云，(b)地面点云生成的DEM.

图7为测区2两种滤波算法对比效果，(a)本发明滤波后地面点局部细节，(b)渐进TIN滤波后的地面点局部细节，c)数据块3影像数据局部效果图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明。

本发明通过借鉴基于坡度滤波的基本原理，充分考虑点云数据组织形式以及地表形态，提出了适合城区分类的自适应坡度滤波算法。基于坡度的滤波算法最早由Vosselman提出来，通过比较相邻两点间的高程变化确定最优滤波函数。对于给定的高差阈值Δh_max(d)，两点间距离的递减，直接降低高程值大的激光点属于地面的可能性。分析相邻两点间高程突变的原因可能由于两激光点位分别位于地面和植被、或房屋之上，因此以比较两点间的高差值的大小来判别拒绝或接受所选的点是基于坡度滤波算法的依据。此外，假设某一测区内地形坡度不会大于30％，考虑到地面点观测值的误差情况，设置5％的置信区间，则滤波函数可定义为：

其中N为激光点集，DEM为地面点，d为两点间距离，若在点p_i的邻域内没有临近点p_j满足式(1)，则p_i就被划分为地面点。式(2)中，σ为地面点的标准差。

本发明提供一种基于自适应坡度的无人机载LiDAR点云滤波方法，并对该滤波方法的滤波质量进行评定。该方法包括两部分内容，一是激光点云数据的预处理，二是滤波过程中的坡度参数的自动更新。本发明能够针对不同类型的城区点云数据进行滤波，并且可以很好的顾及大尺度数据整体的地形起伏和局部细节，为激光点云后续处理及应用提供保障。

机载LiDAR系统获取三维空间中的点云呈随机分布，分别记录了来自裸露地面和地物表面的反射脉冲信号。尤其是城区点云数据，包含不同类型的地物信息，分别是真实地形点、人工地物(建筑物、桥梁等)、自然地物(草地、灌木、高植被等)。一般情况下，影响城区滤波质量的因素可归纳为两个方面：①城市建筑物侧面点云被错误分类为地面点；②将树木靠近地表的多次回波信号误分为地面点。

本发明的方法包括以下步骤：

步骤1：通过无人机载LiDAR获取需要的点云数据；该点云数据为获取的测区内的地形数据。同时，点云数据获取以后，需要对点云数据的噪声进行剔除，即对点云数据进行滤波，剔除高程异常点；本发明中，利用K近邻查询的方式搜索某点周围的激光点，剔除一些高程值异常的激光点。

步骤2：使用虚拟规则格网对点云数据建立索引。由于通过激光获取的激光点集的各点之间是离散分布，没有一定的拓扑关系，致使对于激光点云数据的处理方式和算法都与传统的基于栅格的影像数据有很大不同。目前使用较为广泛的为基于二维分块的格网索引和基于几何的树形索引。

针对具体的数据处理需求采用不同的空间索引，本发明为了后续滤波算法设计，采用虚拟格网的点云组织方式。具体建立过程为：从原始点云入手，通过虚拟规则格网来管理三维点云。结合激光点的平面坐标，与网格建立索引机制，既可以判断某一激光点落在哪个网格，也可以快速查询出一个网格包含哪些脚点。图1则表示了基于虚拟规则格网的索引建立示意图。其中图1(a)、(b)中的圆点代表LiDAR脚点(x，y，z)，长方体块代表虚拟格网。

由图1(a)到图2(b)实现虚拟格网建立的规则是按照下式进行的，其中n_grid代表格网大小，一般大于或者等于点云平均点间距。该类索引的优点便是不仅实现了搜索效率的提高，而且能够保留激光点云的原始信息。

X＝INT(x-x_min)/n_grid

Y＝INT(y-y_min)/n_grid

本步中，也可采用其余的索引方法，例如不规则三角网索引，八叉树索引机制等现有的点云索引方法，此时，需要修改后面的步骤。

步骤3：选择预设尺寸的窗口内的高程最低点作为地面点的初始种子点。

每次迭代循环的最后，都要改变预设窗口大小w_size，即如果选择新的迭代窗口为5倍的上次计算时采用的窗口，则新的迭代窗口＝5*上次迭代窗口w_size。而窗口的移动是根据建立的规则格网进行的，即从数据的整体格网的左上角(x₀，y₀)，依w_size*w_size形成的窗口，按照行—列的方式(从左至右，从上至下)逐一激光点的进行遍历整个激光点云数据。

步骤4：对步骤2中形成的所有格网，利用每个格网中的所有点，使用均方根误差(RMSE)最小原则分别进行平面拟合，同时剔除最小均方根误差大于设定阈值Tp1的平面，保留拟合的平面中保留最小均方误差小于设定阈值的最优平面。其中，拟合平面方程如下式。

ax+by+cz+d＝0

步骤5：计算步骤4中的最优平面和与其最近的地面点之间的距离，以及该地面点到平面中心的坡度，若距离小于设定的距离阈值，且坡度小于设定的坡度阈值，则认为该平面为地面平面，将该地面平面延伸至另一格网，同时判断另一格网内的点到该地面平面的距离和平面中心的坡度，若距离小于设定的距离阈值，且坡度小于设定的坡度阈值，则判断为地面点，依次类推到其它格网，获取地面点集；

如图2图3所示，对步骤(4)中保留的平面，搜索与其最近的地面点，计算地面点到平面点的距离，以及地面点到平面中心的坡度。若距离小于阈值Td1，坡度小于阈值Ts1，则认为该平面作为地面平面，图2中，Point1和point2分别代表在判定条件内的平面点，通过与种子地面点的比较(距离阈值和坡度阈值)来判断是否能成为初始地面点集内的元素。

dist＝distance是平面点point1或者point2到种子地面点的距离。

对位于地面上的平面，计算该平面格网窗口内的点到该平面的距离，若小于距离阈值Td2，则认为该点为平面点；并对当前格网进行区域增长，判断其邻域格网内的点到该平面的距离和平面中心的坡度，若距离小于阈值Td2，坡度小于阈值Ts2，则判为地面点，其标准依然采用图2进行判断，得到初始地面点集。

步骤6：将初始地面点集中的点作为基点，采用反距离加权插值法生成数字地面模型，该数字地面模型的格网尺寸大于步骤2中建立索引时的格网尺寸，例如可设置为原来的格网尺寸的五倍。

步骤7：计算新生的数字地面模型的均值，并从该均值数字地面模型中提取出局部坡度最大值，作为坡度阈值的更新值，重复步骤3～步骤6，直到生成的数字地面模型中的格网尺寸包含测区内最大建筑物的尺寸，获取最终的比较可靠的地面点集。

本发明采用基于地表坡度的滤波思想，从改进算法自动化和尽量保留点云原始信息的角度设计滤波算法，使用时，①鉴于本文实验是大多以城区为基准，其地形相对平坦，为了提高处理效率，采用小尺度的规则格网压缩和组织点云数据。②应对大尺度数据整体的地形起伏同时顾及局部真实的地形细节，采用迭代更新参数的方法，在坡度参数更新的同时，数字地形模型也进行了迭代更新。

本发明还提供对最后生成的滤波结果(地面点集和非地面点集)进行定性和定量评定的方法。

即本发明基于自适应坡度的无人机载LiDAR点云滤波方法系统中，至少应包括以下模块：

点云数据预处理模块：用于海量无人机LiDAR点云的噪声剔除和建立点云索引机制，仅实现了搜索效率的提高，而且能够保留激光点云的原始信息。

自适应坡度滤波模块：考虑到应对大尺度数据整体的地形起伏同时顾及局部真实的地形细节，本发明提出的基于自适应坡度的滤波算法，不仅能够自动更新坡度阈值，而且可以同时更新数字地面模型(DTM)，此外，在应对大尺度数据整体的地形起伏和局部真实的地形细节方面具有较好的优势。

无人机载点云滤波性能评定模块：为了验证本发明在实现基于无人机载海量点云数据的滤波处理方案的可适用性，分别从定性(DEM点集平滑度)和定量(I类误差、II类误差、总体误差和Kappa系数)两个方面进行性能评定，并输出最终评定结果。

与传统无人机载点云滤波算法相比较，本发明的优势体现在：

(1)采用虚拟格网索引，实现点云搜索效率的提高，而且最大程度地保留激光点云的原始信息。

(2)基于自适应坡度的滤波是一个迭代更新参数的过程，实现了在坡度参数更新的同时，数字地形模型也进行了迭代更新，且无需人工设置阈值参数，提高算法运行效率。

(3)针对不同类型的城区点云数据，该滤波算法能够很好的顾及大尺度数据整体的地形起伏和局部细节，为激光点云后续处理以及应用提供保障。

算法验证与分析

测区1数据：该数据来自我国的沙市地区，该测区共包含1,954,659个激光点，X轴向的跨度是904.55米，Y轴向的跨度是789.47米，由于原始数据存在噪声，所以最小高程值是-160.60米，最大高程值是192.33，经过粗差剔除后的最大高程差是38.48米。图4(a)表示的是经过纠正的该测区影像数据，图3(b)是对应的点云数据滤波前构成的数字表面模型(DSM).

测区2数据：该实验数据块选取的是我国湖北老河口地区的一部分，老河口测区的点云数据是采用LSM-Q160激光仪于2013年4月采集的，整个测区共有7条点云条带，共20,141,032个激光点，以及利用DMC-II-140航摄仪获取的12条航带456张影像数据。本测试数据—数据块3的点云个数是443,415，X轴向跨度是545.35米，Y轴向的跨度是505.71米，最小高程值是63.73米，高程最大值是125.64米，图5(a)是数据块3的影像数据，图5(b)则是对应的原始点云DSM。

算法评价—定性：

测区1：如图6所示，从滤波后生成DEM可知，建筑物和树木均被有效地滤除，并且整体效果较为平滑，没有异常的凸起或者凹陷的现象(除了沿河流的两侧的变化)，与实际的地形特征比较吻合。

算法评价—定量：

①混淆矩阵

总分类精度(Overall Accuracy):

变化一致性(Change agreement):

Kappa系数：

其中，n表示原始点云总数，g是滤波后的地面点，o是滤波后的地物点，T_g表示实际地面点(人工或经验得出的样本数据)，T_o则是实际地物点个数。而a表示地面点中被正确判断为地面点的数量，b是地面点中被误判为地物点的数量；c则表示地物点中被误判为地面点的数量，d为地物点中正确判断为地物点的个数。

对于测区2除了采用本发明滤波算法外，还采用传统的渐进TIN滤波算法，图7(a)和(b)表示的两种方法滤波后地面点的细节图，为了更好地进行定性对比，图7(c)则给出了该测区某一部分局部方法的影像数据。表1则是给出了两种方法滤波的混淆矩阵以及相关精度指标，其中样本点中地面点个数是82750，地物点个数是360665。

表1 测区2两种滤波算法的误差统计

由表1可知，两种方法的总体精度都能达到75％以上，即可视为都能得到不错的滤波结果。然而本文的滤波算法是一个迭代更新的过程，不仅会在迭代过程中更新坡度值，还可以更新DTM，能够很好的顾及大尺度数据整体的地形起伏和局部细节，这主要体现在II类误差较之渐近TIN要大幅降低，减少了非地面点误分为地面点的几率，得到比较好的滤波效果，为后续点云数据处理及应用(点云分类、地物识别)提供条件。

以上所述的仅是本发明的优选实施方式，应当指出，对于本领域的技术人员来说，在不脱离本发明整体构思前提下，还可以作出若干改变和改进，这些也应该视为本发明的保护范围。