一种自适应复杂地形结构的点云滤波方法
【技术领域】
[0001] 本发明属于地理空间信息系统技术领域，主要涉及一种自适应复杂地形结构的点 云滤波方法。
【背景技术】
[0002] 随着多平台激光雷达和多角度影像密集匹配技术的快速发展，大量快速可得的点 云数据的高分辨率高精度的地形表面三维重建提供了有效可靠的数据支撑，点云数据处理 已成为国际研究和工业界的热点问题。如何从点云数据中分离地面与非地面点，即点云数 据滤波，是众多后续应用的首要步骤，因此一直是热点前沿问题。
[0003] 在实际生产实践中，由于地表结构的多样性和复杂性，同一套滤波参数难以适应 错综复杂的地形特征，当前国内外点云处理的主流商业软件，如TerraScan、LasTools等， 其自动处理精度难以达到实际生产需求，通常采用半自动滤波加人工交互后处理的方法进 行滤波处理，需要耗费大量的人力、物力；并且由于人为主观因素影响，有时还需要诸如真 正射影像等信息辅助处理，并可能出现不同程度的误差，极大的影响了生产效率。因此，自 动、高效、稳健的点云滤波成为学术界和工业界面临的重大挑战。
[0004] 通常点云滤波需要假设在局部邻域内的最低点为地面点，通过建立临时的地面结 构，通过一定阈值参数对地面点进行加密，并区分地面-非地面点。该方法主要存在缺陷包 括：一方面，地形起伏变化、城市结构形态复杂；另一方面多平台激光雷达和多角度影像密 集匹配产生的点云数据，无论在点云密度还是数据特性上均有较大差异；而现有的滤波算 法通常对同一区域采用一套参数，难以适应错综复杂的地面-非地面物体结构，并且通常 对参数变化较为敏感，导致误分类问题严重。鉴于此，构建对地面起伏变化具有自适应能力 的滤波参数选择方法，成为可能的有效途径之一。另一方面由于点云中的噪声、地面点加密 过程中的误差，在由粗到精的金字塔滤波过程中，难免会受误分类点的影响，如不有效处理 这些误分类点将在渐进式的滤波过程中被不断放大，大幅增加人工交互式后处理工作量。 因此，如何进行点云内插使其具有良好的局部抗噪性是目前亟待解决的问题。
[0005] 目前，现有的主流点云滤波算法可以分为两类：基于内插的滤波方法和基于形态 学运算的滤波方法：
[0006] (1)基于内插的滤波方法
[0007] 基于内插的滤波方法的基本思想是，首先在一较大的窗口内，提取最低点作为地 面种子点，构建不规则三角网（TIN)或数字地表模型（DEM)用来表达地面，然后根据适当阈 值渐进的对地面点进行加密。该方法可十分简洁的扩展为由粗到精的金字塔滤波策略，点 云金字塔中每一层级存储的是对应窗口内的最低值。然而，在不断内插地表模型与迭代滤 波过程中，地面点不可避免包含粗差点，若不有效处理将被不断放大，严重影响最终滤波结 果。目前，Terrascan采用一种渐进三角网加密的滤波策略，LasTools则采用线性预测的滤 波策略，两者均属于此类方法。
[0008] (2)基于形态学运算的滤波方法
[0009] 基于形态学的滤波方法的基本思想来源于数学上的形态学计算一一用形态学开 运算和闭运算从灰度影像上识别目标。建立在最低高程假设之上，非地面点的高程通常高 于周围的地面点，如果将点云数据按高程值转化为灰度影像，由于灰度差异，非地面目标将 被识别出来。对点云数据进行开运算，可以将局部窗口内高程值最小的点识别出来。通过判 断点云与窗口内的最小高程点之间的坡度差异与角度差异作为滤波准则，移动滤波窗口， 实现对全部点云数据的滤波。基于形态学运算的滤波方法的关键点在于选取一个合适的运 算窗口：若选择窗口过小，则只有较小的非地面目标（比如树木、汽车）能被滤除，而对于较 大地物则会在膨胀中再次保留；而若选择窗口过大，则在坡度较大区域，部分地面点也会被 误分为非地面点。该方法虽然效率较高，然而由于对参数选择十分敏感，难以处理复杂地形 点云，实用性不大。
【发明内容】
[0010] 本发明的目的在于克服上述缺陷，提供一种将能显式刻画地形变化特性的弯曲能 量引入到点云数据滤波过程中，使滤波参数随地形变化在局部范围内自适应优化调整；并 将地表一致性约束引入内插算法之中，采用弯曲能量的规则化约束，克服噪声对滤波精度 的影响，并建立一套科学合理的点云自动滤波实施方案。
[0011] 为了实现上述目的，本发明采用的技术方案如下：
[0012] -种自适应复杂地形结构的点云滤波方法，包括以下步骤：
[0013] (1)构建点云金字塔；
[0014] (2)采用格网窗口最大的金字塔层级数据初始化地面点G，并将下一级数据加入 待分类点U，更新待分类点；
[0015] (3)迭代处理点云金字塔的所有层级，处理完所有待分类点，并且保存新增加地面 点
[0016] (31)利用已有地面点，采用具有局部抗噪性的DEM内插算法获取地面DEM，并同时 获取对应区域的弯曲能量；
[0017] (32)根据DEM、弯曲能量以及金字塔尺度信息，采用自适应复杂地形结构的参数 优选方法获取滤波阈值并进行滤波，区分地面点与非地面点；
[0018] (33)重复步骤（31)~（32)，直至无新增地面点；
[0019] (4)存储分类后的地面点与非地面点，其中，地面点标记为2,非地面点标记为1。
[0020] 在构建点云金字塔之前，剔除粗差数据，其具体方式如下：
[0021] 首先，对点云数据构建KD-树的索引结构；然后，对每个点进行半径搜索，计算在5 米半径值内的点的数目，统计其均值μ和标准差σ，将邻域内的点数目小于μ-3 O的点视 为孤立点，进行移除。
[0022] 所述步骤（1)的具体方式如下：
[0023] (11)根据用户输入的最大窗口大小Wmax与最小窗口大小Wmin，均勾获取一系列金 字塔窗口大小w ;
[0024] (12)从最粗层级开始，对原始数据进行规则格网划分，每个格网中记录该格网的 最低点，并从原始点云中移除；
[0025] (13)迭代处理所有层级，建立点云金字塔。
[0026] 所述步骤（31)中具有局部抗噪性的DEM内插算法包括以下步骤：
[0027] (311)在获取所有地面点之后，迭代处理所有地面点，确定所有数据点的三维轴向 外包盒大小，用于后续确定DEM范围、行列数等信息；
[0028] (312)初始化EffiM和弯曲能量栅格结构，并根据金字塔层级信息，计算正则化约束 权重λ ;DEM格网原点、行列数、DEM窗口大小w等信息；该正则化约束权重与地面点中包含 的粗差点比例有关：若粗差比例越大，则权重越大；反之权重越小，直至λ = 〇
[0029] (313)对地面点的二维平面坐标，构建KD-树的空间索引，用于迭代处理DEM格网 的近邻索引，内插对应的格网信息；
[0030] (314)对每一 DEM格网，求取该格网处的二维平面坐标（X，Y)，在KD-树索引中按 自适应半径增长方法搜索至少9个近邻点，并且获取其对应的三维坐标（X，Υ，Ζ)，作为样条 函数的控制点；
[0031] (315)根据控制点三维坐标，更新样条函数矩阵，用于计算样条函数参数，样条函 数的解具有线性方程的闭合解形式，采用QR分解求解样条函数的参数信息；
[0032] (316)利用样条函数参数，计算DEM格网对应二维平面坐标（Χ，Υ)处的高程值ζ并 保存入DEM的栅格结构Z中；
[0033] (317)重复步骤（314)至（316)直至迭代处理所有DEM格网，获取最终的DEM与弯 曲能量栅格图。
[0034] 所述步骤（313)的具体方式如下：根据获取所有地面点的二维平面坐标（X，Υ)，采 用开源算法建立二维地面点数据的KD-树索引；其中KD-树中的索引序号与地面点序号相 同。
[0035] 所述步骤（314)的具体方式如下：
[0036] (3141)计算DEM格网点处的二维平面坐标（X，Υ);
[0037] (3142)获取当前窗口大小w，并且以此初始化一个数组r = [w, 2w, 3w, 4w