1 引言

激光雷达(LiDAR)是近年来发展迅猛的新型遥感技术。LiDAR采用主动扫描的方式获取数据,与传统摄影测量的被动扫描手段相比,具有受天气影响小、可持续扫描的优点。LiDAR具有一定的穿透性,在植被鉴别、建筑物提取方面有着巨大的优势。因此,激光雷达获得越来越多专业人士的认可。

由于地形复杂多样,雷达点云的滤波处理一直是难点。当前滤波算法主要分为4类:1) 数学形态学滤波方法;2) 逐渐加密滤波算法;3) 基于表面滤波方法;4)基于分割滤波方法^[1]。Lindenberger^[2]将数学形态学应用于雷达点云滤波;文献[ 3-5]中尝试了自适应滤波算法,但是仍然存在误差较大的问题;文献[ 6]中提出了区域预测点云滤波算法,通过确定斜率阈值,提高了形态学滤波精度。

相比其他滤波算法,数学形态学滤波算法具有实现效率高、原理简单的优点。准确的阈值参数可大幅提高滤波精度,但形态学算法的难点同样在于阈值参数的选取。

本文针对阈值参数难以确定的问题,提出等值线-形态学滤波算法,利用等值线的闭合性和相邻等值线高度差,选取出形态学滤波算法的关键阈值。

2 等值线法快速获取滤波阈值

2.1 渐进型形态学滤波算法

机载雷达在飞行和扫描中会产生粗差,激光雷达获取的是无序且缺少拓扑关系的离散点云数据,因此利用经验模态分解(EMD)算法去除点云噪声是一个必要的过程^[7-8]。本文采用形态学滤波算法进行点云滤波,需将点云数据划分为规则格网。文献[ 9]中提出虚拟网格法形态学滤波算法,该算法对于格网选择、形态学运算(腐蚀、膨胀、开运算、闭运算)、滤波窗口等都有明确表述,本文不再赘述。

采用一维滤波算子,对滤波窗口内的数据点进行处理。梯度定义为

dh=f-f·g,(1)

式中:d_h为梯度;f为规则化的表面模型;g为结构元素。

梯度为数据点高程减掉数据点开运算后的高程。为了保证每一次滤波的准确性,引入梯度作为检验条件。人工定义高差阈值d_T,取梯度d_h与d_T进行对比,满足条件即为地面点,否则为地物点。渐进滤波算法的实现过程如图1所示。

图 1. 渐进滤波算法实现过程

Fig. 1. Progressive filtering algorithm implementation process

下载图片 查看所有图片

对地形分别采用l₁滤波窗口(l₁窗口为滤波中某一次的滤波窗口大小,l₁遍历完成后窗口增大为l₂,重复l₁的遍历工作)计算窗口内每个点的梯度。在l₁窗口中,P为窗口中的任意一个数据点,如果d_hp,1<d_T,1,则P点为地面点,否则为非地面点。d_hp,1为P点在l₁的梯度,d_T,1为l₁窗口中的高差阈值。如果某窗口内梯度最大值小于阈值,则窗口内全部为地面点。

为解决固定窗口无法很好地保留地形特征的问题,Zhang等^[10]提出渐进型形态学滤波算法,将滤波窗口从小依照一定规律逐步增大,每一个窗口都会对数据进行一次遍历,逐渐增大的窗口使数据点经过数次遍历判断。本文采用线性变化滤波窗口,窗口变化公式为

wk=2kb+1,(2)

式中:w_k为第k次迭代的窗口大小;k为迭代次数,当k最大时,w_k取得最大滤波窗口w_max;b为初始窗口大小。第k次高差阈值的设置为

dT,k=dh0,wk≤2dT,k=C(wk-wk-1)S+dh0,wk>2dT,k≥dhmax,(3)

式中:d_h0为初始高差;C为格网大小;S为斜率;d_hmax为人工定义的最大高差阈值,当逐渐变化的高差阈值超过人为规定的最大阈值时,阈值就会停止增加。为增加高差阈值与窗口阈值的相关性,通常斜率的表达式为

S=dhT,k0.5(wk-wk-1)。(4)

2.2 形态学滤波算法阈值的选取

由渐进性形态学滤波算法可知,该算法需要人工选取阈值,包括最小高差d_h0、最大高差d_hmax、最小窗口大小b、最大窗口大小w_max等参数。准确的阈值可以快速提高滤波算法的准确率,本文采用等高线生成方法判断测区的阈值。

2.2.1 等值线生成流程

文献[ 11]中提出规则格网和不规则三角形等值线生成算法,本文利用区域内离散无序的高程数据点,通过内插方法得到规则格网点的高程值,利用格网点的横坐标值、纵坐标值和高程值生成等值线,等值线生成流程如图2所示。

图 2. 等值线生成流程图

Fig. 2. Contour generation flow chart

下载图片 查看所有图片

2.2.2 等值线法选择最大滤波窗口阈值

雷达点云中的建筑物由于高差和面积较大,一直是滤波算法的难点。建筑物具有明显的特征,建筑边缘大多数都为规则图形,而且建筑边缘的高程值非常接近。基于这一特性,将高程值近似点绘制成等值线。由于建筑物的等值线是闭合的,因此会生成建筑区域内的等值线,等值线所确定的范围即为建筑物的大小。形态学滤波最大窗口阈值需要满足最大建筑物的大小。某测区航拍实景图与等值线示意图如图3~4所示。

图 3. 测区航拍体育馆影像图

Fig. 3. Aerial photography of the survey area

下载图片 查看所有图片

图 4. 等值线确定测区体育馆形状

Fig. 4. Contour to determine the shape of the stadium

下载图片 查看所有图片

本文采用的是一维线性滤波窗口,而闭合等值线常用于确定二维闭合形状。为使一维窗口能够准确确定建筑物大小,本文提出改进方案:利用等值线内所有数据点,计算等值线区域内所有点之间的最大距离,最大距离判定为最大滤波窗口大小。最大距离既保留了建筑物的特性,又能避免二维滤波窗口带来的计算量过大的问题。

Dmax=max(Di,j)=max(xi-xj)2+(yi-yj)2,(5)

式中:D_max为区域内两点之间的最大距离;D_i,j为区域内任意两点间的距离;x_i为i点横坐标;x_j为j点横坐标;y_i为i点纵坐标;y_j为j点纵坐标;i、j为闭合区间任意两点。

激光扫过的高程变化明显的位置既可能被判断为建筑物,也可能被判断为陡坎地形。但建筑物的等高线是一个封闭的图形,而陡坎则通常不是。因此,使用等高线关联有其优势。等高线方法可能会带来阈值误差,但其阈值通常都是偏小的。这是因为:由激光扫描点所构造的等高线与建筑物的外轮廓线并不完全一致,而且始终偏小。为解决该阈值误差,可以增加一个小的常量。事实上,该阈值误差可以忽略,因为建筑物未扫描到的轮廓边界点云是不需要滤波的,而扫描到的建筑物点云则会出现在最大滤波窗口范围内。

2.2.3 等值线法选择最大高差阈值

最大高差是滤波算法中的另一重要阈值。建筑物不同于缓慢变化的山体,建筑梯度变化的陡变性是识别建筑高程的重要依据。等值线的疏密代表地形的平缓与陡峭。建筑底部的激光脚点与顶部的激光脚点虽然会存在横纵坐标的近似,但是高差阈值变化很大。徐旭^[12]在等值线填充的基础上,利用相邻等值线色彩梯度差异,可以快速判断出建筑物的最大高差。建筑区点云等值线填充分布及其雷达点云分布示意如图5~6所示。

图 5. 建筑区点云等值线填充分布

Fig. 5. Point cloud contour filling distribution in the building area

下载图片 查看所有图片

图 6. 建筑区雷达点云分布

Fig. 6. LiDAR point cloud distribution in the building area

下载图片 查看所有图片

2.2.4 坡度角度对于最大高差的检验

对于地面和建筑、地面和大型植被的高差值近似的临界点,高差不足以判断点的属性,因此引入倾斜角作为判断的依据。

tanθ=ΔhΔl,(6)

式中:θ为倾斜角;Δh为相邻点的高差;Δl为相邻点的水平距离。建筑临界点具有陡变性,坡度近似垂直;而高大植被的坡度相对较小。根据倾斜角判断相同高差的地物如图7所示。

2.3 形态学滤波阈值参数的确定

算法步骤如下:

1) 雷达点云去噪。

2) 点云通过虚拟格网法建立规则格网。

3) 对点云数据建立等值线,根据建筑物等值线的闭合性和规律性判断最大建筑的大小,定义最大滤波窗口。

图 7. 根据倾斜角判断相同高差地物

Fig. 7. Tilt angle to determine the same height difference

下载图片 查看所有图片

4) 对等值线进行填充,依据建筑物高差梯度变化和倾斜角定义窗口内的最大高差阈值,确定初始滤波窗口大小。

5) 对滤波窗口的数据进行开运算,判断梯度是否小于阈值:如果满足条件,为地面点;如果不满足条件,为非地面点。

6) 利用(2)式增大滤波窗口,利用(3)式定义不同高差阈值。

7) 滤波窗口最大化,完成滤波。

等值线-形态学滤波算法的流程图如图8所示。

图 8. 等值线-形态学滤波算法流程图

Fig. 8. Isometric line-morphological filtering algorithm flow chart

下载图片 查看所有图片

3 实验结果与分析

雷达点云滤波存在3种误差:1) I类误差,将地面点误分为非地面点;2) II类误差,将非地面点误分为地面点^[13];3) 总误差,两类误差总数与数据点个数的比。实验数据选取国际摄影测量与遥感学会(ISPRS)提供的数据集,用于定量、定性分析滤波精度,通过3组实验进行验证分析。

3.1 实验一:验证形态学滤波精度

选取ISPRS中Samp23数据集用于实验。Samp23有标准的属性,0表示地面点,1表示地物点。等值线法确定测区最大滤波窗口大小和测区雷达点云分布,如图9~10所示。

图 9. 等值线确定窗口阈值

Fig. 9. Contour determination window threshold

下载图片 查看所有图片

图 10. 测区点云分布

Fig. 10. Distribution of point cloud

下载图片 查看所有图片

采用MATLAB软件实现改进型形态学滤波算法,利用等值线生成的最大滤波窗口为30 m,最大高差为15 m,判断地面点与地物点,并统计I、II两类误差。实验结果如表1、图11所示。

通过对表1的精度分析,判断形态学滤波可以较为准确地得到地面点,等值线生成的阈值真实有效。比较滤波前后点云分布,发现滤波算法可以较好地过滤地物,获得地面点。

图 11. Samp23测区滤波前、后的点云分布。(a)滤波前;(b)滤波后

Fig. 11. Samp23 survey area before and after filtering. (a) Before filtering; (b) after filtering

下载图片 查看所有图片

3.2 实验二:验证等值线法确定阈值的滤波精度

为了判断最大窗口阈值和最大高差阈值分别对于滤波结果的影响,引入总误差作为衡量滤波精确度的标准。选取Samp42作为实验样本,样本地形为一大型建筑区。等值线法确定的最佳高差阈值为20 m,最大窗口阈值为30 m。实验二包括窗口阈值不变条件下,高差变化对于滤波结果的影响;和高差阈值不变条件下,窗口变化对滤波结果的影响。实验结果如图12~13所示。

图 12. 不同高差阈值的滤波精度(窗口阈值为30 m)

Fig. 12. Filtering accuracy of different height difference thresholds (window threshold is 30 m)

下载图片 查看所有图片

图 13. 不同窗口阈值的滤波精度(高差阈值为20 m)

Fig. 13. Filtering accuracy of different window thresholds (height difference threshold is 20 m)

下载图片 查看所有图片

分析图12~13可知,等值线方法确定阈值速度快且更加准确,减少了多次阈值选取造成的人力物力损失,有针对性地提取测区范围内的阈值,进一步提高了滤波精度,而过大的阈值误差则会直接影响整个测区的滤波精度。

3.3 实验三:验证等值线-形态学滤波算法的精度

同时采用三角网(TIN)渐进加密算法^[14]对该样本数据进行定性定量验证,对比不同算法的滤波效果。定量实验选取Samp41中测区样本,采用形态学和不规则TIN两种方式分别计算。TIN选取区域内最大尺度阈值为30 m,点到平面的阈值为1 m,倾斜角阈值为6°;形态学滤波最大高差阈值取60 m,最大窗口取20 m。实验结果如表2所示。

表2通过对两种误差数据的定量分析发现:等值线-形态学滤波算法的I、II类误差明显小于TIN滤波算法,说明等值线-形态学滤波算法在这一测区的滤波精度优于TIN滤波算法。

定性实验数据选取FSite8_red2中数据点作为数据集。根据等值线-形态学滤波算法确定最大窗口为20 m,最大高差为30 m;TIN滤波算法设置研究区域最大格网阈值30 m,点到平面距离阈值1 m,最大倾斜角阈值为6°。TIN算法与形态学算法的实验对比结果如图14所示。

图 14. 不同算法在FSite8_red2测区内的滤波结果。(a) TIN算法;(b)形态学算法

Fig. 14. Filtering results of different algorithms in the FSite8_red2 survey area. (a) TIN algorithm; (b) morphology algorithm

下载图片 查看所有图片

从图14可以看出滤波后的地面点分布。相较于TIN滤波算法,形态学算法较好地保留了地貌特征,更精确地还原地貌,而TIN算法由于逐级三角形迭代加密,对于地型易产生I类误差,将树木等较小地物纳入地面点,从而造成较大误差。

4 结论

等值线-形态学滤波算法相较于传统形态学滤波算法,在阈值选取上有更强的目的性,可以更快捷地确定准确的阈值参数,对于不同测区可以选取恰当的滤波阈值,滤波结果的精度更高。但是本文方法仍然没有解决点云滤波的自适应问题,目前的等高线阈值算法确定部分还需要借助人工设置,阈值的自动生成仍有待于进一步研究。