CN105118090A - 一种自适应复杂地形结构的点云滤波方法 - Google Patents

Abstract

本发明公开了一种自适应复杂地形结构的点云滤波方法，解决了现有技术中于地形起伏变化、城市区域地物结构复杂，导致的滤波算法对参数变化敏感，误分类严重的问题。该自适应复杂地形结构的点云滤波方法包括以下步骤：（1）构建点云金字塔；（2）采用最粗层级的金字塔数据初始化地面点G，并将下一级数据加入待分类点U，更新待分类点；（3）迭代处理点云金字塔的所有层级，处理完所有待分类点，并且保存新增加地面点；（32）根据DEM、弯曲能量以及金字塔尺度信息，采用自适应复杂地形结构的参数优选方法获取滤波阈值并进行滤波，区分地面点与非地面点；（33）重复步骤（31）~（32），直至无新增地面点；（4）存储分类后的地面点与非地面点。

Description

一种自适应复杂地形结构的点云滤波方法

技术领域

本发明属于地理空间信息系统技术领域，主要涉及一种自适应复杂地形结构的点云滤波方法。

背景技术

随着多平台激光雷达和多角度影像密集匹配技术的快速发展，大量快速可得的点云数据的高分辨率高精度的地形表面三维重建提供了有效可靠的数据支撑，点云数据处理已成为国际研究和工业界的热点问题。如何从点云数据中分离地面与非地面点，即点云数据滤波，是众多后续应用的首要步骤，因此一直是热点前沿问题。

在实际生产实践中，由于地表结构的多样性和复杂性，同一套滤波参数难以适应错综复杂的地形特征，当前国内外点云处理的主流商业软件，如TerraScan、LasTools等，其自动处理精度难以达到实际生产需求，通常采用半自动滤波加人工交互后处理的方法进行滤波处理,需要耗费大量的人力、物力；并且由于人为主观因素影响，有时还需要诸如真正射影像等信息辅助处理，并可能出现不同程度的误差，极大的影响了生产效率。因此，自动、高效、稳健的点云滤波成为学术界和工业界面临的重大挑战。

通常点云滤波需要假设在局部邻域内的最低点为地面点，通过建立临时的地面结构，通过一定阈值参数对地面点进行加密，并区分地面-非地面点。该方法主要存在缺陷包括：一方面，地形起伏变化、城市结构形态复杂；另一方面多平台激光雷达和多角度影像密集匹配产生的点云数据，无论在点云密度还是数据特性上均有较大差异；而现有的滤波算法通常对同一区域采用一套参数，难以适应错综复杂的地面-非地面物体结构，并且通常对参数变化较为敏感，导致误分类问题严重。鉴于此，构建对地面起伏变化具有自适应能力的滤波参数选择方法，成为可能的有效途径之一。另一方面由于点云中的噪声、地面点加密过程中的误差，在由粗到精的金字塔滤波过程中，难免会受误分类点的影响，如不有效处理这些误分类点将在渐进式的滤波过程中被不断放大，大幅增加人工交互式后处理工作量。因此，如何进行点云内插使其具有良好的局部抗噪性是目前亟待解决的问题。

目前，现有的主流点云滤波算法可以分为两类：基于内插的滤波方法和基于形态学运算的滤波方法：

(1)基于内插的滤波方法

基于内插的滤波方法的基本思想是，首先在一较大的窗口内，提取最低点作为地面种子点，构建不规则三角网(TIN)或数字地表模型(DEM)用来表达地面，然后根据适当阈值渐进的对地面点进行加密。该方法可十分简洁的扩展为由粗到精的金字塔滤波策略，点云金字塔中每一层级存储的是对应窗口内的最低值。然而，在不断内插地表模型与迭代滤波过程中，地面点不可避免包含粗差点，若不有效处理将被不断放大，严重影响最终滤波结果。目前，Terrascan采用一种渐进三角网加密的滤波策略，LasTools则采用线性预测的滤波策略，两者均属于此类方法。

(2)基于形态学运算的滤波方法

基于形态学的滤波方法的基本思想来源于数学上的形态学计算——用形态学开运算和闭运算从灰度影像上识别目标。建立在最低高程假设之上，非地面点的高程通常高于周围的地面点，如果将点云数据按高程值转化为灰度影像，由于灰度差异，非地面目标将被识别出来。对点云数据进行开运算，可以将局部窗口内高程值最小的点识别出来。通过判断点云与窗口内的最小高程点之间的坡度差异与角度差异作为滤波准则，移动滤波窗口，实现对全部点云数据的滤波。基于形态学运算的滤波方法的关键点在于选取一个合适的运算窗口：若选择窗口过小，则只有较小的非地面目标(比如树木、汽车)能被滤除，而对于较大地物则会在膨胀中再次保留；而若选择窗口过大，则在坡度较大区域，部分地面点也会被误分为非地面点。该方法虽然效率较高，然而由于对参数选择十分敏感，难以处理复杂地形点云，实用性不大。

发明内容

本发明的目的在于克服上述缺陷，提供一种将能显式刻画地形变化特性的弯曲能量引入到点云数据滤波过程中，使滤波参数随地形变化在局部范围内自适应优化调整；并将地表一致性约束引入内插算法之中，采用弯曲能量的规则化约束，克服噪声对滤波精度的影响，并建立一套科学合理的点云自动滤波实施方案。

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

一种自适应复杂地形结构的点云滤波方法，包括以下步骤：

(1)构建点云金字塔；

(2)采用格网窗口最大的金字塔层级数据初始化地面点G，并将下一级数据加入待分类点U，更新待分类点；

(3)迭代处理点云金字塔的所有层级，处理完所有待分类点，并且保存新增加地面点

(31)利用已有地面点，采用具有局部抗噪性的DEM内插算法获取地面DEM，并同时获取对应区域的弯曲能量；

(32)根据DEM、弯曲能量以及金字塔尺度信息，采用自适应复杂地形结构的参数优选方法获取滤波阈值并进行滤波，区分地面点与非地面点；

(33)重复步骤(31)～(32)，直至无新增地面点；

(4)存储分类后的地面点与非地面点，其中，地面点标记为2，非地面点标记为1。

在构建点云金字塔之前，剔除粗差数据，其具体方式如下：

首先，对点云数据构建KD-树的索引结构；然后，对每个点进行半径搜索，计算在5米半径值内的点的数目，统计其均值μ和标准差σ，将邻域内的点数目小于μ-3σ的点视为孤立点，进行移除。

所述步骤(1)的具体方式如下：

(11)根据用户输入的最大窗口大小w_max与最小窗口大小w_min，均匀获取一系列金字塔窗口大小w；

(12)从最粗层级开始，对原始数据进行规则格网划分，每个格网中记录该格网的最低点，并从原始点云中移除；

(13)迭代处理所有层级，建立点云金字塔。

所述步骤(31)中具有局部抗噪性的DEM内插算法包括以下步骤：

(311)在获取所有地面点之后，迭代处理所有地面点，确定所有数据点的三维轴向外包盒大小，用于后续确定DEM范围、行列数等信息；

(312)初始化DEM和弯曲能量栅格结构，并根据金字塔层级信息，计算正则化约束权重λ；DEM格网原点、行列数、DEM窗口大小w等信息；该正则化约束权重与地面点中包含的粗差点比例有关：若粗差比例越大，则权重越大；反之权重越小，直至λ＝0

(313)对地面点的二维平面坐标，构建KD-树的空间索引，用于迭代处理DEM格网的近邻索引，内插对应的格网信息；

(314)对每一DEM格网，求取该格网处的二维平面坐标(X,Y)，在KD-树索引中按自适应半径增长方法搜索至少9个近邻点，并且获取其对应的三维坐标(X,Y,Z)，作为样条函数的控制点；

(315)根据控制点三维坐标，更新样条函数矩阵，用于计算样条函数参数，样条函数的解具有线性方程的闭合解形式，采用QR分解求解样条函数的参数信息；

(316)利用样条函数参数，计算DEM格网对应二维平面坐标(X,Y)处的高程值z并保存入DEM的栅格结构Z中；

(317)重复步骤(314)至(316)直至迭代处理所有DEM格网，获取最终的DEM与弯曲能量栅格图。

所述步骤(313)的具体方式如下：根据获取所有地面点的二维平面坐标(X,Y)，采用开源算法建立二维地面点数据的KD-树索引；其中KD-树中的索引序号与地面点序号相同。

所述步骤(314)的具体方式如下：

(3141)计算DEM格网点处的二维平面坐标(X,Y)；

(3142)获取当前窗口大小w，并且以此初始化一个数组r＝[w,2w,3w,4w,5w]，对应半径搜索窗口大小；

(3143)从最小的搜索半径r＝w开始，利用KD-树索引，搜索半径r中包含的点的索引号，并且获取对应的三维坐标(X,Y,Z)；

(3144)将所有三维点分为4个象限，每个象限中的三维点，按与DEM格网点距离大小由小到大排序；

(3145)若每象限点数均大于等于4，则每个象限选取离格网点最近的三维点，作为样条函数控制点；

若其中一个象限点数小于4，则增大搜索半径r＝r_i+1，并且重复步骤(3143)～(3145)，直至r＝5w；

若搜索半径r＝5w，近邻点仍小于9，则将该格网点设置为无效之-9999，并跳过后续步骤(315)和(316)；

(3146)返回搜索到的近邻点的索引号集合，并且获取对应的三维坐标，用于计算样条函数的控制点。

所述步骤(316)的具体步骤如下：

(3161)计算DEM高程值：将DEM格网点坐标(X,Y)和样条函数参数和带入函数Z(X,Y)之中，即求得DEM对应格网处的高程值，并存入DEM对应的行列号之中；

(3162)计算DEM格网处的弯曲能量：根据样条函数参数以及矩阵K，计算样条函数的弯曲能量并存入弯曲能量对应的行列号之中；

(3163)直至计算完所有格网点的信息，并将计算完成的DEM与弯曲能量传出至下一阶段处理。

所述步骤(32)中自适应复杂地形结构的参数优选方法包括以下步骤：

(321)依据金字塔层级信息s，计算当前层级对应的尺度补偿阈值T_s；

(322)对弯曲能量栅格图，统计其最大值与最小值，并采用分段线性内插的方式，将弯曲能量转化为对应的弯曲能量补偿阈值T_b；

(323)对每一待分类点，计算其落入对应层级的格网行列坐标，采用双线性内插方式，获得对应3×3邻域DEM与弯曲能量栅格图上的高程均值z_m和弯曲能量补偿值T_b；

(324)计算自适应滤波阈值，该阈值综合考虑金字塔层级的尺度信息s，弯曲能量信息b，以及初始滤波阈值t，故最终的滤波阈值应是上述三者的函数T＝f(s,b,t)，即T＝T_b+T_s+t；

(325)依据待分类点高程，DEM高程与滤波阈值之间关系判断该点为地面或非地面点；

(326)重复步骤(323)至(325)直到处理完所有待分类点，并且保存新增加地面点。

所述步骤(321)中阈值计算方式如下：

最粗的金字塔层级，补偿阈值为T_s＝T_smax；

最精细的金字塔层级，补偿阈值为0；

中间金字塔层级，补偿阈值根据窗口大小w采用线性内插方法获取，即 $T_s=T_{s\_\max }\frac{w-w_{\min }}{w_{\max }-w_{\min }}.$

所述步骤(322)获取弯曲能量补偿阈值T_b的具体步骤为：

(3221)对所有DEM格网点的弯曲能量从小至大排序；

(3222)将所有弯曲能量均分为10段，并记录分割处的弯曲能量大小[b₀,b₁,·····,b₁₀]；

(3223)设置一最大补偿阈值T_bmax，对所有弯曲能量值b，进行分段内插，若其落入第i段之内，则有

与现有技术相比，本发明具有以下有益效果：

(1)从同时包含地面-非地面点的点云数据中，如何内插出真实地表将直接影响到整体的滤波精度。由于三维点云提供的信息较为有限，只有顾及地面一致性约束，才能一定程度上克服噪声点对DEM内插带来的影响。然而现有通过聚类、平面拟合的约束方法并不具有通用性。因此，本发明依据地面二阶连续平滑的物理意义，直接在内插过程中引入表示平滑程度的弯曲能量正则化约束，在内插过程中顾及地表一致性约束，降低噪声点对DEM内插的影响，整体提高滤波精度。

(2)在得到表达地表形状的DEM之后，需要确定区分地面-非地面点的阈值参数，现有滤波算法通常对同一区域采用相同的参数，而不考虑地形结构信息，因此对参数变化十分敏感。已有研究已经指出，在滤波过程中需要考虑地形结构信息以提高滤波性能的稳定性，现有的方法通常仅是对3×3窗口的DEM值做简单的加权平均处理。而采用本发明，弯曲能量能有效表达局部地形结构信息的，采用弯曲能量使滤波参数自适应的随地形改变，可以显著提高滤波算法的可靠性，降低对滤波参数的敏感程度。

附图说明

图1为本发明各模块之间关系。

图2为本发明算法总体流程图。

图3为顾及弯曲能量正则化约束的局部抗造性内插方法。

图4为DEM与弯曲能量示意图。

图5为顾及弯曲能量的滤波参数自适应优选方法。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。

实施例

本实施例提供了一种自适应复杂地形结构的点云滤波方法，深入挖掘弯曲能量在DEM内插与表达地形起伏变化的特殊作用，在地表二阶连续平滑的假设之上，在DEM内插过程中利用弯曲能量作为正则化约束，克服噪声影像。并且利用弯曲能量显示刻画地形起伏变化，用于滤波参数的自适应优选。在此基础之上，设计一套自动化的滤波流程，采用由粗到精的金字塔滤波策略，提高滤波效率和可靠性。

本实施例中，涉及的各模块关系如图1所示，在由粗到精的滤波过程中，每次迭代包含以下两个主要步骤：1)顾及弯曲能量正则化约束的局部抗噪性DEM内插算法和2)顾及弯曲能量的滤波参数优选方法，本发明的总体流程如图2所示，具体包括以下步骤：

(1)构建点云金字塔；基于地面最低高程假设，在逐级递进的滤波过程中，金字塔中存储的点云数据高程值在局部范围内将逐渐增加；

(2)采用最粗层级的金字塔数据初始化地面点G，并将下一级数据加入待分类点U，更新待分类点；

(3)迭代处理点云金字塔的所有层级，处理完所有待分类点，并且保存新增加地面点，即对每一层级作如下处理：

(31)利用已有地面点，采用具有局部抗噪性的DEM内插算法获取地面DEM，并同时获取对应区域的弯曲能量；

(32)根据DEM、弯曲能量以及金字塔尺度信息，采用自适应复杂地形结构的参数优选方法获取滤波阈值并进行滤波，区分地面点与非地面点；

(33)重复步骤(31)～(32)，直至无新增地面点；

(4)存储分类后的地面点与非地面点，保存为.las格式数据，其中，地面点标记为2，非地面点标记为1。

在构建点云金字塔之前，还需剔除粗差数据，点云粗差通常来自于激光雷达的多回波效应以及密集匹配噪声，考虑到这些粗差点通常游离于点云数据之外，且通常数量稀少，表现为孤立数据，因此，可以考虑采用其邻域内的点数的统计信息进行移除。首先，对点云数据构建KD-树的索引结构，随后对每个点进行半径搜索，计算在某一半径值内的点的数目，统计其均值μ和标准差σ，将邻域内的点数目小于μ-3σ的点视为孤立点，进行移除。

构建点云数据金字塔：根据用户输入的最大窗口大小w_max与最小窗口大小w_min，均匀获取一系列金字塔窗口大小w。从最粗层级开始，对原始数据进行规则格网划分，每个格网中记录该格网的最低点，并从原始点云中移除。迭代处理所有层级，建立完整点云金字塔。

如图3所示，步骤(31)中具有局部抗噪性的DEM内插算法包括以下步骤：

(311)获取所有地面点之后，迭代处理所有地面点，确定所有数据点的三维轴向外包围盒(AxisAlignedBoundingBox，AABB)大小，即三维点云坐标的最大最小值，(Xmin,Xmax,Ymin,Ymax,Zmin,Zmax)。计算方法是，迭代处理所有地面点，对每个点，分别比较该点三维坐标(X,Y,Z)与当前AABB范围，若该点落入AABB之外，则对AABB进行更新；

(312)根据金字塔层级信息，DEM格网原点以及行列数，DEM窗口大小w，正则化约束权重λ等信息，并且初始化DEM和弯曲能量栅格结构。正则化约束权重λ需要与地面点包含粗差比例成正比，随着由粗到精的金字塔滤波的逐级递进，地面点包含噪声比例也随之增大。在最粗的金字塔层级，假设无粗差点，窗口大小为w_max，此时正则化约束权重λ＝0；而在最精细的金字塔层级，此时粗差点最多，窗口大小为w_max，此时正则化约束权重达到最大值λ＝λ_max；中间层级，则采用线性内插方式获取，即初始化DEM和弯曲能量栅格图时，需要利用AABB平面最小值(Xmin,Ymin)，窗口大小w，以及行列数目；

(313)对地面点的二维平面坐标，构建KD-树的空间索引，用于迭代处理DEM格网的近邻索引，内插对应的格网信息；具体方法为，根据获取所有地面点的二维平面坐标(X,Y)，采用开源算法FLANN建立二维地面点数据的KD-树索引。其中KD-树中的索引序号与原始三维点序号相同；

(314)对每一DEM格网，求取该格网处的二维平面坐标(X,Y)，在KD-树索引中按一定规则搜索一定数目的近邻点，并且获取其对应的三维坐标(X,Y,Z)，作为样条函数的控制点；在搜索过程中需要满足一下条件：1)样条函数控制点应分布均匀；2)控制点数目应满足最小要求；3)控制点应选和该格网最相关的三维点，即和格网点最靠近的点。为满足上述要求，具体实时方法如下：

(3141)计算DEM格网点处的二维平面坐标(X,Y)；

(3142)获取当前窗口大小w，并且以此初始化一个数组r＝[w,2w,3w,4w,5w]，对应半径搜索窗口大小；

(3143)从最小的搜索半径r＝w开始，利用KD-树索引，搜索半径r中包含的点的索引号，并且获取对应的三维坐标(X,Y,Z)；

(3144)将所有三维点分为4个象限，每个象限中的三维点，按与DEM格网点距离大小由小到大排序；

(3145)若每象限点数均大于等于4，则每个象限选取离格网点最近的三维点，作为样条函数控制点；

若其中一个象限点数小于4，则增大搜索半径r＝r_i+1，并且重复步骤(3143)～(3145)，直至r＝5w；

若搜索半径r＝5w，近邻点仍小于9，则将该格网点设置为无效之-9999，并跳过后续步骤(315)和(316)；

(3146)返回搜索到的近邻点的索引号集合，并且获取对应的三维坐标，用于计算样条函数的控制点。

(315)传统的DEM内插，通常是通过拟合地面控制点P＝{pi＝(Xi,Yi,Zi)|i＝1,2,3,…,n}，得到一张参数化的曲面Z＝f(X,Y)，在拟合过程中，由于地面控制点中不包含粗差，因此，仅需要保持对地面控制点有较好的拟合度，就可取得较好的效果。即使得数据拟合度ε_data(dataterm)具有较小的值，或完全经过所有控制点，如下式所示

${\mathrm{\ϵ}}_{\mathrm{data}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Z_i-f\left(X_{i,}{Y}_i\right))}^2$

然而，在点云滤波过程中，地面点不可避免的会引入非地面粗差点，因此此时即使ε_data为0，也并不意味着对地面有较好的拟合程度。采用地表一致性约束，可以一定程度上，降低噪声敏感性，提高内插算法的鲁棒性。弯曲能量ε_smooth，如下式所示，

${\mathrm{\ϵ}}_{\mathrm{smooth}}=\∫\∫(f_{\mathrm{xx}}^2(X,Y)+2f_{\mathrm{xy}}^2(X,Y)+2f_{\mathrm{yy}}^2(X,Y))\mathrm{dXdY}$

可以表达一张参数化曲面的二阶连续程度，因此若将ε_smooth作为一个正则化约束引入地面内插过程中，则可隐式的蕴含着地表一致性约束，抵抗噪声点的干扰。即通过优化具有ε＝ε_data+λε_smooth形式的能量函数，进行DEM内插。

给定控制点集P，该优化模型为的解为薄板样条函数，即下式所示，

$Z(X,Y)=a_0+a_{1}X+a_{2}Y+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_{i}R\left(\right|(X_i,Y_i)-(X_j,Y_j)\left|\right)$

为获取该样条函数的参数和，可通过以下形式的线性方程组求解：

$\left[\begin{array}{cc}K& C\\ C& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\→}{p}\\ \stackrel{\→}{a}\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{v}\\ \stackrel{\→}{0}\end{array}\right]$

其中，上述矩阵中的各元素具有以下形式，并且可直接采用QR分解的方式获取最终的参数和

K_ij＝R(|(X_i,Y_i)-(X_j,Y_j)|)+λI_ijα²

R(r)＝r²logr

$\mathrm{\α}=\frac{1}{n^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}|(X_i,Y_i)-(X_j,Y_j)|$

C_i＝[1X_iY_i]

${\stackrel{\→}{v}}_i=Z_i;$

(316)利用样条函数参数，计算DEM格网对应二维平面坐标(X,Y)处的高程值z并保存入DEM的栅格结构Z中；其具体步骤如下：

(3161)计算DEM高程值：将DEM格网点坐标(X,Y)和样条函数参数和带入函数Z(X,Y)之中，即求得DEM对应格网处的高程值，并存入DEM对应的行列号之中；

(3162)计算DEM格网处的弯曲能量：根据样条函数参数以及矩阵K，计算样条函数的弯曲能量并存入弯曲能量对应的行列号之中；

(3163)直至计算完所有格网点的信息，并将计算完成的DEM与弯曲能量传出至下一阶段处理。

(317)重复步骤(314)至(316)直至迭代处理所有DEM格网，获取最终的DEM与弯曲能量栅格图。

如图5所示，步骤(32)中自适应复杂地形结构的参数优选方法包括以下步骤：

(321)依据金字塔层级信息s，计算当前层级对应的尺度补偿阈值Ts；考虑到地形坡度影响，在较粗的金字塔层级，由于窗口w较大，即使平坦区域也可能会存在微小的地形起伏，因此需要根据金字塔层级s和窗口大小w补偿地形坡度起伏影像。具体而言，在最粗的金字塔层级，补偿阈值最大为Ts＝Ts_max；在最精细的金字塔层级，补偿阈值为0；而在中间层级的数据，补偿阈值根据窗口大小w采用线性内插方法获取，即

(322)将弯曲能量栅格图采用分段线性内插转化为对应的弯曲能量补偿阈值T_b。顾及弯曲能量的参数自适应优选方法的基本思想是，在地形起伏较大的地方，如山脊线、地面断裂处等区域，滤波阈值应更大；而在地形起伏较小的区域，如平坦地表等区域，滤波阈值应更小。由于在地形欺负较大的地方，地形变化弯曲能量也相应的会变大，因此滤波阈值应与弯曲能量成正比例变化。因此，其具体获取弯曲能量补偿阈值T_b的具体步骤为；

(3221)对所有DEM格网点的弯曲能量从小至大排序；

(3222)将所有弯曲能量均分为10段，并记录分割处的弯曲能量大小[b₀,b₁,·····,b₁₀]；

(3223)设置一最大补偿阈值T_bmax，对所有弯曲能量值b，进行分段内插，若其落入第i段之内，则有

(323)内插待分类点所处格网的高程均值z_m与弯曲能量补偿值T_b。由于DEM在滤波过程中可能会存在一定噪声，因此不能仅用该格网处的DEM高程值，而考虑采用3×3邻域DEM的高程值z和弯曲能量补偿值T_b。首先，根据待分类点坐标(X,Y,Z)，计算其所对应的格网行列号，r＝(Y-Y_min)/w,c＝(X-X_min)/w。随后，在DEM和补偿阈值中获取对应3×3邻域内的高程均值z_m和补偿阈值T_b；

(324)计算自适应滤波阈值，该阈值综合考虑金字塔层级的尺度信息s，弯曲能量信息b，以及初始滤波阈值t，故最终的滤波阈值应是上述三者的函数T＝f(s,b,t)，即T＝T_b+T_s+t；

(325)依据待分类点高程，DEM高程与滤波阈值之间关系判断该点为地面或非地面点；在3x3邻域内，判断该待分类点和邻域内的9个格网点的相互关系。记N为待分类点高程值Z<z+T的数目，若N≥5，则该待分类点为地面点；否则仍然保留为待分类点，在后续处理中继续判断是否为地面点；

(326)重复步骤(323)至(325)直到处理完所有待分类点，并且保存新增加地面点，将新增地面点传出，供后续步骤处理。

按照上述实施例，便可很好地实现本发明。值得说明的是，基于上述设计原理的前提下，为解决同样的技术问题，即使在本发明所公开的结构基础上做出的一些无实质性的改动或润色，所采用的技术方案的实质仍然与本发明一样，故其也应当在本发明的保护范围内。

Claims (10)

1.一种自适应复杂地形结构的点云滤波方法，其特征在于，包括以下步骤：

(1)构建点云金字塔；

(2)采用格网窗口最大的金字塔层级数据初始化地面点G，并将下一级数据加入待分类点U，更新待分类点；

(3)迭代处理点云金字塔的所有层级，处理完所有待分类点，并且保存新增加地面点

(31)利用已有地面点，采用具有局部抗噪性的DEM内插算法获取地面DEM，并同时获取对应区域的弯曲能量；

(32)根据DEM、弯曲能量以及金字塔尺度信息，采用自适应复杂地形结构的参数优选方法获取滤波阈值并进行滤波，区分地面点与非地面点；

(33)重复步骤(31)～(32)，直至无新增地面点；

(4)存储分类后的地面点与非地面点，其中，地面点标记为2，非地面点标记为1。

2.根据权利要求1所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，在构建点云金字塔之前，剔除粗差数据，其具体方式如下：

首先，对点云数据构建KD-树的索引结构；然后，对每个点进行半径搜索，计算在5米半径内的点的数目，统计其均值μ和标准差σ，将邻域内的点数目小于μ-3σ的点视为孤立点，进行移除。

3.根据权利要求1所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(1)的具体方式如下：

(11)根据用户输入的最大窗口大小w_max与最小窗口大小w_min，均匀获取一系列金字塔窗口大小w；

(12)从最粗层级开始，对原始数据进行规则格网划分，每个格网中记录该格网的最低点，并从原始点云中移除；

(13)迭代处理所有层级，建立点云金字塔。

4.根据权利要求1所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(31)中具有局部抗噪性的DEM内插算法包括以下步骤：

(311)在获取所有地面点之后，迭代处理所有地面点，确定所有数据点的三维轴向外包盒大小，用于后续确定DEM范围、行列数信息；

(312)初始化DEM和弯曲能量栅格结构，并根据金字塔层级信息，计算正则化约束权重λ；

(313)对地面点的二维平面坐标，构建KD-树的空间索引，用于迭代处理DEM格网的近邻索引，内插对应的格网信息；

(314)对每一DEM格网，求取该格网处的二维平面坐标(X,Y)，在KD-树索引中按自适应半径增长方法搜索至少9个近邻点，并且获取其对应的三维坐标(X,Y,Z)，作为样条函数的控制点；

(315)根据控制点三维坐标，更新样条函数矩阵，用于计算样条函数参数，样条函数的解具有线性方程的闭合解形式，采用QR分解求解样条函数的参数信息；

(316)利用样条函数参数，计算DEM格网对应二维平面坐标(X,Y)处的高程值z并保存入DEM的栅格结构Z中；

(317)重复步骤(314)至(316)直至迭代处理所有DEM格网，获取最终的DEM与弯曲能量栅格图。

5.根据权利要求4所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(313)的具体方式如下：根据获取所有地面点的二维平面坐标(X,Y)，采用开源算法建立二维地面点数据的KD-树索引；其中KD-树中的索引序号与地面点序号相同。

6.根据权利要求4所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(314)的具体方式如下：

(3141)计算DEM格网点处的二维平面坐标(X,Y)；

(3142)获取当前窗口大小w，并且以此初始化一个数组r＝[w,2w,3w,4w,5w]，对应半径搜索窗口大小；

(3143)从最小的搜索半径r＝w开始，利用KD-树索引，搜索半径r中包含的点的索引号，并且获取对应的三维坐标(X,Y,Z)；

(3144)将所有三维点分为4个象限，每个象限中的三维点，按与DEM格网点距离大小由小到大排序；

(3145)若每象限点数均大于等于4，则每个象限选取离格网点最近的三维点，作为样条函数控制点；

若其中一个象限点数小于4，则增大搜索半径r＝r_i+1，并且重复步骤(3143)～(3145)，直至r＝5w；

若搜索半径r＝5w，近邻点仍小于9，则将该格网点设置为无效之-9999，并跳过后续步骤(315)和(316)；

(3146)返回搜索到的近邻点的索引号集合，并且获取对应的三维坐标，用于计算样条函数的控制点。

7.根据权利要求4所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(316)的具体步骤如下：

(3161)计算DEM高程值：将DEM格网点坐标(X,Y)和样条函数参数和带入函数Z(X,Y)之中，即求得DEM对应格网处的高程值，并存入DEM对应的行列号之中；

(3162)计算DEM格网处的弯曲能量：根据样条函数参数以及矩阵K，计算样条函数的弯曲能量并存入弯曲能量对应的行列号之中；

(3163)直至计算完所有格网点的信息，并将计算完成的DEM与弯曲能量传出至下一阶段处理。

8.根据权利要求4所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(32)中自适应复杂地形结构的参数优选方法包括以下步骤：

(321)依据金字塔层级信息s，计算当前层级对应的尺度补偿阈值T_s；

(322)对弯曲能量栅格图，统计其最大值与最小值，并采用分段线性内插的方式，将弯曲能量转化为对应的弯曲能量补偿阈值T_b；

(323)对每一待分类点，计算其落入对应层级的格网行列坐标，采用双线性内插方式，获得对应3×3邻域DEM与弯曲能量栅格图上的高程均值z_m和弯曲能量补偿值T_b；

(324)计算自适应滤波阈值：T＝T_b+T_s+t；

(325)依据待分类点高程，DEM高程与滤波阈值之间关系判断该点为地面或非地面点；

(326)重复步骤(323)至(325)直到处理完所有待分类点，并且保存新增加地面点。

9.根据权利要求8所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(321)中阈值计算方式如下：

最粗的金字塔层级，补偿阈值为T_s＝T_smax；

最精细的金字塔层级，补偿阈值为0；

中间金字塔层级，补偿阈值根据窗口大小w采用线性内插方法获取，即

$T_s=T_{s\_\max }\frac{w-w_{\min }}{w_{\max }-w_{\min }}.$

10.根据权利要求8所述的一种自适应复杂地形结构的点云滤波方法，其特征在于，所述步骤(322)获取弯曲能量补偿阈值T_b的具体步骤为：

(3221)对所有DEM格网点的弯曲能量从小至大排序；

(3222)将所有弯曲能量均分为10段，并记录分割处的弯曲能量大小[b₀,b₁,·····,b₁₀]；

(3223)设置一最大补偿阈值T_bmax，对所有弯曲能量值b，进行分段内插，若其落入第i段之内，则有 $T_b=\frac{i+1}{10}{T}_{b\_\max }\frac{b-b_i}{b_{i+1}-b_i}.$