三维激光扫描技术在大范围数字高程模型的高精度实时获取、城市三维模型重建、局部区域的地理信息获取等方面表现出强大优势[1]。原始的点云数据包含大量地物信息，应用于建设数字高程模型时，需要对点云数据进行滤波得到地面点云。目前，点云数据的地面点滤波方法主要包括坡度滤波算法[2]、数学形态学滤波算法[3]、最小二乘曲面拟合滤波算法等[4]。述滤波算法均存在各自的问题，例如坡度滤波阈值的自适应性低，且容易错误的选择地面种子；数学形态学滤波算法存在窗口尺寸不易确定和细节地形方块效应问题；最小二乘曲面拟合滤波算法则对拟合半径依赖性大。本文提出了一种自适应大小的二级网格与曲面结合的滤波方法，根据两种地形点云数据进行试验，对比常用的坡度滤波和形态学滤波，分析试验结果的1类误差、2类误差和总误差。1适应大小的二级网格构建方法1.1虚拟网格坡度滤波虚拟网格是对点云划区分块，将点云区域分成M×N个大小相等的矩形网格，使每个网格都包含离散的激光脚点，网格大小M和N的计算公式：M=xmax-xminS+1N=ymax-yminS+1 （1）式中：xmax和xmin——点云数据中x坐标的最大值和最小值；ymax和ymin——y坐标的最大值和最小值。确定M和N的值后，对每个激光脚点建立索引（X，Y），建立索引的具体计算方法：X=INT(x-xmin)Y=INT(y-ymin) （2）式中：INT——取整数；x和y——激光脚点的坐标；xmin和ymin——点云数据二维平面坐标的最小值。构建虚拟网格后遍历网格内所有的激光脚点，将每个网格内高程最低的激光脚点作为网格内的地面种子点。使用虚拟网格和地面种子点实现坡度滤波算法，设置坡度阈值T，计算虚拟网格内初始地面种子点P0和激光脚点Pi之间的坡度值，将坡度值小于阈值T的点分类为地面点，否则分为非地面点。Pi和P0两点之间的坡度值Gi：Gi=|hi0|×li0 （3）式中：li0——Pi和P0之间的欧式距离；hi0——Pi和P0两点之间的高差。1.2建立自适应大小的二级网格统计坡度滤波后的地面点，根据地面点重新构建虚拟网格。将每个虚拟网格中高程最低点选为初始地面种子点p0，再以p0点作为节点，对虚拟网格进行分割，将虚拟网格初步分割成4个不相等的二级网格。保证每一网格内激光脚点的个数保持在10个左右，因此合并小于10个激光脚点的二级网格，再对二级网格内的脚点建立索引[5]。构建二级网格：（1）选择坡度滤波后的地面点重新构建虚拟网格，网格边长S的大小由以下公式计算：S=mn×10 （4）式中：m——粗提取的地面点云区域的面积；n——地面点云激光脚点的数量。（2）在每个网格内选择初始地面种子点，以种子点为中心将虚拟网格初步分割成四个二级网格。（3）根据虚拟网格内激光脚点与初始地面种子点之间的位置关系建立二级网格的索引，数学关系为：a=xi-x00(yi-y0)≥0 , b=xi-x0≥0yi-y00c=xi-x00(yi-y0)≤0 , d=xi-x0≤0(yi-y0)0 （5）式中：（x0，y0）——地面种子点坐标；（xi，yi）——网格内激光脚点坐标；（a，b，c，d）——虚拟网格内激光脚点所对应的二级网格的索引编号。（4）统计每个二级网格内的点数，4个二级网格中有3个以上网格的点数小于10个，重新合并为一级网格，否则按顺时针方向将小于10个激光脚点的二级网格与相邻的二级网格进行合并，并建立合并后的网格索引，计算公式为：A=yi-y0≥0 , A=a⋃bB=xi-x0≥0 , B=b⋃cC=yi-y00 , C=c⋃dD=yi-y00 , D=d⋃a （6）式中：A，B，C，D——合并后二级网格的编号。1.3基于二级网格拟合二次曲面在每个二级网格中选取6个高程最低的地面点，采用二次曲面拟合的方法拟合曲面[4]：z=a0+a1x+a2y+a3xy+a4x2+a5y2 （7）式中：（x，y）——地面点坐标；z——拟合曲面后的高程；(a0 , a1 , a2 , a3 , a4 , a5)——多项式系数。多项式的系数可由求解多项式的最小残差平方和得到：Emin=∑i=1n[zi-a0+a1x+a2y+a3xy+a4x2+a5y2]2 （8）式中：Emin——最小残差平方和；zi——选取二级网格内的激光脚点的高程。E的最小值是由ai控制，ai对求偏导，将求解公式转换为求极值的过程：∂E∂a0=0 , ∂E∂a1=0 , ∂E∂a2=0∂E∂a3=0 , ∂E∂a4=0 , ∂E∂a5=0 （9）简化为矩阵算法：A=B-1C （10）式中：A——待定系数为[a0 , a1 , a2 , a3 , a4 , a5]T的矩阵；B和C——包含（xi，yi，zi）的矩阵。将每个虚拟网格内粗提取的地面点，带入到公式(10)中计算[a0 , a1 , a2 , a3 , a4 , a5]T的值，即求解出二次曲面的多项式的系数。1.4本文滤波方法本文滤波方法流程如图1所示。10.3969/j.issn.2096-1936.2021.18.025.F001图1算法流程二次曲面由粗提取的地面点拟合而成，必须剔除粗差点减少误差。计算激光脚点拟合前后均方根误差剔除粗差点：σ=vi-v¯T(vi∇-v¯)n-1 （11）式中：νi——激光脚点拟合前后差值；v¯——差值平均值；n——二级网格内的点数；σ——差值的均方根误差。将σ作为误差阈值，若地面激光脚点拟合前后的差值大于2倍的均方根误差，将此脚点判为粗差点进行剔除。剔除粗差点后重新计算二次曲面多项式，并且重新计算曲面点的均方根误差2σ。将坡度滤波得到的非地面点坐标代入对应的二次曲面多项式计算高程值，重新执行式（11），将νi小于2σ的点归类为地面点，否则归类为非地面点。2试验分析本文针对不同的地形区域进行试验，试验1使用平缓的城市区域点云，试验2使用地形起伏区域的点云。滤波结果评价标准使用ISPRS于2003年提出的滤波误差的评判标准，如表1所示。10.3969/j.issn.2096-1936.2021.18.025.T001表1点云滤波误差评判标准原始数据滤波后的数据误差统计地面点非地面点地面点ab1类误差［b/（a+b）×100%］2类误差［c/（c+d）×100%］总误差［（b+c）/n×100%］非地面点cd地势较为平坦的地区，地形坡度的值为0° ~10° ，地势起伏较大的地区的坡度阈值为10° ~30° ，试验1的坡度阈值取为10° ，坡度滤波网格大小为3 m，试验2的坡度阈值取为15° ，网格大小为5 m。2.1地势平缓区域试验试验数据是DublinCity数据集中的一部分，数据编号为T_315500_234500_NE。原始数据集如图2所示。图2点云数据集10.3969/j.issn.2096-1936.2021.18.025.F2a1（a）地面点云10.3969/j.issn.2096-1936.2021.18.025.F2a2（b）非地面点云点云总数12 198 716个，地面点5 957 008个，非地面点个数为6 241 708个。本文滤波方法与虚拟网格坡度滤波和数学形态学滤波相比1类误差减少19.589%、21.293%，2类误差减少15.573%、17.257%，总误差减少17.534%、19.202%。试验1的误差统计如表2所示。10.3969/j.issn.2096-1936.2021.18.025.T002表2试验1的误差统计滤波方法1类误差2类误差总误差虚拟网格坡度滤波31.37629.83530.587本文滤波方法11.78714.26213.053数学形态学滤波33.02631.51932.255%试验结果如图3所示。图3试验1结果10.3969/j.issn.2096-1936.2021.18.025.F3a1（a）虚拟网格坡度滤波10.3969/j.issn.2096-1936.2021.18.025.F3a2（b）数学形态学滤波10.3969/j.issn.2096-1936.2021.18.025.F3a3（c）本文滤波方法2.2地势起伏区域试验原始点云数据如图4所示。图4原始点云数据图10.3969/j.issn.2096-1936.2021.18.025.F4a1（a）地面点云10.3969/j.issn.2096-1936.2021.18.025.F4a2（b）非地面点云总点云数3 123 944个。地面点云和非地面点云通过人工进行分类，地面激光脚点1 225 880个，非地面激光脚点1 898 064个。本文滤波方法与虚拟网格坡度滤波和数学形态学滤波相比1类误差减少22.861%、28.916%，2类误差减少26.942%、23.047%，总误差减少25.351%、25.080%。地面点云结果如图5所示。图5地面点云结果10.3969/j.issn.2096-1936.2021.18.025.F5a1（a）虚拟网格坡度滤波10.3969/j.issn.2096-1936.2021.18.025.F5a2（b）数学形态学滤波10.3969/j.issn.2096-1936.2021.18.025.F5a3（c）本文滤波方法试验2误差统计如表3所示。10.3969/j.issn.2096-1936.2021.18.025.T003表3试验2误差统计滤波方法1类误差2类误差总误差虚拟网格坡度滤波29.34431.57930.712本文滤波方法6.4834.6375.361数学形态学滤波35.39927.68430.711%3结语本文介绍了自适应大小的二级网格与二次曲面相结合的点云的滤波方法，采用拟合二次曲面能够减少地面种子点选取造成的误差。由于三种滤波方法均不能删除建筑物的所有边缘轮廓点，后续可针对此问题进行深入研究。