CN109785261B - 一种基于灰度体元模型的机载lidar三维滤波方法 - Google Patents

Abstract

本发明提供一种基于灰度体元模型的机载LIDAR三维滤波方法，涉及遥感数据处理技术领域。该方法首先读取原始机载LIDAR点云数据，形成原始机载LIDAR点云数据集；将机载LIDAR点云数据规则化为灰度体元模型，其中灰度为体元内激光点的平均强度的离散化表示；基于地面点的局部高程最低特性选取地面种子体元，进而标记与其三维连通且体元值、地形坡度均接近的非0值体元为地面体元。本发明提供的基于灰度体元模型的机载LIDAR三维滤波方法，以三维连通区域构建理论为基础，将强度及地形坡度信息作为辅助信息用于点云数据滤波，可为地面和非地面目标的区分提供更有效的信息，从而提高滤波精度，并扩展基于体元模型的三维滤波方法适用于更复杂的场景。

Description

一种基于灰度体元模型的机载LIDAR三维滤波方法

技术领域

本发明涉及遥感数据处理技术领域，尤其涉及一种基于灰度体元模型的机载LIDAR三维滤波方法。

背景技术

机载激光雷达(Light Detection and Ranging，即LIDAR)数据滤波是目前点云数据处理领域中的关键技术，其滤波的精度直接影响数字高程模型(Digital ElevationModel，DEM)等产品的质量及后续的如建筑物、植被等目标分类的精度。因此，对于LIDAR点云数据滤波技术的研究具有重要意义。

国内外学者围绕点云数据滤波进行了深入研究，提出了非常多的滤波算法，经典的如基于内插、基于坡度、基于数学形态学、基于三维连通集合构建和基于聚类分割的滤波方法等。这些滤波算法所采用的数据结构包括栅格格网、不规则三角网(TIN)、体元和离散点云等。栅格格网和TIN为2.5维数据结构，用其表达机载LIDAR点云数据将导致信息损失并进一步影响基于该类数据结构的滤波结果的完整性。离散点云内部的空间结构及拓扑信息难以利用，导致基于该类数据结构的滤波算法设计困难、效率低下；八叉树的各个节点的尺寸不一，各节点间的邻接关系难以建立，这同样增加了基于该类数据结构的滤波算法设计的难度。而基于体元数据结构的滤波算法可以很好的避免上述方法的缺陷，这是因为：(1)体元结构为真三维数据结构；(2)其内部的体元间隐含有邻接和拓扑信息。因而基于体元结构的滤波算法设计更简单。但是，目前已有的基于体元数据结构的滤波算法均采用二值体元结构(即体元赋值是根据体元内是否含有激光点分别赋值1和0，其中1值体元对应目标，0值体元对应背景)(Liying Wang，Yan Xu，and Yu Li.Aerial LIDAR point cloudvoxelization with its 3D ground filtering application[J].PhotogrammetricEngineering and Remote Sensing，2017，83(2)：95-107)。这种滤波算法认为地面目标会形成三维连通集合，但当其它目标(如低矮植被、建筑物立面)和地面目标相连形成三维连通集合时，该算法则无法有效对二者加以区分。另外，该算法未考虑地形坡度信息，由此导致算法的滤波精度较低。

发明内容

本发明要解决的技术问题是针对上述现有技术的不足，提供一种基于灰度体元模型的机载LIDAR三维滤波方法，解决相连的地面和非地面目标的有效区分问题。

为解决上述技术问题，本发明所采取的技术方案是：一种基于灰度体元模型的机载LIDAR三维滤波方法，包括以下步骤：

步骤1：读取原始机载LIDAR点云数据，形成原始机载LIDAR点云数据集；

步骤2：将原始机载LIDAR点云数据集规则化为灰度体元模型，其中灰度为体元内激光点的平均强度的离散化表示，具体方法为：

步骤2.1：从原始机载LIDAR点云数据集中剔除异常数据，得到剔除异常数据集；

步骤2.1.1：统计原始机载LIDAR点云数据集中各个激光点高程值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.2：确定与真实地形及地物对应的最高高程阈值T_he和最低高程阈值T_le；

步骤2.1.3：针对原始机载LIDAR点云数据集中各个激光点，若其高程值高于最高高程阈值T_he或低于最低高程阈值T_le，则该激光点为高程异常数据，进行剔除，否则保留该激光点，获得剔除高程异常数据集；

步骤2.1.4：统计剔除高程异常数据集中各激光点的强度值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.5：确定与真实地形及地物对应的最高强度阈值T_hi和最低强度阈值T_li；

步骤2.1.6：针对剔除高程异常数据集中各个激光点，若其强度值高于最高强度阈值T_hi或低于最低强度阈值T_hi，则该激光点为强度异常数据，进行剔除，否则保留该激光点，最终获得剔除高程及强度异常数据集；

步骤2.2：将剔除异常数据集规则化为灰度体元模型；

步骤2.2.1：用剔除异常数据集的轴向平行包围盒表示数据集的空间范围；

步骤2.2.2：计算体元分辨率即体元大小，x、y、z方向上的分辨率根据剔除异常数据集中激光点的平均点间距确定；

步骤2.2.3：依据x、y、z方向上的分辨率对剔除异常数据集的轴向平行包围盒进行划分得到三维格网，将每一个三维格网单元称为体元；

步骤2.2.4：将剔除异常数据集中各个激光点映射到三维格网，进而根据体元中包含的激光点的强度属性为各体元赋值，最后将各体元值离散化到{0，…，255}，得到灰度体元模型；

步骤3：基于三维连通区域构建理论，对灰度体元模型中的地面体元进行检测，具体方法为：

步骤3.1：基于地面点的局部高程最低特性，从灰度体元模型中高程最低的非0值体元作为地面种子体元集合；

步骤3.1.1：在水平方向上，利用设定的格网尺寸对灰度体元模型进行分块，并取各块内高程最低的非0值体元为种子体元；

步骤3.1.2：在各块内搜寻与地面种子的高度差小于高度差阈值T_e的非0值体元为地面种子体元，获得加密后的地面种子体元集合V_s，其中，s＝1，2，...；

步骤3.2：标记地面种子体元及与其三维连通且坡度、灰度均接近的体元构成的三维连通区域为地面体元数据集，完成基于灰度体元模型的机载LIDAR三维滤波；

步骤3.2.1：计算灰度体元模型中非0值体元的灰度值的频率直方图，用高斯混合模型拟合灰度频率直方图，确定与地面目标对应的灰度值分布范围；

步骤3.2.2：对任一地面种子体元V_s，遍历灰度体元模型中与当前地面种子体元V_s三维连通、灰度值位于地面目标对应的灰度范围内且局部地形坡度小于坡度阈值T_s的所有未标记体元，并标记为L_g，直至标记完所有地面种子体元V_s的三维连通区域，即地面体元集合；其中，坡度阈值的确定是自适应的，具体方法为：对当前种子体元V_s，检测其空间邻域内是否含有已标记的地面体元，若有，则将已有地面体元之间的最大坡度值确定为该邻域内的地形坡度阈值T_s；否则，将地形坡度阈值设置为90°。

采用上述技术方案所产生的有益效果在于：本发明提供的一种基于灰度体元模型的机载LIDAR三维滤波方法，以三维连通区域构建理论为基础，将强度及地形坡度信息作为辅助信息用于点云数据滤波，可为地面和非地面目标的区分提供更有效的信息，从而提高滤波精度，并扩展基于体元模型的三维滤波方法适用于更复杂的场景。

附图说明

图1为本发明实施例提供的一种基于灰度体元模型的机载LIDAR三维滤波方法的流程图；

图2为本发明实施例提供的原始机载LIDAR点云数据示意图，其中，(a)为Csite1点云数据，(b)为Csite2点云数据，(c)为Csite3点云数据，(d)为Csite4点云数据，(e)为Fsite5点云数据，(f)为Fsite6点云数据，(g)为Fsite7点云数据；

图3为本发明实施例提供的将原始机载LIDAR点云数据集规则化为灰度体元模型的具体流程图；

图4为本发明实施例提供的激光点投影计算凸壳及面积的原理图；

图5为本发明实施例提供的对灰度三维体元数据集进行地面体元检测的流程图；

图6为本发明实施例提供的样本samp41对应的三维体元矩阵V中的非0值体元的灰度频率直方图；

图7为本发明实施例提供的深度优先遍历灰度体元模型中任一种子体元V_s的三维连通区域的流程图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本实施例以某组遥感数据为例，使用本发明的一种基于灰度体元模型的机载LIDAR三维滤波方法对该组数据进行滤波分析。

一种基于灰度体元模型的机载LIDAR三维滤波方法，如图1所示，包括以下步骤：

步骤1：读取原始机载LIDAR点云数据，形成原始机载LIDAR点云数据集。

本实施例中，采用国际摄影测量与遥感协会(International Society forPhotogrammetry and Remote Sensing，ISPRS)第三工作组提供的专门用于滤波算法测试的15个样本数据作为实验数据，以检验方法的有效性和可行性。如图2所示，这些样本数据剪切自7个测试数据，其中包括了不同复杂程度的场景及滤波可能遇到困难的情况，如粗差点的影响、复杂的地物、地物与地面相连、斜坡上或低矮的植被、地形的不连续等，表1列出了本实施例中各个样本数据的特性及基本参数。样本记录了点云的首、末次回波及强度信息(样本数据中并不包含点云数据的强度信息。此处的强度信息是从对应的7个测试数据中提取所得)。

表1样本数据的特性及基本参数

本实施例中设计采用ISPRS提供的参考数据(15个样本数据被准确的分类为地面点和非地面点)作为标准数据评定算法精度。

本实施例中，定义原始机载LIDAR点云数据P＝{p_i(x_i，y_i，z_i)，i＝1，…，n}，其中，i是原始机载LIDAR点云数据的索引，n是原始机载LIDAR点云数据的个数，p_i是第i个原始机载LIDAR点云数据，其坐标为(x_i，y_i，z_i)。

步骤2：将原始机载LIDAR点云数据集规则化为灰度体元模型，如图3所示，具体方法为：

步骤2.1：从原始机载LIDAR点云数据集中剔除异常数据，得到剔除异常数据集。

步骤2.1.1：统计原始机载LIDAR点云数据集中各个激光点高程值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.2：确定与真实地形及地物对应的最高高程阈值T_he和最低高程阈值T_le；

步骤2.1.3：针对原始机载LIDAR点云数据集中各个激光点，若其高程值高于最高高程阈值T_he或低于最低高程阈值T_le，则该激光点为高程异常数据，进行剔除，否则保留该激光点，获得剔除高程异常数据集；

步骤2.1.4：统计剔除高程异常数据集中各激光点的强度值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.5：确定与真实地形及地物对应的最高强度阈值T_hi和最低强度阈值T_li；

步骤2.1.6：针对剔除高程异常数据集中各个激光点，若其强度值高于最高强度阈值T_hi或低于最低强度阈值T_hi，则该激光点为强度异常数据，进行剔除，否则保留该激光点，最终获得剔除高程及强度异常数据集。

本实施例中，最高高程阈值T_he、最低高程阈值T_le、最高强度阈值T_hi和最低强度阈值T_li均为常数，其取值需根据原始机载LIDAR点云数据的空间分布情况确定。

本实施例中，剔除异常数据集记做Q＝{q_i′(x_i′，y_i′，z_i′)，i′＝1，…，t}，其中，i′是剔除异常数据集中激光点的索引，t是剔除异常数据集中激光点的个数，q_i′是剔除异常数据集中第i′个激光点，其坐标为(x_i′，y_i′，z_i′)。

步骤2.2：将剔除异常数据集规则化为灰度体元模型。

步骤2.2.1：用剔除异常数据集的轴向平行包围盒表示数据集的空间范围。

本实施方式中，剔除异常数据集Q的空间范围可由其轴向平行包围盒(Axis-Aligned Bounding Box，AABB)确定，AABB＝{(x，y，z)|x_min≤x≤x_max，y_min≤y≤y_max，z_min≤z≤z_max}，式中，(x_max，y_max，z_max)和(x_min，y_min，z_min)分别代表剔除异常数据集中x、y和z坐标的最大和最小值。

步骤2.2.2：计算体元分辨率即体元大小，x、y、z方向上的分辨率根据剔除异常数据集中激光点的平均点间距确定。

本实施方式中，体元在x、y、z方向上的分辨率Δx、Δy、Δz的计算公式如式(1)所示：

其中，S_xy＝{(x_i′，y_i′)，i′＝1，…，t}Sxz＝{(x_i′，z_i′)，i′＝1，…，t}、S_yz＝{(y_i′，z_i′)，i′＝1，…，t}分别为剔除异常数据集Q在XOY、XOZ及YOZ平面上的投影所得二维点集，C(·)是对应点集的凸壳，A(·)是对应凸壳的面积，如图4所示。

步骤2.2.3：依据x、y、z方向上的分辨率对轴向平行包围盒进行划分得到三维格网，每一个三维格网单元称为体元。

本实施例中，基于体元分辨率(Δx，Δy，Δz)就可以将轴向平行包围盒划分为三维格网，可用三维体元矩阵表示。设V是三维体元矩阵中的体元集合，如式(2)所示：

V＝{v_j(r_j，c_j，l_j)，j＝1，…，m} (2)

其中，j是体元索引；m是体元数；v_j是第j个体元的体元值；(r_j，c_j，l_j)是第j个体元在体元阵列中的坐标，r_j，c_j，l_j分别为坐标的行、列和层号。

步骤2.2.4：将剔除异常数据集中各个激光点映射到三维格网，进而根据体元中包含的激光点的强度属性为各体元赋值，最后将各体元值离散化到{0，…，255}，得到灰度体元模型；

本实施例中，将剔除异常数据集Q中各个激光点映射到3D格网，进而根据体元中包含的激光点的强度属性为各体元赋值。其中，含有激光点的体元赋值为激光点强度值(若含有多个激光点，则赋值高程最低激光点的强度值)、不含有激光点的体元赋值为0，进一步将各体元赋值离散化到{0，…，255}，得到各体元值。由此，得到灰度体元模型，完成对剔除异常数据集的规则化。

步骤3：基于三维连通区域构建理论，对灰度体元模型中的地面体元进行检测，具体流程如图5所示。

步骤3.1：基于地面点的局部高程最低特性，从灰度体元模型中高程最低的非0值体元作为种子体元集合V_s，其中，s＝1，2，...。

步骤3.1.1：在水平方向上，利用设定的格网尺寸对灰度体元模型进行分块，并取各块内高程最低的非0值体元为种子体元。

本实施例中，格网尺寸需根据地形的复杂性确定，如果地形比较复杂，格网尺寸取值较小，反之，可适当增大格网尺寸。但是，最小格网尺寸不能小于原始机载LIDAR点云数据集中最大目标结构(如建筑物)的尺寸，否则，最大结构内部的点将被错判为地面种子体元。

步骤3.1.2：在各块内搜寻与地面种子的高差小于高度差阈值T_e的非0值体元为地面种子，获得加密后的地面种子集合V_s，其中，s＝1，2，...；

本实施方式中，高差阈值T_e为常数，其取值需根据原始机载LIDAR点云数据的空间分布情况确定。

步骤3.2：标记地面种子体元及与其三维连通且坡度、灰度均接近的体元构成的三维连通区域为地面体元数据集，完成基于灰度体元模型的机载LIDAR三维滤波；

步骤3.2.1：计算灰度体元模型中非0值体元的灰度值的频率直方图，用高斯混合模型(Gaussian Mixture Model，GMM)拟合灰度频率直方图，确定与地面目标对应的灰度值分布范围。

本实施例中，以样本samp41为例，统计其对应的V中的非0值体元的灰度值的频率，并以直方图形式显示，如图6所示。由图6可知，其中的灰度分布呈现多峰性。若假设多峰分布为正态多峰分布，则可用高斯混合模型对图1中的直方图进行拟合，得到各个高斯分布的分布特征参数(均值，标准差)。针对本发明的应用目的是滤波得到地面体元，所以与地面对应的高斯分布的分布特征参数μ，σ及分布范围[0，20]被用于确定地面目标的灰度范围，详细方案如下：若令μ-m_l×σ＝0，μ+m_r×σ＝20，则可确定m_l和m_r，由此可确定地面目标的灰度范围，即[μ-mσ，μ+mσ]，其中，乘数m＝max{m_l，m_r}，m取m_l和m_r的最大值是基于高斯分布的对称性确定的。另外，需要注意的是如果灰度范围中存在负数，则将最小值设置为0但不包含0。

步骤3.2.2：对任一地面种子体元V_s，遍历灰度体元模型中与当前地面种子体元V_s三维连通、灰度值位于地面目标对应的灰度范围内且局部地形坡度小于坡度阈值T_s的所有未标记体元，并标记为L_g，直至标记完所有地面种子体元V_s的三维连通区域，即地面体元集合；其中，坡度阈值的确定是自适应的，具体方法为：对当前种子体元V_s，检测其空间邻域内是否含有已标记的地面体元，若有，则将已有地面体元之间的最大坡度值确定为该邻域内的地形坡度阈值T_s；否则，将地形坡度阈值设置为90°。

本实施例中，深度优先遍历灰度体元模型中任一种子体元V_s的三维连通区域的具体方法如图7所示。

本实施例中，基于V_s和其空间邻域内已标记的地面体元计算坡度E的公式见式(3)：

其中，(r_p，c_p，l_p)为V_s坐标；(r_e，c_e，l_e)为已标记的地面体元坐标，t为空间邻域内已标记地面体元个数。

本实施例中，在上述标记过程中应用不同的邻域尺寸会得到不同的滤波结果。最佳邻域尺度将在实验中确定。

本发明提出的滤波方法所得的地面数据是用体元表达的，而参考数据中则是离散的LIDAR激光点的表达形式。为和参考数据做对比进而评价本发明提出的方法精度，本实施例中，首先统计本方法所检测的地面体元内的激光点，然后和参考数据做对比，进而用I类误差(将地面点错分为非地面点比例)、II类误差(将非地面点错分为地面点比例)、总误差(错分的激光点的比例)及Kappa系数等指标定量评价滤波算法精度。

本实施例中，邻域尺度分别为6、18、26、51、56及64时，应用本发明方法对15个测试数据进行滤波，对应的滤波结果的Kappa系数如表2所示。该表中的数据旨在考查不同领域尺寸对滤波结果精度的影响，并由此确定最优邻域尺寸。

表2不同邻域尺寸的滤波方法的Kappa系数

由表2可知，6、18、26、51、56和64邻域的平均Kappa系数分别为71.64％、78.73％、80.81％、84.34％、83.29％和83.28％。这表明：(1)51邻域对应最大的Kappa系数，因此，从Kappa系数指标来看，51邻域是最佳邻域尺寸。(2)邻域尺寸的增加并不意味着滤波精度的必然提高。这是因为：地面目标的信息可以通过灰度三维体元数据集中定义的连通性、灰度和地形坡度强度相似性来传播。以6邻域为例，地面目标的信息只能向其上、下、左、右、前、后6个体元传播，由此导致种子体元只能向其周围平坦地形传播。如果能增加邻域尺寸，如26邻域，其对比6邻域增加的传播方向可能会纳入更多的地面体元，从而提高滤波的准确性。这可以解释为什么26邻域对比6邻域的滤波算法的Kappa系数更高。但是，如果邻域尺度过大，则可能将一些非地面体元错分为地面体元，从而导致滤波准确性的降低。这可以解释为何56邻域对比51邻域的精度反而有所下降。

本实施例中，应用本发明方法以参考数据为标准对15个样本数据的51邻域尺寸下的滤波精度进行的定量评价，滤波结果的精度如表3所示。

表3滤波结果的精度

样本	I类误差(％)	II类误差(％)	总误差(％)	Kappa系数(％)
					samp11	13.49	19.80	16.18	66.86
samp12	3.72	3.54	3.63	92.74
					samp21	1.12	3.47	1.64	95.25
samp22	3.32	5.53	4.01	90.71
					samp23	2.78	5.37	4.00	91.96
samp24	2.32	9.38	4.26	89.21
					samp31	0.94	1.94	1.40	97.18
samp41	2.11	3.99	3.05	93.89
					samp42	1.24	5.99	0.79	98.10
samp51	3.51	15.85	6.20	81.62
					samp52	1.03	31.12	4.19	75.28
samp53	1.58	48.23	3.46	52.93
					samp54	1.91	6.88	4.58	90.83
samp61	0.33	31.62	1.40	76.30
					samp71	2.36	20.28	4.39	77.95

由表3可知：I类误差、II类误差和总误差的误差分布范围为0～14％、1～49％、0～17％；滤波的平均Kappa系数可达84.34％。从而验证了本发明提出的方法的有效性。

本实施例中，应用本发明方法和Wang等(Liying Wang，Yan Xu，and Yu Li.AerialLIDAR point cloud voxelization with its 3D ground filtering application[J].Photogrammetric Engineering and Remote Sensing，2017，83(2)：95-107)提出的基于二值体元模型的三维滤波方法的总误差对比结果如表4所示。

表4本发明方法与基于二值体元模型的三维滤波方法的总误差比较

样本	Wang等方法的总误差(％)	本发明方法的总误差(％)
			samp11	19.49	16.18
samp12	4.02	3.63
			samp21	2.05	1.64
samp22	4.97	4.01
			samp23	5.91	4.00
samp24	6.34	4.26
			samp31	1.58	1.40
samp41	2.17	3.05
			samp42	1.07	0.79
samp51	8.09	6.20
			samp52	4.90	4.19
samp53	3.46	3.46
			samp54	5.62	4.58
samp61	1.93	1.40
			samp71	5.42	4.39
平均	5.13	4.21

由表4可知：本发明方法的总误差明显低于已有的基于二值体元模型的三维滤波算法，即本发明方法精度更高。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (3)

1.一种基于灰度体元模型的机载LIDAR三维滤波方法，其特征在于：包括以下步骤：

步骤1：读取原始机载LIDAR点云数据，形成原始机载LIDAR点云数据集；

步骤2：将原始机载LIDAR点云数据集规则化为灰度体元模型，其中灰度为体元内激光点的平均强度的离散化表示，具体方法为：

步骤2.1：从原始机载LIDAR点云数据集中剔除异常数据，得到剔除异常数据集；

步骤2.2：将剔除异常数据集规则化为灰度体元模型；

步骤3：基于三维连通区域构建理论，对灰度体元模型的地面体元进行检测，具体方法为：

步骤3.1：基于地面点的局部高程最低特性，从灰度体元模型中高程最低的非0值体元作为地面种子体元集合；

步骤3.1.1：在水平方向上，利用设定的格网尺寸对灰度体元模型进行分块，并取各块内高程最低的非0值体元为种子体元；

步骤3.1.2：在各块内搜寻与地面种子的高度差小于高度差阈值T_e的非0值体元为地面种子体元，获得加密后的地面种子体元集合V_s，其中，s＝1,2,…；

步骤3.2：标记地面种子体元及与其三维连通且坡度、灰度均接近的体元构成的三维连通区域为地面体元数据集，完成基于灰度体元模型的机载LIDAR三维滤波；

步骤3.2.1：计算灰度体元模型中非0值体元的灰度值的频率直方图，用高斯混合模型拟合灰度频率直方图，确定与地面目标对应的灰度值分布范围；

步骤3.2.2：对任一地面种子体元V_s，遍历灰度体元模型中与当前地面种子体元V_s三维连通、灰度值位于地面目标对应的灰度范围内且局部地形坡度小于坡度阈值T_s的所有未标记体元，并标记为L_g，直至标记完所有地面种子体元V_s的三维连通区域，即地面体元集合；

步骤3.2.2所述的坡度阈值的确定是自适应的，具体方法为：对当前种子体元V_s，检测其空间邻域内是否含有已标记的地面体元，若有，则将已有地面体元之间的最大坡度值确定为该邻域内的地形坡度阈值T_s；否则，将地形坡度阈值设置为90°。

2.根据权利要求1所述的一种基于灰度体元模型的机载LIDAR三维滤波方法，其特征在于：所述步骤2.1的具体方法为：

步骤2.1.1：统计原始机载LIDAR点云数据集中各个激光点高程值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.2：确定与真实地形及地物对应的最高高程阈值T_he和最低高程阈值T_le；

步骤2.1.3：针对原始机载LIDAR点云数据集中各个激光点，若其高程值高于最高高程阈值T_he或低于最低高程阈值T_le，则该激光点为高程异常数据，进行剔除，否则保留该激光点，获得剔除高程异常数据集；

步骤2.1.4：统计剔除高程异常数据集中各激光点的强度值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.5：确定与真实地形及地物对应的最高强度阈值T_hi和最低强度阈值T_li；

步骤2.1.6：针对剔除高程异常数据集中各个激光点，若其强度值高于最高强度阈值T_hi或低于最低强度阈值T_hi，则该激光点为强度异常数据，进行剔除，否则保留该激光点，最终获得剔除高程及强度异常数据集。

3.根据权利要求1所述的一种基于灰度体元模型的机载LIDAR三维滤波方法，其特征在于：所述步骤2.2的具体方法为：

步骤2.2.1：用剔除异常数据集的轴向平行包围盒表示数据集的空间范围；

步骤2.2.2：计算体元分辨率即体元大小，x、y、z方向上的分辨率根据剔除异常数据集中激光点的平均点间距确定；

步骤2.2.3：依据x、y、z方向上的分辨率对剔除异常数据集的轴向平行包围盒进行划分得到三维格网，将每一个三维格网单元称为体元；

步骤2.2.4：将剔除异常数据集中各个激光点映射到三维格网，进而根据体元中包含的激光点的强度属性为各体元赋值，最后将各体元值离散化到{0,…,255}，得到灰度体元模型。

Images