许多学者都提出了关于机载LiDAR技术的地面滤波方法^[1-2]。点云滤波大致可以分为4类：基于坡度的滤波方法^[3]、基于内插的滤波方法^[4]、基于数学形态学的滤波方法^[5]以及基于聚类的滤波方法^[6]。Sithole等^[7-8]认为聚类分割算法更适合机载LiDAR数据地面滤波，并提出对点云数据进行扫描线分割和基于局部点云特征的分割方法，通过比较分割后相邻集合之间的高度关系进行聚类处理。设地面块低于相邻的地物块。一般来说，在聚类过程中仅仅考虑地面块的拓扑关系和高度关系，可能会引起不合理的聚类结果或者有效信息丢失，因此有必要加入点云的特征数据进行聚类。为了从复杂地形中提取高精度DTM，本文提出一种基于一体化的聚类滤波方法，在基于点云块的聚类过程中，不仅考虑点集合之间的几何关系，而且考虑属性关系。最后通过实验证明了本文算法的有效性和可靠性。

1 一体化聚类滤波

一体化聚类滤波基本思想是在对原始点云数据进行八叉树(Octree)划分的基础上，基于平面拟合进行分割，在一体化聚类的引导下，完成点云聚类分割的粗分类，最后以分割区域代替单个点云数据进行地面三角网渐进加密。

1.1 点云分割和聚类 1.1.1 点云分割

Octree是对几何实体的体元剖分，具有相同的时间和空间复杂度，是等分空间。基于Octree进行点云分割，有助于后续一体化聚类的实现。具体步骤如下：

1) 创建一个Octree对象作为提取点云所用的搜索方法，再创建一个点云索引向量存储实际的点云索引信息，检测到的点云聚类保存在这里。

2) 在分割过程中，离散点云不断向下进行Octree划分，直到所有的节点子集都只包含1个平面，即位于拟合平面上的点云和距离拟合平面一定阈值范围内的点云都归为同一类点云。

1.1.2 点云聚类

由于空间属性和非空间属性是空间数据的固有双重特性，要求空间聚类时同类要素既要在空间上毗邻，又要在非空间属性上具有最大的相似度。这种考虑空间双重属性的思想称为一体化空间聚类。法向量、高斯曲率和平均曲率作为曲面形状的几何代表直接影响着点的非空间属性特征，因此作为点的特征向量进行聚类。

假设三维空间要素集P={p₁, p₂, …, p_n}，n≥2，p_i(x_i, y_i, z_i)中1≤i≤n，非空间属性维数为m。基于机载LiDAR点云数据的特征，通过最小二乘估计得到每个点的特征向量p_i=(x_i, y_i, z_i, a_i, b_i, c_i, K_i, H_i)，其中(a_i, b_i, c_i)和(K_i, H_i)分别是点的法向量和高斯曲率与平均曲率。因此，2个点之间的距离为：

(1)

式中，‖ ‖为范数, 1≤i, j≤n。

一体化空间聚类的具体实现过程如下：

1) 根据八叉树分割结果，得到初始集合k，进而计算聚类中心m₁(0), m₂(0), …, m_k(0)。

2) 根据式(1)计算每个集合和相邻集合的距离。如果距离少于设定的阈值，2个集合进行合并，否则转向下一步。

3) 根据步骤1，更新每个集合的中心，得到更新后的集合中心m₁(t), m₂(t), …, m_k(t)。

4) 重复步骤2，直到集合数量不再发生变化。

1.2 基于块的地面致密化

由于点云数据的复杂性和阈值设置的不完善，尽管一体化聚类方法在大部分区域效果很好，但是还是会有一些地面点不能被正确分类。因此在粗分类的基础上，本文采用基于块的地面三角网渐进加密进行精细迭代判断。

渐进三角网算法利用原始LiDAR点云数据逐渐加密生成三角网，根据一定的判别条件将参与构网的点作为地面点予以保留，而不能参与构网的点当作非地面点予以剔除^[9]。但该方法也存在不足之处：单一的滤波参数不可避免地会造成非地面点与地面点的错误分类。为克服此缺点，本文在基于原始点云的渐进三角网算法基础上，采用分割区域块代替原始单个点云数据进行三角网建立。

1.3 DTM生成

基于一体化聚类滤波生成DTM的流程见图 1。

首先，对原始点云数据进行预处理，剔除异常值。然后，在分析多次回波信息的基础上，选择单次和末次回波信息作为实验数据^[10]，根据一体化聚类滤波区分地面点和非地面点。最后，根据地面点生成DTM。

2 实验与分析 2.1 实验数据

采用ISPRS第三工作组的实验数据和样例数据进行实验，其点云数据的获取方式为Optech ALTM扫描仪，并且记录首次和末次回波信息数据^[11]。选取包含不连续地形特征的实验数据Csite2和样本数据Sample23、Sample24。样例数据中包含了许多不连续现象，比如高频率的地形起伏，中断的陡坡、山脊，间断的地形和不连续的陡坡、沟渠等，详细地形特征参见文献[11]。

2.2 滤波

图 2(a)~(d)为本文提出的滤波方法的滤波过程。首先，采用统计分析技术对原始点云数据进行去噪处理。然后，分析点云的多次回波信息。鉴于本文目的是提取地面点，而点云的多次回波信息反映了不同的地面特征，故首次回波和中间次回波不参与后续的滤波运算。由于选用的ISPRS数据仅提供首次回波和末次回波，故本文只采用末次回波数据进行实验。图 2(a)为对原始点云数据的预处理结果。

基于一体化方法进行空间聚类分割，所得地面聚类结果如图 2(b)所示。聚类合并阈值的设定根据点云密度和粗糙度分析实验得到。剩余点云数据如图 2(c)所示。最后，通过种子地面块建立三角网进行阈值判别，生成地面点，如图 2(d)所示。实验表明，一体化聚类滤波方法可以移除大部分的地物点，种子地面点包含的点云集合比较大，而待分类点所包含的点云集合较零散，符合地物分类特征。

2.3 DTM生成

比较样例数据、本文方法和商业软件Terrasolid所生成的3种DTM，图 3和图 4分别显示了Sample23和Sample24的DSM和由3种方法生成的DTM。从图 3(c)和图 3(d)可以看出，一体化滤波方法生成的DTM很好地保持了地形特征。和图 3(b)相比，尤其是椭圆和方形标识地方，本文方法更接近标准DTM。同样，图 4显示本文方法可以很好地保持了地形特征和地形起伏。

2.4 性能分析

定量评价滤波效果对高质量DTM的生成也是至关重要的。ISPRS发布的滤波方法比较报告将滤波误差分为3类^[8]：Ⅰ类误差是指将地面点误分为地物点的比例；Ⅱ类误差是指将地物点误分为地面点的比例；总误差是指错误分类的地面点和非地面点数量占整个点云的百分比率。利用这些指标可以评价算法的可行性。

表 1显示了样例数据通过一体化聚类滤波所得到的3类误差。图 5显示了15个样例数据通过经典滤波方法得到的总误差。表 1中，a为正确判断的地面个数，b为地面点错判为地物点的个数，c为地物点错判为地面点的个数，d为正确判断的地物点个数。3类误差的计算公式见文献[8]。

表 1和图 5定量显示了一体化聚类滤波的有效性。可以看出，本文提出的方法能达到比较高的精度，总误差低于大部分的经典滤波算法。因此，对于不连续地面，一体化聚类滤波算法可以很好地保持地形特征，但在稀疏植被地区，本文算法还有待于提高。

3 结语

本文引入一种一体化聚类滤波方法进行DTM自动生成，选用ISPRS提供的Csite2数据阐述滤波过程，用样例DTM、商业软件生成的DTM与本文生成的DTM进行分析对比，最后用3类误差对滤波效果进行定量评估。结果表明，本文方法可以在复杂环境中得到高质量的DTM，为进一步的LiDAR应用提供了基础。但本文提出的点云数据滤波算法仍然具有一定的局限性，容易将低密度点云误分为地面点。随着机载LiDAR设备获取地形信息逐渐丰富，获取点云数据的同时还可以获取影像及光谱信息，融合各种信息以提高滤波结果的准确性是下一步研究的方向。