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

赵明君, 刘 超, 高 翔, 李慧慧

(1.安徽理工大学 空间信息与测绘工程学院,安徽 淮南 232001; 2.安徽理工大学 矿山采动灾害空天地协同监测与预警安徽普通高校重点实验室,安徽 淮南 232001)

0 引 言

三维激光扫描技术是一种新型的高时空分辨率、高精度的地理信息数据观测手段[1],所获取的点云数据包括测量区域内所有的地表及地物的三维信息,因此,点云数据后处理算法研究中的关键问题,就是如何从大量的散乱点云数据中快速准确地将地表点云数据与其他地物点云数据区分[2]。点云滤波就是指区分点云中地面点和非地面点的过程[3]。

在点云数据滤波方面,常用的滤波算法有4种,分别为基于数学形态学、基于不规则三角网(triangulated irregular network,TIN)、基于曲面拟合、基于坡度[4-7]。文献[8]提出了基于数学形态学滤波算法,该算法的关键步骤是如何选择合适的滤波窗口;文献[9]提出了一种渐进式形态学滤波方法,虽然取得了较好的点云滤波效果,但是在采用大窗口对点云数据滤波时,仍容易将地形起伏剧烈处的点误判为非地面点,降低了分类精度;文献[10]给出了一种基于TIN的点云滤波算法,该算法是基于二维层次的搜索方法,其计算量和算法的复杂度相对较大;基于曲面拟合算法的核心是通过足够的点拟合出较准确的地形曲面,再通过高差值判断缓冲区内的点是否为地面点[11];文献[12]提出了一种基于坡度的滤波算法,并将其应用到点云滤波过程中,该算法的基本思路是利用地面点与其邻近的非地面点之间的坡度值大于设置阈值来区分非地面点与地面点。

上述4种点云滤波方法中,基于坡度滤波算法的滤波原理较简单,可以适应多种类型的点云数据,但在实际复杂多样的地形中,该算法对整个点云区域按单一的坡度阈值进行滤波是不完全合理的。文献[13]针对该缺点进行了改进,提出了一个基于核函数的滤波方法,增强了坡度滤波算法在陡峭地形区域中的滤波精度;文献[14]中利用数字地面模型(digital terrain model,DTM)进行分层迭代更新得出坡度阈值,得到了坡度阈值随着地形坡度变化而变化的滤波效果;文献[15]提出在陡坎、斜坡等地形变化较大处通过设置坡度增量来区分非地面点与地面点,以减小坡度滤波时的分类误差;文献[16]通过增加不断变化的窗口,利用DTM不断更新来实现坡度阈值的适应性设置,可以顾及点云区域的整体地形起伏和局部细节。这些基于传统滤波算法的改进方法都是基于整体地形变化而更新坡度阈值设置的,没有考虑到局部细节的地形结构;同时多次的迭代计算也增加了算法的复杂度,降低了算法的效率。因此,上述改进方法并未完全实现坡度阈值的自适应设置。

坡度滤波算法需要解决以下2个难点:① 如何根据地表地形变化得出最优坡度滤波阈值;② 在地形变化剧烈的边缘带过度滤波造成的地表点云数据损失问题。为此,本文提出了一种基于多尺度网格的点云自适应坡度滤波算法(简称“本文算法”),通过统计分析局部地形坡度值得出最优坡度阈值,实现了真实地形的坡度阈值自适应设置,从而分离点云数据中的地面点与非地面点。

1 点云自适应坡度滤波算法

本文算法主要内容包括:① 点云数据组织,即进行孤立点去除与网格构建;② 坡度阈值计算与自适应设置,即进行网格属性划分与种子点选取、点云数据的自适应坡度滤波、多尺度网格的点云重复滤波;③ 地面点云数据去噪与数据输出。该算法对复杂地形自适应设置坡度阈值,并通过不同尺度网格重复滤波,尽可能剔除非地面点,保留真实地形。点云自适应坡度滤波算法流程图如图1所示。

图1 点云自适应坡度滤波算法流程图

1.1 多尺度网格构建与孤立点去除

(1) 多尺度网格构建。原始点云数据基本上是散乱无序的,点间无拓扑关系。为了克服算法运算效率低和信息损失的问题,在光探测和测距(Light Detection and Ranging,LiDAR)(简称“激光雷达”)数据处理中引入虚拟网格概念。多尺度网格的示意图如图2所示。

图2a中,圆点表示点云数据,长方体表示相应尺度的网格;图2b所示为在虚拟网格内将三维坐标点投影到X-Y平面,方格的颜色由深到浅表示虚拟网格尺度由小到大;图2c所示为基于点的高程对网格分层时点云的存在形式;图2d所示为包含地面点与植被点时,点云数据在网格中的空间位置。

由图2b可知,可以将点云数据按照其X-Y坐标与网格建立起索引机制,每个点所在的网格或每个网格包含的点均可以快速查询。

图2 多尺度网格示意图

点与网格间的索引关系计算公式为:

(1)

其中:(X,Y)为网格号;(x,y)为点云的平面坐标;(xmin,ymin)为整个数据集的最小平面坐标;m为网格单元的尺度;INT表示对计算结果向下取整。

(2) 孤立点去除。从图2c可以看出,通过在网格内进行分层高程统计,可以得出在不同高程区间内的点云个数,此时若某一层中点云个数不足平坦区域内点云密度的10%,则认为该层的点云为孤立点,并进行删除。此步骤主要是减少孤立点,为后续运算得到较纯净的点云数据。

1.2 坡度阈值自适应设置与点云滤波

1.2.1 网格属性划分和种子点选取

(1) 网格属性划分。从图2d可以看出,该网格包含了地面点云数据和植被点云数据。当网格平面尺度大小一定时,给定网格高度,可以定义空间占比F为:

(2)

其中:hi为网格内第i点高程;hmin′为网格内次低点高程;H为网格高度;n为网格内点云个数。F反映了点云数据在空间存在形式的差异,根据F值可以定义网格的不同语义属性,同时根据网格属性,可以在后续计算中给出不同网格的地形因子,为坡度阈值的计算分配参数。选取1个实验区,网格平面尺度为5 m,H=25 m,计算不同属性网格的F值,结果见表1所列。

表1 实验区不同语义属性的空间占比

(2) 初始地面种子点选取。初始地面种子点的准确选取是基于坡度的滤波算法有效运行的关键因素。原始点云数据经过孤立点去除和点云格网构建后,通常把每个网格内的高程最低点作为地面种子点,但该方法选取的种子点不是最优的地面种子点。为提高所选取的地面种子点的可靠性,对种子点选取方法做如下改进:① 初始地面种子点必须在最大尺度的网格内选取,为保证每个网格内包含地面种子点,网格的最大尺度应尽量大于点云区域内建筑物或树冠的宽度;② 为避免所选取的最低点是局部微地形中位于地面以下的噪声点,对每个网格先按照高程由高到低排序,然后选取高程值从最小值开始位置排列在点数5%的激光点作为地面种子点;③ 如果少量网格包含的点全是非地面点,那么无法选取地面点,此时计算该网格内待拟定的种子点与周围8个邻格内种子点的高差值Δhi(i=1,2,3,…,8),若Δhi均大于网格尺度值的2倍,则认为该种子点为非地面点,予以去除,并删除该网格内的点云数据。上述的3个约束条件的目的是最大概率保证选取的地面种子点为真实地面点。

1.2.2 坡度阈值求解与自适应设置

根据F值划分出网格属性并选出单一网格内的地面种子点后,对于不同的网格属性需要给出相应的地形计算因子σ,σ为网格内点云与种子点间高差的倍数。网格的语义属性仅能粗略地描述地形与地物的空间关系。对于地形复杂的区域,网格内局部地形的变化仍然存在。为更加真实反演网格内的地形变化,在网格内计算出激光点与种子点的平面距离,以种子点为圆心,按平面距离从小到大、等点数间隔逐步分割的原则,将每个网格分割成若干个辐射圈。在每个辐射圈选取高差值为hpt时所对应的激光点为参考点pt,多个参考点pt即为参考点集Pt。

hpt=max{ht≤σumin,umin=minHt}，hpt∈{Ht}

(3)

其中:Ht为辐射圈内点云数据中第t点与该辐射圈内高程最低激光点间的高差;σ为该网格对应的地形计算因子。计算参考点集与种子点间的坡度值si,则其中的最大值为坡度阈值Si。在不同的网格内得出不同的坡度阈值,以解决坡度阈值随地形变化的自适应设置问题。实验区不同语义属性下的σ见表2所列,该σ是通过地面点样本实验得出的。

表2 实验区不同语义属性的地形计算因子

1.2.3 点云多尺度滤波

通过σ得出单个网格的Si后,计算网格内所有点与种子点的坡度值,并与Si进行比较。当坡度值大于Si时,将该点分为非地面点,予以去除。进行滤波操作时,算法按照网格平面尺度由大到小重复滤波。在较大尺度网格进行滤波时,可以去除较大的植被与建筑物点云；在小尺度网格进行滤波时, 可以滤除与地面点接近的地物点,如低矮植被、草丛、低矮地物等。按固定尺度间隔由大到小重复滤波使得较多的地面点数据得以保留。

1.3 地面点去噪与数据输出

通过自适应坡度滤波后的点云中仍然可能存在少量噪声点。这是由于少部分区域的网格属性或σ值很难完全符合真实地形的表达。此类的噪声点多为孤立点,因此设立较小尺度的网格,求取待定点的网格序号及其邻域网格,当邻域网格内不含点云的网格个数大于4时,则认为该点为噪声点,去除该点。最后去噪后的点云数据即为提取出的地面点。

1.4 滤波精度评价参数

为了量化表达基于多尺度网格的点云自适应滤波算法在不同场景中的滤波效果,利用点云数据分类中常用的混淆矩阵来定量评价算法的滤波精度。令T表示正确分类,F表示误分类,P表示地面点,N表示非地面点,则分类结果有以下4种:① TP,即滤波结果中分类正确的地面点数(TP);② FP,即非地面点被误分类为地面点的个数(FP);③ TN,即分类正确的非地面点数(TN);④ FN,即地面点被误分类为非地面点地的个数(FN)。滤波精度的评价参数有Ⅰ类误差e1、Ⅱ类误差e2、总体误差e及Kappa系数。e1、e2、e的计算公式为：

e1=FN/(TP+FN),e2=FP/(FP+TN),

e=(FN+FP)/N总,N总=TP+FP+TN+FN。

Ⅰ类误差反映了提取地面点的完整度,Ⅱ类误差反映了提取地面点的纯度,总体误差反映了算法的可行性。

2 算法实验与结果分析

本文在PCL点云库的基础上采用Matlab进行计算和点云数据滤波成果展示。 在安徽理工大学校内选取多处复杂的地形场景,利用地基站扫描仪和无人机摄影测量获取实验数据,并从定性和定量计算2个角度进行点云滤波实验结果分析。

2.1 单一场景下的点云滤波实验

实验区域为安徽理工大学校园内一处缓坡,所得数据为范围较小的边坡点云数据。该组数据采用海达云HS650脉冲式三维激光扫描仪获得,数据为三维坐标值,无颜色及纹理信息。实验区域的点云平均密度为11 200 点/m2,区域的长度约为18.0 m,宽度约为7.5 m,最大的高程差为5.0 m。实验区主要为缓坡地形,无较大的起伏,最大坡度约为30°,包括的地物为樟树、低矮的阔叶树及少量杂草。 单一场景实验区域的点云数据及滤波结果如图3所示。

图3 单一场景实验点云数据及滤波结果

从图3b可以看出，实验区内的植被点被有效分离,包括与地面较接近的低矮植被也得到滤除,同时整个地面的点云数据得到保留,经过去噪处理后,实现了地面地形的真实还原。 由于采用的数据是由地基站扫描仪得到的,在图3c中出现的裂缝是在扫描过程中被植被遮挡而形成的空洞。从图3d可以看出，在滤除非地面点的同时,尽可能地保证了分离出的非地面点是真实地表中的非地面点,有效减少了地面点的误分类。

2.2 复杂度场景下的点云滤波实验

在安徽理工大学校园内选取3块实验区。实验区的选取标准是尽量表达复杂的自然场景, 包括平地、缓坡、陡坡、河流、道路、建筑群、大小不同的植被、灌木从等,以验证本文算法在复杂地形场景中的滤波效果和算法的实用性。

该3块实验区域采用大疆精灵4RTK无人机飞行器拍摄影像数据,再对影像数据进行三维实景模型构建,最后输出数据格式为LAS的点云数据。点云数据包括点位三维坐标值、点位颜色属性值及点位强度信息等。

2.2.1 多地物的区域滤波实验

该区域为区块1，平均点云密度为260点/m2,区域的长度约为100.0 m,宽度约为56.0 m,最大的高程差为14.0 m。整个区域包含复杂地形和大量植被,是典型的滤波难度大的实验区域。该区域的南边为平地和斜坡,斜坡中部高程值较小,两边地形坡度较大;区域中间部分为宽度约为12 m的河沟,北边为坡度变化小的缓坡;区域内地物种类繁多,有凉亭、假山石、护栏、较高的梧桐树、中等高度的樟树、地表灌木丛、河中的芦苇丛等,同时测区内含有2条完整道路线。区块1的点云数据及滤波结果如图4所示。

图4 区块1的点云数据及滤波结果

算法的初始网格宽度为8.0 m,最低网格宽度为0.5 m。从图4b可以看出，与原始数据对比,地表上的地物基本上得到滤除，包括凉亭、梧桐树、樟树及地表的灌木丛；但是在河中的芦苇丛处仍然存在较多的点云数据,通过实际考察可知，芦苇丛高度与河岸边高度基本一致,在滤波阈值计算时不能很好地分离数据,造成了点云的残留。从图4b还可以看出河沟的边岸线和部分边坡都得到保留,整个地面和坡面数据保留较为完整,无噪声点。图4b中出现的点云空洞是由于数据是由影像数据生成的点云数据, 存在较多地面点被植被遮挡的情况,但是滤波结果对于地形的表达较完整。

2.2.2 植被覆盖度高的区域滤波实验

该区域为区块2，主要为坡度变化明显的斜坡和大量植被,植被有樟树、针叶松、阔叶树、地表灌木、地表草丛等。该区域地形两边为缓坡,中间为凸出的坡地；平均点云密度为960点/m2,区域的长度约为50.0 m,宽度约为40.0 m,最大的高程差为7.8 m。区块2的点云数据及滤波结果如图5所示。

因为区块内的非地面点仅包含植被点,所以以最大树冠直径5.0 m作为初始网格宽度,最低网格宽度为0.3 m,以达到非地面点彻底滤除的效果。 对比图5c与图5d可知,该区域内的非地面点滤除效果良好,并且通过局部地形自适应坡度阈值的设置可以很好地在顾及地形变化的同时剔除非地面点。不同尺度的网格设置可以去除不同大小的植被。点云中的空洞是滤波算法将非地面点剔除和数据源2个方面造成的。

2.2.3 附含建筑物群的区域滤波实验

该区域为区块3，平均点云密度为570点/m2,区域的长度约为80.0 m,宽度约为40.0 m,最大的高程差为22.0 m。该区域地形较简单,主要是缓坡和平坦地面;主要的地物类型是建筑物和地面植被。 区块3的点云数据及滤波结果如图6所示。

从图6b可以看出，点云数据经过滤波后可以得到很完整的地面点云数据,区域中的建筑物和植被均已滤除。在实验中，初始网格宽度设置时需要考虑网格内包含的真实地面点数据,因此初始宽度设置为10.0 m,最低宽度为0.5 m。实验结果表明,建筑物点云很难一次彻底滤除, 对于这种混合地形的点云数据,需要多次计算才能得到真实地面点云。

图6 区块3点云数据及滤波结果

2.3 复杂场景下的点云滤波精度评定

基于TIN的点云滤波、基于二次曲面拟合滤波、传统坡度滤波是常用的3种滤波算法，将本文算法与上述3种滤波算法进行对比。3个区块4种滤波算法的分类结果、滤波误差计算结果和算法耗时对比见表3所列。3个区块的样本点是利用Terrasoild软件和手工交互的方式进行分类后的地面点和非地面点。

表3 3个区块4种滤波算法的滤波误差计算结果和算法耗时对比

点云滤波算法一般要求在控制Ⅱ类误差的基础上尽量降低Ⅰ类误差,即在保证分离出的地面点最大程度上是地面点的同时,尽可能分离出较多的地面点。Kappa系数的大小反映了分离出的地面点位置的准确性程度。总体精度和Kappa系数保证了滤波得到的点云数据的数量和质量。

通过计算结果可以分析得出:

(1) 本文算法与传统坡度滤波算法相比,在总体精度和Kappa系数上都有所提高,表明本文算法分离出的地面点更多,并且与实际地表地形更加符合;但是在计算效率上小于传统的坡度滤波算法,这是由于多尺度网格的设立使得本文算法进行了多次运算,造成了本文算法计算效率降低。同其他2种算法相比,本文算法计算效率较高。

(2) 在同一区块，本文算法总体精度较高,表明本文算法与其他3种滤波算法相比同样具有很好的应用价值,甚至优于其他3种算法;本文算法Kappa系数值较大或者接近于其他3种算法,尤其在区块1的复杂地形中,本文算法的Kappa系数值最高,表明其在地形模拟和地面点提取方面具有很好的效果。

(3) 区块1内4种算法的分类精度和Kappa系数都低于区块2和区块3,这主要是由于区块1的实际地形复杂,地面点和非地面点的分离界限不明确,滤波算法很难从无序点云中提取出地面点;但是Kappa系数值在0.7附近,说明可以获得大体的地形数据,但是局部细节仍有损失。区块2内的地形包括平地、缓坡及斜坡,但是整体地形变化较均匀,含有的低矮植被掩盖了真实地表信息,本文算法的总体精度和计算效率虽然不是最优,但是仍具有较好的滤波效果;区块3内4种算法分类精度和Kappa系数值优于区块1和区块2,但实际实验中该区块由于含有高大建筑物,需要经过多次滤波后才能得出最佳分类效果。

3 结 论

为实现海量LiDAR点云数据中地面点与非地面点的有效分离,本文在传统坡度滤波算法的基础上,提出基于多尺度网格的点云自适应坡度滤波算法。将本文算法在3个复杂场景下进行应用验证，并与3种常用滤波算法的分类精度进行对比,得到以下结论:

(1) 本文提出在构建的虚拟网格内计算点云空间占比,并以此划分网格语义属性,通过地形计算因子计算不同网格的坡度阈值,达到坡度滤波阈值随地形变化自适应设置的目的。

(2) 多级尺度网格构建和坡度阈值的自适应设置在一定程度上解决了坡度滤波算法在边缘和陡坡区域的过滤波问题,有效降低了坡度滤波算法的Ⅱ类误差。

(3) 本文算法总体精度大于85%,Kappa系数值大于0.7,表明该算法在大多数地形中有较好的适用性和分类的准确性。

(4) 本文算法虽然可以自适应设置坡度阈值,但是前期设置的参数需要选取样本值进行实验来获取,增加了计算难度,该问题需进一步研究;对于混合地形区域,本文算法仍有局限性,需要多次实验才能取得良好的滤波效果。