三维激光扫描地形数据获取处理

三维激光扫描地形数据获取处理
惠振阳;吴北平;徐鹏;国继鑫
【摘要】利用三维激光扫描测量技术测得的点云数据获取数字地面模型的关键之处在于将地面点与非地面点进行分离.提出了一种基于高程直方图进行多阈值分割的方法.首先对点云进行离散化,然后进行多阈值分割分层,再对各层中的地面点进行连通,最后利用趋势面拟合法再次对各层地面点进行滤波,最终即可获取纯净的地面点.经试验验证,该方法有效可行.
【期刊名称】《科学技术与工程》
【年(卷),期】2014(014)018
【总页数】6页(P1-5,24)
【关键词】点云数据;高程直方图;单阈值分割;多阈值分割;连通;趋势面拟合
【作者】惠振阳;吴北平;徐鹏;国继鑫
【作者单位】中国地质大学(武汉)信息工程学院,武汉430074;中国地质大学(武汉)信息工程学院,武汉430074;中国地质大学(武汉)信息工程学院,武汉430074;长安大学地质工程与测绘学院,西安710000
【正文语种】中文
【中图分类】P231.5
三维激光扫描测量技术(light detection and ranging，LIDAR)是一种自动立体扫描技术，它是测绘领域继GPS测量技术之后又一项重大突破。

它采用非接触主动
测量方式能够对物体进行立体扫描快速获取三维坐标，突破了传统测量(如全站仪
测量，RTK测量等)单点测量技术的限制，具有速度快、精度高的特点，可以极大
地降低成本，节约时间。

现经常用于获取高精度、高分辨率的数字地面模型。

利用三维激光扫描技术获取的点云数据建立数字地面模型很关键的一步就是要把地面点和非地面点区分开来，从而剔除非地面点获取“纯净”的地面点，这个过程称之为点云滤波。

现如今，许多文献都对点云滤波做出了研究。

主要有以下几类:以
形态学为基础的滤波方法[1—4]，以拟合为基础的滤波方法[5]，以坡度为基础的
滤波方法[6—10]，以扫描线为基础的滤波方法[11—13]。

每种方法都有其局限性，如以形态学为基础的滤波方法滤波结果过分依赖移动窗口大小的选择，对高程变化较大的区域滤波效果不好;以坡度为基础的滤波方法滤波结果则过度依赖坡度阈值
的选择，而且对每个点的K临近进行查询也相当耗时，计算量过大。

本文所提出的方法是基于高程统计来进行点云分类。

首先对点云高程进行离散化，然后按照最大类间方差法(OSTU)获取最佳分割阈值，继而对点云数据进行阈值分割，从而获得地面点云数据。

1 点云高程离散化
由点云数据可获取其最大高程值Hmax以及最小高程值Hmin，通过设定离散化
等级L，可以得到其离散化高程宽度dh，
如果某一点的高程为Hi，则该点所在离散化高程等级为Li，
然后对点云数据按照离散化等级依次进行统计，便可以得到点云的高程统计直方图。

获取高程统计直方图的目的在于判断是否要进行多阈值分割。

通过高程统计直方图可以清楚地看到点云有几类不同高度的目标，如果目标多于两类则需要进行多阈值逐步分割。

2 单阈值分割法
通过观察高程统计直方图，如果不同高度的目标个数不多于两个，便可以利用单阈值分割法将点云分为两类。

进行单阈值分割时可采用最大类间方差法(OTSU)。

最大类间方差法是1979年由日本学者大津(Nobuyuki Otsu)提出的，它是在判决分析最小二乘原理的基础上进一步推导得出的自动选择阈值的二值化方法[14]，所以也简称为大津法(OTSU)。

这种方法主要用于图像分割，它是根据图像灰度特性，将图像分为目标和背景两类。

两类的类间方差越大，说明类的差别越大，错分的概率也就越小。

因此，当两类的类间方差最大时，错分的概率也就最小，此时的分割阈值也就是最佳分割阈值。

借助最大类间分割法在图像分割中的思想，进行点云分割便可以这样做:
设定一阈值将点云分为两类，然后不断调整，当两类点云的方差最大时，此时的阈值就作为分割处理的最佳阈值。

其数学表述如下:
设定阈值K，将点云L分为两类L1、L2。

分别计算L1类的概率ω1和均值μ1以及L2类的概率ω2和均值μ2，则两类点云的数学期望为:
两类的类间方差为:
极大化σ2(k)的过程就是自动确定阈值的过程，最佳阈值Kh为
具体实现步骤如下:
步骤一，初始化分类阈值K(K为点云离散化等级第Li级所对应的高程值)。

步骤二，阈值K将离散化后的点云分为两类，设为类1和类2。

按照式(4)计算此
时类1与类2间的类间方差。

步骤三，遍历点云离散化等级Li(0≤Li≤Lmax)，依次赋予阈值，分别计算该阈值
下的类间方差。

取σ2(K)最大时的阈值K作为整个高程直方图的最佳分割阈值。

此时，便可以将小于阈值K的点归类为地面点，大于阈值K的点归类为非地面点，直接将非地面点进行删除即可。

3 多阈值分割法
如果高程统计直方图表现出多峰现象，即不同高度的目标有好几个，此时如果继续用单阈值分割法就会造成误判。

点云的多阈值分割法可以看成是单阈值分割法的推广，它是一个迭代的过程，对点云逐步进行分割。

其流程图如图1所示。

此流程算法的关键之处在于直方图多峰现象的判断。

所谓的多峰现象，也就是直方图中有两个以上的独立峰。

经试验验证可采用下列三组条件进行独立峰的判断:
1)独立峰要具有一定的宽度，即相邻两波谷的间距≥l。

2)独立峰要具有一定的面积，即相邻两波谷和波峰所围成图形的面积≥s。

3)独立峰要具有一定的峰谷比，即波峰与其相邻两波谷的比值≥δ。

图1 点云分割流程Fig.1 Flow chart of threshold segmentation method
三阈值l、s、δ根据采集数据的实际情况进行设定，可多次设定进行对比，以达到判断多峰现象最佳的目的。

这是个逐步分割的过程，先分割出高程最低的点群并从原始点云中剔除，再在剩余点云中再次迭代分割出高程次低的点群，如此循环便可将原始点云按不同的阈值分割出不同的点群。

4 地面点连通规则
原始点云经过多阈值分割后会分成不同层的点群，一般默认为最下层也就是高程集体表现最低的点群为地面点。

除此之外，地面点还应包括其他各层(尤其是与最下
层相邻的层)点群中由于地面起伏而凸起的点。

将分割后的点云分层进行观察，可以发现每层的点群都会表现出连通区域或大或小的“点云块儿”。

这些块儿状点云主要是由于建筑物、树木以及地面起伏造成的。

显然，应该剔除由于建筑物或者树木造成的块儿状点云，保留由于地面起伏造成的块儿状点云。

如何进行区分，要遵循以下规则:
1)连通区域较小且分散比较开的块儿状点云直接进行剔除，此类块儿状点云是由于一般建筑物或者树木造成的。

2)连通区域较大的块儿状点云可能由两种情况造成，一种是大面积的建筑物(比如
房顶)，另一种则是地面起伏。

前一种情况要剔除，后一种情况要保留。

如何进行
区分，主要依赖于这些块儿状点云与其周边点云高程的落差情况。

如果块儿状点云与其周边点云高差相差较大，则说明该块儿状点云是由于大面积建筑物造成的应该剔除。

反之，则应该保留。

由图2和图3可以看出其区别。

图2 地面起伏形成的块儿状点云Fig.2 Point cloud formed by ground ups and downs in rocks
图3 大面积建筑物形成的块儿状点云Fig.3 Point cloud formed by large area
of building
5 非地面点滤波处理
按照上述规则分别对各层数据进行分析对比，剔除非地面点，将剩余点进行连通即可得到想要的地面点云。

但是此时的点云数据还存在误判性，比如一些低矮的植物，由于其分层后的点群与其周边其他点的高程差也不大，所以很容易将此种类型的点云误分到地面点中去。

为减小这种误判性，可采用趋势面拟合法对连通后的点云进行非地面点滤波处理。

趋势面拟合法的基本思想是利用底层已确定的点云数据建立二次曲面，然后对其余各层中连通的地面点进行高程拟合，最后比较真实值与拟合值的高程差，将高程差超过阈值的点进行删除。

一般地设定二次曲面函数为:
式(6)中，x、y、z为点的三维坐标，a0、a1、a2、a3、a4、a5为待求曲面函数的系数。

首先，通过n个已知的地面点(n≥6)对式(6)建立误差方程，然后按照最小二乘原理求解系数 ai(i=0，1，2，3，4，5)[15]，最后利用求解出来的二次曲面分别对其余各层的地面点进行高程拟合，并计算拟合值与真实值zj的差值Δ，若Δ大于阈值ε则将该点删除掉，否则保留。

6 试验分析
本次试验所采用的仪器是奥地利Riegl公司生产的VZ—400型三维激光扫描仪，其仪器参数如表1所示。

本次试验区域为某教学楼前的一片空地，为便于分析研究选取其中一块儿具有代表性的区域(约为20 m×10 m)，如图4、图5中。

表1 VZ—400型扫描仪仪器参数Table 1 VZ—400 scanner instrument parameter型号 VZ—400类型脉冲测量距离/m 1～500激光级别ⅠFOV(视场角) 360°×150°单点测距精度/mm±5角度分辨率/(°) 0.005最大扫描速率/(点·s-1)12500
图4 扫描区域点云示意图(a)Fig.4 Point cloud of scanning area(a)
图5 扫描区域点云示意图(b)Fig.5 Point cloud of scanning area(b)
从图中可以清楚地看到地面上除了高低不一的独立树还有一些土坎(图中点云中有些空白区域正是由于土坎对扫描脉冲的遮挡所造成的)，如前文所述独立树是要剔除的而土坎则是要保留的。

具体实现步骤如下:
1)对点云高程进行离散化，统计高程直方图。

此区域扫描的点云数量为10091个点，高程范围是22.831～29.213 m。

将离散化等级设置为500级，其高程离散化直方图如图6所示。

表2 分割阈值与点云数量Table 2 Threshold segmentation and the number
of point cloud分层阈值该层点云数量22.831～23.300 m 592723.300～
25.950 m 188725.950～29.213 m 2277
图6 高程离散化直方图Fig.6 Elevation discretization histogram
2)对离散化的点云进行多阈值OTSU分割，计算得出的最佳分割阈值如表2所示。

按照分割阈值对点云进行分层处理，得到的结果如图7～图12所示。

图7 22.831～23.3层的点云数据(a)Fig.7 Point cloud between layer 22.831 and layer 23.3(a)
图8 22.831～23.3层的点云数据(b)Fig.8 Point cloud between layer 22.831 and layer 23.3(b)
3)对除底层外其他各层的点云进行地面点连通判断。

底层点云默认为地面点，第三层由于点云的连通区域较小且该层平均高程与底层点云相差较大故可直接删去。

第二层点云虽然连通区域也不大但由于其极有可能包含地面起伏的点(如土坎等)故需单独进行判断。

分别统计第二层点云以及该层周边点云(这些点存在于第一层中)的高程，计算高差。

设定阈值为1 m，若高差超过1 m，则视为非地面点，否则保留为地面点。

最后，连通底层和第二层中的地面点，结果如图13、图14所示。

图9 23.3～25.95层的点云数据(a)Fig.9 Point cloud between layer 23.3 and layer 25.95(a)
图10 23.3～25.95层的点云数据(b)Fig.10 point cloud between layer 23.3 and layer 25.95(b)
图11 25.95～29.213层的点云数据(a)Fig.11 Point cloud between layer 25.95 and layer 29.213(a)
图12 25.95～29.213层的点云数据(b)Fig.12 Point cloud between layer 25.95
and layer 29.213(b)
4)对上步得到的地面点利用趋势面拟合法再次进行非地面点滤波，可得到最终的地面点云如图15、图16所示。

综上可得每步操作后的点云数据对比如表3所示。

图13 连通后的点云(a)Fig.13 Connected point cloud(a)
图14 连通后的点云(b)Fig.14 Connected point cloud(b)
图15 最终的地面点云(a)Fig.15 Final ground point cloud(a)
图16 最终地面点云(b)Fig.16 Final ground point cloud(b)
7 结论
经试验验证本文方法能够较好地将地面点与非地面点进行分离，并经二次判断能够较好地保留地面起伏点，减小非地面点的误判性。

但本文方法也有一定的弊端，如本文方法比较复杂、计算量大，如果点云数量过大，计算时间会很长。

而且本文进行多阈值分割时需进行多峰判断，多峰判断三条件的阈值选取也直接影响了分割的结果，如果选取不当也会出现点云分割错误的现象，直接影响地面点云滤波的最终结果。

表3 点云数据对比Table 3 Point cloud data contrast类别数量比例10091 100底层地面点 5927 58.74连通后的地面点 6047 59.92趋势面滤波后地面点/%原始点云6010 59.56
参考文献
【相关文献】
1 Roggero M.Airborne laser scanning:clustering in raw data.International Archives of Photogrammetry and Remote Sensing，2001;34(3/W4):227—232
2 Masaharu H，Ohtsubo K.A filtering method of airborne laser scanner data for complex terrain.International Archives of Photogrammetry，Remote Sensing and Spatial Information Science，2002;34(3/B):165—169
3 Zhang K，Chen S，Whitman D，et al.A progressive morphological filter for removing nonground measurements from airborne LiDAR data.IEEE Transactions on Geoscience and Remote Sensing，2003;41(4):872—882
4 Chen Q.Filtering airborne laser scanning data with morphological
methods.Photogrammetric Engineering and Remote Sensing，2007;73(2):175
5 张小红，刘经南.机载激光扫描测高数据滤波.测绘科学，2004;29(6):50—53 Zhang X H，Liu J N.Airborne laser scanning altimetry data filtering.Science of Surveying and Mapping，2004;29(6):50—53
6 Vosselman G.Slope based filtering of laser altimetry data.International Archives of Photogrammetry and Remote Sensing，2000;33(3/W4):203—210
7 Sithole G.Filtering of laser altimetry data using a slope adaptive filter.International Archives of Photogrammetry and Remote Sensing，2001;34(3/W4):203—210
8 Sampath A.Urban modeling based on segmentation and regularization of airborne lidar point clouds.International Archives of Photogrammetry and Remote Sensing，
2004;35:923—941
9 Filin S，Pfeifer N.Segmentation of airborne laser scanning data using a slope adaptive neighborhood.ISPRS Journal of Photogrammetry and Remote Sensing，2006;60:71—80 10 管海燕，张剑清，邓非，等.基于扫描线的城区机载激光扫描数据滤波算法研究.测绘通报，2007;(12):9—13 Guan H Y，Zhang J Q，Deng F，et al.On the filtering algorithm for airborne laser scanning data of urban area based on scanning line.Bulletin of Surveying and Mapping，2007;(12):9—13
11 Sithole G，Vosselman G.Filtering of airborne laser scanner data based on segmented point clouds.International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，2005;36(3):12—14
12 Shan J，Sampath A.Urban DEM generation from raw LiDAR data:a labeling algorithm and its performance.Photogrammetric Engineering and Remote Sensing，
2005;71(2):217—226
13 Meng X，Wang L，Silvancardenas J，et al.A multidirectional ground filtering algorithm for airborne LiDAR.ISPRS Journal of Photogrammetry and Remote Sensing，2009;64(1):117—124
14 Otsu N.A threshold selection method from gray-level histogram.IEEE Transactions on Systems Man Cybernetics，1979;9(1):62—66
15 李丽丹.基于MATLAB的离散数据最小二乘拟合.辽宁工程技术大学学报(自然科学版)，
2011;30:202—204Li L D.Least squares fitting of discrete data based on MATLAB.Journal of Liaoning Technical University(Natural Science)，2011;30:202—204。