改进的渐进加密三角网机载点云滤波方法

王效盖 王 健 刘翔宇 曹 一

(山东科技大学 测绘与空间信息学院, 山东 青岛 266590)

0 引言

机载激光雷达(light detection and ranging,LiDAR)技术是一种高效、快速、大区域采集三维空间点云的主动遥感技术,已广泛应用于数字模型构建、实景三维建设、森林资源管理、地质灾害监测等领域[1]。其中,机载点云滤波是分离地面点和非地面点的关键步骤,也是LiDAR点云数据后续实际应用的基础。目前地面滤波算法主要分为坡度滤波、形态学滤波、分割滤波和插值滤波。

坡度滤波是利用点云之间坡度变化分离地面点的一种滤波方法,这类方法十分依赖坡度阈值的设定,其对连续地形有较好的分割结果,但在陡坡、树木与地面连接处等坡度突变区域容易出现错误滤波现象[2-3]。形态学滤波主要原理是LiDAR点云在形态学开闭运算前后存在高程差异,高程变化小的为地面点,反之为非地面点,其精度依赖窗口大小的设置,需要使用者设定合理的窗口参数,对非专业人员不友好[4]。分割滤波首先根据LiDAR点云某种属性对其进行分割操作,并根据分割属性之间的差异分离地面点和非地面点,在特殊地形中有较好的分割结果,但其精度受分割结果影响较大[5-6]。插值滤波首先需要选择地面种子点,并使用插值方法生成初始参考面,若待分类点与参考面之间满足一定的关系则将其划归为地面点,基于插值的滤波方法原理简单、滤波结果可靠、鲁棒性强,但是其结果容易受种子点构造初始参考面精度的影响[7-8]。Axelsson提出的渐进加密三角网(triangulated irregular network,TIN)滤波方法(progressive TIN densification,PTD)是插值滤波中的经典算法,因其较高的滤波精度和较快的滤波速度而被集成于多个开源或商业软件中,如TerraSolid等[9]。实验发现,PTD算法在对平坦或地形简单的区域有较好的滤波效果,但在复杂林区中因获取初始地面种子点数量变少且可靠性降低,导致在进行迭代运算时误差累积增大,使其滤波效果不太理想。

为此,在PTD算法的基础上,本文提出了一种适用于林区LiDAR点云改进的渐进加密三角网滤波方法(improved progressive TIN densification,IPTD)。IPTD算法在地面种子点选取策略上做出一定改进,借助布料模拟算法和局部薄板样条插值获得大量分布均匀且可靠的地面种子点,用以提高初始TIN的精度,然后对待分类点进行迭代滤波,从而提高算法在林区环境中LiDAR点云的滤波精度。

1 渐进加密三角网滤波

渐进加密三角网滤波算法是Axelsson在2000年提出的一种经典插值滤波方法并且得到广泛的使用[9]。算法流程中选取局部窗口的最低点作为初始地面种子点,以此构建初始三角网,通过计算待定点P到ABC平面的距离d及连线之间的夹角α、β、γ,判断距离和角度是否满足设定阈值,若满足则将P点加入地面点并构建新的三角网,反复迭代至没有新点加入,滤波完成。渐进加密三角网滤波算法原理如图1所示。

图1 渐进加密三角网滤波原理

该算法具体步骤为：

(1)选择地面种子点。对研究区域进行规则格网划分,格网大小设置为区域最大建筑物大小,选取格网中最低点作为初始地面种子点。

(2)构建初始三角网。根据(1)中选择的地面种子点构建一个稀疏的初始三角网。

(3)待定点分类。依次遍历所有待定点,分别计算待定点到三角网平面的距离和夹角,若距离和角度均满足设定阈值,将该点划分为地面点。

(4)更新三角网。使用新加入的地面点对三角网进行更新,生成新的三角网。

(5)重复步骤(3)和(4),迭代至没有新增的地面点,算法结束。

该算法原理简单,效果稳定,对于城市和林区都有较强的适用性。其算法精度主要取决于初始三角网构建的准确性和迭代阈值的设定。在植被密集覆盖且地形崎岖复杂的林区,采用局部最小值获取的种子点可靠性降低,从而影响初始三角网的准确性,最终使得滤波精度较低滤波效果不理想。

2 改进的渐进加密三角网滤波方法

2.1 噪声剔除

原始机载LiDAR数据中噪声来源于扫描仪本身或受环境影响产生的多路径效应。这对使用布料模拟算法获取潜在地面种子点时产生严重影响,导致获取的潜在地面种子点结果不可靠。因此,在进行后续处理前需要剔除原始点云中远离主体点云的离群噪声点。基于空间位置分布的去噪算法(statistical outlier removal,SOR)能够稳健的检测这些离群点。SOR滤波器对每个点进行K邻域统计分析,计算K个邻近点平均距离分布的均值和标准差,将平均距离在给定阈值范围之外的点视为离群点并剔除。如式(1)所示,本文使用12个最近邻计算每个点的均值和中值距离。

Max(d)>Mean(d)+Std(d)

(1)

2.2 获取潜在地面种子点

布料模拟滤波算法是由Zhang等提出的一种新型滤波方法,其基本思想：首先将机载LiDAR点云进行倒置并在测量点上方覆盖一块柔性布料,柔性布料在重力的作用下逐渐黏附在测量点上,最终柔性布料模拟出一个近似的地形表面,如果测量点到柔性布料的距离满足阈值要求则标记为地面点[10]。布料模拟的优点在于：①参数设置简单,仅需较少的参数即可完成滤波;②地形适应性好,可用于连续地形、陡坡和断崖;③运算效率高。PTD算法中只有少数点被选为地面点,导致初始三角网精度较低。相比之下,布料模拟算法可以获取大量分布均匀的潜在地面种子点。布料模拟算法过程如图2所示。

图2 布料模拟算法过程

2.3 获取地面种子点

在应用布料模拟算法滤波后,某些非地面点会错误的包含在潜在地面点中,为了提高原始三角网构建的精度,需要识别并剔除这些非地面点,本文使用薄板样条插值方法来检验并消除非地面点。薄板样条插值是一种空间离散点拟合精度高、鲁棒性强的空间插值方法,在对复杂地形的模拟中表现出较好效果[11]。

对于3维空间中的测量点集(X,Y,Z),薄板样条插值曲面可由式(2)表示。

(2)

(3)

其中,N是用于TPS插值的控制点数;p(x,y)是趋势函数,a是它的系数;wi是与控制点相关的权重;q(ri)是径向基函数。径向基函数定义为

q(r)=r2ln(r)

(4)

r为两点间的欧氏距离。矩阵表达式如下：

(5)

实验发现,非地面点存在局部高程突变的现象,本文认为在邻域点云中高程值异常的点为非地面点,并构造插值曲面从潜在种子点中筛选并去除非地面点。因为求解薄板样条插值函数系数的时间和插值点的数量成正比,所以本文采用局部薄板样条插值代替整体薄板样条插值来提高效率。具体步骤为：①使用规则格网(20×20)组织潜在地面点数据,并选择局部最小值作为插值点;②搜索待估计潜在种子点K邻域的插值点,构造局部薄板样条插值曲面,K设置为12,并将激光脚点坐标代入求解插值高程;③计算真实高程与插值高程的残差,若高程残差大于设定的残差阈值,则认为该点为非地面点并从潜在种子点中剔除;④逐点筛选得到地面种子点。

2.4 改进的渐进加密三角网滤波

传统渐进加密三角网滤波算法选择种子点采用格网区域高程最低点,在林区环境下种子点选取较少且可靠性低,在进行迭代运算时造成误差累积。因此,本文在种子点选取策略上做出改进,使用布料模拟算法和局部薄板样条插值获得大量分布均匀和可靠的地面种子点,解决传统渐进加密三角网滤波在林区环境下种子点稀疏的问题,提高林区下的适应性。改进的渐进加密三角网滤波步骤如下：

(1)种子点选取。首先使用布料模拟算法提取机载LiDAR点云中的潜在地面点,然后通过局部薄板样条插值剔除潜在地面点中的非地面点,获得分布均匀和可靠的地面种子点。

(2)初始三角网构建。在机载LiDAR点云四周设置宽度为m的缓冲区,并以固定间隔n构造辅助点,辅助点高程为最近种子点的高程[12]。通过种子点和辅助点构建初始三角网。

(3)三角网滤波。计算待定点到邻近三角形的距离和角度,若其距离和角度小于设定的距离阈值和角度阈值,则分为地面点。

(4)更新三角网。将步骤(3)地面点与已有种子点重新构建三角网。

(5)迭代运算。重复步骤(3)和(4),迭代运算至没有新的地面点加入。

3 实验及结果分析

采用国际摄影测量与遥感学会(International Society for Photogrammetry and Remote Sensing，ISPRS)标准数据和林区数据进行实验并对滤波结果进行精度分析。本文采用的滤波精度评价指标为Ⅰ类误差(T.Ⅰ)、Ⅱ类误差(T.Ⅱ)、总误差(T.E)和Kappa系数,其中三种误差数值越小和Kappa系数数值越大表示滤波精度越高。计算方法如表1所示,其中a和b分别为正确分类和错误分类的地面点数;d和c分别是正确和错误分类的非地面点数;e是a、b、c和d的总和。

表1 地面滤波精度评价方法

3.1 ISPRS标准数据

采用ISPRS提供的标准数据集中山地部分(samp51-samp71)作为测试数据,验证本文所提方法的有效性。每个样本的参考数据都是人工滤波生成的,样本点被分类为地面点和非地面点。IPTD对ISPRS标准数据集的滤波结果如表2所示。

表2 ISPRS标准数据滤波精度

从表2可以看出：IPTD算法滤波结果Ⅰ类误差为1.23%,Ⅱ类误差为12.29%,样本平均Kappa系数和总误差分别为84.96%和2.16%,说明本文算法有效控制了Ⅰ类误差的产生。samp51、samp61和samp71获得的Ⅰ类误差均小于1%,其中samp51的Ⅰ类误差最小,仅为0.09%,说明本文方法在缓坡滤波中能较较好地保留真实地面点,减小其误分类的可能性。Samp53Ⅱ类误差最大,高达34.06%,其余样本数据都控制在一定范围内,本文算法在断崖处的滤波还有待优化。表中samp51的Kappa系数最高,高达96.04%,samp61的总误差最小,仅为0.85%。samp52和samp53的总误差都较大,samp53的Kappa系数最小,仅为56.06%。图3和图4展示了两个样本的数字地表模型(digital surface model，DSM)、滤波后的数字高程模型(digital elevation model，DEM)及滤波误差的空间分布情况,其中蓝色为Ⅰ类误差,红色为Ⅱ类误差。由图3可见,samp52的数据中存在断崖和陡坡,沿河分布低矮的植被,这些高差突变区域对阈值有一定挑战,判断错误而造成Ⅱ类误差较大。由图4可见,samp53的数据中分布有梯田等落差较大区域,其中误差分布较多的地方地形变化较大,导致Ⅱ类误差较大。

图3 Samp52滤波误差空间分布

图4 Samp53滤波误差空间分布

将本文方法与近年来提出的5种方法进行对比[8,13-16],如表3所示。本文方法在6种方法中总误差平均值最小,较原有滤波方法总误差的平均值减小0.75%,且本文方法在所有样本数据中有3个总误差最小,说明本文方法在滤波精度上具有一定的优势。

表3 本文方法与其他5种方法的总误差对比 单位：%

3.2 林区数据

ISPRS标准数据集点云密度较小且植被特征不够明显,为了进一步验证本文方法在森林地区的适用性,本研究从OpenTopography网站(http://www.opentopography.org/)下载3组高密度林区机载LiDAR数据。3组林区数据具有不同地地形特征：数据Site1位于美国加利福尼亚州内华达山脉南部,为坡度较小和植被覆盖密集的缓坡地形;Site2位于美国犹他州中南部山区,为覆盖分散低矮植被的陡坡地形;Site3位于美国科罗拉多州,地形起伏明显,覆盖较多密集的植被。3组林区数据点云如图5所示。

图5 林区点云数据

采用本文方法和PTD算法对3组样本数据进行滤波处理并进行对比分析,其结果如表4所示。根据表中数据可以看出,本文方法在Ⅰ类误差和Ⅱ类误差上都有大幅的降低,说明本文方法在林区滤波中能够大幅度减少错误分类现象,获得较准确的地面模型。从总误差来看,本文方法相较于PTD算法在Site1区域降低了3.51%,在Site2区域降低8.68%,在Site3区域降低4.02%,说明本文方法在总误差控制上有较好的表现。本文方法Kappa系数最大为94.57,最小为80.14,在Kappa系数有了大幅的提高,其中Site2区域提高仅14.50%。由图6、图7和图8可以看出,本文方法生成的DEM与参考DEM接近。综上所述,相较于传统的PTD算法,本文方法可以应用于林区点云滤波中。

表4 林区数据滤波精度对比

图6 Site1参考数据、PTD和本文方法生成的DEM

图7 Site2生成的DEM

图8 Site3生成的DEM

4 结束语

为了提高渐进加密三角网滤波算法在林区机载点云数据中的滤波精度,提出了一种改进的渐进加密三角网滤波方法。所提方法使用布料模拟滤波和局部薄板样条插值方法代替局部最小值法获取大量均匀可靠的地面种子点,解决了因种子点稀疏导致初始三角网精度较低的问题。使用ISPRS标准数据和林区数据对方法进行验证,并与其他滤波方法进行对比分析,本文方法有较低的总误差和较高的Kappa系数,证明本文方法能有效应用于林区点云滤波中。本文方法在陡坡等复杂地形中滤波误差仍然较大,未来就提高不同地形下林区点云的滤波精度继续展开工作。