一种面向对象的机载激光雷达点云滤波方法
一种面向对象的机载激光雷达点云滤波方法
【技术领域】
[0001] 本发明涉及一种机载激光雷达点云滤波方法，具体地说是一种面向对象的机载激 光雷达点云滤波方法，属于激光雷达点云数据信息提取技术领域。
【背景技术】
[0002] 机载激光雷达设备作业中，激光扫描过程是盲目的，即激光脉冲既可能打在地面 上，又可能打在建筑物、桥梁、电力线、灯塔、车辆等人工地物或者植被上。因此，获取的机载 激光雷达点云数据中既有地面点，又有地物点。点云分类是后续应用的基础。工程应用过程 中，机载激光雷达点云数据分类往往采用了分层分类的思路，即首先把点云数据中的地面 点和地物点（非地面点）分离开、并利用地面点生成数字高程模型(DEM);然后，按照归一化 数字表面高程特征(地物点到DEM的高度特征)可以依次将地物点划分到低矮植被、中等高 度植被、高植被等类别。可见，其中，地面点和非地面点的分离是该类分类方法的基础。从机 载激光雷达点云数据中分离出地形表面激光脚点数据子集的过程，称为滤波。目前，机载激 光雷达点云滤波是机载激光雷达点云数据后处理的关键步骤之一，也是研究的热点和难 点。
[0003] 目前，研究人员已经提出了众多的滤波方法。由于对地面点和非地面点的理解概 念或者思路不同，研究人员从而设计出不同的滤波方法。有些方法可直接应用于不规则分 布的点集，有些方法需将原始点集重采样为规则格网，以利用成熟的图像处理技术;从对点 的处理方式上看，有些方法是考察单点与单点，有些是单点与点集或者点集与点集等。有些 滤波方法通过量测某种不连续测度，如高程差、坡度、到结构表面的最小距离或到参数平面 的最小距离等;各种滤波方法基于各种不同的假设前提，如坡度、平面、表面和分割等;有些 算法的处理方式是一次单步完成，有些是通过迭代处理，这也意味着它们分别在准确性与 速度方面的劣势。目前，比较有代表性的滤波方法有:三角网渐进加密、最小二乘分层稳健 线性内插、基于点云分割的方法、基于数学形态学的方法、基于坡度的方法、基于曲面拟合 的方法等。
[0004] 尽管目前存在众多的滤波方法，但缺乏统一、规范的比较。为此，在2003年，国际摄 影测量与遥感协会(ISPRS)第三委员会的机载激光雷达数据处理小组组织了滤波算法的测 试。测试中，提供了8个数据集，这些数据覆盖了城区、农村区域和森林地区，有八种滤波算 法参加 了测试，算法的作者分别为：Elmqvist、Sohn、Axel sson、Pfeifer、Brovel 1 i、 Roggero、Wack/Wimme;rNSithole/VosselmaruSithole和Vosselman在测试报告中，对比了八 种不同滤波算法在不同场景中的滤波结果，对滤波算法做出了比较详细的分析和比较，旨 在研究现存滤波算法的特点，影响滤波结果的因素和存在的问题。测试结果表明参与比较 的滤波算法，对于不太复杂的场景(坡度平缓，只包含小型建筑物，植被稀疏，能采集到大量 地面点的区域)表现普遍较好，但当处理复杂场景时其可靠性明显下降。综合而言，在已有 的方法中，三角网渐进加密滤波方法取得了最好的效果。
[0005] 在已有的机载激光雷达点云滤波方法中，三角网渐进加密滤波方法原理简洁、所 需参数数量少，且对于多数场景取得了较好的滤波效果。但是，三角网渐进加密滤波方法也 有一个显著的弱点，对于丘陵和山地区域，它很容易将地形陡变处的地面点误认为是非地 面点，如图2所示;对于城市区域，它很容易将车辆下部的地物点误认为是地面点。
[0006] 上述问题导致已有的三角网渐进加密滤波方法在对地形、地物复杂区域(包括:存 在大量的地形断裂现象的高山丘陵区域、车辆密集的城区）的机载激光雷达点云滤波时，存 在滤波精度较低、需要大量人工干预的问题，无法满足用户对快速、高精度、自动滤波的迫 切需求。
【发明内容】
[0007] 针对现有滤波方法的不足，本发明提出了一种面向对象的机载激光雷达点云滤波 方法，其利用对象的特征点替代对象包含的原始点进行面向对象的机载激光雷达点云滤 波，能够有效提尚机载激光雷达点75Γ滤波的精度。
[0008] 本发明解决其技术问题采取的技术方案是:一种面向对象的机载激光雷达点云滤 波方法，其特征是，包括以下步骤：
[0009] 步骤一，加载机载激光雷达点云数据；
[0010] 步骤二，识别机载激光雷达点云中的粗差点，并剔除粗差点；
[0011]步骤三，进行基于光滑表面生长的机激光雷达点云分割以获取对象；
[0012] 步骤四，分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其 所包含的激光雷达点；
[0013] 步骤五，提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与 后续运算；
[0014] 步骤六，进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的 类别。
[0015] 进一步地，所述步骤二具体包括以下步骤：
[0016] (1)初始化预处理，对机载激光雷达点云的类别进行标记，将未分类的点标记为 U-, ?? 丄；
[0017] (2)目视分析机载激光雷达点云的高程分布特征，检查是否存在显著的粗差，如果 该点云数据中存在明显的高位粗差和低位粗差，则进入步骤(3);如果该点云数据中没有明 显的高位粗差和低位粗差，则进入步骤(4);
[0018] (3)建立点云的高程直方图识别机载激光雷达点云显著的粗差；
[0019] (4)建立kd-树空间索引，对于类别号为"Γ的激光雷达点集，利用该点集的水平坐 标信息建立二维的kd-树；
[0020] (5)识别低位粗差和高位粗差，利用kd-tree空间索引寻找粗差点。
[0021] 进一步地，所述高程直方图通过机载激光雷达点云的高程z特征来建立。
[0022] 进一步地，所述步骤三具体包括以下步骤：
[0023 ] (1)估计法向量和残差，点云数据中任意一点及其一定数量的k个三维空间中最邻 近点确定的平面的法向量为该点的法向量，三维空间中一定数量的离散点确定的平面由特 征值法确定，通过特征值法中获取的最小特征值记为离散点对应的残差；
[0024] (2)进行区域生长，将激光雷达点云划分为若干的对象，每个对象包含若干的激光 雷达点；
[0025] (3)赋予对象及相关激光雷达点初始的类别，对于任意一个对象，如果其多回波比 例大于50%，则将该对象包含的所有激光雷达点的类别号标记为地物点。
[0026] 进一步地，所述估计法向量和残差的具体过程包括以下步骤：
[0027] (1)对于类别号为"Γ的激光雷达点集，建立其三维的kd_树空间索引；
[0028] (2)对每一个类别号为"Γ的激光雷达点，逐一利用kd_树空间索引求取其k个最临 近点；
[0029] (3)对于每一点及其k个最临近点，利用特征值法求取其拟合平面的方程，即确定 每一点的法向量Φ i及其残差Μ。
[0030] 进一步地，所述区域生长的具体过程包括以下步骤：
[0031] 1)设置法向量间角度差异阈值α，角度差异阈值α的使用规则如下:对于当前要扩 张的种子点，利用kd_树求取其某一邻接点作为待处理的点，设当前种子点的法向量为h， 待处理的邻接点的法向量为如果待处理的邻接点与种子点之间法向量的点积| |aiaj| 小于阈值cosa，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续判别 其他的未处理的点；如果待处理的邻接点与种子点之间法向量的点积I |^七| I大于等于阈 值cosa，则跳过该待处理的邻接点，处理其他的未处理的点；
[0032] 2)设置给定距离阈值r，如果待处理的邻接点到种子点确定的种子平面的距离小 于给定距离阈值r，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续 判别其他的未处理的点；如果待处理的邻接点到种子点确定的种子平面的距离大于等于给 定距离阈值r则跳过该待处理的邻接点，处理其他的未处理的点；
[0033] 3)检查类别号为"Γ的激光雷达点集包含的点的数量，如果数量为0,则到步骤7)， 否则，将全部类别号为"Γ的激光雷达点的状态记为"未被分割"；从"未被分割"的激光雷达 点集中，寻找出残差λ最小的点，以残差λ最小的点为种子点并将该点压入一个种子点的队 列，开始进行区域生长；
[0034] 4)生长的过程中，取种子点队列中第一个尚未被处理过的种子点，利用kd-tree求 取该种子点的k个最临近点；
[0035] 5)对于每一个临近点，逐一进行判别，如果临近点已经被赋予分割号，则不予以处 理;如果临近点已经在种子点的队列中，则不予以处理;分别按照法向量间角度差异和距离 差异进行相似性的判别，如果邻近点同时满足两个相似性的条件，则将该邻近点加入到种 子点队列中；如果邻近点没有同时满足两个相似性的条件则不予以处理；
[0036] 6)回到步骤4)，检查种子点队列中是否有未被处理过的点，如果有，则进入步骤 5);否则，将种子点队列中的所有点标记一个相同的分割号，且该分割号与已有的分割号不 同，然后，返回步骤3);
[0037] 7)结束分割 。
[0038] 进一步地，所述步骤四具体包括以下步骤：
[0039] (1)统计对象的多回波比例，以点云分割获取的对象为基本组织单元，统计每一个 对象包含的激光雷达点的总数量、多回波激光雷达点的数量以计算多回波的比例；
[0040] (2)根据多回波比例特征判别对象及相关激光雷达点的类别，对于任意一个对象， 如果其多回波比例大于50%，则将该对象包含的所有激光雷达点的类别号标记地物点。
[0041 ] 进一步地，所述步骤五具体包括以下步骤：
[0042] (1)提取对象的轮廓点及轮廓线，采用建立不规则三角网提取对象的轮廓点，依次 连接各轮廓点生成对象的轮廓线；
[0043] (2)提取对象的骨架线及其骨架线附近的点；
[0044] (3)确定对象的特征点，将轮廓点和骨架线附近的点作为对象的特征点。
[0045] 进一步地，所述步骤六具体包括以下步骤：
[0046] (1)建立格网索引，对特征点集建立二维分块格网索引结构；
[0047] (2)选择初始种子点，从二维分块格网索引结构的每个格网单元格中选取格网中 高程最低的特征点，同时将特征点所在对象包含的所有特征点作为地面种子点；
[0048] (3)构建初始种子点的三角网，利用所有的地面种子点构建一个三角网，该三角网 代表该区域的初始数字高程模型；
[0049] (4)判别对象的类别，通过迭代特征点来判别对象的类别；
[0050] (5)更新对象包含原始激光雷达点的类别，将对象的类别赋值给其包含的原始激 光雷达点，从而完成机载激光雷达点云滤波工作。
[0051] 进一步地，所述判别对象的类别是一个迭代的过程，判别对象类别过程中涉及三 个参数:最大角度阈值Θ、最大距离阈值d和最大地形角度阈值t;判别对象类别过程包括以 下步骤：
[0052] 1)迭代次数记为0,对于每一个待判别的类别为"Γ的特征点，找到该特征点落入 的三角形，计算该特征点到三角形构成的平面的距离及该特征点到三角形三点的夹角、并 找出三个夹角中的最大角；如果同时满足:距离小于最大距离阈值d、最大角小于最大角度 阈值Θ，则认为该特征点是地面点，将该点的类别号标记为地面点；
[0053] 2)以对象为基本判别单元，统计该对象中地面点的比例。如果该比例大于50%，将 该对象和相关特征点的类别号标记为地面点；否则，将该对象和相关特征点的类别号标记 为 "1"。
[0054] 3)将获得的地面点插入三角网，同时迭代次数自增1;
[0055] 4)重复上述步骤1)至3)，继续执行直至迭代次数达到5后则停止迭代。
[0056]本发明的有益效果如下：本发明在剔除粗差点后，对机载激光雷达点云进行基于 光滑表面生长的点云分割获取对象，根据对象的多回波比例特征识别潜在的地物对象、并 剔除地物对象及其所包含的激光雷达点，提取对象的特征点、并利用特征点替代对象包含 的原始激光雷达点参与后续运算，以对象作为基本单元参与基于三角网渐进加密的判别以 区分地面对象(地面对象包含的点为地面点）和地物对象(地物对象包含的点为地物点），显 著的提高了滤波的精度。
[0057]与现有技术相比，本发明具有以下特点：
[0058] (1)本发明利用点云分割获取对象、利用对象的特征点替代对象包含的原始点进 行面向对象的机载激光雷达点75Γ滤波，有效提尚了机载激光雷达点75Γ滤波的精度。
[0059] (2)本发明创造性的将粗差剔除、点云分割、多回波分析、特征点提取、基于三角网 渐进加密的对象类别判别等五个步骤进行有机结合起来，形成了一套完整的面向对象的机 载激光雷达点云滤波的技术流程，为机载激光雷达点云滤波技术提供了有效的途径。
[0060] (3)相比现有的机载激光雷达点云滤波方法，本发明的滤波精度更高，对山区、城 区等多种地貌类型具有广泛的适用性。
【附图说明】
[0061 ]图1为本发明的方法流程图；
[0062]图2为现有技术采用的三角网渐进加密滤波方法的局限性分析示意图；
[0063]图3(a)为某一区域的机载激光雷达点云数据示意图；
[0064]图3(b)为机载激光雷达点云数据的高程直方图；
[0065] 图4为本发明采用的点云分割方法的技术流程图；
[0066] 图5(a)为某一对象的初始轮廓点判定的示意图；
[0067] 图5(b)为某一对象的初始轮廓线生成过程的示意图；
[0068] 图6(a)为某对象涉及的LiDAR点的示意图；
[0069] 图6(b)为某对象的特征点的示意图；
[0070]图7(a)为某机载LiDAR点云数据的示意图；
[0071]图7(b)为点云分割的结果示意图；
[0072] 图7(c)为初始的地面种子点的示意图；
[0073] 图7(d)为由种子点构建的三角网的示意图；
[0074]图7(e)为角度和距离判别的示意图；
[0075]图7(f)为镜像过程的示意图；
[0076] 图8为本发明的整体数值计算流程图；
[0077] 图9(a)场景1及Samplell和Samplel2的范围的示意图；
[0078] 图9(b)场景2及Sample21至Sample24的范围的示意图；
[0079]图9(c)场景3及Sample31的范围的示意图；
[0080] 图9(d)场景4及Sample41和Sample42的范围的不意图；
[0081 ] 图9(e)场景5及Sample51至Sample54的范围的示意图；
[0082] 图9(f)场景6及Sample61的范围的示意图；
[0083]图9(g)场景7及Sample71的范围的示意图；
[0084]图10(a)为两种滤波方法在15个样本区的Typel误差分布规律示意图；
[0085]图10(b)为两种滤波方法在15个样本区的Typell误差分布规律示意图；
[0086]图10(c)为两种滤波方法在15个样本区的总误差的分布规律示意图。
【具体实施方式】
[0087]为能清楚说明本方案的技术特点，下面通过【具体实施方式】，并结合其附图，对本发 明进行详细阐述。下文的公开提供了许多不同的实施例或例子用来实现本发明的不同结 构。为了简化本发明的公开，下文中对特定例子的部件和设置进行描述。此外，本发明可以 在不同例子中重复参考数字和/或字母。这种重复是为了简化和清楚的目的，其本身不指示 所讨论各种实施例和/或设置之间的关系。应当注意，在附图中所图示的部件不一定按比例 绘制。本发明省略了对公知组件和处理技术及工艺的描述以避免不必要地限制本发明。
[0088]针对机载激光雷达点云滤波的需要，本发明提供了一种面向对象的机载激光雷达 点云滤波方法，利用对象的特征点替代对象包含的原始点进行面向对象的机载激光雷达点 云滤波，提高了滤波的精度。
[0089] 如图1所示，本发明的一种面向对象的机载激光雷达点云滤波方法，它包括以下步 骤：
[0090] 步骤一，加载机载激光雷达点云数据；
[0091 ]步骤二，识别机载激光雷达点云中的粗差点，并剔除粗差点；
[0092] 步骤三，进行基于光滑表面生长的机激光雷达点云分割以获取对象；
[0093] 步骤四，分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其 所包含的激光雷达点；
[0094] 步骤五，提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与 后续运算；
[0095] 步骤六，进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的 类别。
[0096] 进一步地，所述步骤二具体包括以下步骤：
[0097] (1)预处理
[0098]将待处理机载激光雷达点云涉及所有点的类别号标记为"Γ，"Γ代表了未分类的 点。
[0099] (2)目视分析机载激光雷达点云的高程分布特征
[0100]粗差是影响机载激光雷达点云数据滤波效果的关键因素之一，因为大部分滤波算 法在选择初始地面点时常选取局部最低点为地面点。如果该局部最低点为粗差点，则很容 易将与该粗差点邻近的地面点判别为地物点，会严重地影响一定区域范围内的滤波的效 果。因此，粗差剔除是机载激光雷达数据滤波的前提之一。由于机载激光雷达点云数据中的 部分粗差呈现出显著的高于或者低于其他激光雷达点的现象，如图3(a)所示，对于此类显 著粗差可以通过高程直方图（如图3(b)所示)分析进行剔除。在目视分析机载激光雷达点云 的高程分布特征后，如果发现该点云数据中存在明显的高位粗差和低位粗差，直接进入步 骤(3);如果发现该点云数据中没有明显的高位粗差和低位粗差，此步骤可以省略，直接进 入步骤(4)。
[0101] (3)建立点云的高程直方图识别显著的粗差
[0102] 利用机载激光雷达点云的高程z特征，建立高程直方图，如图3(b)所示。通过观察 高程直方图，人工分别确定一个较小的 高程值ζι?和一个较大的高程值Zhigh。将高程值小于 Z1?和大于Zhlgh的激光雷达点的类别号标记为"7"，"7"代表了粗差点。
[0103] (4)建立kd-树空间索引
[0104] 对于类别号为"Γ的激光雷达点集，利用该点集的水平坐标信息建立二维的kd_ 树。
[0105] (5)识别低位粗差和高位粗差
[0106] 输入经验值k和η，其中k为邻近点的个数、η为高差标准差的倍数。对于任意一个类 别号为"Γ的激光雷达点，利用kd-tree空间索引寻找其最邻近的k个激光雷达点;计算该激 光雷达点与k个邻近激光雷达点的高程值的平均数^[email protected];另外，分别计算该激光雷达点与 k个邻近激光雷达点的高程值之差，计算k个高程值差的标准差σζ;如果该激光雷达点的高 程值与平均数之差的绝对值大于ησ ζ，则将该点类别号标记为"7"，"7"代表了粗差点。
[0107] 进一步地，所述步骤三具体包括以下步骤：
[0108] (1)估计法向量和残差
[0109] 该分割方法的流程图如图4所示。点云数据中，任意一点的法向量定义为该点及其 一定数量的k个三维空间中最邻近点确定的平面的法向量。对于类别号为"2"的激光雷达点 集，任何一点的k个最邻近点可以通过三维kd-树这种空间索引的方法求取。同时，三维空间 中一定数量的离散点确定的平面可以由特征值法(即主成分分析法)确定。同时特征值法中 获取的最小特征值，记为该点对应的残差。法向量和残差估计的具体过程包括：
[0110] 1)对于类别号为"Γ的激光雷达点集，建立其三维的kd-树空间索引；
[0111] 2)对每一个类别号为"Γ的激光雷达点，逐一利用kd-树空间索引求取其k个最临 近点。k的取值由用户输入，默认值为20;
[0112] 3)对于每一点及其k个最临近点，利用特征值法求取其拟合平面的方程，即确定每 一点的法向量Φ i及其残差h。
[0113] (2)进行区域生长
[0114] 区域生长过程与传统的基于图像的区域生长极其相似。两者都包含两个必要的因 素，邻接关系的确定和相似性的度量。其中，邻接关系由kd-树确定，相似性的度量包括邻接 点到种子平面的距离和邻接点与种子点的法向量角度差两个特征。区域生长的结果是，将 激光雷达点云划分为若干的对象，每个对象包含若干的激光雷达点。该区域生长的具体过 程包括：
[0115] 1)输入法向量间角度差异阈值α。角度差异阈值α的使用规则如下:对于当前要扩 张的种子点，利用kd-树求取其某一邻接点作为待处理的点。设当前种子点的法向量为h， 待处理的邻接点的法向量为七。如果待处理的邻接点与种子点之间法向量的点积| |aiaj| 小于阈值cosa，则将此点加入到种子点队列中，并作为新的种子点，继续判别其他的未处理 的点；否则，跳过该点，处理其他的未处理的点。
[0116] 2)输入给定距离阈值r。距离阈值r如果待处理的邻接点到种子点确定的种子平面 的距离小于阈值r，则将此点加入到种子点队列中，并作为新的种子点，继续判别其他的未 处理的点；否则，跳过该点，处理其他的未处理的点。
[0117] 3)检查类别号为"Γ的激光雷达点集包含的点的数量。如果数量为0,则到步骤7)。 否则，将全部类别号为"Γ的激光雷达点的状态记为"未被分割"。从"未被分割"的激光雷达 点集中，寻找出残差λ最小的点。以残差λ最小的点为种子点，将该点压入一个种子点的队 列，开始进行区域生长。
[0118] 4)生长的过程中，取种子点队列中第一个尚未被处理过的种子点，利用kd-tree求 取该种子点的k个最临近点。
[0119] 5)对于每一个临近点，逐一进行判别。如果该点已经被赋予分割号，则不予以处 理。进一步，如果该点已经在种子点的队列中，则不予以处理。接着，分别按照法向量间角度 差异和距离差异进行相似性的判别。如果该邻近点同时满足两个相似性的条件，则将该邻 近点加入到种子点队列中；否则，不予以处理。
[0120] 6)回到步骤4)，检查种子点队列中是否有未被处理过的点。如果有，则进入步骤 5);否则，将种子点队列中的所有点标记一个相同的分割号，且该分割号与已有的分割号不 同。然后，返回步骤3)。
[0121] 7)结束分割。
[0122] (3)赋予对象及相关激光雷达点初始的类别
[0123] 对于任意一个对象，检查其该对象包含的所有激光雷达点的数量。如果该对象包 含的所有激光雷达点的数量小于q(q为经验阈值），则将该对象包含的所有激光雷达点的类 别号标记为"3"，"3"代表了地物点。
[0124] 进一步地，所述步骤四具体包括以下步骤：
[0125] (1)统计对象的多回波比例
[0126] 以点云分割获取的对象为基本组织单元，统计每一个对象包含的激光雷达点的总 数量、多回波激光雷达点的数量以计算多回波的比例。目前，对于多回波机载激光雷达系统 而言，系统记录的回波信息包括单次回波（singular return)和多次回波（multiple returns)，二者的区别在于对同一束激光脉冲是否发生多次反射。即，当激光扫描仪发射的 激光脉冲接触到被测目标时，部分脉冲能量的反射信号会被系统接收并记录，而剩余的脉 冲能量继续传播，当遇到另一目标或原被测目标的另一部分时再次发生反射，直至能量消 耗殆尽。如果激光束只碰到一次目标且回波被接受纪录下来，则称此回波为单次回波;如果 发生多次反射、并使得机载激光扫描系统接收到多个反射信号即多次回波。现有的机载激 光雷达系统都具备记录2-5次(超过5次以上回波的概率几乎为零)多回波信息的能力。对于 多次回波，接收到的第一个回波信号被称为首次回波(f irst return/pu 1 se);接收到的最 后一个回波信号被称为尾次(末次）回波（last pulse);除首次和尾次回波外的中间的回波 信号被顺序编号，称为中间次回波（intermediate pulses)。研究发现，回波信息可以反映 被测目标的类型信息，继而可以辅助激光点云数据的滤波和分类。对于森林地区：单次回波 包含地面激光脚点和植被冠层的激光脚点；首次回波来自茂密且高大的植被冠层或靠近冠 层的枝叶；中间次回波多为高大植被的枝叶或低矮植被;而末次回波多是植被中间层次的 枝叶和地表反射得到的激光脚点。对于城区：单次回波数据主要来源于地表、人工建筑物 (包括立交桥和横跨河流的桥梁）的顶面或墙面、少量植被点；首次回波来源于植被的冠层 和人工建筑物(包括立交桥）的边缘；中间次回波主要来源于植被的枝叶和建筑物的立面； 而末次回波则主要来源于地表，也有部分是来源于复杂的建筑物屋顶和植被低矮层面的枝 叶。上述分析表明，无论对于森林地区还是城市区域，尾次回波和单次回波既可以来自地 面、又可以来自非地面目标;首次回波和中间次回波肯定来源于非地面目标，如植被、建筑 物边缘等。
[0127] (2)判别对象及相关激光雷达点的类别
[0128] 对于任意一个对象，如果其多回波比例大于50%，则将该对象包含的所有激光雷 达点的类别号标记为"3"，"3"代表了地物点。
[0129] 进一步地，所述步骤五具体包括以下步骤：
[0130] (1)提取对象的轮廓点及轮廓线
[0131]对于任意一个类别号标记为"Γ的对象，找出其包含的该对象包含的所有激光雷 达点，并建立这些点的不规则三角网（TIN)，如图5(a)所示。判断轮廓边缘只需分析各三角 形的相邻三角形的个数。通常一个三角形有3个三角形与之相邻，如果个数小于3,则该三角 形位于对象轮廓的边缘，提取相应轮廓点即可获得对象的初始轮廓点。如图5(b)所示，顺时 针依次搜索各初始轮廓点进一步生成对象的轮廓线，步骤如下：
[0132] 1)起始轮廓点确定:搜索所有初始轮廓点中x值最小的激光雷达点，如存在多个x 值相等的点，则取这些X值相等的点中y值最小的点作为起始轮廓点，如图5(b)中"点Γ所 示；
[0133] 2)第2个轮廓点确定：以"点Γ为起点、R为搜索半径，水平向左为起始搜索方向，沿 顺时针方向寻找下一个轮廓点，搜索半径R由点云平均间隔开始逐渐增大，直到在搜索范围 内找到第2个轮廓点（"点2"）；
[0134] 3)顺时针方向搜索下一个轮廓点：以"点2"为起点、R为搜索半径，从"点1"、"点2" 矢量方向的逆时针90°范围开始顺时针搜索下一个轮廓点，搜索半径逐渐增大（图5 (b)中搜 索半径由R1增大到R2)直到找到第3个轮廓点（"点3"）；
[0135] 4)重复类似的步骤3)，直到找到"点Γ则终止搜索，依次连接各轮廓点生成对象的 轮廓线。
[0136] (2)提取对象的骨架线及其骨架线附近的点
[0137] 基于上一步提取的对象轮廓线形成的多边形，提取该多边形的骨架线。在水平面 内，将该对象中距离骨架线在两倍点云平均间隔内、且不属于轮廓点的点标记为该对象的 "骨架线附近的点"。
[0138] (3)确定对象的特征点
[0139] 将对象的轮廓点和"骨架线附近的点"作为该对象的"特征点"。对象的"特征点"将 代替该对象参与后续运算。图6(a)展示了某对象涉及的激光雷达点云，同时 图6(b)展示了 该对象的特征点。
[0140] 进一步地，所述步骤六具体包括以下步骤：
[0141] (1)建立格网索引
[0142] 对于类别号为"Γ的特征点集，首先根据点云数据在x、y方向的坐标范围，构建一 个特征点集的最小包围核。根据给定格网尺寸，即最大建筑物尺寸m，分别对特征点集在X、y 方向划分正方形格网，对特征点集建立二维分块格网索引结构。
[0143] (2)选择初始种子点
[0144] 从每个格网单元格中选取格网中高程最低的特征点，同时将特征点所在对象包含 的所有特征点作为地面种子点、并将该对象和相关特征点的类别号标记为"2"，"2"代表了 地面点。图7(a)展示了某机载激光雷达点云，图7(b)展示了该激光雷达点云的分割效果，同 时图7(c)展示了初始的种子点。
[0145] (3)构建初始种子点的三角网
[0146] 利用所有的地面种子点构建一个三角网，该三角网代表了该区域的初始数字高程 模型DEM。图7(c)和图7(d)分别展示了某机载激光雷达点云数据的初始种子点和这些种子 点构建的三角网。
[0147] (4)判别对象的类别
[0148] 对象类别的判别是一个迭代的过程，该过程涉及三个参数:最大角度阈值Θ、最大 距离阈值d和最大地形角度阈值t。
[0149] 1)迭代次数记为0。对于每一个待判别的类别为"Γ的特征点，找到该点落入的三 角形，计算该点到三角形构成的平面的距离及其该点到三角形三点的夹角、并找出三个夹 角中的最大角，如图7(e)所示。
[0150] 如果同时满足:距离小于最大距离阈值d、最大角小于最大角度阈值Θ，则认为该点 是地面点，将该点的类别号标记为"2"，"2"代表了地面点；否则，还要进一步判别(对当前的 特征点以距离最近的三角形顶点为中心进行镜像，如图7(f)所示。对镜像点进行判别。如果 此镜像点判断为地面点，则仍然判断此点为地面点、将该点的类别号标记为"2" ；否则，不修 改点的类别值）。继续处理下一个类别为"Γ的特征点直至将所有的类别为"Γ的特征点遍 历完毕。
[0151] 2)以对象为基本判别单元，统计该对象中"地面点"的比例。如果该比例大于50%， 将该对象和相关特征点的类别号标记为"2"，"2"代表了地面点；否则，将该对象和相关特征 点的类别号标记为"1"。
[0152] 3)将获得的地面点插入三角网。迭代次数自增1。回到步骤1)继续执行直至迭代次 数达到的5则停止迭代。
[0153] (5)更新对象包含原始激光雷达点的类别
[0154] 将对象的类别赋值给其包含的原始激光雷达点，完成滤波。
[0155] 本发明创造性地提出了一种面向对象的机载激光雷达点云滤波方法，其特色是利 用对象的特征点替代对象包含的原始点进行面向对象的机载激光雷达点云滤波，提高了滤 波的精度。该方法创造性的将粗差剔除、点云分割、多回波分析、特征点提取、基于三角网渐 进加密的对象类别判别等五个步骤连接起来，形成了一套完整的面向对象的机载激光雷达 点云滤波的技术流程。
[0156] 如图8所示，本发明进行面向对象的机载激光雷达点云滤波的数值计算流程如下：
[0157] (I)机载数据。加载机载激光雷达点云数据。
[0158] (II)粗差剔除。在进行滤波之前需要对激光点云数据进行预处理，去除极低局外 点。
[0159] (III)点云分割。对粗差剔除后的激光雷达点云，进行基于光滑表面生长的分割以 获取对象。
[0160] (IV)回波分析。分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对 象及其所包含的激光雷达点。
[0161] (V)提取对象的特征点。以每一个对象为处理单元，取出其该对象涉及的激光雷达 点，并获取每一个对象的边缘点和骨架线附近的点；利用边缘点和骨架线附近的点组成该 对象的特征点集;该特征点集将代替该对象及其包含的原始激光雷达点，参与后续运算。
[0162] (VI)判别。进行格网划分建立二维分块格网索引结构;选择初始种子点；构建初始 种子点的不规则三角网；迭代的通过特征点的判别确定对象的类别，并更新对象包含的原 始激光雷达点的类别。
[0163] 为了验证上述提出的滤波方法的性能，本研究同时采用经典的三角网渐进加密方 法和本发明提出的方法，以ISPRS第三委员会机载激光雷达数据处理小组提供的专门用于 滤波算法测试数据作为实验数据，通过定性和定量分析的手段比较两种算法的性能。
[0164] ISPRS第III小组于2003年组织的滤波算法测试中，提供了一组数据集和参考数 据，用于定量的分析算法的滤波效果（数据来自：111^口 ：//?^^.(301]11]1188；[0113.18口18.0找/ wg3/index.html)。该数据集包含7个场景的点云数据，依次命名为：CStitel、CStite2、 05衍七63、05衍七64、05衍七65、05衍七66、05衍七67，如图9(3)至图9(8)所示。这7个场景的数据 为OptechALTM机载激光扫描仪获取的Vaihingen/Enz和Stuttgart市的机载激光雷达点云 数据。其中的城区数据为Sitel-Site4,点间距为1.0至1.5m;山区为Site5-Site7,点间距为 2.0至3.5m。同时，在7块测试数据中选取了 15个样本数据(具体数据区位如图9(a)至图9(g) 所示），这些样本数据经过人工处理，将激光脚点精确分类为地面点和非地面点这两类。15 个样本数据主要包括了不同地理环境中滤波可能遇到困难的情况，如粗差点的影响、复杂 的地物、地物与地面相连、斜坡上或低矮的植被、地形的不连续等，具体数据特征如表1所 示。利用本发明提出的滤波算法对这7个场景的数据做滤波实验，利用15个参考样本对滤波 结果进行了定性和定量分析，并与已有的三角网渐进加密滤波算法的结果进行了对比分 析。
[0165] ISPRS第三委员会机载激光雷达数据处理小组提出了滤波算法质量评价的两种方 法:定性分析和定量分析。定性分析中，把对滤波后的结果叠加到原始数据上，通过误差图 进行可视化检查、分析和评价。定量分析中，统计了三种类型的误差，即Typel(地面点被误 分为地物点的数量占地面点数量的比例）、Type Π (地物点被误分为地面点的数量占地物点 数量的比例)和总误差(被误分的地面点和和地物点的总数据量占全部脚点数量的比例）， 见表3。
[0166] 表1参考数据的说明
[0167]
[0168] 表2定性分析的等级
[0169]
[0170]~表3定量分析的相应参数及其定义
[0171]
[0172] 注释:表中，a，b，c，d是进行定量分析的基本参数，a代表参考点中是地面点并且滤 波算法也将其分类为地面点的地面点数量，b代表参考点中是地面点但滤波算法将其分类 为地物点的地面点数量，c代表参考点中是地物点而滤波算法将其分类为地面点的地物点 数量，d代表参考点中是地物点且滤波算法也将其分类为地物点的地物点数量;f，g，h，i分 析数据滤波前后的地面点和地物点组成情况，f是参考数据滤波前的地面点的比例，g是参 考数据滤波前的地物点的比例，h是滤波结果中地面点的比例，i是滤波结果中地物点的比 例，；Typel、Type Π 和Total是三类误差;k反映了两类误差的偏离程度。本论文中，j = 0。
[0173] 表4两种滤波方法参数的取值
[0174]
[0176]已有的三角网渐进加密滤波方法需要最大建筑物尺寸、最陡地形倾角、迭代角度、 迭代距离等四个参数。本发明提出的滤波方法，在已有方法所需四个参数基础上增加了三 个参数，即邻接点个数、距离阈值、角度差异阈值。每个场景的实验中，两种方法共有参数的 取值是相同的，本发明所需的额外三个参数的取值需要单独确定。两种方法在7个ISPRS提 供的实验数据的滤波参数的取值情况，见表4。
[0177] 使用表4中列举的参数，利用两种滤波方法对上述7个实验数据进行滤波。同时，从 ISPRS提供的7个场景的滤波结果中，分别选择15个参考样本数据具有相同坐标值的激光脚 点，并利用参考样本的滤波结果评价15个相关区域的滤波效果，得到了三类误差的统计表 5。另外，图10(a)、（b)、（c)分别展示了两种滤波方法的三种类型的误差对比图。
[0178] 表5经典的三角网滤波方法与本发明改进方法的误差对比表(单位：％ )
[0179]
[0180]
[0181]从表5可以看出，已有的三角网渐进加密滤波方法和本发明提出的滤波方法（即表 5以及图10(a)、图10(b)和图10(c)中所述的改进方法)都具有较强的鲁棒性，能广泛的适应 于城区和森林区域各种场景的滤波，且滤波的总误差控制在33%以下，如图10(c)所示。同 时，与经典的滤波方法相比，本发明的滤波方法具有更好的性能，能正确的识别更多的地面 点，尤其是地形陡变处的地面点，其I类误差和总误差显著的低于经典滤波的相应误差。15 个局部区域的滤波效果检验中，除samp21区域之外，其他14个区域中本发明滤波方法的I类 误差显著的低于已 有三角网渐进加密滤波方法的I类误差；除samp21、samp54两个区域之 外，其他13个区域中本发明滤波方法的总误差显著的低于已有三角网渐进加密滤波方法的 总误差，如图10(a)和（c)所示。但是，本发明滤波方法也有一定的局限性，与已有三角网渐 进加密滤波方法相比，它具有更高的II误差。15个局部区域的滤波效果检验中，除samp61之 外，其他14个区域中本发明滤波方法的II类误差高于已有三角网渐进加密滤波方法的II类 误差，如图10(b)所示。整体而言，表5中的统计数据表明，本发明滤波方法比已有三角网渐 进加密滤波方法具有更高的滤波精度。
[0182]以上所述只是本发明的优选实施方式，对于本技术领域的普通技术人员来说，在 不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也被视为本发 明的保护范围。
【主权项】
1. 一种面向对象的机载激光雷达点云滤波方法，其特征是，包括以下步骤： 步骤一，加载机载激光雷达点云数据； 步骤二，识别机载激光雷达点云中的粗差点，并剔除粗差点； 步骤三，进行基于光滑表面生长的机激光雷达点云分割以获取对象； 步骤四，分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其所包 含的激光雷达点； 步骤五，提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与后续 运算； 步骤六，进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的类别。2. 根据权利要求1所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤二具体包括以下步骤： (1) 初始化预处理，对机载激光雷达点云的类别进行标记，将未分类的点标记为"Γ; (2) 目视分析机载激光雷达点云的高程分布特征，检查是否存在显著的粗差，如果该点 云数据中存在明显的高位粗差和低位粗差，则进入步骤(3);如果该点云数据中没有明显的 尚位粗差和低位粗差，则进入步骤(4); (3) 建立点云的高程直方图识别机载激光雷达点云显著的粗差； (4) 建立kd-树空间索引，对于类别号为"Γ的激光雷达点集，利用该点集的水平坐标信 息建立二维的kd-树； (5) 识别低位粗差和高位粗差，利用kd-tree空间索引寻找粗差点。3. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 高程直方图通过机载激光雷达点云的高程z特征来建立。4. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤三具体包括以下步骤： (1) 估计法向量和残差，点云数据中任意一点及其一定数量的k个三维空间中最邻近点 确定的平面的法向量为该点的法向量，三维空间中一定数量的离散点确定的平面由特征值 法确定，通过特征值法中获取的最小特征值记为离散点对应的残差； (2) 进行区域生长，将激光雷达点云划分为若干的对象，每个对象包含若干的激光雷达 占. (3) 赋予对象及相关激光雷达点初始的类别，对于任意一个对象，如果其多回波比例大 于50%，则将该对象包含的所有激光雷达点的类别号标记为地物点。5. 根据权利要求4所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 估计法向量和残差的具体过程包括以下步骤： (1) 对于类别号为"Γ的激光雷达点集，建立其三维的kd-树空间索引； (2) 对每一个类别号为"Γ的激光雷达点，逐一利用kd-树空间索引求取其k个最临近 占. (3) 对于每一点及其k个最临近点，利用特征值法求取其拟合平面的方程，即确定每一 点的法向量Φi及其残差h。6. 根据权利要求4所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 区域生长的具体过程包括以下步骤： 1) 设置法向量间角度差异阈值α，角度差异阈值α的使用规则如下：对于当前要扩张的 种子点，利用kd_树求取其某一邻接点作为待处理的点，设当前种子点的法向量为h，待处 理的邻接点的法向量为aj;如果待处理的邻接点与种子点之间法向量的点积| |aiC[j| |小于 阈值cosa，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续判别其 他的未处理的点；如果待处理的邻接点与种子点之间法向量的点积I|^七|I大于等于阈值 cosa，则跳过该待处理的邻接点，处理其他的未处理的点； 2) 设置给定距离阈值r，如果待处理的邻接点到种子点确定的种子平面的距离小于给 定距离阈值r，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续判别 其他的未处理的点；如果待处理的邻接点到种子点确定的种子平面的距离大于等于给定距 离阈值r则跳过该待处理的邻接点，处理其他的未处理的点； 3) 检查类别号为"Γ的激光雷达点集包含的点的数量，如果数量为0,则到步骤7)，否 贝1J，将全部类别号为"Γ的激光雷达点的状态记为"未被分割"；从"未被分割"的激光雷达点 集中，寻找出残差λ最小的点，以残差λ最小的点为种子点并将该点压入一个种子点的队列， 开始进行区域生长； 4) 生长的过程中，取种子点队列中第一个尚未被处理过的种子点，利用kd-tree求取该 种子点的k个最临近点； 5) 对于每一个临近点，逐一进行判别，如果临近点已经被赋予分割号，则不予以处理； 如果临近点已经在种子点的队列中，则不予以处理;分别按照法向量间角度差异和距离差 异进行相似性的判别，如果邻近点同时满足两个相似性的条件，则将该邻近点加入到种子 点队列中；如果邻近点没有同时满足两个相似性的条件则不予以处理； 6) 回到步骤4)，检查种子点队列中是否有未被处理过的点，如果有，则进入步骤5);否 贝1J，将种子点队列中的所有点标记一个相同的分割号，且该分割号与已有的分割号不同，然 后，返回步骤3); 7) 结束分割。7. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤四具体包括以下步骤： (1) 统计对象的多回波比例，以点云分割获取的对象为基本组织单元，统计每一个对象 包含的激光雷达点的总数量、多回波激光雷达点的数量以计算多回波的比例； (2) 根据多回波比例特征判别对象及相关激光雷达点的类别，对于任意一个对象，如果 其多回波比例大于50%，则将该对象包含的所有激光雷达点的类别号标记地物点。8. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤五具体包括以下步骤： (1) 提取对象的轮廓点及轮廓线，采用建立不规则三角网提取对象的轮廓点，依次连接 各轮廓点生成对象的轮廓线； (2) 提取对象的骨架线及其骨架线附近的点； (3) 确定对象的特征点，将轮廓点和骨架线附近的点作为对象的特征点。9. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤六具体包括以下步骤： (1)建立格网索引，对特征点集建立二维分块格网索引结构； (2) 选择初始种子点，从二维分块格网索引结构的每个格网单元格中选取格网中高程 最低的特征点，同时将特征点所在对象包含的所有特征点作为地面种子点； (3) 构建初始种子点的三角网，利用所有的地面种子点构建一个三角网，该三角网代表 该区域的初始数字高程模型； (4) 判别对象的类别，通过迭代特征点来判别对象的类别； (5) 更新对象包含原始激光雷达点的类别，将对象的类别赋值给其包含的原始激光雷 达点，从而完成机载激光雷达点云滤波工作。10.根据权利要求9所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所 述判别对象的类别是一个迭代的过程，判别对象类别过程中涉及三个参数:最大角度阈值 Θ、最大距离阈值d和最大地形角度阈值t;判别对象类别过程包括以下步骤： 1) 迭代次数记为0,对于每一个待判别的类别为"Γ的特征点，找到该特征点落入的三 角形，计算该特征点到三角形构成的平面的距离及该特征点到三角形三点的夹角、并找出 三个夹角中的最大角；如果同时满足:距离小于最大距离阈值d、最大角小于最大角度阈值 Θ，则认为该特征点是地面点，将该点的类别号标记为地面点； 2) 以对象为基本判别单元，统计该对象中地面点的比例。如果该比例大于50%，将该对 象和相关特征点的类别号标记为地面点；否则，将该对象和相关特征点的类别号标记为 T。 3) 将获得的地面点插入三角网，同时迭代次数自增1; 4) 重复上述步骤1)至3)，继续执行直至迭代次数达到5后则停止迭代。
【专利摘要】本发明公开了一种面向对象的机载激光雷达点云滤波方法，它包括以下步骤：加载机载激光雷达点云数据；识别机载激光雷达点云中的粗差点，并剔除粗差点；进行基于光滑表面生长的机激光雷达点云分割以获取对象；分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其所包含的激光雷达点；提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与后续运算；进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的类别。本发明形成了一套完整的面向对象的机载激光雷达点云滤波的技术流程，有效提高了机载激光雷达点云滤波的精度，对山区、城区等多种地貌类型具有广泛的适用性。
【IPC分类】G06T5/00
【公开号】CN105488770
【申请号】CN201510920332
【发明人】林祥国, 段敏燕, 张继贤
【申请人】中国测绘科学研究院
【公开日】2016年4月13日
【申请日】2015年12月11日
【技术领域】
[0001] 本发明涉及一种机载激光雷达点云滤波方法，具体地说是一种面向对象的机载激 光雷达点云滤波方法，属于激光雷达点云数据信息提取技术领域。
【背景技术】
[0002] 机载激光雷达设备作业中，激光扫描过程是盲目的，即激光脉冲既可能打在地面 上，又可能打在建筑物、桥梁、电力线、灯塔、车辆等人工地物或者植被上。因此，获取的机载 激光雷达点云数据中既有地面点，又有地物点。点云分类是后续应用的基础。工程应用过程 中，机载激光雷达点云数据分类往往采用了分层分类的思路，即首先把点云数据中的地面 点和地物点（非地面点）分离开、并利用地面点生成数字高程模型(DEM);然后，按照归一化 数字表面高程特征(地物点到DEM的高度特征)可以依次将地物点划分到低矮植被、中等高 度植被、高植被等类别。可见，其中，地面点和非地面点的分离是该类分类方法的基础。从机 载激光雷达点云数据中分离出地形表面激光脚点数据子集的过程，称为滤波。目前，机载激 光雷达点云滤波是机载激光雷达点云数据后处理的关键步骤之一，也是研究的热点和难 点。
[0003] 目前，研究人员已经提出了众多的滤波方法。由于对地面点和非地面点的理解概 念或者思路不同，研究人员从而设计出不同的滤波方法。有些方法可直接应用于不规则分 布的点集，有些方法需将原始点集重采样为规则格网，以利用成熟的图像处理技术;从对点 的处理方式上看，有些方法是考察单点与单点，有些是单点与点集或者点集与点集等。有些 滤波方法通过量测某种不连续测度，如高程差、坡度、到结构表面的最小距离或到参数平面 的最小距离等;各种滤波方法基于各种不同的假设前提，如坡度、平面、表面和分割等;有些 算法的处理方式是一次单步完成，有些是通过迭代处理，这也意味着它们分别在准确性与 速度方面的劣势。目前，比较有代表性的滤波方法有:三角网渐进加密、最小二乘分层稳健 线性内插、基于点云分割的方法、基于数学形态学的方法、基于坡度的方法、基于曲面拟合 的方法等。
[0004] 尽管目前存在众多的滤波方法，但缺乏统一、规范的比较。为此，在2003年，国际摄 影测量与遥感协会(ISPRS)第三委员会的机载激光雷达数据处理小组组织了滤波算法的测 试。测试中，提供了8个数据集，这些数据覆盖了城区、农村区域和森林地区，有八种滤波算 法参加 了测试，算法的作者分别为：Elmqvist、Sohn、Axel sson、Pfeifer、Brovel 1 i、 Roggero、Wack/Wimme;rNSithole/VosselmaruSithole和Vosselman在测试报告中，对比了八 种不同滤波算法在不同场景中的滤波结果，对滤波算法做出了比较详细的分析和比较，旨 在研究现存滤波算法的特点，影响滤波结果的因素和存在的问题。测试结果表明参与比较 的滤波算法，对于不太复杂的场景(坡度平缓，只包含小型建筑物，植被稀疏，能采集到大量 地面点的区域)表现普遍较好，但当处理复杂场景时其可靠性明显下降。综合而言，在已有 的方法中，三角网渐进加密滤波方法取得了最好的效果。
[0005] 在已有的机载激光雷达点云滤波方法中，三角网渐进加密滤波方法原理简洁、所 需参数数量少，且对于多数场景取得了较好的滤波效果。但是，三角网渐进加密滤波方法也 有一个显著的弱点，对于丘陵和山地区域，它很容易将地形陡变处的地面点误认为是非地 面点，如图2所示;对于城市区域，它很容易将车辆下部的地物点误认为是地面点。
[0006] 上述问题导致已有的三角网渐进加密滤波方法在对地形、地物复杂区域(包括:存 在大量的地形断裂现象的高山丘陵区域、车辆密集的城区）的机载激光雷达点云滤波时，存 在滤波精度较低、需要大量人工干预的问题，无法满足用户对快速、高精度、自动滤波的迫 切需求。
【发明内容】
[0007] 针对现有滤波方法的不足，本发明提出了一种面向对象的机载激光雷达点云滤波 方法，其利用对象的特征点替代对象包含的原始点进行面向对象的机载激光雷达点云滤 波，能够有效提尚机载激光雷达点75Γ滤波的精度。
[0008] 本发明解决其技术问题采取的技术方案是:一种面向对象的机载激光雷达点云滤 波方法，其特征是，包括以下步骤：
[0009] 步骤一，加载机载激光雷达点云数据；
[0010] 步骤二，识别机载激光雷达点云中的粗差点，并剔除粗差点；
[0011]步骤三，进行基于光滑表面生长的机激光雷达点云分割以获取对象；
[0012] 步骤四，分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其 所包含的激光雷达点；
[0013] 步骤五，提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与 后续运算；
[0014] 步骤六，进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的 类别。
[0015] 进一步地，所述步骤二具体包括以下步骤：
[0016] (1)初始化预处理，对机载激光雷达点云的类别进行标记，将未分类的点标记为 U-, ?? 丄；
[0017] (2)目视分析机载激光雷达点云的高程分布特征，检查是否存在显著的粗差，如果 该点云数据中存在明显的高位粗差和低位粗差，则进入步骤(3);如果该点云数据中没有明 显的高位粗差和低位粗差，则进入步骤(4);
[0018] (3)建立点云的高程直方图识别机载激光雷达点云显著的粗差；
[0019] (4)建立kd-树空间索引，对于类别号为"Γ的激光雷达点集，利用该点集的水平坐 标信息建立二维的kd-树；
[0020] (5)识别低位粗差和高位粗差，利用kd-tree空间索引寻找粗差点。
[0021] 进一步地，所述高程直方图通过机载激光雷达点云的高程z特征来建立。
[0022] 进一步地，所述步骤三具体包括以下步骤：
[0023 ] (1)估计法向量和残差，点云数据中任意一点及其一定数量的k个三维空间中最邻 近点确定的平面的法向量为该点的法向量，三维空间中一定数量的离散点确定的平面由特 征值法确定，通过特征值法中获取的最小特征值记为离散点对应的残差；
[0024] (2)进行区域生长，将激光雷达点云划分为若干的对象，每个对象包含若干的激光 雷达点；
[0025] (3)赋予对象及相关激光雷达点初始的类别，对于任意一个对象，如果其多回波比 例大于50%，则将该对象包含的所有激光雷达点的类别号标记为地物点。
[0026] 进一步地，所述估计法向量和残差的具体过程包括以下步骤：
[0027] (1)对于类别号为"Γ的激光雷达点集，建立其三维的kd_树空间索引；
[0028] (2)对每一个类别号为"Γ的激光雷达点，逐一利用kd_树空间索引求取其k个最临 近点；
[0029] (3)对于每一点及其k个最临近点，利用特征值法求取其拟合平面的方程，即确定 每一点的法向量Φ i及其残差Μ。
[0030] 进一步地，所述区域生长的具体过程包括以下步骤：
[0031] 1)设置法向量间角度差异阈值α，角度差异阈值α的使用规则如下:对于当前要扩 张的种子点，利用kd_树求取其某一邻接点作为待处理的点，设当前种子点的法向量为h， 待处理的邻接点的法向量为如果待处理的邻接点与种子点之间法向量的点积| |aiaj| 小于阈值cosa，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续判别 其他的未处理的点；如果待处理的邻接点与种子点之间法向量的点积I |^七| I大于等于阈 值cosa，则跳过该待处理的邻接点，处理其他的未处理的点；
[0032] 2)设置给定距离阈值r，如果待处理的邻接点到种子点确定的种子平面的距离小 于给定距离阈值r，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续 判别其他的未处理的点；如果待处理的邻接点到种子点确定的种子平面的距离大于等于给 定距离阈值r则跳过该待处理的邻接点，处理其他的未处理的点；
[0033] 3)检查类别号为"Γ的激光雷达点集包含的点的数量，如果数量为0,则到步骤7)， 否则，将全部类别号为"Γ的激光雷达点的状态记为"未被分割"；从"未被分割"的激光雷达 点集中，寻找出残差λ最小的点，以残差λ最小的点为种子点并将该点压入一个种子点的队 列，开始进行区域生长；
[0034] 4)生长的过程中，取种子点队列中第一个尚未被处理过的种子点，利用kd-tree求 取该种子点的k个最临近点；
[0035] 5)对于每一个临近点，逐一进行判别，如果临近点已经被赋予分割号，则不予以处 理;如果临近点已经在种子点的队列中，则不予以处理;分别按照法向量间角度差异和距离 差异进行相似性的判别，如果邻近点同时满足两个相似性的条件，则将该邻近点加入到种 子点队列中；如果邻近点没有同时满足两个相似性的条件则不予以处理；
[0036] 6)回到步骤4)，检查种子点队列中是否有未被处理过的点，如果有，则进入步骤 5);否则，将种子点队列中的所有点标记一个相同的分割号，且该分割号与已有的分割号不 同，然后，返回步骤3);
[0037] 7)结束分割 。
[0038] 进一步地，所述步骤四具体包括以下步骤：
[0039] (1)统计对象的多回波比例，以点云分割获取的对象为基本组织单元，统计每一个 对象包含的激光雷达点的总数量、多回波激光雷达点的数量以计算多回波的比例；
[0040] (2)根据多回波比例特征判别对象及相关激光雷达点的类别，对于任意一个对象， 如果其多回波比例大于50%，则将该对象包含的所有激光雷达点的类别号标记地物点。
[0041 ] 进一步地，所述步骤五具体包括以下步骤：
[0042] (1)提取对象的轮廓点及轮廓线，采用建立不规则三角网提取对象的轮廓点，依次 连接各轮廓点生成对象的轮廓线；
[0043] (2)提取对象的骨架线及其骨架线附近的点；
[0044] (3)确定对象的特征点，将轮廓点和骨架线附近的点作为对象的特征点。
[0045] 进一步地，所述步骤六具体包括以下步骤：
[0046] (1)建立格网索引，对特征点集建立二维分块格网索引结构；
[0047] (2)选择初始种子点，从二维分块格网索引结构的每个格网单元格中选取格网中 高程最低的特征点，同时将特征点所在对象包含的所有特征点作为地面种子点；
[0048] (3)构建初始种子点的三角网，利用所有的地面种子点构建一个三角网，该三角网 代表该区域的初始数字高程模型；
[0049] (4)判别对象的类别，通过迭代特征点来判别对象的类别；
[0050] (5)更新对象包含原始激光雷达点的类别，将对象的类别赋值给其包含的原始激 光雷达点，从而完成机载激光雷达点云滤波工作。
[0051] 进一步地，所述判别对象的类别是一个迭代的过程，判别对象类别过程中涉及三 个参数:最大角度阈值Θ、最大距离阈值d和最大地形角度阈值t;判别对象类别过程包括以 下步骤：
[0052] 1)迭代次数记为0,对于每一个待判别的类别为"Γ的特征点，找到该特征点落入 的三角形，计算该特征点到三角形构成的平面的距离及该特征点到三角形三点的夹角、并 找出三个夹角中的最大角；如果同时满足:距离小于最大距离阈值d、最大角小于最大角度 阈值Θ，则认为该特征点是地面点，将该点的类别号标记为地面点；
[0053] 2)以对象为基本判别单元，统计该对象中地面点的比例。如果该比例大于50%，将 该对象和相关特征点的类别号标记为地面点；否则，将该对象和相关特征点的类别号标记 为 "1"。
[0054] 3)将获得的地面点插入三角网，同时迭代次数自增1;
[0055] 4)重复上述步骤1)至3)，继续执行直至迭代次数达到5后则停止迭代。
[0056]本发明的有益效果如下：本发明在剔除粗差点后，对机载激光雷达点云进行基于 光滑表面生长的点云分割获取对象，根据对象的多回波比例特征识别潜在的地物对象、并 剔除地物对象及其所包含的激光雷达点，提取对象的特征点、并利用特征点替代对象包含 的原始激光雷达点参与后续运算，以对象作为基本单元参与基于三角网渐进加密的判别以 区分地面对象(地面对象包含的点为地面点）和地物对象(地物对象包含的点为地物点），显 著的提高了滤波的精度。
[0057]与现有技术相比，本发明具有以下特点：
[0058] (1)本发明利用点云分割获取对象、利用对象的特征点替代对象包含的原始点进 行面向对象的机载激光雷达点75Γ滤波，有效提尚了机载激光雷达点75Γ滤波的精度。
[0059] (2)本发明创造性的将粗差剔除、点云分割、多回波分析、特征点提取、基于三角网 渐进加密的对象类别判别等五个步骤进行有机结合起来，形成了一套完整的面向对象的机 载激光雷达点云滤波的技术流程，为机载激光雷达点云滤波技术提供了有效的途径。
[0060] (3)相比现有的机载激光雷达点云滤波方法，本发明的滤波精度更高，对山区、城 区等多种地貌类型具有广泛的适用性。
【附图说明】
[0061 ]图1为本发明的方法流程图；
[0062]图2为现有技术采用的三角网渐进加密滤波方法的局限性分析示意图；
[0063]图3(a)为某一区域的机载激光雷达点云数据示意图；
[0064]图3(b)为机载激光雷达点云数据的高程直方图；
[0065] 图4为本发明采用的点云分割方法的技术流程图；
[0066] 图5(a)为某一对象的初始轮廓点判定的示意图；
[0067] 图5(b)为某一对象的初始轮廓线生成过程的示意图；
[0068] 图6(a)为某对象涉及的LiDAR点的示意图；
[0069] 图6(b)为某对象的特征点的示意图；
[0070]图7(a)为某机载LiDAR点云数据的示意图；
[0071]图7(b)为点云分割的结果示意图；
[0072] 图7(c)为初始的地面种子点的示意图；
[0073] 图7(d)为由种子点构建的三角网的示意图；
[0074]图7(e)为角度和距离判别的示意图；
[0075]图7(f)为镜像过程的示意图；
[0076] 图8为本发明的整体数值计算流程图；
[0077] 图9(a)场景1及Samplell和Samplel2的范围的示意图；
[0078] 图9(b)场景2及Sample21至Sample24的范围的示意图；
[0079]图9(c)场景3及Sample31的范围的示意图；
[0080] 图9(d)场景4及Sample41和Sample42的范围的不意图；
[0081 ] 图9(e)场景5及Sample51至Sample54的范围的示意图；
[0082] 图9(f)场景6及Sample61的范围的示意图；
[0083]图9(g)场景7及Sample71的范围的示意图；
[0084]图10(a)为两种滤波方法在15个样本区的Typel误差分布规律示意图；
[0085]图10(b)为两种滤波方法在15个样本区的Typell误差分布规律示意图；
[0086]图10(c)为两种滤波方法在15个样本区的总误差的分布规律示意图。
【具体实施方式】
[0087]为能清楚说明本方案的技术特点，下面通过【具体实施方式】，并结合其附图，对本发 明进行详细阐述。下文的公开提供了许多不同的实施例或例子用来实现本发明的不同结 构。为了简化本发明的公开，下文中对特定例子的部件和设置进行描述。此外，本发明可以 在不同例子中重复参考数字和/或字母。这种重复是为了简化和清楚的目的，其本身不指示 所讨论各种实施例和/或设置之间的关系。应当注意，在附图中所图示的部件不一定按比例 绘制。本发明省略了对公知组件和处理技术及工艺的描述以避免不必要地限制本发明。
[0088]针对机载激光雷达点云滤波的需要，本发明提供了一种面向对象的机载激光雷达 点云滤波方法，利用对象的特征点替代对象包含的原始点进行面向对象的机载激光雷达点 云滤波，提高了滤波的精度。
[0089] 如图1所示，本发明的一种面向对象的机载激光雷达点云滤波方法，它包括以下步 骤：
[0090] 步骤一，加载机载激光雷达点云数据；
[0091 ]步骤二，识别机载激光雷达点云中的粗差点，并剔除粗差点；
[0092] 步骤三，进行基于光滑表面生长的机激光雷达点云分割以获取对象；
[0093] 步骤四，分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其 所包含的激光雷达点；
[0094] 步骤五，提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与 后续运算；
[0095] 步骤六，进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的 类别。
[0096] 进一步地，所述步骤二具体包括以下步骤：
[0097] (1)预处理
[0098]将待处理机载激光雷达点云涉及所有点的类别号标记为"Γ，"Γ代表了未分类的 点。
[0099] (2)目视分析机载激光雷达点云的高程分布特征
[0100]粗差是影响机载激光雷达点云数据滤波效果的关键因素之一，因为大部分滤波算 法在选择初始地面点时常选取局部最低点为地面点。如果该局部最低点为粗差点，则很容 易将与该粗差点邻近的地面点判别为地物点，会严重地影响一定区域范围内的滤波的效 果。因此，粗差剔除是机载激光雷达数据滤波的前提之一。由于机载激光雷达点云数据中的 部分粗差呈现出显著的高于或者低于其他激光雷达点的现象，如图3(a)所示，对于此类显 著粗差可以通过高程直方图（如图3(b)所示)分析进行剔除。在目视分析机载激光雷达点云 的高程分布特征后，如果发现该点云数据中存在明显的高位粗差和低位粗差，直接进入步 骤(3);如果发现该点云数据中没有明显的高位粗差和低位粗差，此步骤可以省略，直接进 入步骤(4)。
[0101] (3)建立点云的高程直方图识别显著的粗差
[0102] 利用机载激光雷达点云的高程z特征，建立高程直方图，如图3(b)所示。通过观察 高程直方图，人工分别确定一个较小的 高程值ζι?和一个较大的高程值Zhigh。将高程值小于 Z1?和大于Zhlgh的激光雷达点的类别号标记为"7"，"7"代表了粗差点。
[0103] (4)建立kd-树空间索引
[0104] 对于类别号为"Γ的激光雷达点集，利用该点集的水平坐标信息建立二维的kd_ 树。
[0105] (5)识别低位粗差和高位粗差
[0106] 输入经验值k和η，其中k为邻近点的个数、η为高差标准差的倍数。对于任意一个类 别号为"Γ的激光雷达点，利用kd-tree空间索引寻找其最邻近的k个激光雷达点;计算该激 光雷达点与k个邻近激光雷达点的高程值的平均数^[email protected];另外，分别计算该激光雷达点与 k个邻近激光雷达点的高程值之差，计算k个高程值差的标准差σζ;如果该激光雷达点的高 程值与平均数之差的绝对值大于ησ ζ，则将该点类别号标记为"7"，"7"代表了粗差点。
[0107] 进一步地，所述步骤三具体包括以下步骤：
[0108] (1)估计法向量和残差
[0109] 该分割方法的流程图如图4所示。点云数据中，任意一点的法向量定义为该点及其 一定数量的k个三维空间中最邻近点确定的平面的法向量。对于类别号为"2"的激光雷达点 集，任何一点的k个最邻近点可以通过三维kd-树这种空间索引的方法求取。同时，三维空间 中一定数量的离散点确定的平面可以由特征值法(即主成分分析法)确定。同时特征值法中 获取的最小特征值，记为该点对应的残差。法向量和残差估计的具体过程包括：
[0110] 1)对于类别号为"Γ的激光雷达点集，建立其三维的kd-树空间索引；
[0111] 2)对每一个类别号为"Γ的激光雷达点，逐一利用kd-树空间索引求取其k个最临 近点。k的取值由用户输入，默认值为20;
[0112] 3)对于每一点及其k个最临近点，利用特征值法求取其拟合平面的方程，即确定每 一点的法向量Φ i及其残差h。
[0113] (2)进行区域生长
[0114] 区域生长过程与传统的基于图像的区域生长极其相似。两者都包含两个必要的因 素，邻接关系的确定和相似性的度量。其中，邻接关系由kd-树确定，相似性的度量包括邻接 点到种子平面的距离和邻接点与种子点的法向量角度差两个特征。区域生长的结果是，将 激光雷达点云划分为若干的对象，每个对象包含若干的激光雷达点。该区域生长的具体过 程包括：
[0115] 1)输入法向量间角度差异阈值α。角度差异阈值α的使用规则如下:对于当前要扩 张的种子点，利用kd-树求取其某一邻接点作为待处理的点。设当前种子点的法向量为h， 待处理的邻接点的法向量为七。如果待处理的邻接点与种子点之间法向量的点积| |aiaj| 小于阈值cosa，则将此点加入到种子点队列中，并作为新的种子点，继续判别其他的未处理 的点；否则，跳过该点，处理其他的未处理的点。
[0116] 2)输入给定距离阈值r。距离阈值r如果待处理的邻接点到种子点确定的种子平面 的距离小于阈值r，则将此点加入到种子点队列中，并作为新的种子点，继续判别其他的未 处理的点；否则，跳过该点，处理其他的未处理的点。
[0117] 3)检查类别号为"Γ的激光雷达点集包含的点的数量。如果数量为0,则到步骤7)。 否则，将全部类别号为"Γ的激光雷达点的状态记为"未被分割"。从"未被分割"的激光雷达 点集中，寻找出残差λ最小的点。以残差λ最小的点为种子点，将该点压入一个种子点的队 列，开始进行区域生长。
[0118] 4)生长的过程中，取种子点队列中第一个尚未被处理过的种子点，利用kd-tree求 取该种子点的k个最临近点。
[0119] 5)对于每一个临近点，逐一进行判别。如果该点已经被赋予分割号，则不予以处 理。进一步，如果该点已经在种子点的队列中，则不予以处理。接着，分别按照法向量间角度 差异和距离差异进行相似性的判别。如果该邻近点同时满足两个相似性的条件，则将该邻 近点加入到种子点队列中；否则，不予以处理。
[0120] 6)回到步骤4)，检查种子点队列中是否有未被处理过的点。如果有，则进入步骤 5);否则，将种子点队列中的所有点标记一个相同的分割号，且该分割号与已有的分割号不 同。然后，返回步骤3)。
[0121] 7)结束分割。
[0122] (3)赋予对象及相关激光雷达点初始的类别
[0123] 对于任意一个对象，检查其该对象包含的所有激光雷达点的数量。如果该对象包 含的所有激光雷达点的数量小于q(q为经验阈值），则将该对象包含的所有激光雷达点的类 别号标记为"3"，"3"代表了地物点。
[0124] 进一步地，所述步骤四具体包括以下步骤：
[0125] (1)统计对象的多回波比例
[0126] 以点云分割获取的对象为基本组织单元，统计每一个对象包含的激光雷达点的总 数量、多回波激光雷达点的数量以计算多回波的比例。目前，对于多回波机载激光雷达系统 而言，系统记录的回波信息包括单次回波（singular return)和多次回波（multiple returns)，二者的区别在于对同一束激光脉冲是否发生多次反射。即，当激光扫描仪发射的 激光脉冲接触到被测目标时，部分脉冲能量的反射信号会被系统接收并记录，而剩余的脉 冲能量继续传播，当遇到另一目标或原被测目标的另一部分时再次发生反射，直至能量消 耗殆尽。如果激光束只碰到一次目标且回波被接受纪录下来，则称此回波为单次回波;如果 发生多次反射、并使得机载激光扫描系统接收到多个反射信号即多次回波。现有的机载激 光雷达系统都具备记录2-5次(超过5次以上回波的概率几乎为零)多回波信息的能力。对于 多次回波，接收到的第一个回波信号被称为首次回波(f irst return/pu 1 se);接收到的最 后一个回波信号被称为尾次(末次）回波（last pulse);除首次和尾次回波外的中间的回波 信号被顺序编号，称为中间次回波（intermediate pulses)。研究发现，回波信息可以反映 被测目标的类型信息，继而可以辅助激光点云数据的滤波和分类。对于森林地区：单次回波 包含地面激光脚点和植被冠层的激光脚点；首次回波来自茂密且高大的植被冠层或靠近冠 层的枝叶；中间次回波多为高大植被的枝叶或低矮植被;而末次回波多是植被中间层次的 枝叶和地表反射得到的激光脚点。对于城区：单次回波数据主要来源于地表、人工建筑物 (包括立交桥和横跨河流的桥梁）的顶面或墙面、少量植被点；首次回波来源于植被的冠层 和人工建筑物(包括立交桥）的边缘；中间次回波主要来源于植被的枝叶和建筑物的立面； 而末次回波则主要来源于地表，也有部分是来源于复杂的建筑物屋顶和植被低矮层面的枝 叶。上述分析表明，无论对于森林地区还是城市区域，尾次回波和单次回波既可以来自地 面、又可以来自非地面目标;首次回波和中间次回波肯定来源于非地面目标，如植被、建筑 物边缘等。
[0127] (2)判别对象及相关激光雷达点的类别
[0128] 对于任意一个对象，如果其多回波比例大于50%，则将该对象包含的所有激光雷 达点的类别号标记为"3"，"3"代表了地物点。
[0129] 进一步地，所述步骤五具体包括以下步骤：
[0130] (1)提取对象的轮廓点及轮廓线
[0131]对于任意一个类别号标记为"Γ的对象，找出其包含的该对象包含的所有激光雷 达点，并建立这些点的不规则三角网（TIN)，如图5(a)所示。判断轮廓边缘只需分析各三角 形的相邻三角形的个数。通常一个三角形有3个三角形与之相邻，如果个数小于3,则该三角 形位于对象轮廓的边缘，提取相应轮廓点即可获得对象的初始轮廓点。如图5(b)所示，顺时 针依次搜索各初始轮廓点进一步生成对象的轮廓线，步骤如下：
[0132] 1)起始轮廓点确定:搜索所有初始轮廓点中x值最小的激光雷达点，如存在多个x 值相等的点，则取这些X值相等的点中y值最小的点作为起始轮廓点，如图5(b)中"点Γ所 示；
[0133] 2)第2个轮廓点确定：以"点Γ为起点、R为搜索半径，水平向左为起始搜索方向，沿 顺时针方向寻找下一个轮廓点，搜索半径R由点云平均间隔开始逐渐增大，直到在搜索范围 内找到第2个轮廓点（"点2"）；
[0134] 3)顺时针方向搜索下一个轮廓点：以"点2"为起点、R为搜索半径，从"点1"、"点2" 矢量方向的逆时针90°范围开始顺时针搜索下一个轮廓点，搜索半径逐渐增大（图5 (b)中搜 索半径由R1增大到R2)直到找到第3个轮廓点（"点3"）；
[0135] 4)重复类似的步骤3)，直到找到"点Γ则终止搜索，依次连接各轮廓点生成对象的 轮廓线。
[0136] (2)提取对象的骨架线及其骨架线附近的点
[0137] 基于上一步提取的对象轮廓线形成的多边形，提取该多边形的骨架线。在水平面 内，将该对象中距离骨架线在两倍点云平均间隔内、且不属于轮廓点的点标记为该对象的 "骨架线附近的点"。
[0138] (3)确定对象的特征点
[0139] 将对象的轮廓点和"骨架线附近的点"作为该对象的"特征点"。对象的"特征点"将 代替该对象参与后续运算。图6(a)展示了某对象涉及的激光雷达点云，同时 图6(b)展示了 该对象的特征点。
[0140] 进一步地，所述步骤六具体包括以下步骤：
[0141] (1)建立格网索引
[0142] 对于类别号为"Γ的特征点集，首先根据点云数据在x、y方向的坐标范围，构建一 个特征点集的最小包围核。根据给定格网尺寸，即最大建筑物尺寸m，分别对特征点集在X、y 方向划分正方形格网，对特征点集建立二维分块格网索引结构。
[0143] (2)选择初始种子点
[0144] 从每个格网单元格中选取格网中高程最低的特征点，同时将特征点所在对象包含 的所有特征点作为地面种子点、并将该对象和相关特征点的类别号标记为"2"，"2"代表了 地面点。图7(a)展示了某机载激光雷达点云，图7(b)展示了该激光雷达点云的分割效果，同 时图7(c)展示了初始的种子点。
[0145] (3)构建初始种子点的三角网
[0146] 利用所有的地面种子点构建一个三角网，该三角网代表了该区域的初始数字高程 模型DEM。图7(c)和图7(d)分别展示了某机载激光雷达点云数据的初始种子点和这些种子 点构建的三角网。
[0147] (4)判别对象的类别
[0148] 对象类别的判别是一个迭代的过程，该过程涉及三个参数:最大角度阈值Θ、最大 距离阈值d和最大地形角度阈值t。
[0149] 1)迭代次数记为0。对于每一个待判别的类别为"Γ的特征点，找到该点落入的三 角形，计算该点到三角形构成的平面的距离及其该点到三角形三点的夹角、并找出三个夹 角中的最大角，如图7(e)所示。
[0150] 如果同时满足:距离小于最大距离阈值d、最大角小于最大角度阈值Θ，则认为该点 是地面点，将该点的类别号标记为"2"，"2"代表了地面点；否则，还要进一步判别(对当前的 特征点以距离最近的三角形顶点为中心进行镜像，如图7(f)所示。对镜像点进行判别。如果 此镜像点判断为地面点，则仍然判断此点为地面点、将该点的类别号标记为"2" ；否则，不修 改点的类别值）。继续处理下一个类别为"Γ的特征点直至将所有的类别为"Γ的特征点遍 历完毕。
[0151] 2)以对象为基本判别单元，统计该对象中"地面点"的比例。如果该比例大于50%， 将该对象和相关特征点的类别号标记为"2"，"2"代表了地面点；否则，将该对象和相关特征 点的类别号标记为"1"。
[0152] 3)将获得的地面点插入三角网。迭代次数自增1。回到步骤1)继续执行直至迭代次 数达到的5则停止迭代。
[0153] (5)更新对象包含原始激光雷达点的类别
[0154] 将对象的类别赋值给其包含的原始激光雷达点，完成滤波。
[0155] 本发明创造性地提出了一种面向对象的机载激光雷达点云滤波方法，其特色是利 用对象的特征点替代对象包含的原始点进行面向对象的机载激光雷达点云滤波，提高了滤 波的精度。该方法创造性的将粗差剔除、点云分割、多回波分析、特征点提取、基于三角网渐 进加密的对象类别判别等五个步骤连接起来，形成了一套完整的面向对象的机载激光雷达 点云滤波的技术流程。
[0156] 如图8所示，本发明进行面向对象的机载激光雷达点云滤波的数值计算流程如下：
[0157] (I)机载数据。加载机载激光雷达点云数据。
[0158] (II)粗差剔除。在进行滤波之前需要对激光点云数据进行预处理，去除极低局外 点。
[0159] (III)点云分割。对粗差剔除后的激光雷达点云，进行基于光滑表面生长的分割以 获取对象。
[0160] (IV)回波分析。分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对 象及其所包含的激光雷达点。
[0161] (V)提取对象的特征点。以每一个对象为处理单元，取出其该对象涉及的激光雷达 点，并获取每一个对象的边缘点和骨架线附近的点；利用边缘点和骨架线附近的点组成该 对象的特征点集;该特征点集将代替该对象及其包含的原始激光雷达点，参与后续运算。
[0162] (VI)判别。进行格网划分建立二维分块格网索引结构;选择初始种子点；构建初始 种子点的不规则三角网；迭代的通过特征点的判别确定对象的类别，并更新对象包含的原 始激光雷达点的类别。
[0163] 为了验证上述提出的滤波方法的性能，本研究同时采用经典的三角网渐进加密方 法和本发明提出的方法，以ISPRS第三委员会机载激光雷达数据处理小组提供的专门用于 滤波算法测试数据作为实验数据，通过定性和定量分析的手段比较两种算法的性能。
[0164] ISPRS第III小组于2003年组织的滤波算法测试中，提供了一组数据集和参考数 据，用于定量的分析算法的滤波效果（数据来自：111^口 ：//?^^.(301]11]1188；[0113.18口18.0找/ wg3/index.html)。该数据集包含7个场景的点云数据，依次命名为：CStitel、CStite2、 05衍七63、05衍七64、05衍七65、05衍七66、05衍七67，如图9(3)至图9(8)所示。这7个场景的数据 为OptechALTM机载激光扫描仪获取的Vaihingen/Enz和Stuttgart市的机载激光雷达点云 数据。其中的城区数据为Sitel-Site4,点间距为1.0至1.5m;山区为Site5-Site7,点间距为 2.0至3.5m。同时，在7块测试数据中选取了 15个样本数据(具体数据区位如图9(a)至图9(g) 所示），这些样本数据经过人工处理，将激光脚点精确分类为地面点和非地面点这两类。15 个样本数据主要包括了不同地理环境中滤波可能遇到困难的情况，如粗差点的影响、复杂 的地物、地物与地面相连、斜坡上或低矮的植被、地形的不连续等，具体数据特征如表1所 示。利用本发明提出的滤波算法对这7个场景的数据做滤波实验，利用15个参考样本对滤波 结果进行了定性和定量分析，并与已有的三角网渐进加密滤波算法的结果进行了对比分 析。
[0165] ISPRS第三委员会机载激光雷达数据处理小组提出了滤波算法质量评价的两种方 法:定性分析和定量分析。定性分析中，把对滤波后的结果叠加到原始数据上，通过误差图 进行可视化检查、分析和评价。定量分析中，统计了三种类型的误差，即Typel(地面点被误 分为地物点的数量占地面点数量的比例）、Type Π (地物点被误分为地面点的数量占地物点 数量的比例)和总误差(被误分的地面点和和地物点的总数据量占全部脚点数量的比例）， 见表3。
[0166] 表1参考数据的说明
[0167]
[0168] 表2定性分析的等级
[0169]
[0170]~表3定量分析的相应参数及其定义
[0171]
[0172] 注释:表中，a，b，c，d是进行定量分析的基本参数，a代表参考点中是地面点并且滤 波算法也将其分类为地面点的地面点数量，b代表参考点中是地面点但滤波算法将其分类 为地物点的地面点数量，c代表参考点中是地物点而滤波算法将其分类为地面点的地物点 数量，d代表参考点中是地物点且滤波算法也将其分类为地物点的地物点数量;f，g，h，i分 析数据滤波前后的地面点和地物点组成情况，f是参考数据滤波前的地面点的比例，g是参 考数据滤波前的地物点的比例，h是滤波结果中地面点的比例，i是滤波结果中地物点的比 例，；Typel、Type Π 和Total是三类误差;k反映了两类误差的偏离程度。本论文中，j = 0。
[0173] 表4两种滤波方法参数的取值
[0174]
[0176]已有的三角网渐进加密滤波方法需要最大建筑物尺寸、最陡地形倾角、迭代角度、 迭代距离等四个参数。本发明提出的滤波方法，在已有方法所需四个参数基础上增加了三 个参数，即邻接点个数、距离阈值、角度差异阈值。每个场景的实验中，两种方法共有参数的 取值是相同的，本发明所需的额外三个参数的取值需要单独确定。两种方法在7个ISPRS提 供的实验数据的滤波参数的取值情况，见表4。
[0177] 使用表4中列举的参数，利用两种滤波方法对上述7个实验数据进行滤波。同时，从 ISPRS提供的7个场景的滤波结果中，分别选择15个参考样本数据具有相同坐标值的激光脚 点，并利用参考样本的滤波结果评价15个相关区域的滤波效果，得到了三类误差的统计表 5。另外，图10(a)、（b)、（c)分别展示了两种滤波方法的三种类型的误差对比图。
[0178] 表5经典的三角网滤波方法与本发明改进方法的误差对比表(单位：％ )
[0179]
[0180]
[0181]从表5可以看出，已有的三角网渐进加密滤波方法和本发明提出的滤波方法（即表 5以及图10(a)、图10(b)和图10(c)中所述的改进方法)都具有较强的鲁棒性，能广泛的适应 于城区和森林区域各种场景的滤波，且滤波的总误差控制在33%以下，如图10(c)所示。同 时，与经典的滤波方法相比，本发明的滤波方法具有更好的性能，能正确的识别更多的地面 点，尤其是地形陡变处的地面点，其I类误差和总误差显著的低于经典滤波的相应误差。15 个局部区域的滤波效果检验中，除samp21区域之外，其他14个区域中本发明滤波方法的I类 误差显著的低于已 有三角网渐进加密滤波方法的I类误差；除samp21、samp54两个区域之 外，其他13个区域中本发明滤波方法的总误差显著的低于已有三角网渐进加密滤波方法的 总误差，如图10(a)和（c)所示。但是，本发明滤波方法也有一定的局限性，与已有三角网渐 进加密滤波方法相比，它具有更高的II误差。15个局部区域的滤波效果检验中，除samp61之 外，其他14个区域中本发明滤波方法的II类误差高于已有三角网渐进加密滤波方法的II类 误差，如图10(b)所示。整体而言，表5中的统计数据表明，本发明滤波方法比已有三角网渐 进加密滤波方法具有更高的滤波精度。
[0182]以上所述只是本发明的优选实施方式，对于本技术领域的普通技术人员来说，在 不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也被视为本发 明的保护范围。
【主权项】
1. 一种面向对象的机载激光雷达点云滤波方法，其特征是，包括以下步骤： 步骤一，加载机载激光雷达点云数据； 步骤二，识别机载激光雷达点云中的粗差点，并剔除粗差点； 步骤三，进行基于光滑表面生长的机激光雷达点云分割以获取对象； 步骤四，分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其所包 含的激光雷达点； 步骤五，提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与后续 运算； 步骤六，进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的类别。2. 根据权利要求1所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤二具体包括以下步骤： (1) 初始化预处理，对机载激光雷达点云的类别进行标记，将未分类的点标记为"Γ; (2) 目视分析机载激光雷达点云的高程分布特征，检查是否存在显著的粗差，如果该点 云数据中存在明显的高位粗差和低位粗差，则进入步骤(3);如果该点云数据中没有明显的 尚位粗差和低位粗差，则进入步骤(4); (3) 建立点云的高程直方图识别机载激光雷达点云显著的粗差； (4) 建立kd-树空间索引，对于类别号为"Γ的激光雷达点集，利用该点集的水平坐标信 息建立二维的kd-树； (5) 识别低位粗差和高位粗差，利用kd-tree空间索引寻找粗差点。3. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 高程直方图通过机载激光雷达点云的高程z特征来建立。4. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤三具体包括以下步骤： (1) 估计法向量和残差，点云数据中任意一点及其一定数量的k个三维空间中最邻近点 确定的平面的法向量为该点的法向量，三维空间中一定数量的离散点确定的平面由特征值 法确定，通过特征值法中获取的最小特征值记为离散点对应的残差； (2) 进行区域生长，将激光雷达点云划分为若干的对象，每个对象包含若干的激光雷达 占. (3) 赋予对象及相关激光雷达点初始的类别，对于任意一个对象，如果其多回波比例大 于50%，则将该对象包含的所有激光雷达点的类别号标记为地物点。5. 根据权利要求4所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 估计法向量和残差的具体过程包括以下步骤： (1) 对于类别号为"Γ的激光雷达点集，建立其三维的kd-树空间索引； (2) 对每一个类别号为"Γ的激光雷达点，逐一利用kd-树空间索引求取其k个最临近 占. (3) 对于每一点及其k个最临近点，利用特征值法求取其拟合平面的方程，即确定每一 点的法向量Φi及其残差h。6. 根据权利要求4所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 区域生长的具体过程包括以下步骤： 1) 设置法向量间角度差异阈值α，角度差异阈值α的使用规则如下：对于当前要扩张的 种子点，利用kd_树求取其某一邻接点作为待处理的点，设当前种子点的法向量为h，待处 理的邻接点的法向量为aj;如果待处理的邻接点与种子点之间法向量的点积| |aiC[j| |小于 阈值cosa，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续判别其 他的未处理的点；如果待处理的邻接点与种子点之间法向量的点积I|^七|I大于等于阈值 cosa，则跳过该待处理的邻接点，处理其他的未处理的点； 2) 设置给定距离阈值r，如果待处理的邻接点到种子点确定的种子平面的距离小于给 定距离阈值r，则将该待处理的邻接点加入到种子点队列中，并作为新的种子点，继续判别 其他的未处理的点；如果待处理的邻接点到种子点确定的种子平面的距离大于等于给定距 离阈值r则跳过该待处理的邻接点，处理其他的未处理的点； 3) 检查类别号为"Γ的激光雷达点集包含的点的数量，如果数量为0,则到步骤7)，否 贝1J，将全部类别号为"Γ的激光雷达点的状态记为"未被分割"；从"未被分割"的激光雷达点 集中，寻找出残差λ最小的点，以残差λ最小的点为种子点并将该点压入一个种子点的队列， 开始进行区域生长； 4) 生长的过程中，取种子点队列中第一个尚未被处理过的种子点，利用kd-tree求取该 种子点的k个最临近点； 5) 对于每一个临近点，逐一进行判别，如果临近点已经被赋予分割号，则不予以处理； 如果临近点已经在种子点的队列中，则不予以处理;分别按照法向量间角度差异和距离差 异进行相似性的判别，如果邻近点同时满足两个相似性的条件，则将该邻近点加入到种子 点队列中；如果邻近点没有同时满足两个相似性的条件则不予以处理； 6) 回到步骤4)，检查种子点队列中是否有未被处理过的点，如果有，则进入步骤5);否 贝1J，将种子点队列中的所有点标记一个相同的分割号，且该分割号与已有的分割号不同，然 后，返回步骤3); 7) 结束分割。7. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤四具体包括以下步骤： (1) 统计对象的多回波比例，以点云分割获取的对象为基本组织单元，统计每一个对象 包含的激光雷达点的总数量、多回波激光雷达点的数量以计算多回波的比例； (2) 根据多回波比例特征判别对象及相关激光雷达点的类别，对于任意一个对象，如果 其多回波比例大于50%，则将该对象包含的所有激光雷达点的类别号标记地物点。8. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤五具体包括以下步骤： (1) 提取对象的轮廓点及轮廓线，采用建立不规则三角网提取对象的轮廓点，依次连接 各轮廓点生成对象的轮廓线； (2) 提取对象的骨架线及其骨架线附近的点； (3) 确定对象的特征点，将轮廓点和骨架线附近的点作为对象的特征点。9. 根据权利要求2所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所述 步骤六具体包括以下步骤： (1)建立格网索引，对特征点集建立二维分块格网索引结构； (2) 选择初始种子点，从二维分块格网索引结构的每个格网单元格中选取格网中高程 最低的特征点，同时将特征点所在对象包含的所有特征点作为地面种子点； (3) 构建初始种子点的三角网，利用所有的地面种子点构建一个三角网，该三角网代表 该区域的初始数字高程模型； (4) 判别对象的类别，通过迭代特征点来判别对象的类别； (5) 更新对象包含原始激光雷达点的类别，将对象的类别赋值给其包含的原始激光雷 达点，从而完成机载激光雷达点云滤波工作。10.根据权利要求9所述的一种面向对象的机载激光雷达点云滤波方法，其特征是，所 述判别对象的类别是一个迭代的过程，判别对象类别过程中涉及三个参数:最大角度阈值 Θ、最大距离阈值d和最大地形角度阈值t;判别对象类别过程包括以下步骤： 1) 迭代次数记为0,对于每一个待判别的类别为"Γ的特征点，找到该特征点落入的三 角形，计算该特征点到三角形构成的平面的距离及该特征点到三角形三点的夹角、并找出 三个夹角中的最大角；如果同时满足:距离小于最大距离阈值d、最大角小于最大角度阈值 Θ，则认为该特征点是地面点，将该点的类别号标记为地面点； 2) 以对象为基本判别单元，统计该对象中地面点的比例。如果该比例大于50%，将该对 象和相关特征点的类别号标记为地面点；否则，将该对象和相关特征点的类别号标记为 T。 3) 将获得的地面点插入三角网，同时迭代次数自增1; 4) 重复上述步骤1)至3)，继续执行直至迭代次数达到5后则停止迭代。
【专利摘要】本发明公开了一种面向对象的机载激光雷达点云滤波方法，它包括以下步骤：加载机载激光雷达点云数据；识别机载激光雷达点云中的粗差点，并剔除粗差点；进行基于光滑表面生长的机激光雷达点云分割以获取对象；分析对象的多回波比例特征识别潜在的地物对象，并剔除地物对象及其所包含的激光雷达点；提取对象的特征点，并利用特征点替代对象包含的原始激光雷达点参与后续运算；进行基于特征点的对象类别判别，并更新对象包含的原始激光雷达点的类别。本发明形成了一套完整的面向对象的机载激光雷达点云滤波的技术流程，有效提高了机载激光雷达点云滤波的精度，对山区、城区等多种地貌类型具有广泛的适用性。
【IPC分类】G06T5/00
【公开号】CN105488770
【申请号】CN201510920332
【发明人】林祥国, 段敏燕, 张继贤
【申请人】中国测绘科学研究院
【公开日】2016年4月13日
【申请日】2015年12月11日
最新回复(0)