CN110335352A - 一种机载激光雷达点云的双基元多分辨率层次滤波方法 - Google Patents

Abstract

本发明公开了一种机载激光雷达点云的双基元多分辨率层次滤波方法。该滤波方法首先基于表面平滑性将原始点云分割为点云块集(块基元)和散点集(点基元)，然后通过改进的多分辨率层次滤波方法判断每个点的属性(地面点或非地面点)，最后根据每个点的属性将点云块分为地面块和非地面块，将散点分为地面点和非地面点。由于本发明方法融入了点云分块算法，然后在判断原始点云中每个点的属性时，引入了M估计进行平面拟合，以替代现有滤波方法中的内插地面参考面，从而提高点云的分类精度；最后利用每个点的属性分别判断点云块和散点的类别，由此能够保留更多的地面点，从而提高原始点云的滤波精度。

Description

一种机载激光雷达点云的双基元多分辨率层次滤波方法

技术领域

本发明涉及一种机载激光雷达点云的双基元多分辨率层次滤波方法。

背景技术

随着激光雷达(Light Detection And Ranging，简称LiDAR)技术的发展，机载激光雷达已逐渐成为收集地面点信息的主流工具。与传统的摄影测量方法相比，机载激光雷达可以穿透植被，能准确获取林地区域地形信息。但机载激光雷达获取的原始点云数据中包含建筑物、桥梁、车辆和植被等各种非地面信息，故在构建数字高程模型(DigitalElevation Model，简称DEM)前需要将这些非地面点滤除。目前常用的滤波方法按照原理可分为基于内插、基于坡度、基于分割和基于形态学的滤波。其中，基于内插的滤波以其较好的鲁棒性被广泛应用于处理地形复杂的森林地区点云。例如，陈传法等提出了一种多分辨率分层机载LiDAR点云滤波方法((Multi-resolution Hierarchical Classification，简称MHC)，该滤波方法在每一层次上通过薄板样条函数(Thin Plate Spline，简称TPS)迭代构造栅格平面，通过预设的阈值迭代完成地面点分类，该滤波方法在地形平坦和场景简单区域可获得较好的滤波效果。然而由于该滤波方法是通过点到TPS内插生成的地面参考面距离来判断该点是否为地面点，因此，该判断方法在茂密植被区域可能会引入部分非地面点，降低原始点云的滤波精度，可见，现有的MHC滤波方法在地形断裂和场景复杂区域滤波精度有待进一步提高。

发明内容

本发明的目的在于提出一种机载激光雷达点云的双基元多分辨率层次滤波方法，以提高原始点云的滤波精度，在尽可能滤除非地面点的同时能有效地保留地形特征点。

本发明为了实现上述目的，采用如下技术方案：

一种机载激光雷达点云的双基元多分辨率层次滤波方法，包括如下步骤：

I.基于表面平滑性将原始点云划分为点云块集和散点集，具体过程如下：

I.1.将原始点云划分为分辨率为Dh的网格，Dh大于测区内最大建筑物尺寸；

I.2.通过局部最小值法选择网格最低点；

I.3.以步骤I.2中网格最低点组成的网格最低点集作为查找数据集，搜索点云每个点在网格最低点集中的r个邻近点，并通过这r个邻近点计算该点的法向量和平面拟合残差；

其中，r为大于或等于5的正整数；

I.4.根据法向量、平面拟合残差和距离限制进行区域增长，将点云分为点云块和散点；

II.通过改进的多分辨率层次滤波方法判断原始点云中每个点的属性，判断过程如下：

II.1.将步骤I.2选择的网格最低点标记为临时地面点；

II.2.标记外层迭代次数iter＝1；

II.3.将原始点云划分为分辨率为h的网格；

II.4.利用局部最小值法从当前已经选择出的临时地面点中选择网格最低点；

II.5.利用TPS插值种子点法插值出没有步骤II.4网格最低点落入的网格中心点高程值；

II.6.将网格中心点和步骤II.4中的临时地面点共同组成参考地面点集；

从参考地面点集中搜索每个网格中心点的n个邻近点，利用这n个邻近点，基于M估计的方法拟合出每个网格对应的地面参考面；

II.7.计算待分类点到m个邻近地面参考面的垂直距离；

II.8.如果待分类点到m个邻近地面参考面的垂直距离中，至少有c个垂直距离小于给定的残差项阈值t，则将待分类点标记为临时地面点，记录新添加的临时地面点数量；

其中，c为大于或等于5的正整数，n、m均为正整数，且n>m>c；

II.9.判断新添加的临时地面点数量是否大于0：

若是，则更新当前的临时地面点数量，并转到步骤II.4，否则，转到步骤II.10；

II.10.令iter＝iter+1，判断更新后的iter值是否小于或等于3：

若是，设定新阈值：h＝0.5h，t＝t+0.1，并转到步骤II.3；若不是，转到步骤III；

III.根据原始点云中每个点的属性，将点云块集中的点云块划分为地面块和非地面块，将散点集中的散点划分为地面点和非地面点。

优选地，步骤II.6中，基于M估计的方法拟合出每个网格对应的地面参考面的过程为：

设地面参考面方程为z＝β₀+β₁x+β₂y，则地面参考面的描述参数为：

β＝[β₀ β₁ β₂]^T，β₀、β₁、β₂为平面方程的系数；

其中，(x_i、y_i)为第i个邻近点的xy坐标，1≤i≤n，n为邻近点数量；

z＝[z₁ z₂ … z_n]^T，z_i表示第i个邻近点的z坐标；

M估计的目标函数表达为：

ρ采用Huber影响函数：

其中，e_i表示第i个邻近点的初始残差，即e_i＝z_i-β₀-β₁x_i-β₂y_i，k为限定常数；

则最终拟合公式表达为：

其中，W为进行拟合的所有邻近点的权重组成的权重矩阵，

权重函数表达为：

其中，u_i为标准化残差：

如上所述，M估计的方法具体步骤为：

①选取最小二乘估计的为迭代初始值，求出每个邻近点初始残差e_i；

②基于公式(3)求出每个邻近点的初始权重

③利用求得代替并求出每个邻近点的新残差：

e_i ^(l)＝z_i-β₀ ^(l)-β₁ ^(l)x_i-β₂ ^(l)y_i；

其中，l为大于或等于1的整数，表示第l次迭代过程中第i个邻近点的残差；β₀ ^(l)、β₁ ^(l)、β₂ ^(l)为第l次迭代过程中平面方程的系数；

④重复步骤②和步骤③，依次迭代计算

当相邻两步的平面系数的差的绝对值的最大值小于预先设定的标准误差时迭代结束，即：

此时，所得到的即为要求的平面方程系数，即所求平面方程为：

z＝β₀ ^(l)+β₁ ^(l)x+β₂ ^(l)y。

优选地，步骤I.3中待分类点的法向量和平面拟合残差的计算过程为：

以r个邻近点计算得到的法向量作为该点的法向量；

以点到r个邻近点拟合平面的垂直距离作为该点的平面拟合残差。

优选地，步骤III具体为：

若点云块中至少一半以上的点标记为临时地面点，则将该点云块中的所有点标记为地面点，同时将该点云块划分为地面块；

否则，将该点云块中的所有点标记为非地面点，同时将该点云块划分为非地面块；

若散点标记为临时地面点，则将该散点划分为地面点，否则该散点划分为非地面点。

本发明具有如下优点：

如上所述，本发明提出了一种机载激光雷达点云的双基元多分辨率层次滤波方法，该滤波方法融入了点云分块算法，即首先将原始点云划分为点云块集(块基元)和散点集(点基元)，然后在判断原始点云中每个点的属性时，引入了M估计进行平面拟合，以替代现有滤波方法中的内插地面参考面，从而提高点云的分类精度；最后利用每个点的属性分别判断点云块和散点的类别，由此能够保留更多的地面点，从而提高原始点云的滤波精度。

附图说明

图1为本发明中机载激光雷达点云的双基元多分辨率层次滤波方法的流程示意图。

图2为本发明实施例中原始点云划分流程示意图。

图3为本发明实施例中改进的区域增长算法中平面拟合示意图。

图4为本发明实施例中改进的区域增长算法与传统区域增长算法分块对比示意图。

图5为本发明实施例中TPS地面种子点插值示意图。

图6为本发明实施例中点到拟合平面和插值平面距离的对比示意图。

图7为本发明实施例的实验中五组数据的数字地表模型示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

结合图1所示，一种机载激光雷达点云的双基元多分辨率层次滤波方法，包括如下步骤：

I.基于表面平滑性将原始点云划分为点云块集(块基元)和散点集(点基元)。

经典的平滑性约束点云分块算法是以空间距离、法向量和平面拟合残差作为限制条件进行区域增长。其中，法向量和平面拟合残差显著影响分块结果，特别是在茂密植被地区。

如果法向量和平面拟合残差计算不准确，则容易出现欠分割或者过分割问题。

鉴于机载激光雷达点云数据为2.5维，在同一平面位置植被点要明显高于地面点，因此本发明实施例在点云块集和散点集的划分过程中采用如下步骤，如图2所示。

I.1.将原始点云划分为分辨率为Dh的网格，其中，Dh大于测区内最大建筑物尺寸。

通过该分辨率的网格选择网格最低点，可以保证选择到的点一定为地面点，从而保证拟合的初始参考面一定是最大程度接近地形的参考面，避免后续引入过多的非地面点。

I.2.通过局部最小值法选择网格最低点。

针对传统的移动窗口局部最小值方法选取的最低点极易受粗差异常值、多路径效应以及激光雷达系统影响，本发明采用改进的扩展局部最小值方法，即：

选择最低点和第二低点作为最低点候选点，然后判断最低点与第二低点之间的高差，若高差小于给定的阈值th(如th＝1m)则将最低点作为最低点，迭代结束；

若高差大于给定的阈值则将最低点作为低位粗差点剔除，并再次选择剩余点中的最低点和第二低点，依次迭代直至最低点与第二低点的高差小于阈值th；

上述改进的扩展局部最小值方法是为了确定选择到的最低点为地面点。

I.3.以步骤I.2中网格最低点组成的网格最低点集作为查找数据集，搜索点云每个点在网格最低点集中的r个邻近点，并通过这r个邻近点计算该点的法向量和平面拟合残差。

计算该点的法向量和平面拟合残差的具体过程如下：

以r个邻近点计算得到的法向量作为该点的法向量；

以待分类点到r个邻近点拟合平面的垂直距离作为该点的平面拟合残差。

本实施例中r为大于5的正整数，例如r＝12。

I.4.根据法向量、平面拟合残差和距离限制进行区域增长，将点云分为点云块和散点。

图3给出了本发明实施例中区域增长算法中平面拟合示意图。

由上述平面拟合示意图可以看出，此方法能够保证杂乱的植被点法向量平行于地面点法向量，同时植被点残差远大于地面点残差，由此不会把植被点误分到地面块中。

如图4(a)为传统的传统区域增长算法分块结果示意图，而图4(b)为本发明实施例中改进的区域增长算法的分块结构示意图。由上述结果对比不难看出：

改进的区域增长算法能够准确的将地面点分为一个块集合，而植被点被分为散点。

本发明针对机载激光雷达点云的数据特征，对区域增长分块算法中的点云法线计算方法进行改进，不再使用传统的邻近点拟合，而是利用预先选择的邻近网格最低点进行平面拟合。由此能够保证各点法线方向和残差的稳定性，避免区域增长过程中引入过多非地面点。

II.通过改进的多分辨率层次滤波方法判断原始点云中每个点的属性。

本发明实施例中整个滤波过程分为3个层次，从低层到高层网格分辨率h和残差阈值t同时增大。在第一层通过从原始点云数据中选择的种子点来拟合地面参考面。

由于种子点个数限制，地面参考面精度不高。为了避免太多非地面点被误添加为种子点，残差阈值t的初值应该足够小，一般取0到1之间的数值。

经过几次迭代后更多的地面点被添加为种子点，下一层拟合生成的地面参考面更加接近真实地面。随着地面参考面的精度提高，适当增大残差阈值t，以尽可能多的添加地面点。

原始点云中每个点的属性具体判断过程如下：

II.1.将步骤I.2选择的网格最低点标记为临时地面点。

II.2.标记外层迭代次数iter＝1。

II.3.将原始点云划分为分辨率为h的网格。

II.4.利用局部最小值法从当前已经选择出的临时地面点中选择网格最低点。

该步骤II.4中采用的局部最小值法，也采用上文步骤I.2中提到的改进的扩展局部最小值方法，以确定选择到的最低点为临时地面点，具体过程不再详细赘述。

在该步骤II.4中再次选择网格最低点的原因有以下两个方面：

1.若不再次进行网格最低点选择，则通过网格中心点邻近的临时地面点进行地面参考面拟合时，可能会由于各个邻近点之间的距离过小而造成拟合计算过程中的矩阵奇异；

2.由于初始几次迭代的地面参考面与真实地形存在差异，故通过该参考面选择的地面点中可能会含有部分非地面点，若直接用网格中心的邻近点进行参考面拟合，便会由于非地面点的存在而影响参考面精度。通常非地面点会高于地面点，而通过再次选择网格最低点便可以尽可能的避免非地面点对参考面拟合的影响。

II.5.利用TPS插值种子点法插值出没有步骤II.4网格最低点落入的网格中心点高程值。

在建筑物区域或地面点密度不高时，部分网格的n个邻近点通常远离网格中心。因此，若直接用这些邻近点进行平面拟合，会使拟合平面精度较低。

为了提高精度，本发明在没有种子点的网格借助TPS插值种子点的方法，如图5所示。

II.6.将网格中心点和步骤II.4中临时地面点共同组成参考地面点集。

从参考地面点集中搜索每个网格中心点的n个邻近点，利用这n个邻近点，基于M估计的方法拟合出每个网格对应的地面参考面。

本发明实施例中基于M估计的平面拟合替代内插地面参考面的原因在于：

传统的多分辨率层次滤波方法通过点到TPS内插生成的地面参考面距离来判断该点是否为地面点，该判断方法在茂密植被区域可能会引入部分非地面点。

如图6所示，当非地面点P₁位于网格边缘时，即使它到内插网格平面的距离d_dem1很小，但到真实地面的距离d_plane1却很大。同样，当斜坡上的地面点P₂位于网格边缘时，它到内插网格平面的距离d_dem2较大，但到真实地面的距离d_plane2却很小。

因此，本发明以M估计拟合出的平面代替内插地面参考面，以提高点云分类精度。

由于初始地面种子点不可避免的存在非地面点(即粗差)，为了尽可能的减少这些粗差对拟合平面的影响，同时考虑计算的时间复杂度和鲁棒性，本发明实施例优选采用具有高抗差性的M估计进行平面拟合。M估计稳健回归是通过迭代加权最小二乘估计回归系数，即根据前一次回归残差确定样本权重，以此减小异常值对最终结果的影响。

基于M估计的方法拟合出每个网格对应的地面参考面的过程为：

设地面参考面方程为z＝β₀+β₁x+β₂y，则地面参考面的描述参数为：

β＝[β₀ β₁ β₂]^T，其中，β₀、β₁、β₂为平面方程的系数。

其中，(x_i、y_i)为第i个邻近点的xy坐标，1≤i≤n，n为邻近点数量。

z＝[z₁ z₂ … z_n]^T，z_i表示第i个邻近点的z坐标。

M估计的目标函数表达为：

ρ采用Huber影响函数：

其中，e_i表示第i个邻近点的初始残差，即e_i＝z_i-β₀-β₁x_i-β₂y_i。

k为限定常数，本实施例中取值k＝2.5。

则最终拟合公式表达为：

其中，W为进行拟合的所有邻近点的权重组成的权重矩阵，

权重函数表达为：

其中，u_i为标准化残差：

如上所述，M估计的方法具体步骤为：

①选取最小二乘估计的为迭代初始值，求出每个邻近点初始残差e_i。

②基于公式(3)求出每个邻近点的初始权重

③利用求得代替并求出每个邻近点的新残差：

e_i ^(l)＝z_i-β₀ ^(l)-β₁ ^(l)x_i-β₂ ^(l)y_i。

其中，l为大于或等于1的整数，e_i ^(l)表示第l次迭代过程中第i个邻近点的残差；β₀ ^(l)、β₁ ^(l)、β₂ ^(l)为第l次迭代过程中平面方程的系数。

④重复步骤②和步骤③，依次迭代计算

当相邻两步的平面系数的差的绝对值的最大值小于预先设定的标准误差时迭代结束，即：

此时，所得到的即为要求的平面方程系数，即所求平面方程为：

z＝β₀ ^(l)+β₁ ^(l)x+β₂ ^(l)y。

通过当前网格中心点邻近的n个种子点拟合当前网格的地面参考面。

本实施例中n为大于5的整数，例如n＝12。

本发明将传统多分辨率层次滤波方法中用插值网格平面来进行点云类别判断的方法更改为用拟合平面的方法，由于拟合平面带有方向特征，因此能够更好的贴合真实地面。

II.7.计算待分类点到m个邻近地面参考面的垂直距离，例如，m＝9。

这m个邻近地面参考面分别为当前待分类点垂直落入的网格平面，垂直落入网格的上、右上、右、右下、下、左下、左、左上网格平面。

II.8.如果待分类点到m个邻近地面参考面的垂直距离中，至少有c个垂直距离小于给定的残差项阈值t，则将待分类点标记为临时地面点，并记录新添加的临时地面点数量。

其中，c为大于或等于5的正整数，n、m均为正整数，且n>m>c。

II.9.判断新添加的临时地面点数量是否大于0：

若是，则更新当前的临时地面点数量，并转到步骤II.4，否则，转到步骤II.10。

临时地面点的具体更新过程为：利用本次迭代过程中当前已经选择出的临时地面点与新添加的临时地面点的和作为下一次迭代过程中当前已经选择出的临时地面点。

在地形起伏区域，当在同一网格分辨率下只进行一次迭代时，由于临时地面点数量较少，拟合出的地面参考面不能很好地贴合真实地形，故只能添加部分离当前地面参考面较近的地面点。而通过迭代，由于每次迭代都会添加一些新的临时地面点，此时再拟合出的地面参考面将更贴合真实地形，故可以添加更多的临时地面点。通过不断迭代，直到当前分辨率下不能添加新的地面点为止，即此分辨率下地面参考面已极大程度的接近真实地形。此时须通过缩小网格尺寸(增大网格分辨率)来进一步迭代拟合，即进入下一层iter的迭代过程。

II.10.令iter＝iter+1，判断更新后的iter值是否小于或等于3：

若是，设定新阈值：h＝0.5h，t＝t+0.1，并转到步骤II.3；若不是，转到步骤III。

III.根据原始点云中每个点的属性，将点云块集中的点云块划分为地面块和非地面块，将散点集中的散点划分为地面点和非地面点，具体过程如下：

若点云块中至少一半以上的点标记为临时地面点，则将该点云块中的所有点标记为地面点，同时将该点云块划分为地面块；

否则，将该点云块中的所有点标记为非地面点，同时将该点云块划分为非地面块；

若散点标记为临时地面点，则将该散点划分为地面点，否则该散点划分为非地面点。

由于本发明方法融入了点云分块算法，然后，在判断原始点云中每个点的属性时，引入了M估计进行平面拟合，以提高点云的分类精度；最后，利用每个点的属性分别判断点云块和散点的类别，由此能够保留更多的地面点，从而提高原始点云的滤波精度。

下面对本发明提出的双基元多分辨率层次滤波方法的有效性进行验证。

本发明从OpenTopography网站下载了五组1:1000比例尺标准图幅范围的山地地形数据对算法进行验证。其中，图7a(数据1)和图7b(数据2)主要为覆盖茂密植被的斜坡地形，图7c(数据3)、图7d(数据4)和图7e(数据5)为包含建筑物和植被的斜坡地形。

以上五组数据的平均点间距均为1.25m。每组数据都在其原始滤波结果的基础上结合影像进行人工识别，从而保证标准地面点和非地面点的精确分类。除了本发明方法外，经典多分辨率层次滤波算法(MHC)以及近些年提出的新滤波方法也用于数据滤波，如布料模拟(CSF)、最大局部斜率(MLS)、自适应TIN(ATIN)和渐进形态学(Morph)。

以上各个方法滤波结果主要借助I类、II类和总误差以及kappa系数进行精度评价。

首先，利用本发明方法处理以上五组数据，所用的参数如表1所示。除了初始分辨率外，本发明方法对其它参数取值并不敏感，表明滤波结果具有较高稳定性。

表1本发明方法处理数据时用的参数

本发明方法对以上五组数据处理精度如表2所示。通过表2不难发现，本发明方法的总误差小于5％，Kappa系数大于88％，且其中4组数据的Kappa系数大于92％。因此，本发明滤方法在覆盖茂密植被的斜坡区域具有较高的精度，且具有很好的稳健性。

表2本发明方法滤波误差和Kappa系数(％)

本发明方法与其它六种滤波方法滤波结果如表3所示。由表3能够看出，本发明方法滤波精度明显优于其他滤波方法，其平均总误差为3.55％，平均Kappa系数为92.62％。较经典的多分辨率层次滤波方法(MHC)总误差提高约57％，Kappa系数提高约10％。

由表3不难看出，本发明方法在五组数据中，数据3的处理精度最高，Kappa系数为95.37％，数据4的处理精度最低，Kappa系数为88.88％，本发明方法具有很高的稳健性。

表3 DMHF与其它6种经典滤波方法的总误差和Kappa系数(％)

综上，本发明方法在精度和稳健性上较其它经典滤波方法均有所提高。

为了进一步比较各种滤波方法的效果，本发明还对数据3处理结果进行分析。其中，数据3包含典型的山脊、山谷地形，且在山坡上分布有明显的植被和建筑物。

通过对各种方法处理结果分析发现：

本发明方法分成地面点的非地面点最少，且分类出的真实地面点最多。

而MHC算法、ATIN算法和Morph算法在斜坡地形上保留了过多的非地面点，CSF算法、Morph算法和MLS算法在斜坡区域损失了过多的地面点。

综上，本发明方法在尽可能滤除非地面点的同时能有效地保留地形特征点。

此外，本发明还对各种算法生成的DEM进行了分析。

通过对比表明，本发明方法生成的DEM表面平滑、地形细节完整，且跟真实DEM具有最高的相似度。其余五种算法生成的DEM表面均有多余的尖刺，其中：

MHC算法和Morph算法仅在斜坡区域存在尖刺，CSF算法、ATIN算法和MLS算法在平整区域也出现尖刺，而且CSF算法在斜坡区域还损失了过多的地形细节(斜坡部分有明显的过度平滑现象)。由此可见，本发明滤波后地面点生成的DEM很好的还原真实地形。

通过以上结果分析表明，本发明方法可有效滤除原始点云中的建筑物、植被和其他非地面点，能较好的适应不同地形环境，在覆盖茂密植被的斜坡区域具有突出滤波优势。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (4)

1.一种机载激光雷达点云的双基元多分辨率层次滤波方法，其特征在于，

包括如下步骤：

I.基于表面平滑性将原始点云划分为点云块集和散点集，具体过程如下：

I.1.将原始点云划分为分辨率为Dh的网格，Dh大于测区内最大建筑物尺寸；

I.2.通过局部最小值法选择网格最低点；

I.3.以步骤I.2中网格最低点组成的网格最低点集作为查找数据集，搜索点云每个点在网格最低点集中的r个邻近点，并通过这r个邻近点计算该点的法向量和平面拟合残差；

其中，r为大于或等于5的正整数；

I.4.根据法向量、平面拟合残差和距离限制进行区域增长，将点云分为点云块和散点；

II.通过改进的多分辨率层次滤波方法判断原始点云中每个点的属性，判断过程如下：

II.1.将步骤I.2选择的网格最低点标记为临时地面点；

II.2.标记外层迭代次数iter＝1；

II.3.将原始点云划分为分辨率为h的网格；

II.4.利用局部最小值法从当前已经选择出的临时地面点中选择网格最低点；

II.5.利用TPS插值种子点法插值出没有步骤II.4网格最低点落入的网格中心点高程值；

II.6.将网格中心点和步骤II.4中的临时地面点共同组成参考地面点集；

从所述参考地面点集中搜索每个网格中心点的n个邻近点，利用这n个邻近点，基于M估计的方法拟合出每个网格对应的地面参考面；

II.7.计算待分类点到m个邻近地面参考面的垂直距离；

II.8.如果待分类点到m个邻近地面参考面的垂直距离中，至少有c个垂直距离小于给定的残差项阈值t，则将所述待分类点标记为临时地面点，记录新添加的临时地面点数量；

其中，c为大于或等于5的正整数，n、m均为正整数，且n>m>c；

II.9.判断新添加的临时地面点数量是否大于0：

若是，则更新当前的临时地面点数量，并转到步骤II.4，否则，转到步骤II.10；

II.10.令iter＝iter+1，判断更新后的iter值是否小于或等于3：

若是，设定新阈值：h＝0.5h，t＝t+0.1，并转到步骤II.3；若不是，转到步骤III；

III.根据原始点云中每个点的属性，将所述点云块集中的点云块划分为地面块和非地面块，将所述散点集中的散点划分为地面点和非地面点。

2.根据权利要求1所述的双基元多分辨率层次滤波方法，其特征在于，

所述步骤II.6中，基于M估计的方法拟合出每个网格对应的地面参考面的过程为：

设地面参考面方程为z＝β₀+β₁x+β₂y，则地面参考面的描述参数为：

β＝[β₀ β₁ β₂]^T，β₀、β₁、β₂为平面方程的系数；

其中，(x_i、y_i)为第i个邻近点的xy坐标，1≤i≤n，n为邻近点数量；

z＝[z₁ z₂…z_n]^T，z_i表示第i个邻近点的z坐标；

M估计的目标函数表达为：

ρ采用Huber影响函数：

其中，e_i表示第i个邻近点的初始残差，即e_i＝z_i-β₀-β₁x_i-β₂y_i，k为限定常数；

则最终拟合公式表达为：

其中，W为进行拟合的所有邻近点的权重组成的权重矩阵，

权重函数表达为：

其中，u_i为标准化残差：

如上所述，M估计的方法具体步骤为：

①选取最小二乘估计的为迭代初始值，求出每个邻近点初始残差e_i；

②基于公式(3)求出每个邻近点的初始权重

③利用求得代替并求出每个邻近点的新残差：

e_i ^(l)＝z_i-β₀ ^(l)-β₁ ^(l)x_i-β₂ ^(l)y_i；

其中，l为大于或等于1的整数，e_i ^(l)表示第l次迭代过程中第i个邻近点的残差；β₀ ^(l)、β₁ ^(l)、β₂ ^(l)为第l次迭代过程中平面方程的系数；

④重复步骤②和步骤③，依次迭代计算

当相邻两步的平面系数的差的绝对值的最大值小于预先设定的标准误差时迭代结束，即：

此时，所得到的即为要求的平面方程系数，即所求平面方程为：

z＝β₀ ^(l)+β₁ ^(l)x+β₂ ^(l)y。

3.根据权利要求1所述的双基元多分辨率层次滤波方法，其特征在于，

所述步骤I.3中，待分类点的法向量和平面拟合残差的计算过程为：

以r个邻近点计算得到的法向量作为该点的法向量；

以点到r个邻近点拟合平面的垂直距离作为该点的平面拟合残差。

4.根据权利要求1所述的双基元多分辨率层次滤波方法，其特征在于，

所述步骤III具体为：

若点云块中至少一半以上的点标记为临时地面点，则将该点云块中的所有点标记为地面点，同时将该点云块划分为地面块；

否则，将该点云块中的所有点标记为非地面点，同时将该点云块划分为非地面块；

若散点标记为临时地面点，则将该散点划分为地面点，否则该散点划分为非地面点。