基于主动学习的机载LiDAR点云滤波方法与流程

本发明涉及地理空间信息系统技术领域，特别是涉及一种基于主动学习的机载lidar点云滤波方法。

背景技术：

机载lidar(lightdetectionandranging，激光雷达)技术近年来发展十分迅速，现已成为获取地球空间信息的重要的测量方式。机载lidar系统主要由全球定位系统(gps)、激光扫描仪(ls)以及惯性导航系统(ins)组成。该系统能够主动地向地面发射激光脉冲，并接收来自地面目标物体的回波信息。相较于传统的被动遥感技术，机载lidar技术不受外界环境的影响，能够24小时全天候地进行数据采集。此外，机载lidar系统发射的激光脉冲能够穿透植被打到地面，因此能够有效解决传统光学影像数据易受树木遮挡的影响。现如今，机载lidar技术已广泛应用于众多的地球空间信息领域，例如三维城市建模、森林参数估测、电力线提取等。

实现大多数的点云后处理应用往往需要首先将地形点与地物点进行有效分离，进而建立数字地面模型(dtm)，此过程通常称之为点云滤波。近年来，针对机载lidar点云滤波具有较多的研究。

此外，随着人工智能方法的发展，有相关技术将机器学习方法或者深度学习方法应用到点云滤波中，但现有技术中，此类监督学习的方法最大的问题在于需要大量的样本标记，过多的样本标记将是耗时且费力的，而且滤波精度不够理想。

技术实现要素：

本发明的目的在于解决现有技术需要大量的样本标记，滤波精度不够理想的问题，提出一种基于主动学习的机载lidar点云滤波方法。

一种基于主动学习的机载lidar点云滤波方法，包括以下步骤：

s1，获取点云数据并进行低位噪声点剔除；

s2，采用多尺度形态学运算自动获取并标记训练样本集

s3，对训练样本集进行特征提取并建立svm模型；

s4，采用训练模型对候选样本集进行分类，分为候选地面点集u^g和候选非地面点集u^ng，将神谕(oracle)设置为候选点集至拟合曲面距离的s型函数，每次迭代分别从候选地面点集和候选非地面点集中各选取q个点加入训练样本集t^k中并更新训练模型，一直迭代，直到候选地面点集和候选非地面点集中的点云个数均不再大于q为止，最后将最新训练模型的分类作为点云滤波结果；

s5，滤波优化。

根据本发明提供的基于主动学习的机载lidar点云滤波方法，首先采用多尺度形态学运算自动获取并标记训练样本，然后，采用主动学习的方式通过不断迭代更新支持向量机(svm)训练模型来获得精度越来越高的滤波结果，实现了在无需人工样本标记前提下点云数据的自动分类，有效解决了传统监督学习方法需要过多样本标记的问题，实验结果表明，该方法能够在复杂地形区域在无人为干预条件下获得不错的滤波结果，实验结果表明，相较于其它滤波方法，本发明能够获得最小的滤波总误差，在不同的地形环境下均能获得良好的滤波效果，本发明提供的方法其平均ⅰ类误差和ⅱ类误差较为平衡，说明方法在有效去除非地面点的同时又能保护有效地形点，使得该方法能够适应复杂的地形环境，并能取得较高的滤波精度。

另外，根据本发明上述的基于主动学习的机载lidar点云滤波方法，还可以具有如下附加的技术特征：

进一步地，所述步骤s1具体包括以下步骤：

s11，将三维点云数据data(x,y,z)转换为二维栅格数据dsm(i,j)，保留转换映射关系r；

s12，对dsm(i,j)进行均值滤波得到对比dsm(i,j)和之间各个栅格特征值的变化量，将dsm(i,j)中变化量大于阈值栅格的特征值替换为均值滤波的结果，公式表示如下：

式中，abs(·)为取绝对值，t1为栅格特征值变化阈值；

s13，按照数据转换映射关系r对dsm(i,j)进行逆变换，得到

s14，对比各个点的观测高程值z和滤波后的高程值将差值大于阈值的点判定为噪声点并进行剔除。

进一步地，所述步骤s2具体包括：

分别采用大窗口50×50和小窗口3×3进行形态学滤波，将大窗口形态学滤波结果中的地面点作为初始训练正样本，而将小窗口形态学滤波结果中的地物点作为初始训练负样本。

进一步地，所述步骤s3包括：

将几何特征向量作为训练模型的输入值，几何特征向量主要由两部分组成，一部分是基于局部结构协方差张量计算得到，另一部分则基于点云的高程值计算得到，具体包括以下步骤：

遍历各个训练样本点，通过构建kd树找到各点的k个临近点集合kn(p)；

利用这k个临近点集合计算各样本点的局部结构协方差张量：

式中，为临近点集合kn(p)的中心点，公式计算如下：

协方差covp为正定矩阵，计算得到该矩阵的三个特征值λ0≥λ1≥λ2≥0，以及对应的特征向量e0、e1和e2，利用这三个特征值和三个特征向量分别计算得到：不变性、本征熵、各向异性、平面性、线性、点性、曲面变化性、垂直性，公式表示如下：

各向异性：(λ0-λ2)/λ0(4)

平面性：(λ1-λ2)/λ0(5)

线性：(λ0-λ1)/λ0(6)

点性：λ2/λ0(7)

曲面变化性：λ2(8)

基于点云的高程值可以计算另一部分的特征向量，公式表示如下：

高程变化幅度：zmax{kn(p)}-zmin{kn(p)}(9)

低点度：zp-zmin{kn(p)}(10)

高点度：zmax{kn(p)}-zp(11)。

进一步地，所述步骤s4具体包括：

将神谕设置为各待定点到拟合曲面距离的sigmoid函数，公式表示如下：

式中，p为待定点，zp为该点的观测高程值，为采用径向基函数进行曲面拟合，f(p)为点p到拟合曲面的距离；

然后按照主动学习算法流程，每次迭代分别从地面点集{g}^k中选出q个s(f)最小的点作为地面点，从地物点集{ng}^k中选出q个s(f)最大的点作为地物点，加入到训练样本集合；

更新训练模型，一直迭代直到地面点集和地物点集中点的个数不再大于q为止。

进一步地，所述步骤s5具体包括以下步骤：

对滤波结果进行格网剖分；

获取各个格网内的最低点作为地面种子点；

利用这些地面种子点采用rbf函数对点云中的各个点(pi(xi,yi,zi))，i＝1,2,…,n)进行高程差值拟合计算，获取各个点的拟合高程值

根据拟合曲面计算获取各个点所对应的横纵方向的坡度变化值依次遍历点云数据，将满足公式(13)的点判定为非地面点并进行移除；

式中，zi为各个点的观测值，为对应的拟合高程值，te为阈值常量，表示平均相邻地面点至拟合曲面的距离。

附图说明

本发明实施例的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为本发明实施例提供的基于主动学习的机载lidar点云滤波方法的流程图；

图2为神谕设置示意图；

图3为三组实验数据的地形特征示意图，其中，(a)为sample11；(b)为sample12；(c)为sample21；

图4为sample11数据的滤波结果对比图；

图5为sample12数据的滤波结果对比图；

图6为sample21数据的滤波结果对比图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供的基于主动学习的机载lidar点云滤波方法，包括以下步骤s1～s4，请同时结合图1，对本实施例的方法进行详细说明：

s1，获取点云数据并进行低位噪声点剔除；

其中，由于受到仪器自身或者外界环境的影响，所获取的点云数据中往往包含有噪声点。噪声点存在，尤其是低位噪声点会对点云接下来的处理操作带来干扰。如许多滤波方法往往假定局部最低点为地面点，低位噪声点的存在会对滤波结果带来严重误判。此外，噪声点的存在也会对特征提取带来干扰，影响模型建立精度。因此，需要首先将低位噪声点进行去除。

具体的，步骤s1包括以下步骤：

s11，将三维点云数据data(x,y,z)转换为二维栅格数据dsm(i,j)，保留转换映射关系r；

s12，对dsm(i,j)进行均值滤波得到对比dsm(i,j)和之间各个栅格特征值的变化量，将dsm(i,j)中变化量大于阈值栅格的特征值替换为均值滤波的结果，公式表示如下：

式中，abs(·)为取绝对值，t1为栅格特征值变化阈值；

s13，按照数据转换映射关系r对dsm(i,j)进行逆变换，得到

s14，对比各个点的观测高程值z和滤波后的高程值将差值大于阈值的点判定为噪声点并进行剔除。

具体流程可以参见表1

表1点云去噪算法流程

s2，采用多尺度形态学运算自动获取并标记训练样本集；

其中，传统的监督学习方法虽然能够获得不错的滤波精度，但需要大量的样本标记，不仅耗时费力而且算法自动化程度较低。如何在无人为干预情况下实现非监督学习，自动获取并标记正、负训练样本(即地面点和地物点)是研究的重点。本实施例采用多尺度形态学运算实现非监督学习。

采用不同尺度的滤波窗口进行形态学开运算，可以获得不同的滤波结果。当采用大尺寸的滤波窗口时，虽然能够有效地滤除大型建筑物，但无法保护地形细节(如小的地形凸起被滤除)。当采用小尺寸的滤波窗口时，虽然小的地形凸起得到了保护，但部分建筑物并未滤除。从另一个角度分析，若采用大窗口进行形态学滤波，虽然部分地形凸起会被滤除，但未被滤除区域的地面点却是准确的。若采用小窗口进行形态学滤波，虽然部分建筑物未被滤波，但滤除的建筑物部分可以保证是真实的地物点。

因此，具体在本实施例中，分别采用大窗口(50×50)和小窗口(3×3)进行形态学滤波，将大窗口形态学滤波结果中的地面点作为初始训练正样本，而将小窗口形态学滤波结果中的地物点作为初始训练负样本。

s3，对训练样本集进行特征提取并建立svm模型；

其中，步骤s3包括：

将几何特征向量作为训练模型的输入值，几何特征向量主要由两部分组成，一部分是基于局部结构协方差张量计算得到，另一部分则基于点云的高程值计算得到，具体包括以下步骤：

遍历各个训练样本点，通过构建kd树找到各点的k个临近点集合kn(p)；

利用这k个临近点集合计算各样本点的局部结构协方差张量：

式中，为临近点集合kn(p)的中心点，公式计算如下：

协方差covp为正定矩阵，计算得到该矩阵的三个特征值λ0≥λ1≥λ2≥0，以及对应的特征向量e0、e1和e2，利用这三个特征值和三个特征向量分别计算得到：不变性、本征熵、各向异性、平面性、线性、点性、曲面变化性、垂直性，公式表示如下：

各向异性：(λ0-λ2)/λ0(4)

平面性：(λ1-λ2)/λ0(5)

线性：(λ0-λ1)/λ0(6)

点性：λ2/λ0(7)

曲面变化性：λ2(8)

基于点云的高程值可以计算另一部分的特征向量，公式表示如下：

高程变化幅度：zmax{kn(p)}-zmin{kn(p)}(9)

低点度：zp-zmin{kn(p)}(10)

高点度：zmax{kn(p)}-zp(11)。

s4，采用训练模型对候选样本集进行分类，分为候选地面点集和候选非地面点集，将神谕设置为候选点集至拟合曲面距离的s型函数，每次迭代分别从候选地面点集和候选非地面点集中各选取q个点加入训练样本集中并更新训练模型，一直迭代，直到候选地面点集和候选非地面点集中的点云个数均不再大于q为止，最后将最新训练模型的分类作为点云滤波结果；

其中，在主动学习(activelearning，al)中，学习器可以“主动地”向学习器之外的某个神谕(oracle)进行查询来获得训练例的标记，然后再将这些有标记的示例作为训练例来进行监督学习。因此，需要研究正确、恰当的“神谕”设置方法，以获得越来越精确的点云滤波结果。

因此，在本实施例中，将神谕设置为各待定点到拟合曲面距离的sigmoid函数，公式表示如下：

式中，p为待定点，zp为该点的观测高程值，为采用径向基函数进行曲面拟合，如图2曲线所示，f(p)为点p到拟合曲面的距离；

然后按照表2的主动学习算法流程，每次迭代分别从地面点集{g}^k中选出q个s(f)最小的点作为地面点，从地物点集{ng}^k中选出q个s(f)最大的点作为地物点，加入到训练样本集合；

更新训练模型，一直迭代直到地面点集和地物点集中点的个数不再大于q为止。

表2主动学习下点云滤波算法流程

s5，滤波优化。

其中，通过计算点基元的特征向量来建立svm模型，进而实现对点云数据的分类，但基于点基元获取特征向量往往会存在误差，因此分类的结果有可能会出现误判，本实施例采用基于坡度的方法来对分类结果进行优化。

为进一步优化滤波结果，具体的，步骤s5包括以下步骤：

对滤波结果进行格网剖分；

获取各个格网内的最低点作为地面种子点；

利用这些地面种子点采用rbf函数对点云中的各个点(pi(xi,yi,zi))，i＝1,2,…,n)进行高程差值拟合计算，获取各个点的拟合高程值

根据拟合曲面计算获取各个点所对应的横纵方向的坡度变化值依次遍历点云数据，将满足公式(13)的点判定为非地面点并进行移除；

式中，zi为各个点的观测值，为对应的拟合高程值，te为阈值常量，表示平均相邻地面点至拟合曲面的距离。

为了验证本发明的效果，本实施例采用国际摄影测量与遥感学会(isprs)提供的专门用于检验滤波效果的三组数据(sample11、sample12、sample21)进行实验分析(https://www.itc.nl/isprs/wgiii-3/filtertest/)。这三组实验数据由optechaltm三维激光扫描仪获取，点间距为1-1.5m。这三组实验数据包含有不同的地形、地貌特征，因此有利于检验滤波方法在不同环境下的滤波效果。如图3(a)所示，sample11区域地形坡度变化较大，存在有较为密集的植被，并且斜坡上建有房屋。sample12区域主要滤波难点在于存在有复杂的建筑物(图3(b))，有利于检测滤波方法在主城区的滤波效果。对于大多数滤波算法，与地形相连的地物(如桥梁)一般较难被剔除，通过对sample21(图3(c))进行实验，能检测本实施例的方法是否能够有效滤除此类地物。

图4、图5和图6分别为三组实验数据滤波前后的对比图。(a)为各个样本数据生成的数字表面模型(dsm)，(b)为由人工选取的准确的地面点生成的数字地面模型(dtm)，(c)为由本实施例的方法获取的地面点生成的数字地面模型(dtm)。从图中可以看出本实施例的滤波结果十分接近准确的滤波结果，在三种不同的复杂地形环境下均能获得不错的滤波效果。sample11中的低矮植被，斜坡上的房屋，sample12中的复杂建筑物以及sample21中与地形相连的桥梁都得到了有效的去除。由此可以看出，本实施例的方法能够适应多种复杂的地形环境，滤波方法具有较强的鲁棒性。

为了更加客观的评价本实施例的滤波方法，选用ⅰ类误差(t1)、ⅱ类误差(t2)和总误差(total)进行定量评价。ⅰ类误差又称为拒真误差，指的是地面点误判为地物点的比例；ⅱ类误差又称为纳伪误差，指的是地物点误判为地面点的比例；总误差值得所有误判点所占的比例。如表3所示，建立滤波结果的交叉矩阵，并按照公式(14)-(16)计算三类误差。

表3交叉矩阵

t1＝b/(a+b)(14)

t2＝c/(c+d)(15)

total＝(a+b)/(a+b+c+d)(16)

本实施例的方法的滤波结果以及近年来一些代表性滤波方法的三类误差对比结果如表4-6所示。jahromi等(2011)将ann应用于点云滤波中。mongus和zalik(2012)通过进行薄板样条(tps)来实现点云滤波。zhang和lin(2013)采用分割点云的方式来改进传统的渐进加密不规则三角网(ptd)算法。li等(2014)提出一种改进的高帽滤波方法以保护地形细节。hui等(2016)将传统的曲面拟合滤波法和形态学滤波法进行结合，获得了不错的滤波效果。zhang等(2016)提出一种布料模拟滤波算法，能够在较少的参数调节下获得不错的滤波结果。hu和yuan(2016)将cnn应用于点云滤波中取得了不错的滤波结果。需要指出的是，表4给出的是该方法在采用十组训练数据建立模型的分类结果。如果采用1.7亿个样本数据进行训练，能够获得更小的滤波误差。li等(2017)提出了一种不需要指定最大滤波窗口的形态学滤波改进方法，提升了算法的适用性。ni等(2018)采用图割的方法来实现点云滤波。rizaldy将全卷积网络(fcn)应用于点云滤波中，能够在少量训练样本下获得较高的滤波结果。

从表4可以看出，相较于其他10中滤波方法，本实施例的方法能够获得最小的平均总误差(5.51％)。由此可见，本实施例的方法能够获得较高的滤波精度。此外，本实施例的方法在三种样本数据上均能获得较小的总误差。由此可以得出，本实施例的方法对不同的地形环境具有较强的适应性。从表5和表6可以看出，本实施例的滤波方法的平均ⅰ类误差(6.68％)和平均ⅱ类误差(4.84％)较为接近，表明本实施例的方法既能有效去除非地面点又能有效地保护地形细节信息不被破坏。此外，sample11的ⅰ类误差和ⅱ类误差都相对较大，这主要是因为sample11的地形坡度变化较大，部分地形凸起区域被误判为非地面点而被剔除，形成ⅰ类误差；部分建筑物屋顶与地形相连，被误判为地面点，形成ⅱ类误差。在sample12区域，部分低矮地物未被正确剔除，使得ⅱ类误差较大。sample21区域地形较为平坦，只要少部分点被错分，滤波结果十分接近正确结果，滤波总误差也最小(1.23％)。

由此可以得出，本实施例的方法在地形平坦区域滤波效果较好，而在地形坡度变化较大区域滤波效果较差。此特点也与其他的大多数滤波方法相符。

表4总误差对比结果

表5ⅰ类误差对比结果

表6ⅱ类误差对比结果

上述实验结果表明，本实施例的方法相较于其他十种近年来代表性的滤波方法，能够获得最小的滤波总误差。此外，三组实验数据的总误差都相对较小，表明本实施例的方法在不同的地形环境下均能获得良好的滤波效果。本实施例的方法的平均ⅰ类误差和ⅱ类误差较为平衡，表明本实施例的方法在有效去除非地面点的同时又能保护有效地形点。

综上，根据本发明提供的基于主动学习的机载lidar点云滤波方法，首先采用多尺度形态学运算自动获取并标记训练样本，然后，采用主动学习的方式通过不断迭代更新支持向量机(svm)训练模型来获得精度越来越高的滤波结果，实现了在无需人工样本标记前提下点云数据的自动分类，有效解决了传统监督学习方法需要过多样本标记的问题，实验结果表明，该方法能够在复杂地形区域在无人为干预条件下获得不错的滤波结果，实验结果表明，相较于其它滤波方法，本发明能够获得最小的滤波总误差，在不同的地形环境下均能获得良好的滤波效果，本发明提供的方法其平均ⅰ类误差和ⅱ类误差较为平衡，说明方法在有效去除非地面点的同时又能保护有效地形点，使得该方法能够适应复杂的地形环境，并能取得较高的滤波精度。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。