1 引言

随着GPS和INS(inertial navigation system)组合定位精度的提高，使用LiDAR数据制作DTM(digital terrain models)已经变得非常普及。与传统的摄影测量方式相比，LiDAR数据采集对天气、季节以及时段要求较小，并且可以快速获取地表三维空间信息^{[1, 2]}。原始获取的LiDAR点云数据包括地面点和非地面点(可能包含少量噪声点)，如何高精度地分离地面点与非地面点是DTMs制作的关键技术。十多年来，人们将重点放在对单一滤波算法的研究上，具有代表性的算法大概分为以下几种：

(1) 以形态学为基础的滤波算法。该类算法通过定义一系列形态学算子进行地面点和非地面点分类^{[3, 4, 5, 6, 7]}，较为适合城区点云滤波，但其结果过分依赖移动窗口尺寸，对非平坦地面效果不理想。

(2) 以坡度为基础滤波算法。该算法以通过计算当前点与邻近点之间的高程差值和坡度值，并以此为依据判别地面点与非地面。文献[8]最早提出坡度的计算方法，该方法简洁易行，但在实践过程中滤波结果过度依赖所设置的坡度阈值，并且每个点的K邻近查询也相当耗时。

(3) 以点云平差为基础的线性预测滤波算法。该算法由文献[9]提出，并在Inpho公司(http://www.inpho.de)进行产品化。其基本思想是通过计算每个高程点与预测趋势面的残差，估计每个点的内插权重，通过权重的改变自动分离地面点与非地面点。线性预测滤波算法计算复杂度高，滤波效果受所选择权函数的影响较大，适合处理密林地区。

(4) 以内插不规则三角网为基础的渐进加密算法。该算法根据一定的地面光滑阈值条件自动剔除非地面点，某个版本已在TerrainScan (http://www.terrasolid.fi)中产品化^[10]，整体适应性较好，但为了达到理想效果，需要针对不同类型区域频繁调整阈值参数。

(5) 以扫描线为基础的滤波算法。该类算法多是通过比较扫描线上相邻点间的高程差异或坡度值而实现，优点是处理效率高，其中文献[11]设计了一种基于扫描线的坡度滤波方法，通过对原始LiDAR点云进行分割而进一步分离地面与非地面点，而文献[12]则提出基于一维双向扫描线标记(one-dimensional and bi-directional labeling,OBL)方法。为了克服一维扫描线邻接关系考虑不足这一缺点，文献[13]提出一种多方向扫描线滤波(multi-direction ground filtering,MGF)方法，从水平和垂直两个方向将一维扫描线扩展到二维扫描线，并取得一定效果。总体来看，大部分扫描线滤波方法只能选择有限个扫描方向，无法充分顾及足够真实的地面坡度信息，可能会在一定程度上影响滤波质量。

(6) 其他算法。文献[14, 15]先后使用两次回波之差进行激光脚点分类，尽管该方法一定程度上发挥了LiDAR点云自身优势，但仅仅使用回波信息进行点云分类难以适合大部分地形。文献[16]提出一种基于多分辨率方向预测的点云滤波方法，该方法在不同分辨率上以邻近点高程差作为滤波依据，其效果与基于坡度滤波方法较为接近。

ISPRS(International Society for Photogrammetry and Remote Sensing)组织学者对各类滤波算法进行比较研究，认为大多数滤波算法针对特定地形可以取得良好结果，没有哪一种算法适合各种地形条件，并指出未来的滤波算法需要融合多元数据分析、分类识别等辅助手段^[17]。影响城区滤波质量的两个主要因素是：① 将建筑物侧面离散点误分为地面点；② 将树干近地端的多次回波信号误分为地面点。为了克服上述两个因素影响，本文算法以内插三角网为滤波基础，采用面向对象分类方法对原始LiDAR点云进行地面点与非地面点初步分类，将分类结果作为先验知识自适应调整地面点构网条件，从而提高地面点与非地面点分离的可靠性。最后针对真实城区数据进行滤波试验以验证效果。

针对城区点云特点，采用一种基于知识引导的三角网渐进滤波方法：首先对点云进行规则格网化(其中格网间隔取点云平均间距)，并对栅格数据进行面向对象分割及地面点与非地面点粗分类处理；其次在已分类的地面对象中选取一定数量种子点，构建初始地面三角网，针对剩余地面点采用设定阈值渐进加密地面三角网；最后对非地面点数据进行强阈值内插构网，最终参与构网的点判定为地面点，未参与构网的点判定为非地面点。基于知识的点云滤波流程如图 1所示。

图 1 基于知识的点云滤波流程图 Fig. 1 Flowchart of urban point cloud filtering based on knowledge

图选项

数据集是徕卡ALS50激光雷达系统所采集的原始LiDAR点云，绝对地理位置位于河南省某地区。系统数据采集基本参数信息如下：航高约800 m；航飞速度约70 km/h；最大扫描角度绝对值约24°；最大回波记录次数4次；点阵密度约4点/m²，平均点距约0.5 m；点云平面误差小于1.0 m，高程误差小于0.5 m，总点数目384 293个。原始点云主体地物包含高层建筑、低矮建筑、树木、道路等，在三维环境下按高程渲染后的顶视图如图 2所示,白色方框标注区域为典型树木，对应三维透视图在右上方，树冠到地表均有不同密度点；椭圆框标注区域为典型建筑，对应三维透视图在右下方，建筑物侧边存在大量点。

图 2 点云按高程分层三维渲染后顶视图 Fig. 2 Top-view of point cloud rendering by different elevation

图选项

从分类学角度讨论，点云滤波是一个对点云进行地面点与非地面点精细分类的过程。但是单一的滤波算法很难得到理想结果^[17]，本文使用由粗到精的滤波思路：首先采用面向对象的分类手段将点云粗略分为地面点与非地面点两类，并以此作为知识指导后续精滤波处理。面向对象分类又分为面向对象分割及动态聚类两个子过程。

面向对象分类技术的基础是生成高质量初始图斑对象。本文采用一种拓扑启发式影像分割算法生成初始对象。拓扑启发式影像分割是一个自下而上，逐步合并的过程：分割由单个像元开始，通过启发式搜索方式寻找局部最优分割区域对，合并区域并维护拓扑网络的动态更新，迭代上述过程直到分割结束。分割结果如图 3所示，其中图斑以拓扑网络形式进行表达。

图 3 面向对象分割结果图 Fig. 3 Result of object-oriented segmentation

图选项

Otsu方法又称最大类间方差法^[19]，是在判决分析最小二乘法原理的基础上推导得出的自动选择阈值的二值化方法，基本思想是将图像用直方图分割成两组，当被分割成的两组方差最大时，此灰度值就作为二值化处理的阈值，数学表述如下：

设定阈值k，将灰度级L分为两组，C₀、C₁分别代表背景和目标；C₀概率为ω₀=ω(k)，C₁概率为ω₁=1－ω(k)；C₀均值为u₀，C₁均值为u₁，则两组的数学期望为

u=ω₀u₀+ω₁u₁

σ²(k)=ω₀(u₀-u)²+ω₁(u₁-u)²

从定性角度讨论，如果分类结果中类内方差越小，而类间方差越大，则认为分类结果是稳定的。在一定的特征空间中，Otsu方法使得分类后的两组数据类间方差最大化，在动态阈值选择上具有良好的适应性。依据上述原理设计迭代Otsu聚类算法对地面对象和非地面对象初步聚类。其中，分类特征空间为二维向量，即对象均值和标准差，二者算术和作为分类特征值。聚类算法步骤如下：

步骤1 标记所有对象为地面区域，计算所有地面区域的分类特征值；

步骤2 按特征值数值大小将地面对象进行由小到大排序；

步骤3 依次遍历对象列表，逐一计算类间方差数值并记录；

步骤4 选择类间方差最大值所对应对象的分类特征值作为分割阈值，将低于该分割阈值的对象标记为地面区域，高于该阈值的对象标记为非地面区域，并将相邻的非地面区域进行合并；

步骤5 重复步骤2到步骤4,当σ(k)小于设定阈值σ₀时，结束迭代。

σ₀的实际意义是分割所得到的两类间的最小差别，也表征着被分割地形的起伏程度。由于城区地形起伏较小，一般设为2.0 m。迭代Otsu方法聚类合并后，结果如图 4所示。图 4(a)为聚类合并后的拓扑网络效果图，其中大部分地物对象得到正确合并；图 4(b)是对地物对象进行填充后的效果图，可看出地面与非地面对象得到正确分离。

图 4 迭代Otsu方法聚类效果图 Fig. 4 Result of iteration Otsu clustering

图选项

上述聚类结果为一幅二维类别索引图，将点云数据进行逐一投影可获得各点的类别属性。依据各点类别属性进行基于知识的精细分类处理，其关键技术包括地面点判据与滤波策略两个部分。

文献[10]最先提出采用两类判据进行地面点与非地面点分类，将同时满足两个判据的点标记为地面点，并参与构网，不满足条件的点标记为非地面点，且不参与构网。假设当前点为P，待内插三角面元为平面ABC，P到ABC的距离为H，直线PA、PB、PC与平面的夹角分别为α、β、γ，如图 5，设定距离阈值为H_max，角度阈值为θ_max，则两类判据定义如下：

判据1 P点到平面ABC的距离小于设定阈值，即H＜H_max；

判据2 直线PA、PB、PC与平面ABC最大夹角小于设定阈值，即max(α,β,γ)＜θ_max。

图 5 判据计算示意图 Fig. 5 Schematic diagram of criterion calculating

图选项

基于知识的滤波方法仍沿用上述两类判据，并采用一定滤波策略进行判据阈值自适应调整，滤波策略描述如下：

步骤1 在每个地面对象中，分别选择区域内最低点作为种子点并构建初始地面三角网；

步骤2 对分类类别为地面点的数据进行优先构网，判据阈值为设定阈值；

步骤3 对分类类别为非地面点的数据进行构网，判据阈值为较强阈值，其取值应与扫描系统点位精度以及滤波精度相关，一般取H₀＜0.5 m、θ₀＜3°；

步骤4 输出参与构网的点云数据，即地面点数据。

基于知识的自适应滤波过程如图 6所示。图 6(a)中黄色圆状点为初始地面种子点提取结果，种子点基本上呈均匀分布；图 6(b)中黄色三角网为种子点所构建的初始地面三角网；图 6(c)为自适应三角网滤波结束后的地面点(非黑色区域)，其中设定阈值为H_max=2.5 m、θ_max=10°，强阈值为H₀＜0.5 m、θ₀＜3°；图 6(d)为与图 2中相对应的两处局部地物滤波后结果，即图(c)中白色框标注区域。

图 6 自适应滤波过程及结果 Fig. 6 Procedure and results of self-adaption filtering

图选项

以自主研发的海量点云编辑平台(SurfacePeak)进行全手工处理(包括断面编辑与真三维编辑)，处理后所得到的分类数据作为标准滤波结果，如图 7(a)；以传统渐进三角网算法滤波结果作为对比数据，且滤波参数同样为H_max=2.5 m、θ_max=10°，如图 7(b)所示。经比较可见，图 7(a)点云色彩分布均匀，地物边界滤波效果良好；而图 7(b)建筑物边界处存在少量残留斑块，呈红、黄色亮斑。

图 7 辅助数据效果图 Fig. 7 Results of reference data

图选项

文献[17]提出采用两类误差统计方法进行滤波算法质量的定量评价。两类误差统计可看做是分类混淆矩阵的特殊形式，其中第1类误差(Type Ⅰ)是地面点被错误分为地物点的百分比率；第2类误差(Type Ⅱ)是地物点被错误分为地面点的百分比率。基于知识的滤波算法与传统算法两类误差统计结果如表 1。

地物点被误分为地面点(Type Ⅱ)是影响DTM生产精度的主要因素，因此滤波分类的基本原则是控制第1类误差的同时，尽量降低第2类误差。依据统计结果可看出，本文算法第1类误差较传统算法降低5个百分点，且第2类误差得到较大改善。

点云滤波的一个重要目的是生产高精度数字高程模型(DEM)，两类误差统计是一种对算法质量评价的有效工具，但不能直接评价DEM生产精度。本文定义以下评价机制，对由滤波结果内插所获得DEM的生产精度进行评价：

假设DEM水平向长度为a，垂直向长度为b，二维空间R²={0＜x＜a;0＜y＜b}，则

(1) 在R²中随机选取一定数量坐标串，如Samp{(x₀,y₀),(x₂,y₂)，…，(x_n,y_n),n∈N}；

(2) 在待评价数据、比较数据及参考数据分别进行高程内插，对应高程序列为

(3) 统计待评价数据与参考数据差值序列Z_LF=|Z_L-Z_F|以及比较数据与参考数据差值序列Z_CF=|Z_C-Z_F| 对应各项精度指标，如最大残差、最小残差、残差均值以及残差中误差等。

不同滤波结果内插后DEM以及随机抽取1024个样本点的分布如图 8所示。

图 8 不同算法内插DEM与随机样本点分布 Fig. 8 DEM of different algorithms and distribution of random sample points

图选项

各项精度指标统计结果如表 2所示。

无论是基于内插DEM高程图像的定性评价，还是从DEM生产精度统计的定量评价，都表明在相同参数阈值条件下，本文算法较传统算法具有更好的滤波质量。

本文详细讨论了一种基于知识引导的三角网渐进滤波方法，对该滤波方法的各项关键技术进行深入探讨。针对较大数据量(约500 m×400 m)典型城区LiDAR点云数据进行试验，并采用多种判定手段进行质量评价，结果表明基于知识引导的滤波方法确实可以进一步提高滤波质量。