基于数学形态学的LiDAR数据滤波新方法

罗伊萍，姜 挺，王 鑫，张 锐，罗 胜

(信息工程大学测绘学院，河南郑州 450052)

基于数学形态学的LiDAR数据滤波新方法

罗伊萍，姜 挺，王 鑫，张 锐，罗 胜

(信息工程大学测绘学院，河南郑州 450052)

数学形态学滤波是从激光雷达点云数据中识别地面点、创建数字高程模型的一种重要方法。在分析现有滤波方法的优劣性以及数学形态学滤波方法存在的问题基础上，提出一种新的具有一定自适应性的数学形态学滤波算法。该方法通过分析LiDAR数据的特点，利用形态学算子提取原始点云数据的空白区域，进行精确的重采样生成DSM;然后将DSM多尺度形态学滤波的结果作为初始的DEM，分析由于局部地形突变导致可能存在误分类区域;最后提取该区域边缘外的地面点作为新的种子点，利用最小二乘平面拟合的方法进行区域生长搜索误分类地面点，从而提高分类精度。试验结果表明该方法能够有效识别地面点和地物点，并且保留地形的细节信息。

激光雷达;滤波;数学形态学滤波;数字高程模型

一、引 言

机载激光雷达(light deteation and ranging，LiDAR)技术自20世纪80年代引入摄影测量领域以来，受到广大研究者的极大关注。LiDAR系统利用机载激光雷达测距系统和GPS/IMU能直接获取地面点的三维坐标，形成离散的、不规则的三维点云数据。为了从离散的点云数据中获取数字高程模型(DEM)，必须对点的属性进行区分，将原始的点云数据分为地面点集和非地面点集，这个过程被称为滤波。

目前，已有不少国内外学者对激光雷达点云数据的滤波方法进行了深入的研究，主要的滤波算法可分为以下三类:基于表面的内插滤波算法［1-2］、基于区域的滤波算法［3-5］和基于约束曲面的滤波算法［6-7］。基于表面的内插滤波算法的核心思想是通过一个较粗的起始DEM，逐步从备选数据点筛选并内插加密DEM达到分类的目的。因此，滤波结果受到初始DEM影响较大并且误差会随着迭代的过程积累。而基于Snake样条曲面等约束曲面的滤波方法认为地面是一个连续且平缓变化的表面，这样就过于强调地形的平缓变化而忽略了地形的复杂性［8］。

数学形态学滤波方法属于基于区域的滤波算法，与基于表面的内插滤波算法相反，它是一种自下而上、从局部出发扩展到全局的滤波方法。借鉴栅格图像的处理思想，根据激光点云数据生成的深度图像，利用形态学开运算剔除高于地面的点，通过逐步增大滤波窗口，得到逼近地形的一个表面。然而在滤波过程中，仅考虑窗口内的地形特征，容易受到局部地形的影响，例如道路旁存在沟渠，当滤波窗口尺寸大于道路宽度时会导致路面点被分为非地面点。另外，滤波过程中的参数需要用户根据不同地形设置，自适应性不强［3-4］。针对数学形态学滤波算法存在的问题，本文提出一种新的多尺度数学形态学滤波策略用于LiDAR数据滤波。

二、多尺度数学形态学滤波方法

数学形态学的语言和理论是基于集合运算原理提取图像中的特征。腐蚀和膨胀运算是形态学图像处理的基础，通常用于“减少”(腐蚀)或“增大”(膨胀)图像中特征形状的尺寸。基于灰度图像的腐蚀和膨胀运算是在结构元素定义的邻域内选择图像像素值和结构元素相作用后的最小或最大像素值。

在LiDAR点云数据生成的规则化数字表面模型(DSM)中，腐蚀和膨胀运算可以定义为:

腐蚀

膨胀

其中，f为规则化的DSM;g为结构元素;Z(i，j)为腐蚀或膨胀运算后规则化的DSM中第i行第j列的高程值;w为结构元素的窗口。

腐蚀和膨胀运算组合后，形成开运算和闭运算用于LiDAR点云数据滤波，开运算和闭运算可以定义为:

开运算

闭运算

式中，开运算通过先腐蚀后膨胀的方法，先将比结构元素尺寸小的树木点等非地面点从原始点云数据中移除掉，然后通过膨胀运算恢复被腐蚀掉的建筑物边缘形状。然而，利用一个固定窗口大小的结构元素很难移除各种尺寸的非地面物体，如果窗口过小，则大部分地面点被保留，只有小的非地面目标，比如车或者树木被移除，而城区的大型建筑物不能被移除。因此必须利用多尺度的滤波窗口进行迭代计算，在迭代过程中逐步增大滤波窗口的尺寸。

每次迭代开运算时，滤波的窗口尺寸wk按指数增长，wk=2bk+1，式中 k为迭代的次数，k=1，2，…，M;b为初始的窗口尺寸大小，一般为2。为了保证建筑物都被去除掉，最后一次迭代的窗口尺寸wk必须大于试验区域中最大建筑物的面积，一般迭代次数为5或6次。为了确保在滤波的过程中保留细节的地形特征，设置每次迭代高差阈值dhT，k

式中，wk为第k次滤波的窗口大小;c为DSM的格网间距;dh0和dhmax分别为最小和最大高差阈值;s为地形坡度参数。从式(5)可以看出，当滤波窗口增大时，高差阈值dhT，k随之增大，增幅大小由地形坡度s决定。如果DSM中某个格网点开运算前后的高程值之差小于本次迭代的高差阈值dhT，k，则认为是地面点，否则认为是非地面点。

三、DEM提取过程

1.点云数据格网化

1)在点云数据的xy平面范围内，选择合适的格网间距构成m×n平面格网。对于每个离散的激光脚点，根据其平面坐标将其分配到相应的格网中。为了使尽可能多的激光点用于初始的DSM，一般选取比点间距稍小的格网间距值。由于点云数据分布的不均匀性，如图1(a)所示，一个格网中可能会没有或者有一个或多个激光点。如果一个格网内有多个激光点落入，高程值越低的点作为地面点的概率越大，因此选取最低高程值作为格网的取值。记录下每个格网中激光点的原始索引号，用于最后的精度评定。

2)在初始的DSM中，由于选取了格网内所有激光脚点的最低值，因此由激光扫描系统误差和随机误差引起的少量无意义低值粗差点也被保留下来。考虑到粗差点与相邻点之间的高差形成脉冲波峰，在初始表面模型中的5×5邻域内，如果中心点高程值与24个邻域点高程值的差值均大于高差阈值3 m，那么判定该点为噪声点并且该格网中没有激光点。

以m×n的二值标记矩阵记录下对应格网中是否有激光点，矩阵值为0则没有激光点，如图1(c)所示。

2.DSM重采样

DSM重采样是针对格网化后的DSM中的空白格网进行数据内插。由于激光雷达的视场角比传统摄影测量相机视场角小，飞行时容易造成扫描漏洞，如图1(c)的顶部有细长的扫描漏洞。另外，水体对1 064～1 600 nm波段的激光有吸收效应，也会造成数据空白区域，如图1(c)的左侧中部区域。如果存在高层建筑物，由于建造物的遮挡也会产生数据空白区域。由于这种数据的不均衡性，导致数据内插的搜索半径的大小很难确定，尤其是存在大量湖泊河流的区域，并且对建筑物遮挡的空白区域内插生成的DSM扩大了建筑物的面积。因此，要先针对这种大面积的数据空白区域进行数据重采样，然后以较小的搜索半径对剩余的空白格网内插，以提高重采样的效率和精度。

首先，利用3×3的形态学窗口对二值标记矩阵作闭运算，去掉由于激光点分布不均匀造成的稀疏空白格网。其次，为了消除水域上的零散激光点，采用小面积消除法，求得标记矩阵内所有值为1的四连通区域，其个数小于40的连通区域认为是由水体反射的激光点造成的，将它们赋值为0。然后，二值矩阵取反，再次使用小面积消除法消除零散的无激光点格网，得到图1(d)。此时仅保留了面积大于40个格网数的水域和因飞行扫描漏洞造成的数据空白区域，为了求得这些区域的边缘，将标记矩阵用3×3的形态学窗口作膨胀运算，再减去原标记矩阵即得到如图1(e)所示的边缘。一般认为水域是比较平缓的高程一致性区域，因此将水域边缘的最低值作为水域的高程取值。扫描漏洞和建筑物遮挡区域也作相同处理。最后，由合适的搜索半径对其余零散的空白格网进行最邻近值内插，得到重采样的DSM，如图1(b)所示。

图1 点云数据重采样

3.多尺度数学形态学滤波

在点云数据规则化和重采样后，本文采用上述的基于坡度和高差阈值的多尺度数学形态学滤波方法。针对城市区域地形起伏较小以及建筑物的最大面积一般不会超过120 m×120 m的实际情况，滤波的参数选择如表1所示。

表1 滤波参数设置

4.质量控制

由于数学形态学滤波仅考虑窗口内的地形特征，在局部地形有突变的情况下，会导致相当数目的地面点分类为非地面点，从而影响分类的精度。在这种情况下由于地面的连续性，误分类区域在一定范围内与另外的地面区域相连接，通过提取相连接地面区域的边缘点作为种子点，通过区域生长法搜索误分类区域中的地面点，从而改进分类精度。

将多尺度数学形态学滤波后的结果作为初始的DEM，并且从中提取非地面点集的二值化标记图像。用小面积消除法剔除标记图像中小面积的树木点，再用3×3的窗口对标记图像进行膨胀运算再减去标记图像，得到大面积非地面点区域的外边缘点，即相连接地面区域的边缘点。从这些边缘点中，根据它们3×3邻域的最大和最小高程值之差，剔除建筑物的外边缘点，保留有可能存在误分类地面点集的边缘点作为种子点。取种子点5×5邻域中的地面点作最小二乘平面拟合，以拟合平面作为初始地面，搜索标记图像中与初始地面的高差小于高差阈值的点，作为误分类的地面点。

四、试验与分析

1.试验结果

试验采用ISPRS网站提供的点云数据，这些点云数据是在2003年由ISPRS Commission III提供给广大学者的，用于比较不同滤波算法的试验。试验选取城市区域的四块测试数据Site1-Site4的末次回波激光点作为原始数据，点间距为1.5 m，重采样的格网间距为1 m。试验区域散布着大、中、小型的不规则形状的建筑物，同时有桥梁、隧道、铁路、汽车以及与建筑物相邻的树木。在四块测试数据中选取了九个样本数据，对样本数据进行了手工分类，将激光脚点精确分类为地面点和非地面点两类。利用分类后的样本数据，可以对滤波算法进行分类误差的定量分析。

图2(a)给出了应用多尺度数学形态学算法对重采样后Site2的DSM进行滤波的结果，其中空白地区为分类后的非地面点。从图2(a)可以看出，在A、C处由于地下通道与旁边的路面存在高程突变，而在B处也由于地面不连贯变化，导致部分地面点明显被误分为非地面点。对滤波结果进行质量控制，搜索误分类的地面点后得到最终分类结果，如图2(c)所示。

图3分别给出了Site1、Site3和Site4的分类后地面点重采样生成的DEM。试验结果表明该算法能够有效地识别地面和地物点，地面细节信息得到保留;树木被有效识别，所有建筑物被滤除掉。

图2 Site2试验结果

图3 地面点重采样的DEM效果图

2.试验分析

文献［9］对不同滤波方法从算法上进行了试验比较，表2给出了本文滤波方法与文献［9］和文献［4］中所涉及算法的总误差对比。

表2 样本的总误差对比

由表2可以看出，在陡坡地带的样本11、复杂建筑物密集的地形高程变化较大的样本23以及铁路等高频率的地形起伏样本42地区，Axelsson的内插滤波方法的精度稍微高于本文和文献［4］的形态学方法，因为从粗到细的处理方式能避免大地形起伏导致距离计算值较大而引起的误判。而在地形平缓、建筑物密集的样本12、样本21、样本31以及样本41地区，本文方法明显优于内插方法和文献［4］的形态学方法，尤其是数据样本21和样本41，这是由于多尺度形态学滤波方法采取由局部出发扩展到全局的策略，在平缓地区不容易导致误判，尤其是在有数据空白和一个大型建筑物的样本41，Axelsson的内插滤波方法对数据缺失比较敏感，而文献［4］的形态学方法以高程较低值填补数据空白后，则会将高于数据空白地区的周边地面点分类为非地面点，同时该方法的滤波窗口最大值小于建筑物的面积，导致部分建筑物点被分类成地面点，而本文以形态学滤波的结果作为初始DEM，通过进一步搜索误分类的地面点明显提高了分类的精度。

五、结束语

本文提出一种新的数学形态学滤波策略，该方法对点云数据进行合理的重采样，并用多尺度数学形态学滤波后作为初始的DEM，通过搜索误分类地面点的方式提高分类精度。该方法解决了数学形态学滤波算法中存在的两个问题:局部地形不连续对分类的影响以及坡度或高差阈值的人工选取。试验结果表明该方法能够有效识别地面点和地物点，并且保留地形的细节信息。

［1］KRAUS K，PFEIFER N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1998，53(4):193-203.

［2］PFEIFER N，STADLER P.Derivation of Digital Terrain Models in the SCOP++Environment［C］∥Proceedings of OEEPE Workshop on Airborne Laser Scanning and Interferometric SAR for Detailed Digital Elevation Models.Stockholm:［s.n.］，2001.

［3］ZHANG K，CHEN S C，WHITMAN D，et al.A Progressive Morphological Filter for Removing Nonground Measurement from LiDAR Data［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41(4):872-882.

［4］CHEN QGONG P，BALDOCCHI D，et al.Filtering Airborne Laser Scanning Data with Morphological Methods［J］.Photogrammetric Engineering and Remote Sensing，2007，73(2):175-185.

［5］张熠斌，隋立春，曲佳，等.基于数学形态学算法的机载 LiDAR点云数据快速滤波［J］.测绘通报，2009(5):16-18.

［6］ELMQVIST M JUNGERT E，LANTZ F，et al.Terrain Modelling and Analysis Using Laser Scanner Data［J］.International Archives of Photogrammetry and Remote Sensing，2001，34(3/W4):219-224.

［7］BROVELLI M A，CANNATA M.Digital Terrain Model Reconstruction in Urban Areas from Airborne Laser Scanning Data:the Method and an Example of the Town of Pavia(Northern Italy)［J］.Computers ＆ Geosciences，2004，30(4):325-331.

［8］黄先锋，李卉，王潇，等.机载LiDAR数据滤波方法评述［J］.测绘学报，2009，38(5):465-469.

［9］SITHOLE G，VOSSELMAN G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds［J］.ISPRS Journal of Photogrammetry and Remote Sensing， 2004，59(1-2):85-101.

［10］穆超，余洁，许磊，等.基于高分辨率遥感影像的DSM建筑物点的提取研究［J］.武汉大学学报:信息科学版，2009，34(4):414-417.

A New Filtering Method for LiDAR Data Based on Mathematic Morphological Approach

LUO Yiping，JIANG Ting，WANG Xin，ZHANG Rui，LUO Sheng

0494-0911(2011)03-0015-05

P237

B

2010-02-23

罗伊萍(1982—)，女，湖南永州人，博士生，研究方向为LiDAR数据处理和航天遥感工程。