激光雷达技术，在很多基础经济建设如农林业^[1]、水力电力勘察^[2]、道路设计^[3]及城市规划^[4]中已经被普遍使用，而点云滤波分类和特征提取是激光雷达点云数据处理的重要部分，也是后续各类数字产品生产和应用的基础，因此很多学者开展了相关的滤波方法研究，归纳起来可分为以下几类：基于坡度的滤波方法^[5-9]、基于数学形态学的滤波算法^[10-14]及拟合法^[15-20]等。

基于坡度的滤波方法能够较好地保留地形信息，但是坡度阈值的自适应性有待提高；基于数学形态学的滤波算法存在的主要问题是窗口尺寸的人工选取和细节地形的方块效应，拟合法则对拟合半径的依赖性大。在拟合法的基础上，许多学者提出了移动曲面拟合的滤波算法，文献[21]提出了移动曲面拟合滤波算法，该算法简单明了，自适应性强，运算速度快，但是固定的窗口尺寸和阈值大小会影响滤波效果。文献[22]提出的改进移动窗口曲面拟合法可以自适应地调整网格尺寸大小，同时克服滤波过程中高差阈值难以确定的不足，但是该算法使用6个格网中的最低点进行拟合，容易出现方块效应。文献[23]提出的多级移动曲面拟合滤波算法，建立区域格网索引，设置自适应阈值。该方法普适性较好，效率高，解决了粗差点、阈值自动设置等问题，但预处理需要人为设定参数，并且未考虑每个分块的邻域信息。

在综合分析上述算法优点的基础上，本文提出多级移动曲面拟合的自适应阈值点云滤波算法。该算法增加了格网邻域并实现格网大小、邻域大小和阈值的自动设置，有效提高了滤波效率，降低了滤波误差，提高算法的自适应性。

1 算法原理与流程

传统移动曲面拟合的滤波算法是利用2×3的移动窗口，通过寻找每个窗口的最低点计算出一个粗略的拟合曲面，然后判断这6个移动窗口中的点与拟合曲面的高程差，将所有高程差超过给定高差阈值的点剔除，得到地面点。传统算法仅用6个点拟合曲面会导致拟合的曲面未更好地考虑周围地形的影响，与真实的地形相差比较大，从而影响滤波效果，而且在滤波过程中，对于不同的地形采用固定窗口尺寸和高差阈值，其滤波效果会更差。

针对上述问题，本文考虑到地形的连续性，增加了窗口邻域，认为每个窗口的拟合曲面均是基于邻域窗口的最低点获得，最小邻域为3×3(点数均大于6)。为了解决格网大小的选取问题，采用网格尺寸按倍数逐级减小的方法，在改变网格尺寸的同时自动设置不同的窗口邻域和高差阈值，具体流程如图 1所示。

1.1 数据预处理

预处理主要是剔除异常点(噪声点)，本文通过对点云数据生成频率直方图，设定阈值以去除首尾两端的极大值和极小值，实现去噪的目的。

1.2 格网分割及索引

点云格网化也称点云分割，即用规则的格网将点云数据全部覆盖，对于不同的地形，采用不同的网格大小，将点云按平面坐标分配到相应格网中，格网的划分及属性如图 2所示。

将点云数据进行编号，按照定义的格网大小划分，为了保证每个格网中的最低点为真实的地面点，首先设置格网大小大于最大建筑物的尺寸以去除具有较大尺寸的地物点，然后按倍数逐级减小格网大小，直到剔除多数地物点，拟合出最佳地形表面为止。其中每个格网有4种不同的属性：

(1) 点集属性。计算每个点对应的格网，并记录点号，从而得到落在每个格网中点的ID及点数。

(2) 最低点ID。用邻域格网中的最低点进行曲面拟合，因此记录每个格网的最低点标号，可以提高计算效率。

(3) 边界坐标。为了保证计算结果的准确性，划分格网时记录格网的边界信息，在拟合曲面时，将所有邻域格网的坐标数据减去邻域格网中左下角坐标。

(4) 标号属性。用二维矩阵的行列号对每个格网进行标号，在后面的计算中根据该标号进行点云资料的快速存储，提高滤波、插值的计算速度。

1.3 建立拟合曲面

假设地形表面是一个复杂的空间曲面，该曲面的局部面元可以用二次曲面来逼近，如式(1)所示

真实的地面点得到的拟合曲面反映了地面的起伏，将每个格网中的最低点作为真实地面点，利用每一个格网的邻域窗口中的最低点拟合曲面，然后计算该格网中点的拟合高程值与真实高程值的差值，通过高差阈值进行筛选，每一个格网都用不同的曲面方程来表达，并且该曲面方程具有一定的连续性。该过程的关键是格网邻域大小的自动设置以及拟合曲面参数(a₀, a₁, a₂, a₃, a₄, a₅)的求解。

1.3.1 格网邻域大小的设置

用于拟合地形表面的二次曲面共有6个参数，许多研究都是基于6个格网中的最低点拟合曲面，然而由于数据分布的不均匀性，有些格网中没有点，因此最低点的个数小于或等于参数的个数，得到的拟合曲面与真实的地形相差较大，滤波效果降低。为了保证拟合曲面的连续性，需要考虑每一个格网周围格网中点云数据的分布特征。为了保证该格网在中间位置，一般设置邻域大小为奇数如3×3、5×5，如图 3所示。邻域3×3即利用9个格网中的最低点拟合一个曲面，用此曲面判断格网(i, j)中的点云数据。

本文邻域大小的选择根据地形坡度设置，其中平均坡度阈值(5°，10°)的选取建立在大量试验基础之上。首先计算每个格网对于3×3邻域格网的平均坡度值，即每个格网的最低值和周围8邻域格网最低值的差值与距离比值的平均值，如式(2)所示

式中，Z_min (i, j)是第i行第j列格网的最低点高程值。

如果存在格网的坡度平均值超过0.176(tan10°)，则邻域大小设置为3×3，否则计算每个格网对于5×5邻域格网的平均坡度值，如果存在格网的坡度平均值超过0.087(tan5°)则设置格网大小为5×5，否则设置为7×7。

1.3.2 曲面参数的求解

本研究利用邻域窗口的方法增加了求解参数的点数，即方程的个数大于待求解参数个数，利用广义逆矩阵及最小二乘法求解曲面方程的参数。

1.4 高差阈值的确定

在改变网格尺寸和设置不同的窗口邻域大小时，为了有效地将地面点和非地面点分开，需要设置不同的高差阈值。

根据文献[24]提出的一种假设：大量离散的地面点云数据在自然状态下呈正态分布，而非地面点数据则影响了正态分布的状态。因此，相对于地面点，非地面点为误差点，其标准偏差大于3σ。由此确定总体数据的高差阈值，然而每个格网包含的地物特征不同，不能直接对整块数据选取同样的阈值。考虑到地面点和地物点具有分层现象，因此可以采用谱系聚类的思想对拟合高程差值进行分类，根据不同的分类数，确定高差阈值^[25]。然而有些格网中的地面点与地物点分层现象不明显，因此将两种方法进行综合，具体步骤如下：

(1) 总体高差阈值的确定。假设总数据有N个点，计算每个点云与其对应拟合曲面方程的拟合高程差值，并将N个拟合高差数据从小到大顺序排列，设x₁ < x₂ < … < x_N，计算相邻数据间距离l_{i(1≤i≤N-1)}和其期望值μ_l及标准方差σ，则总体阈值为x₁+3σ。

(2) 每个格网高差阈值的确定。以其中一个格网为例，假设该格网有n个点，按照与总体高差阈值相同的步骤计算出相邻两个数据间距离，并按照从大到小顺序排列成l₁≥l₂≥l₃≥ …≥l_n-1利用类间距离l₁，将格网中的全部数据A分为A₁和A₂两类，两类的直径(同一类中最大数据和最小数据的差)分别为Q₁和Q₂，两类的中心(同一类中最大数据与最小数据的算术平均数)分别为O₁和O₂，若这些参数满足下面其中一个条件，则分为以下两类，否则不进行分类：① l₁≥max{Q₁，Q₂}，即类间距离相对较大，A₁和A₂为不同的类别；② l₁＜max{Q₁，Q₂}，l₁>min{Q₁，Q₂}, 且(O₁+O₂)/2∈A₁或A₂，即两类中心的算术平均值落在其中一类之中，说明这一类数据集中于两类中心的算术平均值附近，而另一类则相反，两类数据特征相异。

高差阈值的取值与分类数相关，其取值规定如下，设A中相邻两个数据间距离最大的两个值为l₁和l₂, 且l₁≥l₂：①当A不分类时，数据间的分层现象不明显，高差阈值设为总体高差阈值；②当A分成两类，且两类的最大直径为Q^*时，当Q^*≤l₁≤2Q^*时，高差阈值取值为l₂处的拟合高差值，其他情况下，高差阈值取值为l₁处的拟合高差值；

根据每个格网本身的实际情况，在多级滤波的情况下自动设置每个格网的高差阈值，最大可能避免了以往手动设置同一高差阈值而产生的过度滤波或是不完全滤波问题。

2 试验结果及分析 2.1 试验数据及预处理

为了验证本算法的精度，试验数据采用ISPRS公开发布的激光雷达数据集。该数据为Optech ALTM系统采集，由8块测试数据组成。这些测试数据包含代表性的地形地物特征，ISPRS提供了8块测试数据的15个样本数据，并且对样本数据进行了手工分类，将激光脚点精确分类为地面点和非地面点。

对样本数据集进行多级滤波，首先根据地物的最大尺寸设置第1次滤波的初始格网大小，之后滤波的格网大小按倍级减小，15个样本集的数据特征、初始格网大小设置见表 1。

2.2 滤波结果

滤波的前提条件是每个格网的最低点为地面点，然而有的样本数据存在异常点，因此在滤波之前需去除异常点，并将其作为地物点，之后再进行多级滤波。

为了体现多级滤波的处理过程，以sample 12为例，展示每级滤波的结果，其初始格网大小为20，异常点个数为81，共进行3级滤波，每级滤波的地面点和非地面点结果分别如图 4和图 5所示，图 4中白色的标记点为第1类分类错误的点，图 5中白色的标记点为第2类分类错误的点。

从图 4和图 5可以看出，第1次滤波有大量的非地面点被分为地面点，而非地面点基本无地面点，说明第1次滤波窗口尺寸及高差阈值设置大，去除了具有较大尺寸的高大建筑物；第2次滤波地面点中的非地面点明显减少，一些低矮的建筑物及建筑物之间的植被被剔除，而有少部分的地面点被当作非地面点；第3次滤波地面点中的非地面点进一步减少，而一部分地面点被当作非地面点，采用分级滤波可以逐步滤除地物点，保持地面细节信息，达到算法设计的要求。

2.3 滤波结果分析与评价

为了对滤波结果进行定量评价，使用ISPRS在2003年提出的滤波误差的评判标准。第1类误差为将地面点分类为非地面点的误差，第2类误差为将非地面点分类为地面点的误差，第3类为总误差即第1类和第2类误差点数占总点数的比例，如表 2所示。其中e和f表示参考数据中地面点和非地面点的个数，g和h分别是滤波后数据中地面点和非地面点的个数，n为扫描数据点总个数，b/e×100%为第1类误差，c/f×100%为第2类误差，(b+c)/n×100%为总误差。

根据滤波误差定义计算15个样本的滤波误差，如表 3所示。结果表明，文中滤波方法对于sample 12、sample 21、sample 42、sample 54、sample 61有相对好的结果，总误差均在4.00%以下，其中sample 61中包含陡坡，说明该方法对于连续的山地区域可以达到较好的滤波效果。而sample 11、sample 23、sample 24、sample 41的滤波结果较差，根据这些样本数据的数据特征可知，sample 11、sample 23、sample 24、sample 41均存在不连续的地形。文中滤波方法是通过拟合地形的曲面方程，判断点云到拟合曲面的距离是否在高差阈值范围内，对于不连续的地形(如sample 11中的梯田)拟合出来的曲面与实际的地形相差会很大，因此滤波效果相对较差。

为了进一步分析滤波算法的准确性，将总误差与ISPRS公布的8大经典滤波算法^[26-33]进行比较，结果见表 4。从误差平均值结果中可以看出文中滤波方法精确度优于7种滤波算法，其中sample 12和sample 61的效果最好。本文算法在地形连续的城市区域和山地区域可以取得较好的效果。

整体来说，文献[26]提出的滤波算法效果比本文方法的效果好，但是其滤波计算量大，需要反复构建三角网并进行高程、角度的计算，速度较慢。其他7种算法的整体效果没有本文滤波效果好，其中文献[28]、文献[27]、文献[30]、文献[32]提出的滤波算法在山地区域滤波效果较差，而本文滤波算法考虑了地形的连续性信息，因此在地形连续的山地区域可以达到较好的滤波效果。文献[28]算法在无地形起伏区域(sample 21、sample 31、sample 42)效果较好，因为该算法多次迭代计算拟合地形高程，充分考虑了地形和地物间的联系。文献[31]、文献[33]滤波算法对地形突变地区滤波难度较大，且其滤波算法阈值的确定需要根据先验知识判断，不同的地形需设置不同的坡度阈值，坡度阈值的自适应性较差，而本文算法是基于自适应阈值的，根据数据本身的特征选取合适的阈值，因此能达到较好的效果。

3 结论

本文采用基于多级移动曲面拟合的自适应阈值点云滤波方法在传统滤波算法的基础上增加了窗口邻域，根据区域范围中最大建筑物尺寸设置初始格网大小，并采用多级滤波的方法，在逐级改变网格尺寸的同时, 根据坡度自动设置窗口邻域，根据拟合高差值自适应性设置高差阈值。试验表明本文算法精度高，具有一定的自适应性。与ISPRS公布的8大滤波算法的总误差对比分析，比其中7种滤波算法的效果好。总体来说，本文算法在地形复杂且连续变化的区域能够获得很好的滤波效果，不足之处在于对不连续地形的区域滤波效果有待进一步提高。