基于多尺度的城区机载点云滤波方法研究

钟良，汤璇，管海燕，邬建伟

1.长江委长江空间信息技术工程有限公司，武汉430010

2.武汉军械士官学校计算机系，武汉430075

3.加拿大滑铁卢大学地理与环境学院，加拿大滑铁卢，N2L2R7

4.武汉大学遥感信息工程学院，武汉430079

1 前言

机载Lidar点云数据为三维空间中呈现随机分布的点云。在点云中，有些点是真实地形点，有些则是人工地物（比如建筑物、桥梁、塔、车辆等）或者是自然地物（如树木、灌木等）。从激光扫描点云中区分地形点子集与非地形点子集（包含人工地物与自然地物）称之为滤波。目前，自动区别地形点和地物点一直是研究的热点，特别是针对场景中含有很多各种复杂地形地貌特征情况。目前学术界已经提出了很多滤波方法：数学形态学[1-4]、线性预测[5-6]、渐进加密[7-8]和分割[9-11]，其他滤波方法还有全波形[12]、机器学习[13-14]、扫描线[15-16]。

以上算法都是基于对地面点和非地面点的理解概念或者思路来设计的，从算法设计角度来说可分为直接应用于不规则分布的点集，或者将原始点集重采样为规则格网再利用成熟的图像处理技术进行处理；从算法对激光点的几何关系的分析上可分为考察单点与单点之间的关系，单点与点集或者点集与点集之间的关系。从算法对分析点云数据特征中某种不连续测度上可分为高程差、坡度、到结构表面的最小距离或到参数平面的最小距离等。总的来说上述的各种滤波算法都基于各种不同的假设前提，如坡度、平面、表面和分割等；有些算法的处理方式是一次单步完成，有些是通过迭代处理，这也意味着它们分别在准确性与速度方面的存在着或多或少的劣势[17]。

考虑到对复杂城区点云数据进行滤波、生成DTM，针对城区建筑物偏多，局部地形较为平坦的特性，本文提出了一种以离散度分析为基础的平面点提取与基于强度方差的区域增长相结合的组合滤波算法。算法分为三步：首先计算每个点的局部空间分布特性，将激光点大致分为平面点，边缘点和离散点；其次将平面点作为分析对象跟踪出一组联通的区域，利用开运算和闭运算去除零散的小区域，基于强度方差测度将其分类为地面，非地面以及未确定；最后后向搜索从边缘点以及未确定点中加密地形点。

2 复杂城区激光点云滤波算法

本算法分为三个步骤：基于离散度的初始激光点粗分类，基于区域的精分类，以及后向地形点加密。算法首先分析并计算每个点与其周围一定邻域内激光点之间的几何特征值关系，将点云粗分类为平面点、边缘点和离散点；然后对平面点进行区域跟踪，利用高程以及强度方差等将平面点分类为地面点、建筑物点以及为确定类别点；最后对地面点构建不规则三角网，反向分析确定类别点以及边缘点用以加密地面点集。

2.1 基于离散度的初始激光点粗分类

激光点云空间离散度可以通过计算每个点与其周围一定邻域范围内所有点的空间几何分布来获取。一般树木、小型地物（如汽车等）等在空间各个方向上都较为分散，而裸露地面以及建筑物在空间的离散程度可以认为沿二维表面（不一定水平）分布，这种离散性可以通过离散矩阵的特征值来分析。具体方法：首先构建激光点云的二维K-D树，搜索激光点云中每个点邻域内的全部相邻点，建立该激光点在空间上的3×3离散矩阵，见公式（1）：

其中，Sj为第j个点的3×3离散矩阵，n为第j个点邻域内相邻点数目，vi为第j个点的第i相邻点的空间坐标vi=(xi，yi，zi)，为vi的转置矩阵，M为激光点云点数。

将离散矩阵Sj作奇异值分解，可获取该点矩阵的三个特征值λ0、λ1、λ2，并将特征值从大到小排列，e0、e1、e2为三个特征值对应的特征矢量[18]。则公式（1）可以表示为：

根据特征值，一般设定三个类别：（1）λ2＜＜λ1≈λ0，则该点被标记为平面类（图1（a））；（2）λ2≈λ1＜＜λ0，则该点被标记为边缘类（图1（b））；（3）λ0≈λ1≈λ2，并且都足够大，则标记为空间离散类（图1（c））。

图1 离散度计算

根据以上的描述，无法确定阈值来判定这三个类别。因此公式（3）可表示为：

在(Slinear，Splanar，SSphere)中，最大的值就为该Sj所属的类别。一般大部分房屋顶由平面组成，因此呈现出平面特征，但是其中部分由于屋顶结构复杂、形状多变，因此会有些屋顶出现少量离散的点，且地面上小物体（如小汽车等）也会造成少量离散点；基本上裸露地面会出现大量的平面点特性；而由于城市中树木稀疏，因此各向离散性能够很好地表现出来，这样就可以去除大部分植被点以及建筑物边缘点等。经过离散度计算，平面点作为下步处理的输入数据来获取地形点。

2.2 基于区域的精分类

经过离散度计算后，Lidar点集中保留有地面点和建筑物点，以及其他零散的激光点，为此需要对区域进行精分类以将地面点从建筑物点中分离出来。由于地面是由连续联通的面片构成，因此本文不是处理单个的点，而是将区域作为分析单元。如图2所示为本阶段的流程图。

图2 基于区域分类滤波流程图

（1）建立Delaunay三角网数据结构。其目的是用来寻找临近点，形成一个个联通区域。

（2）开运算与闭运算的运用。由于地形和建筑物也是具有一定面积的连续平面，开运算和闭运算可用于去除那些激光点很少的小区域。一般这个点数量阈值Tφ需要根据点密度来参考。在本研究中，点数量阈值Tφ设为20，如果区域内点的数量小于20个点，这个区域将被归并到边缘点中。

（3）计算全局高程阈值和局部高程阈值。全局高程阈值确定HG是根据输入数据的高程直方图分析获取的。如图3所示，对建筑物和裸露地面，在高程直方图中会出现多个波峰，一般第一个波谷就是区别二者的全局高程阈值（图中大约为150m）。局部高程阈值的确定HL是通过计算每个待分类区域的最大范围W以及中心点坐标。然后基于中心点，在原始点云获取一个2W的区域。则2W区域内高程平均值被作为局部高程阈值HL。

图3 全局高程阈值的确定

（4）如果待判定区域，其平局高程均大于全局阈值和局部阈值，则该区域被认为是建筑物平面；区域平均高程均小于全局和局部阈值，则被分类为地面点；剩下的不符合这两个判断条件的，将进行进一步的分析。

（5）根据上一个步骤所获取的地面点和建筑物类别点，再分别计算已分类地面点和建筑物点的强度方差，其方差计算公式如下：

通过以上三个步骤，平面点基本分类为地面点，建筑物点和待定点。地面点则被用于建立Delaunay三角网作为初始地形。待定点和前面分类为边缘点将被带入下一个步骤，进行后向地面点加密获取。因为边缘点一般是处于断裂线附近，其中有部分地面点。

2.3 后向地形点加密

对于任意一个Lidar点如果同时满足两个地形判定条件，则将其作为地形点。假设N1～N3为已判定的地形点，Lp为待判定激光点，d为待判定激光点Lp到N1～N3所构成的三角面片的距离，(α，β，λ)分别为Lp到N1～N3构成的三角面片的夹角。给定两个阈值：距离阈值dmax和角度阈值θmax，则两类判据分别定义为[7]：

其中，如果待定激光点Lp同时满足这两个判据，则Lp被认定为地形点。

图4 判据计算示意图

将初始地面点集构建不规则三角网，从边缘点集中需找潜在的地面点加入到地面点集中。遍历不规则三角网，计算落在三角形中每个激光点到三角形的距离以及三个夹角，如果满足给定的高程差阈值以及角度阈值则认为该点是地面点。在TIN中所有三角形遍历结束后，将获得的地面点重新加入不规则三角网。重复迭代执行，直至满足给定的迭代最大次数或者迭代搜索获取的地面点小于给定最少地面点，迭代结束。

3 实验结果与分析

实验1 ISPRS测试数据

本文从ISPRS网站提供的实验数据中选择了四个城区测试数据集中的几个典型场景（数据来自http：//www.commission3.isprs.org/wg3/index.html），如图5所示。其中点云的水平分辨率均为0.67个点/m2（点间距为1.0～1.5m）。

图5 ISPRS sample数据

对激光扫描数据进行滤波的主要目标是生成DTM，因此滤波结果中存在的两种误差将影响DTM的生成质量。第一种误差称为I型误差，指的是将地形点错误的分类为非地形点；第二种误差称为I型误差，指的是将非地形点错误的分类为地形点。I型误差造成精度上的损失，II型误差将导致DTM质量的恶化。因此I型误差和II型误差也成为评价滤波算法的主要指标。根据ISPRS提供的测试数据，表1为对本文滤波结果的I型误差和II型误差。

根据表1，其滤波结果还是相当满意的，尤其是Type II误差较小，基本上保持了地形的特征。数据Sample11、Sam ple22和Sample24的I型误差较高，产生的主要原因是地形的坡度变化比较大。由于算法设计主要是针对城市地区，所有滤波的假设条件就使得算法对这类地形具有一定的缺陷。Sam ple22、Sample23的II型误差相对较高，主要原因是一些建筑物和植被小于算法设定的阈值，因此将一部分点错误的分类为地面点了。

实验2 三个不同国家地区数据实验

图6 三个不同国家地区数据实验

表1 滤波数据的I型误差和II型误差分析表

实验2选取了三块数据进行实验，数据一是山西太原部分城区激光点云，如图6（a）所示；数据二是加拿大多伦多城区数据，如图6（b）所示；数据三是德国Toposys公司提供的Mannheim城区数据（点平均密度为4 points/m²），如图6（c）所示。这三块原始点云数据是根据高程进行分色显示的，因此从图中可以看出这三块地区地势较为平坦，是典型的城区地貌。图6（d）、（e）和（f）是三角网加密迭代滤波实验结果。由于没有实地的DEM数据做参考，难以定量地衡量其提取的DTM的质量，但是总体感觉上提取的质量还是不错的，基本上去除了非地形点，保留了地形特征。图6（f）显示的是将Mannheim城区激光点云滤波实验结果投影到其配准好的航空影像上，图中绿色像素是投到二维影像上的三维激光点。从实验结果可以看出，基本上将非地形点去除掉，保留了基本地形点。

4 结语

通过上述的实验比较与分析可知，本文提出的多尺度从粗到细的滤波算法不仅考虑了地形的连续性，也最大顾及了整个扫描区域地形的特征，其算法性能也能保证DTM的精度和质量。这种基于点与其周围领域内其他点空间关系来计算离散度能够很好地反映点的特性，可以去除大量的离散点和边缘点；同样通过基于区域的滤波精分类，很大程度上减少了候选地形点中的非地形点，提高DTM的质量，减少了II型误差。由于边缘点中包含地面点和非地面点，后向加密地形点可以发现更多的地面点，从而提高滤波的质量。实验证明利用多尺度的滤波算法能够在大面积、大数据量、复杂城区得到运用，为下一步Lidar数据的应用打下了很好的基础。

[1]Lindenberger J.Laser-Profilmenssungen zur topographischen Gelandeaufnahme[D].Munich：Deutsche Geodätische Kommission，1993.

[2]Kilian J，Haala N.Capture and evaluation of airborne laser scanner data[J].International Archives of Photogrammetry and Remote Sensing，1996，31：383-388.

[3]Vosselman G.Slope based filtering of laser altimetry data[C]//Proceedings of International A rchives of Photogrammetry，Remote Sensing and Spatial Information Sciences，WG III/3，Amsterdam，2000.

[4]Silván-Cárdenas J L，Wang L.A multi-resolution approach for filtering Lidar altimetry data[J].ISPRS Journal of Photogrammetry and Remote Sensing，2000，61（1）：11-22.

[5]K raus K，Pfeifer N.Determination of terrain models in wooded areas with airborne laserscanner data[J].ISPRS Journal of Photogrammetry and Remote Sensing，1998，53：193-203.

[6]Kobler A，Pfeifer N.Repetitive interpolation：a robust algorithm for DTM generation from Aerial laser scanner data in forested terrain[J].Remote Sensing of Environment，2007，108：9-23.

[7]Axelsson P.DEM generation from laser scanner data using adaptive TIN models[C]//Proceedings of International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences XXXIII（B4/1），1999：110-117.

[8]Sohn G，Dowman I.Terrain surface reconstruction by the use of tetrahedron model with the MDL criterion[C]//Proceedings of the Photogrammetric Computer Vision Symposium，2002.

[9]Voegtle T，Steinle E.On the quality of object classification and automated building modelling based on laserscanning data[C]//Proceedings of IAPRSIS XXXIV 3W 13，2003：149-155.

[10]Tóvári D，Pfeifer N.Segmentation based robust interpolation-a new approach to laser data filtering[C]//Proceedings of ISPRSWG III/3，III/4，V/3 Workshop“Laser Scanning 2005”，Enschede，the Netherlands，2005.

[11]Forlani G.Complete classification of raw LIDAR data and 3D reconstruction of building[J].Pattern Analysis and Applications，2006，8（4）：357-374.

[12]Doneus M，Briese C.Digital terrain model ling for archaeological interpretation within forested areas using full waveform laserscanning[C]//Proceeding of the 7th International Symposium on Virtual Reality，A rchaeology and Cultural Heritage VAST，2006.

[13]Lodha S，K reps E J，Helmbold D，et al.Aerial lidar data classification using Support Vector Machines（SVM）[C]//Proceedings of the 3rd International Symposium on 3D Data Processing，Visualization，and Transmission（3DPVT’06），2006：567-574.

[14]Lu W L，Murphy K P，Little J J，et al.A hybrid conditional random field for estimating the underlying ground surface from airborne lidar data[J].IEEE Transactions on Geoscience and Remote Sensing，2009，47（8）：2913-2922.

[15]Shan J，Sampath A.Urban DEM generation from raw Lidar data：a labeling algorithm and its performance[J].Photogrammetric Engineering&Remote Sensing，2005，71（2）：217-226.

[16]M eng X L，Wang L.A multi-directional ground filtering algorithm for airborne lidar[J].ISPRS Journal of Photogrammetry and Remote Sensing，2009，64：117-124.

[17]Sithole G.Experimental comparison of filtering algorithm s for bare-earth extraction from airborne laser scanning point clouds[J].ISPRS Journal of Photogrammetry and Remote Sensing，2004，59（1/2）：85-101.

[18]Poullis C，You S.Delineation and geometric modeling of road networks[J].ISPRS Journal of Photogrammetry and Remote Sensing，2010，65：165-181.