, 5w]，对 应半径搜索窗口大小；
[0038] (3143)从最小的搜索半径r = w开始，利用KD-树索引，搜索半径r中包含的点的 索引号，并且获取对应的三维坐标（X，Y，Z);
[0039] (3144)将所有三维点分为4个象限，每个象限中的三维点，按与DEM格网点距离大 小由小到大排序；
[0040] (3145)若每象限点数均大于等于4,则每个象限选取离格网点最近的三维点，作 为样条函数控制点；
[0041] 若其中一个象限点数小于4,则增大搜索半径r = Γι+1，并且重复步骤（3143)~ (3145)，直至 r = 5w ;
[0042] 若搜索半径r = 5w，近邻点仍小于9,则将该格网点设置为无效之-9999,并跳过后 续步骤（315)和（316);
[0043] (3146)返回搜索到的近邻点的索引号集合，并且获取对应的三维坐标，用于计算 样条函数的控制点。
[0044] 所述步骤（316)的具体步骤如下：
[0045] (3161)计算DEM高程值：将DEM格网点坐标（X，Y)和样条函数参数#和5，带入 函数Z (Χ，Υ)之中，即求得EffiM对应格网处的高程值，并存入DHM对应的行列号之中；
[0046] (3162)计算DEM格网处的弯曲能量：根据样条函数参数J5以及矩阵K，计算样条 函数的弯曲能量6 = 并存入弯曲能量对应的行列号之中；
[0047] (3163)直至计算完所有格网点的信息，并将计算完成的DEM与弯曲能量传出至下 一阶段处理。
[0048] 所述步骤（32)中自适应复杂地形结构的参数优选方法包括以下步骤：
[0049] (321)依据金字塔层级信息s，计算当前层级对应的尺度补偿阈值Ts;
[0050] (322)对弯曲能量栅格图，统计其最大值与最小值，并采用分段线性内插的方式， 将弯曲能量转化为对应的弯曲能量补偿阈值Tb;
[0051] (323)对每一待分类点，计算其落入对应层级的格网行列坐标，采用双线性内插方 式，获得对应3 X 3邻域DEM与弯曲能量栅格图上的高程均值Zni和弯曲能量补偿值T b;
[0052] (324)计算自适应滤波阈值，该阈值综合考虑金字塔层级的尺度信息s，弯曲能量 信息b，以及初始滤波阈值t，故最终的滤波阈值应是上述三者的函数T = f (s，b，t)，即T = Tb+Ts+t ;
[0053] (325)依据待分类点高程，DEM高程与滤波阈值之间关系判断该点为地面或非地 面点；
[0054] (326)重复步骤（323)至（325)直到处理完所有待分类点，并且保存新增加地面 点。
[0055] 所述步骤（321)中阈值计算方式如下：
[0056] 最粗的金字塔层级，补偿阈值为Ts= T smax ;
[0057] 最精细的金字塔层级，补偿阈值为0 ;
[0058] 中间金字塔层级，补偿阈值根据窗口大小w采用线性内插方法获取，即
[0059] 所述步骤（322)获取弯曲能量补偿阈值Tb的具体步骤为：
[0060] (3221)对所有DEM格网点的弯曲能量从小至大排序；
[0061] (3222)将所有弯曲能量均分为10段，并记录分割处的弯曲能量大小
[b0, bu....., b10]；
[0062] (3223)设置一最大补偿阈值Tbmax，对所有弯曲能量值b，进行分段内插，若其落入 第i段之内，贝1J肩
[0063] 与现有技术相比，本发明具有以下有益效果：
[0064] (1)从同时包含地面-非地面点的点云数据中，如何内插出真实地表将直接影响 到整体的滤波精度。由于三维点云提供的信息较为有限，只有顾及地面一致性约束，才能一 定程度上克服噪声点对DEM内插带来的影响。然而现有通过聚类、平面拟合的约束方法并 不具有通用性。因此，本发明依据地面二阶连续平滑的物理意义，直接在内插过程中引入表 示平滑程度的弯曲能量正则化约束，在内插过程中顾及地表一致性约束，降低噪声点对DEM 内插的影响，整体提高滤波精度。
[0065] (2)在得到表达地表形状的DEM之后，需要确定区分地面-非地面点的阈值参数， 现有滤波算法通常对同一区域采用相同的参数，而不考虑地形结构信息，因此对参数变化 十分敏感。已有研究已经指出，在滤波过程中需要考虑地形结构信息以提高滤波性能的稳 定性，现有的方法通常仅是对3 X 3窗口的DEM值做简单的加权平均处理。而采用本发明，弯 曲能量能有效表达局部地形结构信息的，采用弯曲能量使滤波参数自适应的随地形改变， 可以显著提高滤波算法的可靠性，降低对滤波参数的敏感程度。
【附图说明】
[0066] 图1为本发明各模块之间关系。
[0067] 图2为本发明算法总体流程图。
[0068] 图3为顾及弯曲能量正则化约束的局部抗造性内插方法。
[0069] 图4为DEM与弯曲能量示意图。
[0070] 图5为顾及弯曲能量的滤波参数自适应优选方法。
【具体实施方式】
[0071] 下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于 下列实施例。
[0072] 实施例
[0073] 本实施例提供了一种自适应复杂地形结构的点云滤波方法，深入挖掘弯曲能量在 DEM内插与表达地形起伏变化的特殊作用，在地表二阶连续平滑的假设之上，在DEM内插过 程中利用弯曲能量作为正则化约束，克服噪声影像。并且利用弯曲能量显示刻画地形起伏 变化，用于滤波参数的自适应优选。在此基础之上，设计一套自动化的滤波流程，采用由粗 到精的金字塔滤波策略，提高滤波效率和可靠性。
[0074] 本实施例中，涉及的各模块关系如图1所示，在由粗到精的滤波过程中，每次迭代 包含以下两个主要步骤：1)顾及弯曲能量正则化约束的局部抗噪性DEM内插算法和2)顾 及弯曲能量的滤波参数优选方法，本发明的总体流程如图2所示，具体包括以下步骤：
[0075] (1)构建点云金字塔；基于地面最低高程假设，在逐级递进的滤波过程中，金字塔 中存储的点云数据高程值在局部范围内将逐渐增加；
[0076] (2)采用最粗层级的金字塔数据初始化地面点G，并将下一级数据加入待分类点 U，更新待分类点；
[0077] (3)迭代处理点云金字塔的所有层级，处理完所有待分类点，并且保存新增加地面 点，即对每一层级作如下处理：
[0078] (31)利用已有地面点，采用具有局部抗噪性的DEM内插算法获取地面DEM，并同时 获取对应区域的弯曲能量；
[0079] (32)根据DEM、弯曲能量以及金字塔尺度信息，采用自适应复杂地形结构的参数 优选方法获取滤波阈值并进行滤波，区分地面点与非地面点；
[0080] (33)重复步骤（31)~（32)，直至无新增地面点；
[0081] (4)存储分类后的地面点与非地面点，保存为.Ias格式数据，其中，地面点标记为 2,非地面点标记为1。
[0082] 在构建点云金字塔之前，还需剔除粗差数据，点云粗差通常来自于激光雷达的多 回波效应以及密集匹配噪声，考虑到这些粗差点通常游离于点云数据之外，且通常数量稀 少，表现为孤立数据，因此，可以考虑采用其邻域内的点数的统计信息进行移除。首先，对点 云数据构建KD-树的索引结构，随后对每个点进行半径搜索，计算在某一半径值内的点的 数目，统计其均值μ和标准差σ，将邻域内的点数目小于μ-3 σ的点视为孤立点，进行移 除。
[0083] 构建点云数据金字塔：根据用户输入的最大窗口大小Wmax与最小窗口大小w min，均 匀获取一系列金字塔窗口大小w。从最粗层级开始，对原始数据进行规则格网划分，每个格 网中记录该格网的最低点，并从原始点云中移除。迭代处理所有层级，建立完整点云金字 塔。
[0084] 如图3所示，步骤（31)中具有局部抗噪性的DEM内插算法包括以下步骤：
[0085] (311)获取所有地面点之后，迭代处理所有地