崔 浩,高 飞,余 敏,叶周润

(合肥工业大学 土木与水利工程学院,安徽 合肥 230009)

在工程建设、道路设计、滑坡监测、土地覆盖分类、森林管理和道路设计中,高分辨率数字地形模型(digital terrain model,DTM)有着至关重要的作用。光探测与测距技术(light detection and ranging,LiDAR)又称激光雷达，是一种有效的收集上空三维点云的方法,并广泛用于生成DTM。为了获取高精度DTM,需要将地面和非地面测量数据进行分离,这是一个滤波过程。目前地面滤波算法主要有基于坡度的方法、基于数学形态的方法和基于地表的方法等。

基于坡度的方法通常假设在一个邻近区域内,地形的坡度变化通常是渐进的,而建筑物或树木与地面之间的坡度变化很大。基于这个假设文献 [1]开发一种基于坡度的滤波算法,利用LiDAR点与其相邻点之间坡度的差值来进行滤波。该方法在复杂地形如陡坡、断裂地形处容易出现滤波错误,具有一定的局限性[2]。

基于数学形态的方法是通过数学形态学来剔除非地面LiDAR点。对于这种滤波方法,选择最佳窗口大小是至关重要的[3-6]。较小尺寸的窗口可以有效地滤除较小的物体,但会保留较大的建筑物在地面点中,而大的窗口尺寸往往会滤除平滑的地形细节,如山峰、山脊和悬崖,故通常需要利用研究区域的额外先验知识来确定合适的窗口大小。

基于地表的方法主要是建立地形模型,然后通过迭代运算近似得到原始地形,其中最经典的就是文献[7]提出的渐进式不规则三角网(triangular irregular network,TIN)算法。基于TIN滤波的可靠性和准确性已经被大量的研究所证实,但其在处理不连续地形时滤波效果较差[8-10]。

针对上述问题,许多研究者提出改进的方法。文献[11]利用数字形态学算法提取种子点;文献[12]利用均方差、点密度和小格网高程统计数据来选取种子点;文献[13]采用邻域卷积和进行种子点的判断与虚拟种子点的构造;文献[14]采用多核并行处理技术加快三角网渐进加密滤波方法的算法效率。

单一的算法往往都具有一定的局限性,而绝大多数的改进算法都是在TIN算法基础上进行的,因此本文提出一种结合布料模拟滤波(cloth simulation filtering,CSF)[15]与TIN的点云滤波算法。CSF算法的主要优点在于:① 可以应用于陡坡、间断区域和狭长地形;② 相比于其他滤波算法要设置复杂的滤波参数,CSF算法参数设置较为简便;③ 其他算法针对不同的地形要对算法参数进行调试,找出最优的地形值,而CSF算法的布料模拟率可以针对不同地形进行选择,自动进行处理[16-17]。但是CSF算法最后的阈值设置仅仅依赖于点云数据与格网数据之间的距离,因此存在一定的局限性。本文将CSF算法与TIN优势互补,通过CSF获取初始种子点,然后利用初始种子点整体构建格网,通过判断该点与三角网之间的角度和距离是否小于一定的阈值来获取地面点,将该地面点与之前的三角网再次整体构网,最终获取完整地面点。

1 本文算法原理

本文算法主要包括点云粗差剔除、CSF粗滤波和改进的TIN精滤波3个步骤。

1.1 点云的粗差提取

为减弱粗差对地面点提取的影响,本文利用CloudCompare软件对机载LiDAR点云进行粗差探测和剔除。

1.2 布料模拟粗滤波

传统的CSF中,往往是通过判断点云数据中的点与格网粒子之间的距离是否小于阈值,从而确定其是否为地面点,但是这种方法通常会导致误判或者漏判,并不能作为判断地面点的准确依据。本文首先采用CSF对粗差剔除后的点云进行滤波,滤除建筑物、构筑物等较大的地物,然后进行种子点的选取,将其作为三角网滤波的初始种子点,解决TIN初始种子点选取的难题,提高算法的精度和自适应性。

CSF是基于表面的点云滤波算法,主要是对于物理过程的模拟。假设一块足够柔软的布被放在地面上,这块布因重力作用而下落,可以黏着在物体表面上,地形表面空洞的点最终会由于黏着在空洞周围的非空洞点而拉回到地面,最终,布料表面的形状就是数字地表模型(digital surface model,DSM)。但是,如果将地形上下翻转,并且该布料具有一定的硬度,那么最终布料的形状就是DTM,以此为基础将原始点划分为地面和非地面部分。

本文中CSF的主要步骤如下:

(1) 将机载点云高程数据上下颠倒,如图1所示。

图1 点云数据倒置

(2) 设置布料网格分辨率,从而得到格网粒子,但是格网中所有点的高程坐标在反转点云的最高点之上。

(3) 把所有点云数据与格网数据投影到一个水平面,并为每个粒子找到其最邻近点(corresponding point,CP),记录其投影前的高程(intersection height value,IHV)。

(4) 对于所有可移动的格网粒子,粒子的位置移动取决于外部和内部2个方面的驱动因素,表达式为:

(1)

其中:m为粒子的质量,其值通常设置为1;X为时刻t粒子的位置;Fext(X,t)为格网粒子受到的外力(重力或者遇到障碍物的碰撞力);Fint(X,t)为粒子在位置X和时刻t的内力(相互连接的力)。内部力和外部力随着时间t变化而变化。

粒子受到外部驱动因素影响而产生的位移表达式为:

(2)

其中:Δt为时间步长;G为常数;X(t+Δt)为下一步点的位置;X(t)为当前点位置;X(t-Δt)为上一步点的位置。若时间步长和初始位置都确定,则可以计算粒子的当前位置。

将粒子当前位置与该格网粒子对应CP点的IHV进行对比,若格网粒子高程等于或者低于IHV,则该粒子的高程被设置为IHV并且为不可移点。

(5) 对于格网粒子,计算每个粒子受到内部驱动因素所产生的位移,若2个相邻点是可移动点,则其中一点向上移动,另一点向下移动;若2个相邻点中的一点是不可移动点,另一点是可移动点,则可移动点的移动距离与布料的硬度相关,布料硬度越大,移动的距离越小,如图2所示。

图2 CSF内力作用

不可移动点与可移动点之间的距离为d,当布料硬度为1时,可移动点的距离上移d/2;当布料硬度为2时,可移动点的距离上移3d/4;当布料硬度为3时,可移动点的距离上移7d/8。

位移量公式如下:

(3)

其中:di为粒子的位移量;当粒子为可移动时,b=1,不可移动时b=0;P0为待定的移动点;Pi为P0的相邻粒子;n为点在垂直方向上的单位向量[0 0 1]T。

(6) 重复进行步骤(4)、步骤(5)计算,当迭代次数达到设置的最大迭代次数或者所有粒子的移动距离都小于阈值时,停止迭代。

(7) 辨别地面点与非地面点的方法为:若LiDAR点云与格网粒子之间的距离小于预先设置的阈值,则认为其是地面点;反之则为地物点。阈值可以选择一个较小的值,从而为以后的TIN滤波提供一定数量的种子点。

1.3 改进的TIN精滤波

传统的TIN滤波算法中,初始种子点的选取主要依靠最大建筑物的尺寸及其范围,通常会导致种子点的密度较小,初始种子点构成的TIN不能充分表示地貌,因此,本文采用CSF获取格网粒子作为构建初始TIN的种子点,解决传统TIN点云滤波算法中初始种子点密度过小的问题,增加算法的鲁棒性和自适应性。

本文改进的TIN精滤波主要步骤如下:

(1) 首先对CSF后的点云进行分析,设置网格分辨率,对其进行格网划分,并且选取每个格网中的最低点为初始种子点。

(2) 利用点云整体四周的边界点为辅助点,边界点的高程一般为最近的种子点的高程,然后用种子点与辅助点构造初始的Delaunay三角网。

(3) 遍历所有三角形和点,确定每个点在哪个三角形的投影中,并且该三角形与该点构成一个四面体,如图3所示。

图3 TIN滤波原理示意图

其中点P为三角形上方的任意离散点,点P到点A的连线与三角形平面构成∠a,点P到点B的连线与三角形平面构成∠b, 点P到点C的连线与三角形平面构成∠c,点P到三角形平面的距离为d′,如果d′、∠a、∠b、∠c都小于某一设定的阈值,那么该点为地面点,反之为地物点,角度阈值一般设置为6°～50°,每次迭代过程中可以适当减小阈值,最终阈值一般为2°～10°。

(4) 取滤波后的地面点和种子点重新构建TIN。

(5) 重复步骤(3)、步骤(4),进行迭代,直到没有新的地面点加入或达到最大迭代次数。

2 实验结果和分析

为了验证该算法的有效性,选取国际摄影测量与遥感学会(International Society for Photogrammetry and Remote Sensing,ISPRS)网站(https://www.isprs.org/)的3组滤波测试样本数据进行实验,原始点云数据如图4所示。第1组为Samp11点云数据(图4a),覆盖的区域地形起伏较大,包含较多复杂建筑物,还有路灯、汽车等小型移动地物,左上方有一个山坡,整个区域包括37 922个离散点,区域高差为107.31 m;第2组为Samp31点云数据(图4b),覆盖区域为城市区域,地势较为平缓,包含大量建筑物,区域内包含28 836个离散点,区域高差为35.49 m;第3组为Samp54点云数据(图4c),地势较为平缓,包含大量植被及房屋,区域内包含8 596个离散点,区域高差为32.99 m。

图4 3组原始点云数据

本文算法采用点云库(Point Cloud Library,PCL)C++语言编程对上述3组点云数据进行滤波实验,并且将滤波后的点云数据与ISPRS官网生成的点云数据(参考数据)进行对比,如图5～图7所示。从图5～图7可以看出,本文算法对各种场景的点云数据有较好的滤波效果,可较好地滤除建筑物、构筑物和树木等。

图5 Samp11滤波结果

图6 Samp31滤波结果

图7 Samp54滤波结果

根据ISPRS规定的Ⅰ类误差、Ⅱ类误差评价标准对实验结果进行精度评定,为了便于对比分析,同时采用CSF算法、传统TIN算法对3组实验数据进行滤波实验并进行精度评价,Ⅰ类误差、Ⅱ类误差、总误差、E1和E2结果见表1所列。表1中:Ⅰ类误差为将原本是地面点的点划分入地物点集;Ⅱ类误差为将非地面点错误地划分入地面点集;总误差为Ⅰ类误差和Ⅱ类误差占总点数的比例;E1为本文算法与CSF算法相比的误差变化量;E2为本文算法与TIN算法相比的误差变化量。

从Ⅰ类误差和Ⅱ类误差看,本文算法在Samp11中的Ⅰ类误差较低,说明该算法对于坡度较大的地区,地物点的精度较高,而在Samp31和Samp54中,本文算法Ⅰ类误差与CSF算法相比偏高，但与TIN算法相比较小,说明其在地势坡度较低区域的地物点精度较低,因此本文算法在平坦地势上低矮地物的激光脚点存在欠滤波的现象,对于此类地物该算法的滤波效果有待提高。

与传统TIN算法相比,本文Ⅰ类误差降低幅度最大,在坡度较大的区域Samp11中降低50%以上,在建筑物密集区域Samp31、Samp54中降低幅度较小;与CSF算法相比,本文算法在Samp11中Ⅰ类误差降低7.35%。

对Ⅱ类误差分析可知,相比于CSF算法,本文算法在降低Ⅱ类误差方面具有一定的优势,在Samp31中降低1.45%;相比于TIN算法,本文算法Ⅱ类误差值增大,滤波精度降低。

从总误差看,本文算法在Samp11中较传统TIN算法最大降低26.51%,说明该算法与传统TIN算法相比具有更高的精度;与CSF算法相比,本文算法总误差在Samp11中增加0.47%,在Samp31中增加1.68%,在Samp54中增加5.68%。

综上所述,本文算法在坡度较大的区域可以显著地降低Ⅰ类误差,并将Ⅱ类误差控制在一定的范围内,验证了该算法具有较好的普适性与有效性,能够获得更接近真实地形的高精度数字高程模型(digital elevation model,DEM)。

表1 3种算法滤波结果误差对比 单位:%

3 结 论

本文提出一种CSF 与TIN相结合的地面点滤波算法,该算法在坡度较大的区域中可以提高地面点提取精度。与传统TIN滤波方法相比,本文算法在CSF的结果上适当获取种子点,不仅可以避免地形因素的干扰,而且在种子点选取上可以不受建筑物尺寸影响,解决了种子点选取敏感问题和阈值选取问题。与CSF算法相比,本文算法可以在较大坡度的区域提高地面点提取精度,但是在坡度较低的区域内,提取点云点的精度相对降低,但可以保持在一定的范围内。本文算法需针对坡度较低的区域进行研究,进而保证其在坡度较低区域中地面点提取精度。