一种基于虚拟种子点的机载LiDAR地面点云滤波方法

摘要

本发明公开了一种基于虚拟种子点的机载LiDAR地面点云滤波方法，属于机载LiDAR点云数据分类领域。本发明的步骤为：首先，将机载LiDAR点云格网化，利用格网内点云的分布维度和高程分布直方图构建地面虚拟种子点；其次，使用多尺度形态学对虚拟种子点评估，剔除不属于地面的虚拟种子点，并且对点云进行初步滤波；最后，利用评估后的虚拟种子点构建初始TIN网，迭代加密TIN网完成LiDAR数据的地面点云滤波。本发明既降低了滤波算法的时间耗费，又保证了滤波算法的滤波精度，并且适用范围广。

说明书

技术领域

本发明涉及一种机载LiDAR点云数据分类方法，特别是涉及一种基于虚拟种子点的机载LiDAR地面点云滤波方法。

背景技术

机载激光雷达技术(Light Detection And Ranging，简称LiDAR)作为一种快速发展的新型探测手段，可以直接获取目标的三维坐标，具有主动性强、探测精度高、作业周期短等优势。该技术可用于数字地面模型提取，建筑物建模，植被参数估计，电力线提取等诸多领域。在具体应用过程中，最重要的数据处理步骤就是地面点与非地面点的分离，即LiDAR数据的地面点云滤波(以下简称点云滤波)。LiDAR数据量大，人工分类耗时耗力，故亟需自动化的点云滤波方法。

LiDAR数据的地面点云滤波的主要目的是对地面点进行特征分析，结合其与非地面之间的区别，构造相应的分类函数，把未知类别的点映射到所构造的地面类别或非地面类别中，从而得到地面点集和非地面点集。目前，点云滤波主要有基于坡度，基于表面、基于区块、基于聚类的四种模型。2000年，Axelsson在期刊《International Archives of Photogrammetry and Remote Sensing》中发表“DEM Generation from Laser Scanner Data Using Adaptive TIN Models”一文，研究了渐进三角网加密滤波算法，其基本原理是先建立粗糙的不规则三角网，然后自上而下地对粗糙不规则三角网进行精细化，得到精细的三角网模型，参与构建精细的不规则三角网的点被分为地面点，其他点被分为非地面点；2003年，Zhang等在期刊《IEEE Transactions on Geoscience and Remote Sensing》41卷第4期发表“A Progressive Morphological Filter for Removing Nonground Measurements From Airborne LIDAR Data”一文，研究了多尺度数学形态学滤波算法，其基本原理是地面点低于其领域内的非地面点，然后自下而上增加领域的范围，得到不同尺度下非地面点集，最后不属于非地面点集的点被归为地面点集。然而在很多时候，单独使用自上而下的渐进三角网加密滤波算法效率低，耗时长；单独使用自下而上的多尺度形态学滤波算法精度低，效果差。这是因为地面通常比较复杂，难以用简单的模型表示。

发明内容

1.发明要解决的技术问题

针对现有技术中，渐进三角网耗时比较大，难以快速处理点云数据，但是滤波精度高，多尺度形态学当地形比较复杂时，滤波精度低，但是运行效率高的问题。本发明提供一种基于虚拟种子点的机载LiDAR地面点云滤波方法，通过构建虚拟种子点将渐进三角网和多尺度形态学结合起来，既降低了滤波算法的时间耗费，又保证了滤波算法的滤波精度，并且适用范围广。

2.技术方案：

原理：点云滤波是根据地面的特征构造一个地面模型，把机载LiDAR系统获取的点映射到相应的地面类和非地面类中，从而得到机载LiDAR的地面分类结果。本发明通过点云格网化，分析格网内点云维数特征与高程直方图信息，构建地面虚拟种子点，并利用多尺度形态学评估初始虚拟种子点，最后由虚拟种子点生成初始TIN网，并加密分类出地面点集和非地面点集。

一种基于虚拟种子点的机载LiDAR地面点云滤波方法，包括以下步骤(具体的步骤流程图见图1)：

步骤1：通过机载LiDAR系统，获取描述地表的LiDAR点云数据，用二维格网在XY平面上划分LiDAR点云数据，得到格网集合GridSet；

步骤2：按照公式①～⑤计算GirdSet中的每个格网内点云三维分布维数特征参数D₁、D₂、D₃，D₁、D₂、D₃分别表示一维、二维、三维的特征参数，将D₂>0.5的格网归为呈二维分布的格网集合GirderSet2，将D₂≤0.5的格网归为非二维分布的格网集合GirdSet1，取格网集合GridSet2中的每个格网内点云在X、Y方向的平均值为虚拟种子点在X、Y方向的坐标值x、y，Z方向的高程分布直方图的峰值所对应的高程为虚拟种子点在Z方向的坐标值z，所有虚拟种子点构成虚拟种子点集合VP；

S＝{X,Y,Z}-mean({X,Y,Z}) 公式①

λ₁，λ₂，λ₃＝eigenvalue(S^T×S)公式②

式中，{X,Y,Z}表示格网内三维点云坐标的集合所形成的矩阵，mean()表示求解矩阵列的平均值，λ₁，λ₂，λ₃为特征值，eigenvalue()表示求解矩阵的特征值，S^T表示S的转置，D₁、D₂、D₃表示一维、二维、三维的特征参数；

步骤3：按照公式⑥～⑧，以初始窗口W_min对虚拟种子点集合VP进行形态学运算并比较，将Δh>M_threhold的点归为非地面点，然后按照公式⑨～⑩增加窗口尺寸和阈值，进行再一轮的形态学运算分类，重复执行此过程，直至窗口的尺寸大于最大建筑物的尺寸；

Δh＝Z_p-d(e(p))>

W_k＝W_min×2^k公式⑨

M_threhold＝S×W_k公式⑩

式中，(x_p,y_p,z_p)是三维点云的坐标，w是窗口，即点的在XY平面上的领域，d(p)表示数学形态学膨胀运算，e(p)表示数学形态学腐蚀运算，d(e(p))表示数学形态学的开运算(即先腐蚀在膨胀)，Z_p表示点p的高程值，W_k和M_threhold表示第k次循环使用的窗口尺寸和阈值，S表示地面的坡度；

将集合VP中点分为两类，不属于地面点的种子点DS和真正属于地面的虚拟种子点GS，将DS剔除，保留真正属于地面的虚拟种子点GS构成集合VP′；

步骤4：利用虚拟种子点集合VP′中的点构建初始TIN网，然后，根据初始TIN网模型，迭代地将待判定点集和GirdSet1中每个格网内满足判定条件的点加密到TIN网中，并更新TIN网，直至没有满足条件的点，加入TIN网的点被分为地面点，其他点被分为非地面点；

步骤5：归并步骤3和步骤4的分类结果得到地面点集和非地面点集(如图2f)。

步骤3中，所述初始窗口W_min＝5m。

步骤3中，将集合VP中点分为两类，即利用集合VP中点的类别对格网集合GridSet2中点进行分类，具体方法如下：

DS所在格网内与DS高程差大于点云测量误差的点被分为待判定点集，其他点被分为非地面点，将GS所在格网内与GS高程差小于点云测量误差的点被分为地面点，其他点被分为非地面点。

步骤4中判定条件为：点到三角面的距离d_h，以及点到三角面各顶点的最大角度值θ，都小于设定的阈值。

3.有益效果

相比现有技术，本发明将渐进三角网加密和多尺度形态学方法结合起来，合理的建立了未知分类点和类别之间的映射，提高了目标识别的性能。具体有益效果如下：

第一，引进了虚拟种子点的选取。与已有的地面点云滤波相比，本发明在利用机载LiDAR点云数据本身的特点添加了部分虚拟地面种子点，有效地避免了地面低值点对滤波结果的影响。

第二，改进了渐进三角网的优化。本发明利用多尺度形态学对点云数据完成了初步分类，减少了需要迭代进行分类的点云数量，有效地减少了迭代加密的耗时。

第三，实现了多层次的点云滤波。本发明步骤2中引入了点云分布维度，对格网进行了分类，构建了虚拟种子点，为多层次的点云滤波做了准备工作，利用多尺度形态学完成了初始滤波，利用渐进三角网加密完成了第二次滤波。

综上，本发明通过构建虚拟种子点结合了渐进三角网加密与多尺度形态学方法，能够用多尺度形态学完成初步滤波并引导后续渐进三角网加密，实践证明，该方法能够有效地进行地面点云滤波，得到比较满意的结果，并具有高效率、高精度的特点。

附图说明

图1为基于虚拟种子点的机载LiDAR地面点云滤波方法流程图；

图2为基于虚拟种子点的机载LiDAR地面点云滤波方法各步骤示意图。

具体实施方式

本实施例采用图2a所示LiDAR点云数据作为地面点云滤波识别的初始点云数据，该数据采用Optech ALTM Gemini系统于2009年11月采集获得，点云密度约为3.62点/m²。实施例区域的图像大小为68.16m×74.59m。本实例具体实施采用MATLAB和TerraSolid实现。具体的流程图如图1，实施步骤如下：

步骤1：通过机载LiDAR系统，获取描述地表的LiDAR点云数据，用二维格网在XY平面上划分LiDAR点云数据，得到格网集合GridSet(如图2b所示)；

即以MATLAB创建四维数组Dataset来读取LiDAR点云数据，初始化点的类别为3(表示待判定点)。根据LiDAR点云的范围，以3m的单元格网划分LiDAR点云，创建二维胞元数组保存所在格网内点的索引(即：点在Dataset中的行号)。

核心代码为：

％％读取LiDAR点云数据，初始化所有点的类别为3

Dataset＝dlmread('D:Lidar_DataShiyituRo_Shiyitu.txt')；

Dataset(:,4)＝3；

％以偏移量转换数据，便于运算处理

Min_d＝[min(Dataset(:,1))min(Dataset(:,2))min(Dataset(:,3))0]；

Trandataset＝bsxfun(@minus,Dataset,Min_d)；

X_max＝max(Dataset(:,1))；X_min＝min(Dataset(:,1))；

Y_max＝max(Dataset(:,2))；Y_min＝min(Dataset(:,2))；

X_range＝X_max-X_min；Y_range＝Y_max-Y_min；

％创建胞元数组Grid

Cell_row＝ceil(X_range/Grid_size)；

Cell_col＝ceil(Y_range/Grid_size)；

Grid＝cell([Cell_row,Cell_col])；

％％点云数据的格网划分并将其索引存储在胞元数组Grid中

for i＝1:size(Trandataset,1)

i_row＝ceil(Trandataset(i,1)/Grid_size)；

i_col＝ceil(Trandataset(i,2)/Grid_size)；

if i_row＝＝0

i_row＝1；

end

if i_col＝＝0

i_col＝1；

end

if isempty(Grid{i_row,i_col})

Grid{i_row,i_col}＝i；

else

Grid{i_row,i_col}＝[Grid{i_row,i_col},i]；

end

end

步骤2：按照公式①～⑤计算GirdSet中的每个格网内点云三维分布维数特征参数D₁、D₂、D₃，D₁、D₂、D₃分别表示一维、二维、三维的特征参数，将D₂>0.5的格网归为二维分布的格网集合GirdSet2(如图2d中的浅灰色和深灰色区域)，格网内的点云类别设为2。将D₂≤0.5的格网归为非二维分布的格网集合GirdSet1(如图2d中的黑色区域)，取格网集合GridSet2中的每个格网内点云在X、Y方向的平均值为虚拟种子点在X、Y方向的坐标值x、y，Z方向的高程分布直方图的峰值所对应的高程为虚拟种子点在Z方向的坐标值z，所有虚拟种子点构成虚拟种子点集合VP(如图2c)。

S＝{X,Y,Z}-mean({X,Y,Z})公式①

λ₁，λ₂，λ₃＝eigenvalue(S^T×S)>

式中，{X,Y,Z}表示格网内三维点云坐标的集合所形成的矩阵，mean()表示求解矩阵列的平均值，λ₁，λ₂，λ₃为特征值，eigenvalue()表示求解矩阵的特征值，S^T表示S的转置，D₁、D₂、D₃表示一维、二维、三维的特征参数；

核心代码为：

步骤3：按照公式⑥～⑧，以初始窗口W_min＝5m对虚拟种子点集合VP进行形态学运算并比较，将Δh>M_threhold的点归为非地面点，然后按照公式⑨～⑩增加窗口尺寸和阈值，进行再一轮的形态学运算分类，重复执行此过程，直至窗口尺寸大于最大建筑物的尺寸，本实例中最大建筑的尺寸为25m，所以要重复三次，其窗口分别为10m、20m、40m。本实例的坡度S取0.1。

Δh＝Z_p-d(e(p))公式⑧

W_k＝W_min×2^k公式⑨

M_threhold＝S×W_k公式⑩

式中，(x_p,y_p,z_p)是三维点云的坐标，w是窗口，即点的在XY平面上的领域，d(p)表示数学形态学膨胀运算，e(p)表示数学形态学腐蚀运算，d(e(p))表示数学形态学的开运算(即先腐蚀在膨胀)，Z_p表示点p的高程值，W_k和M_threhold表示第k次循环使用的窗口尺寸和阈值，S表示地面的坡度；

核心代码为：

将集合VP中点分为两类，不属于地面点的种子点DS和真正属于地面的虚拟种子点GS，将DS剔除，保留真正属于地面的虚拟种子点GS构成集合VP′(如图2e中的点)。利用集合VP中点的类别对格网集合GridSet2中点进行分类，具体方法如下：

DS所在格网内与DS高程差大于点云测量误差的点被分为待判定点集，其他点被分为非地面点，将GS所在格网内与GS高程差小于点云测量误差的点被分为地面点，其他点被分为非地面点。本实例中点云的垂直方向的测量误差约为0.2m。

步骤4：利用虚拟种子点集合VP′中的点构建初始TIN网(如图2e)，然后，根据初始TIN网模型，迭代地将待判定点集和GirdSet1中每个格网内满足判定条件的点加密到TIN网中，并更新TIN网，直至没有满足条件的点，加入TIN网的点被分为地面点，其他点被分为非地面点，判定条件为：点到三角面的距离d_h，以及点到三角面各顶点的最大角度值θ，都小于设定的阈值。本实例阈值的取1.4m和6°。

步骤5：归并步骤3和步骤4的分类结果得到地面点集和非地面点集(如图2f)。

为了测试本发明的性能，采用了ISPRS Commission III Working Group III/3提供的14个代表不同地形地物特征的典型样本数据进行了实验。并以这14个样本数据的ISPRS地面点云的真实分类结果为依据建立混淆矩阵，计算用以评价分类精度的Ⅰ类误差、Ⅱ类误差、总体误差和Kappa系数。计算结果见表1。

表1滤波精度评价表(单位：％)

实验表明，本发明的滤波性能较好，实验的平均总体误差为8.72％、Kappa系数为91.27％，鲁棒性较好，可以用于多种复杂地形。并且需要渐进三角网滤波的数据量较少，有效提升了其运行效率。