CN103745436B - 基于区域预测的LiDAR点云数据形态学滤波方法 - Google Patents

Abstract

本发明涉及基于区域预测的LiDAR点云数据形态学滤波方法。方法包括以下步骤：步骤一，构建滤波初始表面模型；步骤二，去除初始表面模型的粗差点，得到滤波表面模型；步骤三，对滤波表面模型进行分块：步骤四，滤波处理。本发明的方法克服了现有技术中形态学滤波不能很好对复杂场景滤波的缺点,以及使用恒定的人工设置的地形参数进行滤波处理而导致的自适应不强的问题，利用了分块原理和使用各个分块区域预测地形参数s能够根据地形的起伏情况调整高差阈值等优点，从而能够自适应地对点云数据进行滤波，最终获得较好的滤波结果。

Description

基于区域预测的LiDAR点云数据形态学滤波方法

技术领域

本发明属于机载LiDAR数据处理技术领域，更进一步涉及机载LiDAR数据滤波技术领域中基于区域预测的形态学点云数据滤波算法。可应用于机载LiDAR扫描地面得到的点云数据，可以有效的将点云数据滤波并提取DEM。

背景技术

机载LiDAR是一种主动式遥感技术，能够扫描大面积区域直接获取地面三维信息，其数据具有高精度、高密度、离散不规则等特点。

LiDAR点云数据包含有真实地面、建筑物和植被等多种信息，对点云数据进行滤波得到真实地面点，形成数字高程模型(DEM)具有重要意义。目前许多点云滤波算法采用基于形态学的滤波方法。

赵明波等人在文献“赵明波，何峻，田军生，付强，‘基于改进的渐进多尺度数学形态学的激光雷达数据滤波方法’[J]，光学学报，33(3)，292-301(2013)”中提出了改进的基于渐进形态学的激光雷达数据滤波算法。该方法的通过改进形态学开运算来处理大面积空白区域的点云数据，改进的开运算具体是腐蚀运算时把当前网格邻域内的空白区域高程值用最大高程值替代，同理，膨胀运算则是用最小高程值替代。该方法首先将原始激光雷达点云数据进行格网化并剔除粗差，然后使用窗口尺寸混合增长的方式和改进的开运算对格网进行迭代的形态学滤波，最终得到原始点云数据的分类。此算法中需人工设置大量参数，自适应性不强。

隋立春等人在文献“隋立春，张熠斌，柳艳等，‘基于改进的数学形态学算法的LiDAR点云数据滤波’[J]，测绘学报，39(4)，390-396(2010)” 中提出了改进的形态学LiDAR数据滤波方法。此方法对形态学算法进行扩展和改进，提出了针对不同地形特点的自适应滤波算法，具体的是在数学形态学“开”算子的基础上,提出了增加一个“带宽”参数用于点云数据滤波的方法，最后对再山包地形和山坡地形进行了判别和处理。该方法提出的动态滤波结合“带宽”参数的改进的数学形态学算法可以有效地用于LiDAR点云数据的滤波，并在一定程度上避免和降低了地面存在山包时错误滤波的后果。但是此算法没有对ISPRS提供的测试数据进行实验，无法与其他算法比较，并且算法在不同地形、不同高度和不同形状山包情况下还要进行进一步试验和完善。

综上所述，基于形态学的点云滤波方法在点云滤波方面虽然获得了较好的效果，但是一般来说，基于形态学的点云滤波方法中开运算后的高差阈值的确定需要许多参数，这些参数需根据不同地形设置，并且假设地形坡度恒定，针对复杂场景这是不合理的，导致算法自适应性不强，给对不同地形的区域滤波带来困难。

发明内容

本发明的目的是针对现有技术基于形态学的点云数据滤波方法中需人工设置滤波阈值导致自适应不强的缺点，提出了基于区域预测的形态学LiDAR点云数据滤波算法。

在本发明中，对点云数据的基于区域预测的滤波原理包括了四部分：点云数据格网化、去除粗差点、滤波准则和地面点的确定。

本发明充分利用形态学滤波算法自下而上、从局部出发扩展到全局的特点和速度快、效率较高的优点，同时利用标准差代表地形的起伏情况的特点，使用分块区域的标准差来预测地形坡度参数的取值以进行滤波。最终得到较 好的滤波结果，又能最大限度地保留实验区域的地形特征。

为此，本发明提供的基于区域预测的LiDAR点云数据形态学滤波方法包括以下步骤：

步骤一，构建滤波初始表面模型；

步骤二，去除初始表面模型的粗差点，得到滤波表面模型；

步骤三，对滤波表面模型进行分块：

按适当尺寸对滤波表面模型进行分块，其中分块n记为R_n，n＝1,2,3,…,N，N为分块的总数；该分块R_n中包含M个点云数据，其中点云数据m记为r_mn，m＝1,2,3,…,M；

步骤四，滤波处理：

对各分块分别进行滤波处理，其中对分块R_n，n＝1,2,3,…,N的处理过程如下：

从最小结构元素窗口尺寸开始执行以下操作，之后结构元素窗口的半径每增加一米，执行一次以下操作直至最大结构元素窗口尺寸：

分别对分块R_n中所有点云数据进行开运算，得到各点云数据开运算后的高程值，

分别计算各点云数据的高程值与其自身的开运算后的高程值的差值，差值大于e_n的点云数据标记为非地面点，且将各点云数据的高程值更新为开运算后的高程值，其中：e_n＝s_n×w×c，w为当前结构元素窗口的半径，c为格网边长，且1米≤c≤5米，s_n＝a×σ_n+b，σ_n为分块R_n中所有点云数据高程值的标准差，0.005≤a≤0.02，0.005≤b≤0.02；

滤波处理后未标记的点云数据为地面点。

进一步，本发明的基于区域预测的LiDAR点云数据形态学滤波方法还包括：

将步骤四确定的地面点插值得到数字高程模型，

接着分别确定所有点云数据的初始高程值与数字高程模型之间的垂直距离，其中垂直距离小于0.5米的点云数据为最终地面点。

优选的，步骤二中去除初始表面模型的粗差点是将所有点云数据的高程值取负值，然后使用高差阈值为5米、结构元素窗口半径为1米的形态学滤波来标记低位粗差点为非地面点。

优选的，步骤三中按尺寸大小为（20m×20m）～（50m×50m）对滤波表面模型进行分块。

优选的，步骤四中最小结构元素窗口半径为1米，最大结构元素窗口半径为5～30米。

优选的，对于分块时不能被整分的滤波表面模型，被分块后剩余的格网，采用如下方法进行滤波处理：

从最小结构元素窗口尺寸开始执行以下操作，之后结构元素窗口的半径每增加一米，执行一次以下操作直至最大结构元素窗口尺寸：

分别对剩余格网中所有点云数据进行开运算，得到各点云数据开运算后的高程值，

分别计算各点云数据的高程值与其自身的开运算后的高程值的差值，差 值大于e的点云数据标记为非地面点，且将各点云数据的高程值更新为开运算后的高程值，其中：e＝0.3×w×c，w为当前结构元素窗口的半径，c为格网边长，且1米≤c≤5米，滤波处理后未标记的点云数据为地面点。

与现有技术相比，本发明具有以下优点：

第一，本发明在对点云数据进行形态学滤波之前，结合先验知识判断实验区域的地貌特征，对规则化DSM按一定的尺寸进行了分块。由于实验区域一般地形复杂多变，将整个区域分块能够针对分块中的相似地形进行滤波，以便有效的滤去地物点保留地面点。克服了现有技术中形态学滤波不能很好对复杂场景滤波的缺点,使得使用本发明中的方法可以获得较好的滤波结果。

第二，本发明在对点云数据进行形态学滤波处理时，使用了各个分块区域预测的地形参数来确定高差阈值。由于标准差代表了地形的起伏情况，所以以标准差对地形参数进行预测从而得到高差阈值可以自适应地对数据滤波。克服了现有技术中的地形参数需要人工设定并为恒定，自适应不强的问题，使得使用本发明中的方法可以根据地形自适应地对点云数据进行滤波。

本发明应用于点云数据滤波处理中，能够有效去除地物点并保留地面点，尤其能有效减小第II类误差，同时维持较低的第I类误差和总误差。

附图说明

图1为实施例1中对sample12滤波的仿真图；

图2为实施例2中对sample31滤波的仿真图；

图3为实施例3中对sample53滤波的仿真图；

图4为实施例4中对sample54滤波的仿真图。

具体实施方式

本发明的构建初始表面模型可采用本领域常规的方法，具体包括对点云数据格网化：根据点云数据横坐标x和纵坐标y的最大值和最小值确定点云 数据的xy平面范围，并选择合适的格网尺寸构成规则格网，使得每个格网中至少包含一个点云数据，合适的格网尺寸一般为（1m×1m）～（5m×5m），且格网为正方形；如果在一个格网中没有点云数据，使用插值来得到缺失的数据点；如果一个格网中有多个点云数据，则选取格网中最低的点作为格网的取值。

本发明的去除初始表面模型的粗差点的方法可采用本领域传统的方法，优选以下方法：粗差点分为高位粗差点和低位粗差点，高位粗差点在滤波处理过程中可以被滤掉，而低位粗差点不能，需对低位粗差点进行单独处理，处理方法是将所有点云数据的高程值取负值，然后使用高差阈值为5米、结构元素窗口半径为1米的形态学滤波来标记低位粗差点为非地面点。

本发明中对滤波表面模型进行分块时，通过观察实验区域的地貌特征，对滤波表面模型按适合的尺寸进行分块，尽量保证块内的地形变化有相似性，并且分块面积不能小于最大建筑物面积以避免某个分块落入最大建筑物内，分块时的合适尺寸（20m×20m）～（50m×50m）。

本发明的滤波使用结构元素窗口尺寸线性增长的方式来迭代地进行滤波，初始结构元素窗口半径为1m，最大结构元素窗口尺寸大于实验区域中最大建筑物的尺寸，一般最大结构元素窗口半径为大于5～30m。

本发明所述的点云数据的高程值是指该点云数据的Z坐标值，对于点云数据r_mn来说，其高程值为Z值。各点云数据的初始高程值为其自身的Z坐标值。在滤波处理过程中各点云数据的最初高程值为其自身的Z坐标值，之 后，每进行一次开运算，其高程值便更新一次。

本发明所述的插值为文献参见文献：Thomas J.Pingel,Keith C.Clarke,WilliamA.McBride.An improved simple morphological filter for the terrainclassification of airborne LIDAR data[J].ISPRS Journal of Photogrammetry andRemote Sensing,2013(77):21-30.中公开的spring-metaphor inpainting technique。

本发明所述点云数据开运算后的高程值计算方法是首先对点云数据开运算指先对当前分块中所有点云数据进行腐蚀运算，再对分块中所有点云数据进行膨胀运算。具体包括：

首先对当前分块中所有点云进行腐蚀运算：

遍历每个点云数据，

设当前点云数据为r_mn，

为点云数据r_mn的坐标，为点云数据r_mn被腐蚀后的高程值，(x,y)为当前以r_mn为圆心的结构元素窗口内的点云数据坐标，Z(x,y)为点云数据(x,y)的高程值；

待对当前块中的所有点云数据进行完腐蚀运算后，分别对各点云数据进行膨胀运算：

Z₁(x,y)为当前以r_mn为圆心的结构元素窗口内的点云数据(x,y)被腐蚀后的高程值，为点云数据r_mn被膨胀后的高程值；

即为点云数据r_mn开运算后的高程值。

以下是发明人提供的具体实施例，以对本发明的技术方案作进一步解释说明。

实施例1：

该实施例中采用本发明的方法对ISPRS（国际摄影测量与遥感学会，International Society for Photogrammetry and Remote Sensing）提供的测试样本sample12进行滤波的仿真。

仿真条件是在MATLAB2010a软件下进行。

参照图1，对测试样本sample12，包含52119个点的点云数据进行仿真实验。图1(a)是初始表面模型（DSM），反应实验区域的原始地形；图1(b)显示样本数据中的地面点云分布，所有空白为应当滤去的地物点；图1(c)是本算法滤波后地面点点云的分布结果；图1(d)是最终生成的数字高程模型（DEM）。

从图1(a)可以看出，samp12的数据特征为山坡上的建筑物和植被的混合，其中的建筑物形状多为长矩形。从图1(c)和图1(b)可以看出滤波结果保持了地形的原有结构，将人工建筑物和植被滤去并很好的保留了地面点。

实施例2：

该实施例中采用本发明的方法对ISPRS提供的测试样本sample31进行滤波的仿真。

仿真2仿真条件是在MATLAB2010a软件下进行。

参照图2，对测试样本sample31，包含28862个点的点云数据进行仿真 实验。

图2(a)是DSM，反应实验区域的原始地形；图2(b)显示样本数据中的地面点点云分布，所有空白为应当滤去的地物点；图2(c)是本发明滤波后地面点点云的分布结果；图2(d)是最终生成的DEM。

从图2(a)可以看出，sample31主要是不连续的地形和低点，包含复杂建筑物、植被等目标。从图2(c)和图2(b)可以看出提出的滤波算法的效果较好。

实施例3：

该实施例中采用本发明的方法对ISPRS提供的测试样本sample53进行滤波的仿真。

仿真3仿真条件是在MATLAB2010a软件下进行。

参照图3，对测试样本sample53，包含34378个点的点云数据进行仿真实验。图3(a)是DSM，反应实验区域的原始地形；图3(b)显示样本数据中的地面点点云分布，所有空白为应当滤去的地物点；图3(c)是本算法滤波后地面点点云的分布结果；图3(d)是最终生成的DEM。

从图3(a)可以看出，samp53的数据主要为间断的地形，建筑物较少，但是其数据结构比较复杂，不连续的陡坡分层断开，并且陡坡上分布了建筑物。从图3(c)和图3(b)可以看出建筑物可以很好的滤除，断层结构上存在较少的误判，有过度滤波的现象。

实施例4：

该实施例中采用本发明的方法对ISPRS提供的测试样本sample54进行滤波的仿真。

仿真4仿真条件是在MATLAB2010a软件下进行。

参照图4，对测试样本sample54，包含8608个点的点云数据进行仿真实验。图4(a)是DSM，反应实验区域的原始地形；图4(b)显示样本数据中的地面点点云分布，所有空白为应当滤去的地物点；图4(c)是本算法滤波后地面点点云的分布结果；图4(d)是最终生成的DEM。

从图4(a)可以看出，sample54的数据表现的是斜面上的低分辨率的建筑物，建筑物为较密集的房屋，房屋周围的植被也很多。从图4(c)和图4(b)可以看出本算法可以很好的将地物点滤去，保留地面点，但插值形成的DEM右边边缘有地物保留。

使用ISPRS小组提出的误差评判标准，对四个实施例的仿真结果计算三类误差。结果如表1所示。

ISPRS测试报告中提供了八种滤波算法结果的误差，其中滤波方法分别包含Elmqvist滤波方法、Sohn滤波方法、Axelsson滤波方法、Pfeifer滤波方法、Brovelli滤波方法、Roggero滤波方法、Wack滤波方法、Sithole滤波方法。将这8种滤波方法的误差与本方法的误差进行比较，如表1所示。

从表1的统计结果看出，与八种经典滤波算法相比，本发明算法满足点云数据滤波严格控制第II类误差的要求，并且在控制第II类误差较小的基础上，同时使得第I类误差和总误差较小，因此本发明的滤波方法在滤除各种地物的同时能较好地保留地面点。

Claims (4)

1.一种基于区域预测的LiDAR点云数据形态学滤波方法，其特征在于，方法包括以下步骤：

步骤一，构建滤波初始表面模型；

步骤二，去除初始表面模型的粗差点，得到滤波表面模型；

步骤三，对滤波表面模型进行分块：

按适当尺寸对滤波表面模型进行分块，其中分块n记为R_n，n＝1,2,3,…,N，N为分块的总数；该分块R_n中包含M个点云数据，其中点云数据m记为r_mn，m＝1,2,3,…,M；

步骤四，滤波处理：

对各分块分别进行滤波处理，其中对分块R_n，n＝1,2,3,…,N的处理过程如下：

从最小结构元素窗口半径开始执行以下操作，之后结构元素窗口的半径每增加一米，执行一次以下操作直至最大结构元素窗口半径：

分别对分块R_n中所有点云数据进行开运算，得到各点云数据开运算后的高程值，

分别计算各点云数据的高程值与其自身的开运算后的高程值的差值，差值大于e_n的点云数据标记为非地面点，且将各点云数据的高程值更新为开运算后的高程值，其中：e_n＝s_n×w×c，w为当前结构元素窗口的半径，c为格网边长，且1米≤c≤5米，s_n＝a×σ_n+b，σ_n为分块R_n中所有点云数据高程值的标准差，0.005≤a≤0.02，0.005≤b≤0.02；

滤波处理后未标记的点云数据为地面点；

对于分块时不能被整分的滤波表面模型，被分块后剩余的格网，采用如下方法进行滤波处理：

从最小结构元素窗口半径开始执行以下操作，之后结构元素窗口的半径每增加一米，执行一次以下操作直至最大结构元素窗口半径：

分别对剩余格网中所有点云数据进行开运算，得到各点云数据开运算后的高程值，

分别计算各点云数据的高程值与其自身的开运算后的高程值的差值，差值大于e的点云数据标记为非地面点，且将各点云数据的高程值更新为开运算后的高程值，其中：e＝0.3×w×c，w为当前结构元素窗口的半径，c为格网边长，且1米≤c≤5米，滤波处理后未标记的点云数据为地面点。

2.如权利要求1所述的基于区域预测的LiDAR点云数据形态学滤波方法，其特征在于，方法还包括：

将步骤四确定的地面点插值得到数字高程模型，

接着分别确定所有点云数据的初始高程值与数字高程模型之间的垂直距离，其中垂直距离小于0.5米的点云数据为最终地面点；

所述所有点云数据的初始高程值为步骤四所述各点云数据开运算后的高程值。

3.如权利要求1所述的基于区域预测的LiDAR点云数据形态学滤波方法，其特征在于，步骤三中按尺寸大小为(20m×20m)～(50m×50m)对滤波表面模型进行分块。

4.如权利要求1所述的基于区域预测的LiDAR点云数据形态学滤波方法，其特征在于，步骤四中最小结构元素窗口半径为1米，最大结构元素窗口半径为5～30米。