CN102663237B - 基于网格分块与移动最小二乘的点云数据全自动滤波方法 - Google Patents

Abstract

一种基于网格分块与移动最小二乘的点云数据全自动滤波方法：对激光点云数据进行分块处理，得到多个网格；采用动态四叉树方法对各网格建立相应的索引；判断各网格内激光点云数据密度是否大于密度阈值，对激光点云数据密度大于密度阈值的网格进行抽稀处理；对所有网格内激光点云数据用移动最小二乘法拟合数字高程模型，并以数字高程模型作为参考面；求取所有网格内每个激光点分别到参考面的距离，将距离大于距离阈值的激光点进行删除，剩余的激光点云数据保留；对保留的激光点云数据，返回执行直至所有网格内激光点到参考面的距离小于或等于当前的距离阈值时停止。

Description

基于网格分块与移动最小二乘的点云数据全自动滤波方法

技术领域

本发明属于海量点云数据处理领域，特别是涉及了基于网格分块与移动最小二乘的全自动滤波方法。

背景技术

由于激光能在短时间能获得地物三维坐标信息，并且数据量极大，如何快速从海量激光点云数据中提取有用的信息是目前研究的热点和难点，激光点云数据的滤波处理就显得尤为重要。国内外许多文献都对点云滤波进行了讨论和研究，并且提出了许多滤波算法，包括基于数学形态学的滤波算法，基于坡度的滤波算法，基于TIN的渐进加密算法和基于数据分割的滤波算法等等，也取得了一些研究成果，其中还存在一些问题尚未解决，比如点云滤波算法的效率低以及精度不能满足要求。Vosselman等使用Delaunay三角网组织数据，并利用灰度级形态学中的腐蚀运算根据坡度去除非地面点。该算法首先使用训练样区确定滤波函数，然后根据每个激光点与周围点的距离和滤波函数判断该点是否为地面点。该算法的缺点是：要找到合适的训练样区以获取滤波函数，用来处理其他各种情况是比较困难的。Roggero和Sithole分别对Vosselman的方法进行了改进，通过构造一些新的复杂的滤波函数来改进算法的适应性。与Vosselman的方法不同，Zhang提出一种渐进式形态学滤波的方法。该方法首先将原始激光点云数据插值为规则格网，再使用设定的初始结构元素进行形态学开运算，并计算每点经过开运算后的高度值与原先高度值的差，如果该高度差大于一定的阈值，就为非地面点。逐渐增大结构元素的尺寸和相应的阈值，反复进行上面的运算去除非地面点，直到结构元素的尺寸大于最大的建筑物。这样反复进行开运算不仅需要很大的计算量，而且容易误删除地面点。虽然可以对开运算前后的高度差设定阈值来判断是否为地面点，但此阈值的设定也是一个难点。

不管是机载Lidar数据还是地面Lidar数据大部分是基于激光点云中高程突变信息进行滤波，假定点云中高程低的点为地面点，高程较高的点为非地面点，由于系统误差存在，这种情况未必是对的。一些算法还将离散的观测值内插成格网，内插后数据与原始数据之间存在误差，处理得到的结果与实际情况也存在误差，造成精度损失。目前滤波算法需要大量的人工干预，属于半自动或者人工操作，耗时耗力。半自动算法一般是基于对高程值的统计分析，这些都是目前滤波算法中的一些待解决的问题。另外一些滤波算法适用范围有限，比如一些针对机载Lidar数据，一些只能用于地面Lidar数据滤波，还有一些只适合地形起伏变化小的地方。从上述问题可以看出提出一种简单、快速、适用范围广、效率高的点云滤波算法是非常必要的。

由于激光点云数据量大，并且点云数据的不规则，散乱，复杂等性质决定了点云数据处理工作的复杂困难。

发明内容

针对激光点云数据的特点，本文提出了一种基于网格分块与移动最小二乘的点云数据全自动滤波方法。

本发明提供的技术方案是一种基于网格分块与移动最小二乘的点云数据全自动滤波方法，包括以下步骤：

步骤1，对激光点云数据进行分块处理，得到多个初始网格；

步骤2，判断各初始网格内激光点云数据密度是否大于预设的密度阈值，对激光点云数据密度大于预设的密度阈值的网格进行抽稀处理；所述抽稀处理的方法为，根据激光点云数据的密度大小进一步分割初始网格，保留每个分割后网格内高程较低的激光点云数据；

步骤3，基于所有分割后网格内保留的激光点云数据，用移动最小二乘法拟合数字高程模型，并以数字高程模型作为参考面；

步骤4，求取所有分割后网格内每个激光点分别到步骤3所得参考面的距离，将距离大于距离阈值的激光点进行删除，剩余的激光点云数据保留；

步骤5，对保留的激光点云数据，返回执行步骤1～步骤4，每次执行步骤1时，分块处理的初始网格尺寸小于上一次，每次执行步骤4时，距离阈值小于上一次；直至所有分割后网格内激光点到参考面的距离小于或等于当前的距离阈值时停止。

而且，步骤1对激光点云数据进行分块处理后，采用四叉树方法对所得各初始网格建立相应的索引。

而且，所述四叉树方法为经典四叉树方法或动态四叉树方法。

而且，步骤3采用精简移动最小二乘法拟合数字高程模型，包括以下分步骤，

步骤3.1，对于每个分割后网格，找出步骤2所保留的激光点云数据中的每个最低点周围的K个邻近点，并判断最低点与周围K个邻近点分别的高差，如果最近点与最近邻近点的高差大于最近邻近点与其他K-1个邻近点之间最小距离N倍，则认为最近点不是地面点，反之则认为是地面点；

步骤3.2，通过将所有分割后网格内地面点的坐标值作为已知量，代入以下形函数公式拟合DEM，

$N^k=w(x-x_i)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{p_i\left(x\right)p\left(x\right)}{(p_i\left(x\right),p_i\left(x\right))}$

其中，w(x-x_i)为权重函数，p_i (x)为点集x和权函数w_i所构成集合的正交函数集，点集x由坐标点x_i构成，i＝1，2，...n。

本发明提出了先将点云数据进行网格分块，保证点云数据的原始性，减少单次数据处理量，然后对分块数据建立空间索引，可以快速对点云进行管理与应用，提高点云数据处理的效率。最后针对每块点云数据进行基于移动最小二乘法点云滤波处理，得到每个分块内的点云滤波结果，然后将这些分块滤波数据合并，根据设定阈值对整个数据进行检查，得到最终滤波结果。因此，本发明技术方案能够实现无人工干预的全自动点云数据滤波处理，并且精度、效率都比较高。

附图说明

图1是本发明实施例的点云滤波处理流程图。

图2是本发明实施例的点云分块示意图。

图3是本发明实施例的动态四叉树空间索引示意图。

图4是本发明实施例的原始的激光点云数据示意图。

图5是本发明实施例的地表的激光点云数据示意图。

具体实施方式

本发明技术方案可由本领域技术人员采用计算机软件技术实现自动运行流程。以下结合附图和实施例详细说明本发明技术方案。参见图1，本发明实施例实现过程包括以下步骤：

步骤1，对激光点云数据进行分块处理，得到多个初始网格，如图2所示。分块目的为了降低数据复杂度，减小每次数据处理的运算量。

实施例利用内存映射技术，读取海量的激光点云数据，判断海量点云数据范围，然后进行分块。由于windows32的限制，进程的虚拟内存不能达到4GB只有3.7G左右，由于点云数据量比较大，使内映射技术能够突破这个限制，可以用来打开超过4GB的文件。

基于网格的点云数据分块技术，分块的大小直接影响到数据处理层次及深度，相应地影响算法的效率。分块越小，分割越细，效率就越低，其合并的区域相对增大，数据的压缩比就越高；反之，效率就越高，而压缩比相对降低。一般情况下，分块大小受到仪器的最小采样间距的限制，低于此间距的划分就显得毫无意义。此外，不同距离上的仪器采样间距是不一样的，被测目标一般都分布在一定的纵深距离内。为数据处理的方便，建议数据处理以过最大纵深距离点且与仪器的视场中心轴垂直的平面作为投影面，将所有采样点投影到该平面，然后对所有采样点基于该平面进行等格网分割。最小格网大小的选择应是最小采样间距的整数倍，具体数值的确定取决于被测对象的复杂度、仪器的最小采样间距以及期望的数据压缩比。

为便于数据管理以及提高数据处理效率，实施例采用动态四叉树方法对各网格建立相应的索引。对所有存在点的网格与点云之间建立索引关系后，在后续步骤对点云进行处理和过滤时，只需要考虑存在点的网格即可。具体实施时，也可采用其他索引方法，例如经典四叉树方法。

为了高效的管理和存储分块网格及网格内的激光点云数据，需要对所有存在点的网格与点云之间建立索引关系，在对点云进行处理和过滤时，只需要考虑存在点的网格即可。四叉树结构的生成和维护相对比较简单，且当空间数据对象分布比较均匀时，基于四叉树的空间索引可以获得比较高的空间数据插入和查询效率，因此四叉树是空间数据库中常用的索引之一。由于激光点云的分布不均匀和边界形状极其不规则，为了克服常规四叉树空间索引结构中的问题，本发明实施例采用已有的空间动态四叉树的方法对分块网格建立索引关系。为了便于实施参考起见，进一步提供了详细说明：

其算法要点为，在开始建立四叉树时，不需要事先确定工作区域的范围，只须把要插入空间数据库的第1个空间对象的MBR(minimum bounding rectangle)最小外包矩形，中心作为四叉树的顶点，随着数据处理工作的进行，对作业空间进行分解。

为便于实施参考起见，提供具体算法如下：

(1)根据第1个插入的空间对象，确定四叉树空间的中心点，如图3所示点A；点A的4个叶子节点的某一个方向是开放的，即叶子节点中没有其它子节点或者子节点个数小于四个，是不饱和的；

(2)计算空间对象的MBR；

(3)搜索出所有包含该空间对象的叶子节点，参见图3；

(4)判断这些叶子节点所包含的空间对象数是否超出阈值，如果节点中的对象数超出阈值，则分2种情况进行处理：

1)如果是开放边界的节点，则要首先计算该节点内的包含所有空间对象的MBR，然后把该节点分解成4个新的叶子节点；

2)否则，如果该节点是非开放节点，则直接把该节点分成4个新的叶子节点；

(5)重新计算新生成叶子节点中所包含的空间对象；

(6)对新生成叶子节点中所包含的空间对象重复(1)～(5)，直到处理完所有激光点云数据。

经过以上流程，图3中得到的四叉树结构为：A为根节点，BCDE为A的叶子节点，FGHI为B的叶子节点，JKLM为E的叶子节点，NOPQ为K的叶子节点，除了BEK其它叶子节点均为开放的叶子节点。

步骤2，判断各初始网格内激光点云数据密度是否大于密度阈值，对激光点云数据密度大于密度阈值的网格进行抽稀处理。由于点云数据不规则散乱，即使进行了分块处理仍然不能保证每个块内的点云数据分布均匀，仍然会出现某些块呢的点云数据量过大，可以对各网格内激光点云数据密度比较大的数据进行抽稀处理，如果分块后点云数据密度较小，则这一步可以省略掉。预设的密度阈值为固定值。根据不同的地形情况，阈值的取值也不同，例如山地阈值可以取大些，而平地的阈值可以设置小一点

建议抽稀方法为，根据激光点云数据的密度大小进一步分割初始网格，保留每个分割后所得的网格内高程较低的激光点云数据。设某个分割后网格内激光点云数据的最低高差为X1，只保留高程在[X1，X2]之间的激光点云数据。X2＞X1，X2和X1之间的差值可以预先设定。

步骤3，对所有分割后网格内激光点云数据用移动最小二乘法拟合数字高程模型(DEM)，并以数字高程模型作为参考面。实施例找到网格内地面点后，则用移动最小二乘法拟合DEM，并以此DEM数据作为参考面。对于每个分割后网格，找出步骤2所保留的激光点云数据中的每个最低点周围的K个邻近点，并判断最低点与周围K个邻近点分别的高差，如果最近点与最近邻近点的高差大于最近邻近点与其他K-1个邻近点之间最小距离N倍，则认为最近点不是地面点，反之则认为是地面点。具体实施时，K和N可以根据具体情况设置。

本发明提出采用移动最小二乘法，提供一种较高次数的多项式逼近方式对散乱点云进行曲面拟合，并且要求拟合函数在各个节点处的误差的平方和最小，能够保证比较高的精度。用该算法拟合的曲面比较平滑，与实际曲面相近似。但是移动最小二乘算法比较复杂，运算效率不高，如果点云数据量比较大，则处理比较困难，因此需要对该方法进行精简，在保证精度前提下提高运算效率。实施例设计的精简移动最小二乘法采用带权的正交函数作为基函数，在求解系数矩阵的时候就可以只考虑对角线的元素，不用求逆矩阵，减少了运算量提高运算效率，同时也提高了精度，适合于数据量比较大的点云数据拟合。

为了便于实施参考起见，本发明进一步提供了实施例的步骤4详细说明：

最小二乘拟合函数：

$f\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\left(x\right){\mathrm{\α}}_i\left(x\right)=p^T\left(x\right)\mathrm{\α}\left(x\right)---\left(1\right)$

其中α_i(x)为第i项待求系数(i＝1，2，...n)，是坐标X函数，p_i(x)为第i项基函数(i＝1，2，...n)，它是一个k阶完备的多项式，n是基函数的项数。p^T(x)为n个基函数之和，α(x)为n个待求系数之和，坐标点x_i，i＝1，2，...n，构成点集x。

待求系数α_i(x)通过有权重的最小二乘拟合得到：

$J=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}w(x-x_i){[f\left(x\right)-f\left(x_i\right)]}^2$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}w(x-x_i){[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\left(x_i\right){\mathrm{\α}}_i\left(x\right)-f\left(x_i\right)]}^2$ (2)

其中：w(x-x_i)为第i项权重函数。J表示待求系数α_i(x)的拟合结果、f(x_i)是拟合函数在坐标点x_i处的函数值。

式(2)可写为：

J＝(Pα-f)^T W(x)(Pα-f)(3)

式中，P为基函数矩阵，f为函数值矩阵，W(x)为权重函数矩阵，α为待求系数α_i(x)构成的行矩阵：

f^T＝[f(x₁)f(x₂)...f(x_n)](4)

$P=\left[\begin{array}{cccc}{p}_1\left(x_1\right)& p_2\left(x_1\right)& L& p_n\left(x_1\right)\\ p_1\left(x_2\right)& p_2\left(x_2\right)& L& p_n\left(x_2\right)\\ M& M& O& M\\ p_1\left(x_n\right)& p_2\left(x_n\right)& L& p_n\left(x_n\right)\end{array}\right]---\left(5\right)$

$W\left(x\right)=\left[\begin{array}{cccc}w(x-x_1)& 0& L& 0\\ 0& w(x-x_2)& L& 0\\ M& M& O& M\\ 0& 0& L& w(x-x_n)\end{array}\right]---\left(6\right)$

待求系数α(x)通过求取目前函数J的极值获取，见式(7)：

$\frac{\∂J}{\∂\mathrm{\α}}=A\left(x\right)\mathrm{\α}\left(x\right)-B\left(x\right)f=0---\left(7\right)$

$A\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}w(x-x_i)p\left(x_i\right)p^T\left(x_i\right)$

B(x)＝[w(x-x₁)p(x₁)，w(x-x₂)p(x₂)，...，w(x-x_n)p(x_n)]

式中w(x-x_i)代表权重函数，p_i(x)为第i项基函数，A(x)为n项权函数与对应基函数乘积之和，B(x)为n项权函数与对应基函数乘积的行矩阵。

而这时如果对于点集x和权函数w_i构成的集合{w_i}(i＝1，2，...n)若存在一组函数p_i(x)(i＝1，2，...n)满足

$(p_k,p_i)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{w}_i{p}_k\left(x_i\right)p_i\left(x_i\right)=\left\{\begin{array}{c}0,(k=i)\\ A_k,(k\≠i)\end{array}\right.(k,i=\mathrm{1,2}...n)---\left(8\right)$

A_k为所得到的函数值。

则称p_i(x)(i＝1，2，...n)为点集x和权{w_i}的正交函数集，那么移动最小二乘可变为：

$\left[\begin{array}{cccc}{p}_1\left(x_1\right)& 0& L& 0\\ 0& p_2\left(x_2\right)& L& 0\\ M& M& O& M\\ 0& 0& L& p_n\left(x_n\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\left(x\right)\\ {\mathrm{\α}}_2\left(x\right)\\ M\\ {\mathrm{\α}}_n\left(x\right)\end{array}\right]=\left[\begin{array}{c}(p_1\left(x\right),f\left(x_1\right))\\ (p_2\left(x\right),f\left(x_2\right))\\ M\\ (p_n\left(x\right),f\left(x_n\right))\end{array}\right]---\left(9\right)$

那么很容易就解出α_i(x)， ${\mathrm{\α}}_i\left(x\right)=\frac{(p_i\left(x\right),f\left(x_i\right))}{(p_i\left(x\right),p_i\left(x\right))}.$

$f\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\left(x\right)\·\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}w(x-x_I)p_i\left(x\right)}{(p_i\left(x\right),p_i\left(x\right))}=N^{k}f---\left(10\right)$

其中形函数N^k为：

$N^k=w(x-x_i)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{p_i\left(x\right)p\left(x\right)}{(p_i\left(x\right),p_i\left(x\right))}---\left(11\right)$

从中可以看出，很容易就求得待定系数α_i(x)，避免了求矩阵的逆，提高了运算效率，并且避免了求解病态方程，同时又提高了精度。

将所有网格内的点云数据的坐标值作为已知量，代入形函数公式(11)，拟合得到DEM。

步骤4，求取所有分割后网格内每个激光点分别到步骤3所得参考面的距离，将距离大于距离阈值的激光点进行删除，剩余的激光点云数据保留。同时将删除的数据另外保存，作为数据备份。

步骤5，对保留的激光点云数据，返回执行步骤1～步骤4，每次执行步骤1时，分块的初始网格尺寸小于上一次，每次执行步骤4时，距离阈值小于上一次；直至所有分割后网格内激光点到参考面的距离小于或等于当前的距离阈值时停止，输出本次在步骤4保留的激光点云数据即可。根据经验，后一次初始网格分块大小和距离判别值可取为前一次的1/2～2/3。依次过滤3-5次，即可获取地表的激光点云数据。如图4所示原始激光点云数据，经本发明提供技术方案处理后，获取地表的激光点云数据如图5所示。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (2)

1.一种基于网格分块与移动最小二乘的点云数据全自动滤波方法，其特征在于，包括以下步骤： 

步骤1，对激光点云数据进行分块处理，得到多个初始网格； 

步骤2，判断各初始网格内激光点云数据密度是否大于预设的密度阈值，对激光点云数据密度大于预设的密度阈值的网格进行抽稀处理；所述抽稀处理的方法为，根据激光点云数据的密度大小进一步分割初始网格，保留每个分割后网格内高程较低的激光点云数据，实现方式为，设某个分割后网格内激光点云数据的最低高差为X1，只保留高程在[X1,X2]之间的激光点云数据，X2＞X1，X2和X1之间的差值预先设定； 

步骤3，基于所有分割后网格内保留的激光点云数据，用移动最小二乘法拟合数字高程模型，并以数字高程模型作为参考面；进一步包括采用精简移动最小二乘法拟合数字高程模型，包括以下分步骤， 

步骤3.1，对于每个分割后网格，找出步骤2所保留的激光点云数据中的每个最低点周围的K个邻近点，并判断最低点与周围K个邻近点分别的高差，如果最近点与最近邻近点的高差大于最近邻近点与其他K-1个邻近点之间最小距离N倍，则认为最近点不是地面点，反之则认为是地面点； 

步骤3.2，通过将所有分割后网格内地面点的坐标值作为已知量，代入以下形函数公式拟合DEM， 

其中，w(x-x_i)为权重函数，p_i(x)为点集x和权函数w_i所构成集合的正交函数集，点集x由坐标点x_i构成，i＝1,2,…n； 

步骤4，求取所有分割后网格内每个激光点分别到步骤3所得参考面的距离，将距离大于距离阈值的激光点进行删除，剩余的激光点云数据保留； 

步骤5，对保留的激光点云数据，返回执行步骤1～步骤4，每次执行步骤1时，分块处理的初始网格尺寸小于上一次，每次执行步骤4时，距离阈值小于上一次；直至所有分割后网格内激光点到参考面的距离小于或等于当前的距离阈值时停止。 

2.如权利要求1所述基于网格分块与移动最小二乘的点云数据全自动滤波方法，其特征在于： 步骤1对激光点云数据进行分块处理后，采用四叉树方法对所得各初始网格建立相应的索引。

3.如权利要求2所述基于网格分块与移动最小二乘的点云数据全自动滤波方法，其特征在于：所述四叉树方法为经典四叉树方法或动态四叉树方法。