1 引言

激光雷达测量技术是从20世纪中后期发展起来的高新技术,大多采用主动扫描方式,广泛应用于电力线提取与重建^[1]、道路设计^[2]、林业^[3]、环境监测^[4]等。点云滤波是雷达数据处理的核心步骤之一,主要任务是分离地面点与非地面点,其结果直接影响特征提取、点云重建等后续工作的精度,因此很多学者对滤波算法展开了研究。目前较为成熟的滤波算法有以下几类:基于表面的滤波算法^[5-6]、基于坡度的滤波算法^[7]、基于形态学的滤波算法^[8-9]、基于曲面拟合的滤波算法^[10]。基于表面的滤波算法能有效地去除低矮物体,具有较好的地形特征。但该方法需要反复构建三角格网,通过角度、距离来判断点和三角网的归属问题,计算量大,耗时长。基于坡度的滤波方法的效率高、计算简单,但对于坡度阈值过于敏感,且单一坡度阈值无法满足自适应问题。基于形态学的滤波算法难以确定窗口尺寸,窗口过大、过小分别会导致滤波不足与过度滤波,窗口边界存在板块化问题。基于曲面拟合的滤波算法简洁明了,运行速度快,自适应强,但对拟合点要求高,对拟合范围依赖性强。因此如何选取最优种子点生成拟合曲面、设置自适应高差阈值是决定滤波效果的关键。

文献[ 10]提出了经典的移动曲面算法,但是该方法是直接选取未经筛选的6个最低点进行曲面拟合,拟合效果不佳,同时固定的高差阈值无法自适应真实的起伏地形。文献[ 11]通过置信区间的方法对选取的种子点进行判断,若种子点在置信区间内,则视为无粗差点,其原理是在自然状态下地面点服从正态分布。虽然该方法能使整体上的种子点满足要求,但是拟合曲面的生成都是基于局部面元,拟合曲面和真实地形的相似度提升不明显。文献[ 12]采用多级滤波策略,通过聚类方法对高差值进行分类,根据不同分类数确定自适应高差阈值。该方法有较好的滤波效果,但是通过频率直方图选取的拟合种子点仍然存在粗差问题。文献[ 13]通过不断扩大计算窗口逐步过滤地物点,但是通过直接剔除异常点的方式无法去除低粗差点,存在的粗差种子点在后续迭代运算中会严重影响滤波精度。

针对移动曲面算法存在的种子点粗差问题,本文提出一种多级种子点优化策略。首先确定无粗差一级种子点,然后利用一级种子点对备选种子点进行筛选、生长,最后使用无粗差种子点进行曲面拟合并采用顾及地形的自适应高差阈值来解决采用唯一经验高差阈值时存在的滤波效果不佳与自适应差问题。

2 算法原理与流程

经典移动曲面滤波算法利用6个移动窗口寻找6个最低点,然后利用这6个点进行曲面拟合,判断6个窗口内所有点的拟合高程和真实高程的差值与高差阈值的大小,所有高差小于阈值的点则视为地面点。针对粗差种子点和固定高差阈值问题,本文将曲面拟合范围定为3×3窗口(有3个约束条件),然后对3×3窗口内的种子点进行筛选判断。利用筛选的种子点进行曲面拟合,采用顾及地形起伏的自适应高差阈值剔除非地面点,具体流程如图1所示。

图 1. 算法流程图

Fig. 1. Flow chart of the algorithm

下载图片 查看所有图片

2.1 数据预处理

本文算法对噪声点不敏感,将通过直接剔除异常值的方式进行预处理。

2.2 建立格网索引

首先对点云数据进行格网化,用规则格网将点云数据全部划分,每个格网内有一定数量的点,并建立格网索引。在图2(a)中,黑点表示散乱分布的数据点,虚线框表示每个一级格网,格网的大小与地形相关,取地形内的最大建筑面积^[14],并将每个一级格网分成3×3个二级格网如图2(b)所示。

I=ceilingxmax-xminl,(1)J=ceilingymax-yminl,(2)

式中:I,J为数据点的最大行列号;l为格网步长(取最大建筑边长);ceiling(·)为向上取整函数;x_max,x_min,y_max,y_min分别为点云数据坐标的最大值与最小值。

图 2. 格网索引示意图。(a)一级格网;(b)二级格网

Fig. 2. Grid index diagram. (a) First-level grid; (b) second-level grid

下载图片 查看所有图片

2.3 种子点筛选与判断

本文使用KD树^[15]进行海量三维点云的邻域搜索。KD树是一个二叉树,其索引策略是通过不断地将某一空间递归划分成两个子空间来实现搜索,如图3所示。根节点包含了整个空间,使用二维KD树将整个空间视为一个矩形,通过对矩形的划分实现了点云搜索。

图 3. KD树各节点与平面子区域的对应关系

Fig. 3. Correspondence between the nodes of the KD tree and the plan sub region

下载图片 查看所有图片

图 4. 种子点筛选判断示意图。(a)种子点分布示意图;(b)坡度筛选二级种子点

Fig. 4. Screening and judgement diagram of seed point. (a) Seed point distribution; (b) slope screening for second-level seed point

下载图片 查看所有图片

选取一级格网内最低点作为种子点(p点),通过KD树求以p点为圆心的最邻近k(本文k=7)个点的空间索引,通过最小二乘法求这k个点的特征曲面。求出p点到特征曲面的距离d,若d<d_average(d_average为所有拟合点到拟合曲面的平均距离),则p点为无粗差种子点;否则取一级格网中第二低点为种子点,并使用上述方法判断种子点是否为无粗差种子点,重复操作直到选出无粗差种子点。其原理是在更小的局部面元内,所有的地面点都满足局部地形特征,而粗差点明显不满足局部地形特征。

一级种子点所在格网必定是9个二级格网中的一个,选取其他8个二级格网的最低点作为二级种子点,通过判断二级种子点和一级种子点的坡度来确定二级种子点是否为粗差点,如图4所示。图4(a)中的空心点是一级种子点,其余8个黑点则是二级种子点。分别计算黑点和空心点之间的坡度θ,如图4(b)所示。若θ<θ_t(θ_t为坡度阈值)则视为种子点,否则视为非种子点。由于坡度与高程和位置有关,对于二级格网最低点不满足条件的情况,继续判断该二级格网内第二低点与第三低点是否满足条件。如果三个点都不满足条件,则视该二级格网无种子点。

θ=arctanz2-z1(x2-x1)2+(y2-y1)2,(3)

式中:z₁,z₂分别为p₁,p₂的高程值;(x₁,y₁),(x₂,y₂)分别为p₁,p₂的坐标值。坡度阈值θ_t则根据具体地形而定。参考不同地形的坡度,阈值一般为30°~45°。

2.4 种子点生长

若经过判断、筛选后的一级格网内种子点数少于6个,则无法进行曲面拟合。本文以一级种子点p为参考点,采用基于法向量和距离特征的表面生长算法进行种子点生长^[16-17]。首先求出点p的法向量,然后求p点及其邻近点的特征平面。某一点的法向量定义为以该点及其一定数量的n个最邻近点确定的平面法向量,如图5所示。如果其他点到特征平面的距离s、其他点与种子点的法向量夹角α均在给定阈值内,则将该点视为生长的种子点。

图 5. 点的法向量与特征平面

Fig. 5. Points normals and feature plane

下载图片 查看所有图片

参考文献[ 16-17],本文根据多次实验,设置距离阈值s_t=0.1,夹角阈值α_t=10°,n=6。

2.5 曲面拟合

假设地表是一个复杂空间曲面,则其局部面元可用二次曲面表示,即

f(x,y)=a0+a1x2+a2y2+a3xy+a4x+a5y。(4)

通过求解拟合曲面参数a₀,a₁,a₂,a₃,a₄,a₅就可得到曲面方程。利用筛选、生长的种子点对一级格网进行曲面拟合,并利用每个一级格网的拟合曲面计算一级格网内所有点的拟合高程值,求出拟合高程值和真实高程值的高差,并判断高差与阈值H_t的大小。如果差值小于H_t则将该点视为地面点,否则视为非地面点。由于每个一级格网的地形起伏不一样,采用传统的固定经验阈值时滤波效果不佳,因此使用顾及地形起伏的自适应高差阈值进行滤波。

Ht=(Hmax-Hmin)×R,(5)

式中:H_max为一级格网内的最大高程值;H_min为一级格网内的最小高程值;R根据多次实验,取0.1^[13]。R过大、过小分别会导致滤波不足与过度滤波的问题。

3 分析与讨论

3.1 实验数据

为了对本次实验进行精确分析,使用ISPRS (国际摄影测量与遥感协会)公布的国际滤波标准雷达点云数据集。数据集是由Optech ALTM系统采集,包含8个测区,分为4个城镇测区和4个乡村测区,含有15个样本。其样本属性和初始参数见表1。

3.2 实验分析

ISPRS于2003年提出了点云滤波误差标准。总共有三类误差:Ⅰ类误差是将地面点分为非地面点误差;Ⅱ类误差是将非地面点分为地面点误差;Ⅲ类误差是分类错误点占总点数的比例,称总误差。假设a为滤波后正确分类的地面点点数,b为滤波后地面点被误分为非地面点点数,e为所有真实地面点点数,c为滤波后非地面点被误分为地面点点数,d为滤波后正确分类的非地面点点数,f为所有真实非地面点点数,g为滤波后所有地面点点数,h为滤波后所有非地面点点数,n为样本点总数,则:Ⅰ类误差为: be×100%;Ⅱ类误差为: cf×100%;总误差为: b+cn×100%。

图6为Sample 54滤波后的分布情况,圆点代表正确分类点,三角形代表Ⅰ类误差点,叉代表Ⅱ类误差点。其中误差点大多分布在边界、高程突变处,具有不连续性。其原因是拟合曲面对于此处真实地形的相似度不够,此处的少数个别点不适用于阈值条件,存在误分类的情况。但是整体上滤波效果较好,Sample 54的Ⅰ类误差、Ⅱ类误差、总误差分别为4.16%、5.92%、5.11%。

图 6. 滤波后点云分布图

Fig. 6. View of points distribution after filtering

下载图片 查看所有图片

表2给出本文算法与经典移动曲面算法滤波结果的3类误差统计表。图7为两种算法滤波结果的3类误差的对比图。

通过表2可以看出,两种算法都能取得较好的滤波效果,特别是对于规则建筑物区域,这是因为高差阈值对于规则建筑的判别能力强。经典算法Ⅰ类误差、Ⅱ类误差、总误差的平均值都低于17.62%,优化算法的3类误差的平均值更低,低于12.95%。但是两种算法对于地形起伏较大且不连续的区域滤波效果较差,例如Sample 41是高密度的不连续区域,经典移动曲面算法对该样本的Ⅰ类误差高达62.35%。其原因是经典移动曲面算法直接选取最低6个点作为种子点进行曲面拟合,种子点中存在粗差较大的点且地形不连续,使得拟合曲面与真实地形的差异较大。本文通过种子点的优化策略能有效地降低Ⅰ类误差,虽然Ⅱ类误差增加了2.69%,但是总误差降低了24.2%,滤波效果提升明显。

图7表明本文算法能有效减少Ⅰ类误差,但是对于Ⅱ类误差的精度提高不大。同时本文需要设置更多的阈值参数,计算复杂。但是对于简单样本, Sample 21、Sample 42、Sample 51、Sample 61、Sample 71的总误差提升不明显,分别只减少了1.56%、0.26%、1.97%、1.48%、1.81。但对15个样本的Ⅰ类误差、Ⅱ类误差、总误差的平均值分别降低了7.30%、4.67%、5.57%,有一定的改进效果。

本文将总误差与ISPRS公布的8种算法进行比较,以进行算法验证^[18],如图8所示。根据总误差的平均值,本文算法可以取得较好的滤波效果,但是整体上本文算法的精度略低于Axelsson这一经典不规则三角格网渐进加密滤波算法。但是不规则三角格网渐进加密滤波算法的计算量大,需要反复更新三角格网并计算高程和角度,还需判断镜像点,计算速度比较慢。Pfeifer的分层插值滤波算法效果整体上与本文相似,且对不同地形的适应性强,所有样本总误差均没有超过20%,但本文对关键的高差阈值采用自适应阈值,进一步提高了算法的适应性。其他算法滤波效果较差,适应性差,不同样本的总误差起伏较大。

图 7. 两种算法的3类误差对比。(a) Ⅰ类误差;(b) Ⅱ类误差;(c)总误差

Fig. 7. Comparison of the three types of errors for the two algorithms. (a) Error Ⅰ; (b) error Ⅱ; (c) total error

下载图片 查看所有图片

图 8. 与其他8种滤波算法总误差对比

Fig. 8. Comparison with total errors of other eight filtering algorithms

下载图片 查看所有图片

4 结论

本文在经典移动曲面的基础上提出一种多级种子点的优化策略,基于筛选的种子点拟合曲面,获得的曲面更接近于真实地形。同时采用顾及地形起伏的自适应高差阈值,能有效地判别地面点与非地面,在一定程度上提高了算法的自适应性。通过与ISPRS公布的算法和经典移动曲面算法进行误差对比分析,证明本文算法精度较高,能在连续的地形处取得较好的滤波效果,但是用来优化种子点的坡度阈值仍是根据地形设置的经验阈值。