基于点云空洞修复和TPS模型的数学形态学机载LIDAR滤波

胡永杰，姚晓伟，胡倩伟,韩凤纳

(1.新疆兵团勘测设计院(集团)有限责任公司，新疆 乌鲁木齐 830002；2.东华理工大学 测绘工程学院，江西 南昌 330013；3.北京四维远见信息技术有限公司，北京 100070)



基于点云空洞修复和TPS模型的数学形态学机载LIDAR滤波

胡永杰1,2，姚晓伟1，胡倩伟3,韩凤纳1

(1.新疆兵团勘测设计院(集团)有限责任公司，新疆 乌鲁木齐 830002；2.东华理工大学 测绘工程学院，江西 南昌 330013；3.北京四维远见信息技术有限公司，北京 100070)

摘要：在分析现有滤波算法的基础上，结合数据驱动和模型驱动式算子各自的优点，提出基于点云空洞修复和TPS变形模型的数学形态学机载LIDAR点云滤波，该方法首先提取和修复由水域造成的大面积点云空洞，采用多尺度形态学开算子作用于修复的数据，得到近似裸露地表面；然后利用2D空间的TPS变形模型，以近似地表面为基础，插值原始点云，根据插值与原始点云高程的差值大小去识别地面点和非地面点。通过定量分析，验证该方法不仅有较高的滤波精度，而且也能较好的保留裸露地表的细节特征，同时该方法有助于辅助人工处理，提高数据处理的质量。

关键词：空洞修复；多尺度形态学开算子；2D空间TPS模型；机载LIDAR滤波

机载LIDAR技术是一种主动式遥感技术，它能够快速获取大面积区域内高精度的空间三维地形信息，已经在地球空间信息学科的许多领域得到广泛应用[1-2]。机载LIDAR获取的原始数据，其结构是离散的，以不规则方式分布在地形表面，因此称为点云。以点云为基础进行滤波分类，是该系统数据后处理重要的组成部分，也是进行后续高层次数据处理的基础。而点云数据滤波及其配套的质量控制是生成DEM最关键也是最耗时的一步，约占后处理流程60%～80%的时间,因此找到一种滤波精度较高的方法对提高DEM的质量是至关重要的，同时近几年DEM一个重要的应用领域就是虚拟地理环境(Geographic Virtual Environments)，要提高人在虚拟环境的沉浸感，对地形起伏的仿真度提出更高要求[3]。这就要求滤波算法不仅要尽量减少分类误差，同时也要尽可能保留地形特征的细节信息[4]。

目前滤波算法大致可以分为两类：一类是数据驱动式的滤波算法，例如基于数学形态学的方法[5]、基于TIN结构的分层方法[6]；另一类是模型驱动的滤波算法，例如基于分层稳健线性估计的滤波算法[7]、基于样条插值和区域增长的滤波算法[8]。这两类算法大多都从单一的高差或坡度作为分类的依据，对地形情况复杂的场景很难获得理想的结果，另一方面这两类算法几乎都没有考虑到大面积的数据空洞对滤波保留地形特征的影响。

针对这两类滤波算法所存在的这些不足，本文提出基于形态学理论和TPS变形模型的机载LIDAR点云滤波，该算法对点云的滤波不仅综合考虑高差和坡度的影响，而且也减少大面积的空洞对滤波保留地形特征的影响。通过试验证明该方法不仅有较高的滤波精度，而且也能较好地保留地形的细节特征信息,同时该方法也有助于辅助人工处理，提高数据处理的质量。

1算法原理及流程

1.1近似裸露地表的获取

1.1.1水体边界的提取及填充

首先需要将离散的点云栅格化，根据三维坐标的最大最小值，确定栅格图的范围，每个栅格值的大小取1 m范围内最近点的高程值。由于水体对激光的吸收，导致传感器无法接收到回波信号，从而使部分区域点云过于稀疏,这可能会导致有一部分栅格无值。通常水体的高程在一个相邻的区域内高程最低，因此可以用水体边界高程值最低的点填充水体空洞区域[9]。首先将栅格图生成二值图f(x)，这里引入形态学梯度算子检测水体空洞区域的外边界，算子如式(1)所示。

(1)

其中g称为结构窗口，窗口大小根据式(2)确定。

(2)

图1　水域边界提取

盲目的用高程最低值填充点云空洞会造成局部地形区域的不连续性，破坏地形特征[10]，因此水体造成的大面积空洞，采用水体边界高程最低值的点填充,最终得到以点云高程值为基础的栅格表面SF。

1.1.2获取近似裸露地表面

离散点云形成栅格图的时候，采用最近点填充栅格的方式，这必然会导致栅格图像中有植被、房屋等非地面栅格单元，因此必须剔除这些非地面栅格，这里采用多尺度形态学开算子，该方法能够较好的保留地形起伏的裸露地表特征[11]。其过程如下：

1)将栅格表面SF赋给SL，形成多尺度形态学窗口向量w={w1,w2,…,|wi≤wmax}，其中

wi=2bik+1.

(3)

式中:k=0.5,bi=i(i=0,1,2…)。

2)计算栅格滤波阈值et。从w中取出一个值

et=swic.

(4)

式中：s为测区的最大坡度估值；c为栅格大小。

3)对表面SL执行形态学开算子得到表面ST。开算子定义为先对目标图像f通过结构窗口g被腐蚀，然后在腐蚀的目标图像上进行膨胀操作，其数学过程为

(5)

其结构窗口g大小为wi。

4)标示地面栅格。当ST-SL>et,在SL上对应栅格标示为地面栅格。

5)执行SL=ST,重复2)到5)的步骤直到wi等于wmax为止。

通过以上的步骤，最终可以得到机载LIDAR扫描区的一个近似的表面S，该表面位于该区域DEM与DSM表面之间，因此可以作为点云滤波的参考表面。

1.2点云滤波

1.2.12D空间的TPS变形模型

薄板样条函数TPS模型模拟一个薄金属板在点约束下的弯曲变形，用一个简单的代数式表示弯曲的能量变化，属于非线性变换方法。假定有n观测值z(xi,yi)，则观测值可由式(6)表示[12]。

(6)

其中i=1,2,…,n，f(xi,yi)是观测值z(xi,yi)的估值函数，ε(xi,yi)是估计误差。

若使估计函数最大限度的接近观测值，则需要使用最小化函数实现，算式为

(7)

其中平滑参数λ>0,本文中λ=1，其中If为估计函数的平滑项，定义为

(8)

f(xi,yi)定义为

(9)

(10)

1.2.2滤波过程

TPS相对于其它的插值方法具有反映高程异常变化的物理特性，具有光滑、连续、弹性好的特点，而且该方法不需要测量点呈规则分布[13],因此该方法较适合离散点云的插值，它能插值出光滑的表面，有利于提高插值精度[14]。滤波的过程如下：

2)求每个离散点Pi(xi,yi,zi)对应地表面S上的近似坡度，计算滤波阈值hyi。坡度的算式为

(11)

其中，Fx(i,j)，Fy(i,j)为水平垂直梯度，水平梯度算式为

(12)

式中：i,j为第i个点对应栅格S的行列号，需要注意的是：第1列上的水平梯度计算方法为第2列与第1列的差值，最后一列为最后两列的差值。垂直梯度与水平梯度计算方法类似。

(13)

式中：hy为固定的一个阈值；sc为缩放因子(一般取值为：0.0到2.5)。

2试验结果与分析

实验数据是由国际摄影测量与遥感协会网(http://www.itc.nl/isprswgIII-3/filte-rtest)提供的8种场景LIDAR数据，专门用于测试滤波算法，其中城市和乡村地区的数据各占 4 景;这些数据包含平原、植被、建筑物、公路、铁路、桥梁、电线、水域、陡坡等不同的地形地貌特征。

本文算法用MATLAB实现，定量分析采用ISPRS提出的I类误差、II类误差和总误差作为算法精度的评定准则，参数和实验结果如表1所示。

Terrasolid开发人员将上述实验数据在Terracan中半自动化地进行了处理，滤波结果达到该处理的最佳状态。这里将该软件的处理结果与本方法进行比较，比较结果如图2所示。

图2　My Method 与TerraScan的比较

I类误差和II类误差反映算法的适应性，总误差反眏算法的可行性[15],图2中该方法与Terracan软件中的算法滤波结果比较，一类误差和总误差明显优于Terracan软件的滤波，这反映出该算法能够较好的保留裸露地表上的点，从机载LIDAR点云中能够提取质量较高的DEM；总误差较低能总体反映算法在实际生产应用中具有较高的应用价值和优越性，II类误差虽然较Terracan软件的滤波结果高，但是平均误差也能控制在7.6%以下，这也表明该方法能够适应不同复杂程度的地形环境。表1中该方法II类误差较I类误差稍大，这一特点与当下大部分滤波算法刚好相反(即I类误差大于II类误差)，经验表明，与Ⅰ类误差相比，Ⅱ类误差的点云大部分有着高于地面点云的高程，所以在可视化环境下，人工干预时能更容易地识别和修复Ⅱ类误差带来的分类错误，这能较好地提高数据处理效率和数据质量。因此该算法无论是自动化处理还是辅助人工处理都能较好提高滤波精度和数据质量。

如图3所示，S12和S52分别为城市和山区的浮雕图，图3(b)为人工逐点分类的结果图，可认为是真实的裸露地表。对比图3(c)、(b)和(a)可以明显地看到，本文方法无论在城区还是山区都能够较好地保持原始裸露地貌，而且对于地形起伏较大的山区也能保留地形起伏的轮廓细节特征；图中红色箭头所指示的陡坎、斜坡区域内的植被几乎都被剔除，较好地保持了原始的地形地貌，没有出现过度剔除地形特征的现象。这种情况在虚拟现实中能够改善场景的真实性，为客户提供较满意的虚拟沉浸感。

3结论

本文结合数据驱动和模型驱动式算子各自的优点，提出基于点云空洞修复和TPS变形模型的数学形态学机载LIDAR点云滤波，通过试验验证该方法相对于商用软件Terracan，滤波精度较高。I类误差平均为2.19%，II类误差平均为7.54%，总误差平均为3.06%，因此该算法不仅提高滤波精度，而且该算法能够适应不同复杂程度的地形环境，适应性强;I类误差小于II类误差，说明该方法能较好地保存原始裸露地表上的点，用该方法提取的DEM用于三维可视化和虚拟现实中，可以明显提高场景真实感和沉浸感，所以在可视化环境下，人工干预时能更容易地识别和修复Ⅱ类误差带来的分类错误，这有助于辅助人工处理提高工作效率，提升数据处理的质量。当然该方法有些不足，比如在试验区S53中二类误差较大，S53为人为改造出现的间断地形，因此该方法应进一步改进，以降低在间断地形上的二类误差。

图3　S12和S52浮雕图

参考文献：

[1]吴芳,张宗贵,郭兆成,等.基于机载LiDAR点云滤波的矿区DEM 构建方法[J].国土资源遥感,2015,27(1):62-67.

[2]肖春蕾,郭兆成,张宗贵,等.利用机载LiDAR数据提取与分析地裂缝[J].国土资源遥感,2014,26(4):111-118.

[3]PINGEL T,CLARKE K C.An improved simple morphological filter for the terrain classification of airborne LIDAR data[J].ISPRS Journal of Photogrammetry and Remote Sensing,2013,77:21-30.

[4]张小红.机载激光雷达测量技术理论与方法[M].武汉:武汉大学出版社，2007:93-95.

[5]沈晶,刘纪平,林祥国. 用形态学重建进行机载激光雷达数据滤波的方法研究 [J]. 武汉大学学报(信息科学版), 2011, 36(2)：167-170.

[6]亢晓琛，刘纪平，林祥国．多核处理器的机载激光雷达点云并行三角网渐进加密滤波方法[J]．测绘学报,2013,42(3):331-336.

[7]PFEIFER N, STADLER P, BRIESE C. Derivation of Digital Terrain Models in the SCOP Environment [C]. OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models. Stockholm:Sweden, 2001:1-12.

[8]成晓倩，赵红强.基于区域生长的LIDAR点云数据滤波[J].国土资源遥感，2008，20(4)：6-8.

[9]QI Chen,PENG Gong,DENNLS B,et al.Filtering Airborne Laser Scanning Data With Morphological Methods[J].Photogrammetric Engineering & Remote Sensing,2007(73):175-185.

[10] 刘文.基于数学形态学原理和TerraScan的Lidar点云数据分类方法研究[D].青岛：国家海洋局第一研究所，2008:28-45.

[11] 隋立春，张熠斌，柳艳，等. 基于改进的数学形态学算法的LiDAR点云数据滤波[J].测绘学报，2010,39(4)：390-396.

[12] BOER E P J, MDEBEURS K,HARTKAMP A D.Kriging and thin plate splines for mapping climate variables[J].International Journal of Applied Earth Observation and Geoinformation,2001,3(2):146-154.

[13] 杜国军，贾良文. 薄板样条函数在空间数据插值中的应用[J]. 计算机工程与应用,2009,45(36):238-240.

[14] MONGUS D, ZLIK B. Parameter-free ground filtering of LiDAR data for automatic DTM generation[J]. ISPRS Journal of Photogrammetry and Remote Sensing,2011,67:1-12.

[15] 吴丛丛,卢小平. 基于TIN的LiDAR数据滤波算法研究[J].测绘通报，2013(3):32-34.

[责任编辑：张德福]

Filtering of airborne LiDAR based on mathematical morphological of the repairing point cloud hole and TPS model

HU Yongjie1,2,YAO Xiaowei1,HU Qianwei3,HAN Fengna1

(1.XPCC Surveying & Designing Institute(Group) Co.,Ltd., Urumqi 830002,China;2.Institute of Surveying and Mapping, East China Institute of Technology, Nanchang 330013,China;3.Beijing Geo-VisionTech.Co.,Ltd., Beijing 100070,China)

Abstract：Based on the basic analysis of existing filtering algorithm,combiningrespective advantages of data-driven and model-driven,a filtering method is proposed based on mathematical morphological of the repairing point cloud hole and TPS model.This method first extracts and repaires the point cloud,using multi-scale morphological opening operator acting on the restoration of data to form an approximate bare ground surface.Then by the TPS deformation model of 2D space based on approximate surface,interpolation of the original point cloud,the difference between the original point cloud elevation interpolation,ground point and non-ground points are identified.In this paper,through quantitative analysis shows that the method not only has high filtering precision but also better retains the minutiae of bare surface.Meanwhile,the method can also help auxiliary artificial processing and improve the quality of data processing.

Key words：repairing holes;multi-scale morphological opening operator;2D space TPS model;airborne LiDAR filtering

DOI:10.19349/j.cnki.issn1006-7949.2016.07.006

收稿日期：2015-08-20

作者简介：胡永杰(1986-)，男，助理工程师.

中图分类号：P237.3

文献标识码：A

文章编号：1006-7949(2016)07-0028-05