面点，确定所有数据点的三维轴 向外包围盒（Axis Aligned Bounding Box，AABB)大小，即三维点云坐标的最大最小值， (Xmin, Xmax, Ymin, Ymax, Zmin, Zmax)。计算方法是，迭代处理所有地面点，对每个点，分别比 较该点三维坐标（X，Y，Z)与当前AABB范围，若该点落入AABB之外，则对AABB进行更新；
[0086] (312)根据金字塔层级信息，DEM格网原点以及行列数，DEM窗口大小w，正则化约 束权重λ等信息，并且初始化DEM和弯曲能量栅格结构。正则化约束权重λ需要与地面 点包含粗差比例成正比，随着由粗到精的金字塔滤波的逐级递进，地面点包含噪声比例也 随之增大。在最粗的金字塔层级，假设无粗差点，窗口大小为《_，此时正则化约束权重λ =0 ;而在最精细的金字塔层级，此时粗差点最多，窗口大小为《_，此时正则化约束权重达 到最大值λ = λ_;中间层级，则采用线性内插方式获取，g卩
初始 化DEM和弯曲能量栅格图时，需要利用AABB平面最小值（Xmin，Ymin)，窗口大小w，以及行 列数目；
[0087] (313)对地面点的二维平面坐标，构建KD-树的空间索弓丨，用于迭代处理DEM格 网的近邻索引，内插对应的格网信息；具体方法为，根据获取所有地面点的二维平面坐标 (X，Y)，采用开源算法FLANN建立二维地面点数据的KD-树索引。其中KD-树中的索引序号 与原始三维点序号相同；
[0088] (314)对每一 DEM格网，求取该格网处的二维平面坐标（X，Y)，在KD-树索引中按 一定规则搜索一定数目的近邻点，并且获取其对应的三维坐标（X，Υ，Ζ)，作为样条函数的控 制点；在搜索过程中需要满足一下条件：1)样条函数控制点应分布均匀；2)控制点数目应 满足最小要求；3)控制点应选和该格网最相关的三维点，即和格网点最靠近的点。为满足 上述要求，具体实时方法如下：
[0089] (3141)计算DEM格网点处的二维平面坐标（X，Υ);
[0090] (3142)获取当前窗口大小w，并且以此初始化一个数组r = [w，2w，3w，4w，5w]，对 应半径搜索窗口大小；
[0091] (3143)从最小的搜索半径r = W开始，利用KD-树索引，搜索半径r中包含的点的 索引号，并且获取对应的三维坐标（X，Y，Z);
[0092] (3144)将所有三维点分为4个象限，每个象限中的三维点，按与DEM格网点距离大 小由小到大排序；
[0093] (3145)若每象限点数均大于等于4,则每个象限选取离格网点最近的三维点，作 为样条函数控制点；
[0094] 若其中一个象限点数小于4,则增大搜索半径r = Γι+1，并且重复步骤（3143)~ (3145)，直至 r = 5w ;
[0095] 若搜索半径r = 5w，近邻点仍小于9,则将该格网点设置为无效之-9999,并跳过后 续步骤（315)和（316);
[0096] (3146)返回搜索到的近邻点的索引号集合，并且获取对应的三维坐标，用于计算 样条函数的控制点。
[0097] (315)传统的DEM内插，通常是通过拟合地面控制点P = {pi = (Xi，Yi，Zi) I i = 1，2, 3,…，η}，得到一张参数化的曲面Z = f (X，Y)，在拟合过程中，由于地面控制点中不包 含粗差，因此，仅需要保持对地面控制点有较好的拟合度，就可取得较好的效果。即使得数 据拟合度edata(data term)具有较小的值，或完全经过所有控制点，如下式所示
[0099] 然而，在点云滤波过程中，地面点不可避免的会引入非地面粗差点，因此此时即使 ε ^3为〇,也并不意味着对地面有较好的拟合程度。采用地表一致性约束，可以一定程度 上，降低噪声敏感性，提高内插算法的鲁棒性。弯曲能量es_th，如下式所示，
[0101] 可以表达一张参数化曲面的二阶连续程度，因此若将ε s_th作为一个正则化约束 引入地面内插过程中，则可隐式的蕴含着地表一致性约束，抵抗噪声点的干扰。即通过优化 具有ε = edata+X es_th形式的能量函数，进行DEM内插。
[0102] 给定控制点集P，该优化模型为的解为薄板样条函数，即下式所示，
[0104] 为获取该样条函数的参数罗和歹，，可通过以下形式的线性方程组求解：
[0106] 其中，上述矩阵中的各元素具有以下形式，并且可直接采用QR分解的方式获取最 终的参数声和α &
CN 105118090 A I兄明 8/9 页
[0112] (316)利用样条函数参数，计算DEM格网对应二维平面坐标（X，Y)处的高程值z并 保存入DEM的栅格结构Z中；其具体步骤如下：
[0113] (3161)计算DEM高程值：将DEM格网点坐标（X，Υ)和样条函数参数彦:和足，带入 函数Z (Χ，Υ)之中，即求得EffiM对应格网处的高程值，并存入DHM对应的行列号之中；
[0114] (3162)计算DEM格网处的弯曲能量：根据样条函数参数#以及矩阵Κ，计算样条函 数的弯曲能量二户'/φ 5并存入弯曲能量对应的行列号之中；
[0115] (3163)直至计算完所有格网点的信息，并将计算完成的DEM与弯曲能量传出至下 一阶段处理。
[0116] (317)重复步骤（314)至（316)直至迭代处理所有DEM格网，获取最终的DEM与弯 曲能量栅格图。
[0117] 如图5所示，步骤（32)中自适应复杂地形结构的参数优选方法包括以下步骤：
[0118] (321)依据金字塔层级信息s，计算当前层级对应的尺度补偿阈值Ts ;考虑到 地形坡度影响，在较粗的金字塔层级，由于窗口 w较大，即使平坦区域也可能会存在微小 的地形起伏，因此需要根据金字塔层级S和窗口大小W补偿地形坡度起伏影像。具体 而言，在最粗的金字塔层级，补偿阈值最大为Ts = Ts_max ;在最精细的金字塔层级，补 偿阈值为〇;而在中间层级的数据，补偿阈值根据窗口大小w采用线性内插方法获取，BP
[0119] (322)将弯曲能量栅格图采用分段线性内插转化为对应的弯曲能量补偿阈值Tb。 顾及弯曲能量的参数自适应优选方法的基本思想是，在地形起伏较大的地方，如山脊线、地 面断裂处等区域，滤波阈值应更大；而在地形起伏较小的区域，如平坦地表等区域，滤波阈 值应更小。由于在地形欺负较大的地方，地形变化弯曲能量也相应的会变大，因此滤波阈值 应与弯曲能量成正比例变化。因此，其具体获取弯曲能量补偿阈值Tb的具体步骤为；
[0120] (3221)对所有DEM格网点的弯曲能量从小至大排序；
[0121] (3222)将所有弯曲能量均分为10段，并记录分割处的弯曲能量大小
[b0, bu....., b10]；
[0122] (3223)设置一最大补偿阈值Tbmax，对所有弯曲能量值b，进行分段内插，若其落入 第i段之内，贝1J有
[0123] (323)内插待分类点所处格网的高程均值Zni与弯曲能量补偿值T b。由于DEM在滤 波过程中可能会存在一定噪声，因此不能仅用该格网处的DEM高程值，而考虑采用3X3邻 域DEM的高程值z和弯曲能量补偿值Tb。首先，根据待分类点坐标（X，Y，Z)，计算其所对应 的格网行列号，r= (Y-Ynin)/w，C= (X-Xnin)/￥。随后，在DEM和补偿阈值中获取对应3X3 邻域内的高程均值zn和补偿阈值T b;
[0124] (324)计算自适应滤波阈值，该阈值综合考虑金字塔层级的尺度信息s，弯曲能量 信息b，以及初始滤波阈值t，故最终的滤波阈值应是上述三者的函数T = f (s，b，t)，即T = Tb+Ts+t ;
[0125] (325)依据待分类点高程，DEM高程与滤波阈值之间关系判断该点为地面或非地 面点；在3x3邻域内，判断该待分类点和邻域内的9个格网点的相互关系。记N为待分类点 高程值Z〈z+T的数目，若N多5,则该待分类点为地面点；否则仍然保留为待分类点，在后续 处理中继续