一种基于体元的机载lidar点云三维滤波方法
【专利摘要】针对现有研究所采用的数据组织形式均不利于发挥机载LIDAR点云真3D特性的局限性，本发明提出一种基于体元的机载LIDAR点云三维滤波方法，属于遥感数据处理领域，方法包括：读取原始LIDAR点云数据；异常数据辨识及剔除；二值体元数据结构模型建立；基于体元的3D滤波。该方法将离散的点云规则化为二值3D体数据，根据体元中是否包含有激光点分别为体元赋值1和0，包含激光点赋值1，无激光点赋值0，1和0分别代表目标和背景，进而选择局部高程最低作为地面种子体元并标记其邻域目标体元为地面体元，本发明有助于3D数据表达及处理。
【专利说明】
一种基于体元的机载LIDAR点云三维滤波方法
技术领域
[0001 ]本发明属于遥感数据处理领域，具体涉及一种基于体元的机载LIDAR点云三维滤 波方法。
【背景技术】
[0002] 机载激光雷达(LIDAR，即Light Detection And Ranging)点云数据滤波质量直接 影响着后续产品的质量及该技术在生产应用中的成败，是LIDAR技术研究的重点与难点。国 内外学者围绕点云数据滤波进行了深入研究，提出了非常多的滤波算法，经典的如基于内 插、基于坡度、基于数学形态学、基于方向扫描和基于聚类分割的滤波方法等。这些滤波算 法所采用的数据结构包括基于规则格网的、基于不规则三角网（TIN)的、基于剖面或一维扫 描线的和基于离散点的方法等。其中，扫描线为ID数据表达方式;栅格格网和TIN为2.5D的 数据表达方式，二者在建模连续表面(如裸露地面等)时表现极佳，但不能用于3D表达，其在 相同的水平坐标(X，Y)对应多个高程(或Z)值；点云可以是真3D表达，但其拓扑和邻域信息 很难利用。可见，现有方法所采用的数据组织均不利于发挥机载LIDAR真3D的技术优势。
【发明内容】
[0003] 针对现有技术的不足，本发明提出采用二值体元数据结构表达机载LIDAR点云数 据，并在此基础上提出了一种基于体元的三维滤波算法。二值体元数据结构模型是将离散 非规则分布的激光点集表示为以二值体元为基本单元的规则空间数据结构，是一种更为简 单的真3D数据表达方式，且其内部体元间隐含有易于被利用的拓扑和邻域关系。二值体元 模型的构建就是将包含场景目标的空间划分成三维体元网格，并将点云映射到三维体元网 格，进而根据体元中有无激光点为各体元赋值1或0。
[0004] 本发明的技术方案是：
[0005] -种基于体元的机载LIDAR点云的三维滤波方法，包括以下步骤：
[0006] 步骤1:读取原始LIDAR点云数据；
[0007] 步骤2:异常数据辨识及剔除；
[0008] 步骤3:二值体元数据结构模型建立；
[0009] 步骤4:基于体元的3D滤波。
[0010]所述步骤2的具体步骤如下：
[0011] 步骤2.1:对LIDAR点云建立三维格网索引；
[0012] 步骤2.1.1:确定数据集的轴向包围盒；
[0013 ]步骤2.1.2:依据激光点云的密度确定格网尺寸，进而将轴向包围盒划分为三维格 网；
[0014] 步骤2.1.3:将各激光点映射到相应的格网单元，获取其对应的格网索引；
[0015] 步骤2.2:根据激光点所在的格网单元间的拓扑关系筛选候选点；
[0016] 步骤2.3:利用局部距离离群因子判定候选点的异常程度；
[0017]步骤2.4:删除各个异常数据点，得到异常数据剔除后的数据集。
[0018]所述步骤3的具体步骤如下：
[0019]步骤3.1:体元结构构建，将包含场景的三维空间划分成三维体元网格，其依赖于 场景所包络的三维空间范围和体元大小。
[0020]步骤3.1.1:用异常数据剔除后的数据集的轴向包围盒表示三维空间范围；
[0021 ]步骤3.1.2:计算空间分辨率(即体元大小）；
[0022]步骤3.1.3:对轴向包围盒进行划分得到3D体元网格；
[0023]步骤3.2:体元赋值，将点云映射到三维体元网格，进而根据体元中是否包含有激 光点分别为体元赋值1和〇，包含激光点赋值1，无激光点赋值〇，1和〇分别代表目标和背景。
[0024] 所述步骤4的具体步骤如下：
[0025] 步骤4.1:地面种子体元搜寻；
[0026] 步骤4.2:3D滤波，标记地面种子体元的邻域体元且体元值为1的体元，并采用深度 优先策略遍历同属地面的所有体元；
[0027]步骤4.2.1 :初始化，建立空栈；
[0028] 步骤4.2.2:种子体元入栈并标记为地面体元；
[0029]步骤4.2.3:如果栈中为空，停止;否则，进入步骤4.2.4;
[0030]步骤4.2.4:从栈顶弹出一个元素，获取其26邻域内体元值为1且未标记的体元，标 记为地面体元并存入栈中；
[0031] 步骤4.2.5:转至步骤4.2.3。
[0032]本发明针对机载LIDAR点云数据提出一种基于体元数据结构的3D滤波方法。该方 法将离散的点云规则化为二值3D体数据，进而选择局部高程最低作为地面种子体元并标记 其邻域目标体元为地面体元。以ISPRS网站提供的15个实测滤波测试样本数据为实验数据 评定本发明所提出算法的精度，结果表明，算法的平均kappa系数为78.62% (其中，城区和 山区的分别为87.36%和65.51%)，算法精度和已有的滤波算法精度相当，本发明有助于3D 数据表达及处理。
【附图说明】
[0033]图1为本发明【具体实施方式】中基于体元的3D滤波方法流程图；
[0034]图2为本发明【具体实施方式】中步骤2的具体流程图；
[0035]图3为本发明【具体实施方式】中步骤2中的平均点间距计算原理图；
[0036]图4为本发明【具体实施方式】中步骤3的具体流程图；
[0037]图5为本发明【具体实施方式】中步骤3中的空间分辨率计算原理示意图；
[0038]图6为本发明【具体实施方式】中步骤3中的体元坐标系示意图；
[0039]图7为本发明【具体实施方式】中步骤4的具体流程图；
[0040]图8为本发明【具体实施方式】中步骤4中的3维体元阵列分块示意图；
[0041]图9为本发明【具体实施方式】中的Samp31样本滤波结果图，其中（a)为Samp31点云 图，（b)为异常数据剔除后点云图，（c)为二值体元数据结构模型图，（d)为3D滤波结果图。
[0042]图10为本发明【具体实施方式】中应用本发明方法和现有的Alexsson、Pfeifer及 Sohn提出方法的精度对比图。
【具体实施方式】
[0043]下面结合附图对本发明的【具体实施方式】做详细说明。
[0044] 一种基于体元的机载LIDAR点云的三维滤波方法，如图1所示，包括以下步骤：
[0045]步骤1:读取原始LIDAR点云数据；
[0046]本实施方式中，定义LIDAR数据集为P = {pi(xi，zi，zi)，i = l，. . .，n}，其中，i是点 的索引，η是激光点的个数，pi是第i个激光点、其坐标为(Xi，yi，Zi)。
[0047]步骤2:异常数据辨识及剔除；
[0048] 如图2所示，具体步骤如下：
[0049] 步骤2.1:对LIDAR点云建立三维格网索引；
[0050]步骤2.1.1:计算数据集P的轴向包围盒；
[0051 ] 轴向包围盒可由包围盒的左下角坐标（Xmin，ymin，Zmin )和右上角坐标（Xmax，ymax， Zmax)确定，其中，（Xmax，ymax，Zmax)和(Xmin，ymin，Zmin)是P中X、y和Z坐标的最大和最小值，（Xmax， Ymax 7 Zmax) -HlElX {(Xi，yi，Zi)，i - I，···，Il}，( Xmin，ymin，Zmin ) -IIlin{(Xi，yi，Zi)，i - I，···，Il};
[0052] 步骤2.1.2:将轴向包围盒划分为三维格网，其中，格网边长取平均点间距；
[0053]
(1)
[0054] 兵甲，Sxy= Uxi，yi)，i = l，…，η}，见图3，C( ·)是点集Sxy的凸壳，A( ·)是凸壳C的 面积，由此可得M X N X L个立方体格网单元。
[0055]
(2)
[0056] 步骤2.1.3:将各激光点映射到各个格网单元，获取各激光点的格网索引；
[0057] 将各激光点m映射到各个格网单元，
[0058.
(3)
[0059] 其中，（u，v，w)代表pi所在的格网单元的索引。
[0060] 步骤2.2:筛选候选异常数据；
[0061] 若某个激光点所在格网单元的26邻域内激光点个数小于一定的阈值Tn，即为候选 点;循环直至所有激光点处理完毕。
[0062]在本实施方式中，Tn的取值需根据原始点云的空间分布情况确定(比如4)。
[0063]步骤2.3:利用局部距离离群因子（Local Distance Outlier Factor，LD0F)对候 选点进行异常判定；
[0064] pi的k邻域集（不包括p〇Ni由距pi最近的k个点构成，记Ni = {gj(xj，yj，zj)，j = 1，~^}，则
[0065]
[0066] 点间的距离，d(gj，gy)为幻，gy两点间的距离，LDOF描述了激光点与其邻域的偏离程度，数 值越大表示偏离程度越高，反之越低，可预设阈值T_F进行判定。
[0067]在本实施方式中，k为常数，其值大于各异常点群中包含的激光点个数的最大值， Tldof的取值需根据原始点云的空间分布情况确定。
[0068]步骤2.4:删除各个异常数据点；
[0069] 记异常数据剔除后的LIDAR数据集为Q= {qi(xi，yi，zi)，i = l，…，ηι}，其中，i是点 的索引，m是激光点的个数，qi是第i个激光点，其坐标为(Xi，yi，Zi)。
[0070] 步骤3:建立二值体元数据结构模型；
[0071]如图4所示，具体步骤如下：
[0072]步骤3.1:将包含场景的三维空间划分成三维体元网格；
[0073]步骤3.1.1:计算场景所包络的三维空间范围；
[0074]用数据集Q的轴向包围盒表示三维空间范围，轴向包围盒的确定参见步骤2.1.1;
[0075] 步骤3.1.2:计算空间分辨率；
[0076] 空间分辨率依据激光点密度确定，
[0077]
⑶
[0078] 其中，（Δ X，Δ y，Δ z)分别为体元模型沿x、y及z轴的分辨率；Sxy = {(Xi，yi)，i = 1，···，n}，SXZ = {(xi，zi)，i = l，· · ·，n}，Syz = {(yi，zi)，i = l，· · ·，n}，见图5，C( ·)是点集Sxy 的凸壳，A( ·)是凸壳C的面积；
[0079] 步骤3.1.3:对轴向包围盒进行划分得到3D体元阵列；
[0080]设V是3D体元阵列中的体元集合，
[0081] V={vj(rj,cj,lj), j = l,··· ,m}, (6)
[0082] 其中，j是体元的索引，m是体元的个数，Vj是第j个体元的体元值，（rj，cj，lj)是第j 个体元在体元阵列中的坐标(行、列和层号，见图6);
[0083] 步骤3.2:将{qi，i = 1，…，m}映射到V，建立二值3D体元阵列，
[0084] .(7)
[0085]步骤4:基于体元的3D滤波；
[0086] 如图7所示，具体步骤如下：
[0087] 步骤4.1:搜寻地面种子体元；
[0088]在水平方向上，利用设定的格网尺寸对3D体元阵列进行分块，如图8所示，并取各 块内高程最低且体元值为1的体元为种子体元。
[0089] 本实施方式中，格网尺寸需根据地形的复杂性确定，如果地形比较复杂，格网尺寸 取值较小，反之，可适当增大格网尺寸。但是，最小格网尺寸不能小于点云中最大结构（如建 筑物)的尺寸，否则，最大结构内部的点将被错判为地面种子体元；
[0090] 步骤4.2:标记地面种子体元的邻域且体元值为1的体元为地面体元；
[0091]步骤4.2.1:初始化，建立空栈；
[0092] 步骤4.2.2:种子体元入栈并标记为地面体元；
[0093]步骤4.2.3:如果栈中为空，停止;否则，进入步骤4.2.4;
[0094]步骤4.2.4:从栈顶弹出一个元素，获取其26邻域内体元值为1且未标记的体元，标 记为地面体元并存入栈中；
[0095] 步骤4.2.5:转至步骤4.2.3。
[0096] 本发明可以在CPU为Core(TM)2P86002.40GHz、内存4GB、Windows 7旗舰版系统上 使用MATLAB7.11.0软件编程实现仿真。
[0097]本实施方式中设计采用ISPRS第三工作组提供的专门用于滤波算法测试的15个样 本数据作为实验数据以检验算法的稳定性和可行性，这些样本数据包括了不同地理环境中 滤波可能遇到困难的情况，如粗差点的影响、复杂的地物、地物与地面相连、斜坡上或低矮 的植被、地形的不连续等，表1列出了本实施方式中各个样本数据的特性及基本参数。
[0098] 本实施方式中设计采用ISPRS提供的参考数据（15个样本数据被准确的分类为地 面点和非地面点)作为标准数据评定算法精度。
[0099] 表1样本数据的特性及基本参数
[0101] 图9为Samp31样本数据滤波结果图，其中（a)为Samp31点云图，（b)为异常数据剔除 后点云图，（c)为二值体元数据结构模型图，（d)为3D滤波结果图。
[0102] 表2为本实施方式中，以ISPRS提供的参考数据为标准对本发明方法对15个样本数 据3D滤波结果进行的定量评价。可以看出，基于体元的3D滤波方法平均kappa系数为 78 · 62%，其中，城区的平均kappa系数达到了87 · 36%，山区的平均kappa系数相对较低，为 65.51 %，总体精度较高。另外，需注意的是，由于本发明的滤波结果为地面体元，而标准数 据中是将原始点云分成地面点和非地面点，所以必须统计本发明的滤波所得地面体元内包 含的地面点，以便和标准数据做对比进而得出本发明的滤波精度。
[0103] 表2本发明方法滤波结果定量评价
[0105] 图10为本实施方式中应用本发明方法和现有的Alexsson、Pfeifer及Sohn提出的 方法（已发表的算法中精度最高的，且Al exs son所提出的算法被应用到商业软件 TerrasoI id中）的精度对比，可以看出，本发明的滤波精度和Alexsson所提出的算法精度相 当。从而验证了本发明使用基于体元的3D滤波算法的有效性。
[0106] 以上所述，仅为本发明中最基础的【具体实施方式】，但本发明的保护范围并不局限 于此，任何本技术领域人士在本发明所揭露的技术范围内，例如改变体元模型的空间分辨 率等可理解到的替换，都应涵盖在本发明的包含范围之内。因此，本发明的保护范围应该以 权利要求书的保护范围为准。
【主权项】
1. 一种基于体元的机载LIDAR点云三维滤波方法，其特征在于，包括以下步骤： 步骤1:读取原始LIDAR点云数据； 步骤2:异常数据辨识及剔除； 步骤3:二值体元数据结构模型建立； 步骤4:基于体元的3D滤波。2. 根据权利要求1所述的基于体元的机载LIDAR点云三维滤波方法，其特征在于，所述 步骤2的具体步骤如下： 步骤2.1:对LIDAR点云建立三维格网索引； 步骤2.1.1:确定数据集的轴向包围盒； 步骤2.1.2:依据激光点云的密度确定格网尺寸，进而将轴向包围盒划分为三维格网； 步骤2.1.3:将各激光点映射到相应的格网单元，获取其对应的格网索引； 步骤2.2:根据激光点所在的格网单元间的拓扑关系筛选候选点； 步骤2.3:利用局部距离离群因子判定候选点的异常程度； 步骤2.4:删除各个异常数据点，得到异常数据剔除后的数据集。3. 根据权利要求1所述的基于体元的机载LIDAR点云三维滤波方法，其特征在于，所述 步骤3的具体步骤如下： 步骤3.1:体元结构构建，将包含场景目标的三维空间划分成三维体元网格，三维体元 网格用3D体元阵列表示，其依赖于场景目标所包络的三维空间范围和体元大小； 步骤3.1.1:用异常数据剔除后的数据集的轴向包围盒表示三维空间范围； 步骤3.1.2:计算空间分辨率； 步骤3.1.3:对轴向包围盒进行划分得到3D体元网格； 步骤3.2:体元赋值，将点云映射到三维体元网格，进而根据体元中是否包含有激光点 分别为体元赋值，包含激光点赋值1，无激光点赋值0，1和0分别代表目标和背景。4. 根据权利要求1所述的基于体元的机载LIDAR点云三维滤波方法，其特征在于，所述 步骤4的具体步骤如下： 步骤4.1:地面种子体元搜寻； 步骤4.2:3D滤波，标记地面种子体元的邻域体元且体元值为1的体元，并采用深度优先 策略遍历同属地面的所有体元； 步骤4.2.1:初始化，建立空栈； 步骤4.2.2:种子体元入栈并标记为地面体元； 步骤4.2.3:如果栈中为空，停止;否则，进入步骤4.2.4; 步骤4.2.4:从栈顶弹出一个元素，获取其26邻域内体元值为1且未标记的体元，标记为 地面体元并存入栈中； 步骤4.2.5:转至步骤4.2.3。
【文档编号】G06T15/00GK105844602SQ201610201234
【公开日】2016年8月10日
【申请日】2016年4月1日
【发明人】王丽英, 王春艳, 赵泉华
【申请人】辽宁工程技术大学