TLS三维点云聚类滤波算法应用研究

段利媛，赵超英，张 勤，牛玉芬,张静

(长安大学 地质工程与测绘学院，陕西 西安 710054)

TLS三维点云聚类滤波算法应用研究

段利媛，赵超英，张 勤，牛玉芬,张静

(长安大学 地质工程与测绘学院，陕西 西安 710054)

现有地面三维激光扫描点云数据滤波算法较少，针对地形复杂区域的点云滤波效果更是不甚理想，因此对二维聚类算法进行改进，提出三维点云聚类滤波算法，并对其在地形复杂区域的TLS数据滤波中的应用进行研究。以重庆鸡冠岭危岩体的TLS数据为例，分别采用曲率平滑滤波方法和文中提出的点云聚类滤波方法处理，并对两种方法处理过的数据进行形变量计算和分析。实验证明，针对植被覆盖茂密、地形复杂的山体，该方法的点云滤波效果较好，且处理速度有较大提升，能为点云后期形变量计算提供较好的基础。

TLS；滤波；三维聚类算法；曲率平滑滤波；形变量

随着计算机技术、光电技术与传感器技术的发展，获得高精度、高空间分辨率的三维地形信息不再像过去那样遥不可及，地面三维激光扫描(Terrestrial Laser Scanner, TLS)技术凭借其受观测环境影响小，可以实现非接触、快速、高精度获取研究区域逼真三维形态的优势迅速崛起，被广泛应用在地形测绘[1]、滑坡调查与变形监测[2]、危岩体位移监测[3]、地面沉降监测[4]等领域。三维激光扫描技术在降低了外业工作强度的情况下，将主要的作业内容放在内业数据处理上，数据处理流程主要包括：点云滤波、DEM生成、DEM配准。为了得到准确的形变监测结果，在对点云数据进行后处理分析前必须先滤除不属于地形信息的点( 如植被、电线、建筑物等)，所以点云滤波一直是影响三维激光扫描数据精度的重要因素，也是近年来学者们研究的热点问题。

目前针对TLS数据的点云滤波方法较少，主要有平滑滤波法、移动最小二乘曲面法[5]、双边滤波法等，这些方法原理简单、较易实现，计算速度快，却容易出现过度光顺的现象，造成细部特征损失，模型扭曲变形。针对机载、车载三维激光扫描(也称LiDAR，Light Detection And Ranging)的点云滤波方法较多，也相对成熟，TLS数据滤波也可以借鉴使用，主要有：坡度或法向滤波算法[6-7]、移动窗口法[8]、迭代线性最小二乘滤波[9]、形态学及其改进方法[10-11]等。这些算法依赖于扫描数据的空间位置信息，所处理点云的尺度较大，并且要求其空间坐标系为平面—高程系统，对自由坐标系下的地面三维激光扫描数据处理并不完全适用，还存在计算量大、计算速度慢，对于低矮草丛、灌木的处理较困难等问题，特别是在研究区域地形复杂的情况下。基于此，本文提出一种基于R语言的以点云距离为判别特征的三维点云聚类滤波算法，该算法在原有的二维聚类算法[12-14]上进行改进，有效避免二维聚类算法使用中进行点云投影时不可避免的点云空间结构信息的损失，充分利用TLS点云数据的高维度和高密度特点。为此，本文选取重庆武隆鸡冠岭危岩体为研究对象，进行点云滤波处理，与常用曲率平滑滤波方法进行对比分析，并对滤波处理后的数据进行形变量计算与分析。该危岩体与地面呈近90°分布，岩体表面风化严重，凹凸不平，并常年覆有茂密植被，大大增加该危岩体形变监测的难度，这些情况能够更好地验证本文所提出滤波算法对植被茂密的复杂地形点云数据滤波处理的效果。

1 点云聚类滤波算法原理

该算法主要包括建立模型和模型参数估计两部分。

1.1 建立模型

将研究区域的特征点与噪声点看作叠加在一起且分别服从齐次泊松分布[12]，首先计算每个数据点到其周围K个最近点的欧氏距离DK，由于高密度区点的DK值较低密度区小，所以特征点的DK值比噪声点的DK值小，由此做基于距离值的空间聚类及特征点提取。

二维泊松过程理论认为，在二维空间中，

(1)

只考虑一个齐次泊松过程时，可以得出距离DK的分布，对于x∈[0,)，

(2)

式(2)指假设每个数据点都有一个以该点为中心的半径为x的圆。如果DK比x大，这个圆里肯定有0,1,…,K-1个点，因此可以得到密度算式为

(3)

这是一个转换的伽马随机变量γ～Γ(K,λπ)，其中γ=(DK)2，所以DK～Γ1/2(K,λπ)，即一个广义伽马分布。在二维齐次泊松过程中，强度(过程速率)λ的极大似然估计量为

(4)

其中，di是DK的计算值。

当DK为混合分布时：将特征点和噪声点的分布模型做简单叠加，假设DK的分布(大约)是两个分布组合的泊松过程，DK的模型为

(5)

假设每个数据点都有一个δi∈{0,1}，如果第i个数据点在特征聚类群中，δi=1，反之(第i个点在噪声点聚类群中)则δi=0。因此每个数据点都有一个DK的计算值di和一个未知量δi，以此将数据点聚类到特征点群和噪声点群中。

1.2 模型参数估计

用期望最大化(Expectation Maximization，EM)算法估计DK模型参数λf，λn，p。EM算法可以更详细地分为E(Expectation)算法和M(Maximization)算法。E算法估计缺失值δi为

(6)

M算法估计λf，λn，p为

(7)

(8)

(9)

根据在混合数据点中每个点的观测值DK是否有更高的密度这一标准，再用上述估计值，就可以对点云数据进行聚类分析。

1.3 三维点云聚类滤波算法原理

在n维过程，可以用n维超球体的分布结果代替圆的分布。

n维超球体的体积为

(10)

然后用Vnxn代替πx2，重复上面的过程，得出结果：γ～Γ(k,λVn) ,其中γ=(DK)n。

在三维空间中，

(11)

2 工程实例

2.1 研究区域概况

鸡冠岭危岩体位于四川省涪陵区武隆县兴顺乡核桃村乌江边滩峡左岸。本文扫描区域为紧邻1994年鸡冠岭岩崩段的新生危岩体[15-16]，分布标高为540～740 m，长约300 m，宽约70 m，分布面积约 2.1×104m2，方量约为 45.06×104m3，如图1所示。乌江在此段切割强烈，地形陡峭，且该岩体地处桐麻湾背斜核部地层中[15,17]，由灰岩和页岩夹煤层组成， 具有上硬下软的地质结构特征。由于背斜核部受到强烈的褶皱作用，应力集中，导致上覆层状岩体呈现出弯曲变形的特征。且此段危岩下部岩溶十分发育，对上部岩体的稳定性非常不利，极易形成大型崩塌。

项目组使用徕卡HDS8800扫描仪于2012年10月、2013年8月对该区域进行扫描，扫描范围如图1所示，本文选取图1中两块矩形区域(下文称面积较小区域为区域一，面积较大区域为区域二)验证所设计滤波算法对植被茂密的复杂地形点云数据滤波处理效果。所选研究区域植被茂密(高大树木与低矮植被并存)，且山体表面凹凸不平，能有效验证滤波算法针对地形复杂多植被区域的滤波效果。

2.2 数据处理

对区域一和区域二的两期点云数据进行处理，主要步骤包括：滤波、DEM建模、模型配准及形变量计算。本文主要工作为滤波及形变量计算两步，其中滤波方法分为两种：曲率平滑滤波方法、三维点云聚类滤波方法。后期形变量使用Geomagic软件获取。

图1 鸡冠岭危岩体

2.2.1 滤波结果与分析

区域一和区域二的滤波结果分别如图2、图3所示。图中矩形区域为主要噪声点所在处，处理过程中主要比较矩形内的点云滤波效果。图2(a)为区域一原始点云，点云数量为200,278，图2(b)为经人工和曲率平滑滤波处理的点云，点云数量为184,338，从图2(b)可以看出，人工和曲率平滑滤波只是滤除了明显的噪声点(高大树木的树冠部分)，对于矮小植被和高大树木剩余的树干部分则无法进一步处理，在地形复杂的part1和part4两处尤其明显。图2(c)为经三维点云聚类算法处理的点云，K=8时滤波效果最好，此时点云数量为148,923。用EM算法估计区域一DK模型的拟合，结果如表1所示，得出模型参数的估计值λf=308.15，λn=71.41，p=0.94。图3(a)为区域二原始点云，点云数量为581,601，图3(b)为经人工和曲率平滑滤波处理的点云，点云数量为547,003，图中part3和part4两部分残余噪声点明显，其中part4部分皆是矮小的草木，且地形精细复杂，常用曲率平滑滤波方法对此几乎失效，图3(c)为经三维点云聚类算法处理的点云，点云数量为436,767，part4部分噪声点全部被滤除，且完整保存该处的山体特征。估计区域二DK模型的拟合，结果如表2所示，模型参数估计值λf=580.13，λn=76.47，p=0.93。

注：a：原始点云；b:人工、曲率平滑滤波方法处理点云;c:三维点云聚类算法处理点云(K=8)图2 点云滤波结果对比(区域一)

注：a:原始点云;b:人工、曲率平滑滤波处理点云;c:三维点云聚类算法处理点云(K=10)图3 点云滤波结果对比(区域二)

表1 EM估计模型拟合的迭代计算(区域一)

表2 EM估计模型拟合的迭代计算(区域二)

图4为两区域的噪声探测结果，经试验，K=8时区域一的滤波效果最好，K=10时区域二的滤波效果最好。从图4中可以看出，探测出的噪声点主要分布于矩形中，与噪声的实际分布一致，零星分布的其余红色点则是由于其扫描点云质量不好，太过稀疏引起的噪声点误判，所以外业扫描时得到质量较好的点云数据是非常关键的。

2.2.2 形变量结果与分析

对滤波后的点云进行DEM建模、模型配准后，对两期DEM模型进行3D误差分析，计算山体形变量，结果见图4，图中白色部分为两期数据都缺失部分，灰色部分为因一期数据缺失而无法进行比较的部分(2013年8月数据局部点云过于稀疏，质量不好)。从图可以看出：区域一使用曲率平滑滤波方法处理后的最大形变量为±0.070 m，平均形变量为-0.015～0.015 m，标准偏差(见式(12))为0.019 m；使用三维点云聚类算法处理后的最大形变量降为±0.040 m，平均形变量为-0.014～0.014 m，标准偏差为0.017 m。区域二使用曲率平滑滤波方法处理后的最大形变量为±0.2 m，平均形变量为-0.025～0.024 m，标准偏差为0.037 m，使用三维点云聚类算法处理后的最大形变量为±0.06 m，平均形变量为-0.014～0.014 m，标准偏差为0.018 m。两种处理方法的形变量结果中，最大形变量相差较大，平均形变量、标准偏差值相差较小，这是因为曲率平滑滤波无法将植被茂密区域的噪声点全部滤除，其中仍然存在较明显噪声点(见图2(b)、图3(b))，故图4(a)中该区域的形变量值较大，同时影响整个研究区域的平均形变量和标准偏差值，反映到图5中则表现为：两种方法处理后的形变量分布大体相同，只在矩形部分差异较大。从形变量结果中可以看出，三维点云聚类滤波算法针对地形复杂区域的TLS数据滤波效果较好，可以为后期形变量计算提供更好的基础，在TLS点云数据滤波处理中具有明显优势。模型配准时，由于2013年8月数据点云过于稀疏，导致该无标靶数据配准误差较大，为1.4 cm，对形变量结果影响较大，故针对复杂地形的无标靶点云数据的配准亦是今后努力的方向。

(12)

图4 噪声探测

注：(a-1)：曲率平滑滤波处理后的区域一形变量；(b-1)：三维点云聚类算法处理后的区域一形变量；(a-2)：曲率平滑滤波处理后的区域二形变量；(b-2)：三维点云聚类算法处理后的区域二形变量图5 研究区域形变量

3 结 论

针对现有滤波算法在实现地形复杂、植被茂密区域的TLS数据滤波处理方面的不足提出改进的三维点云聚类滤波算法，更充分地利用TLS数据的高维度、高密度特点。以鸡冠岭危岩体数据为例，通过对其进行滤波处理，并与常用曲率平滑滤波方法进行对比分析，证明三维点云聚类算法对多植被的复杂地形点云数据进行滤波处理的优越性，通过计算与分析研究区域的形变量，验证该算法对植被茂密的复杂地形点云数据滤波的适用性。在此基础上对该算法的优点及不足之处做出以下结论：

1)三维点云聚类算法能够对地形复杂、多植被区域的TLS数据进行快速滤波，且能较好保存区域的细部特征，滤波效果较好。

2)该算法中K值的选择没有可依据的准则，需按地形的复杂程度试验确定，但一般情况下K≥3。

3)该算法对外业扫描数据质量要求较高，点云过于稀疏时可能会引起部分误判。

[1] ARDIZZONE F, CARDINALI M, GALLI M, et al. Identification andMapping of Recent Rainfall-Induced Landslides Using Elevation Data Collected by Airborne Lidar[J]. Natural Hazards and Earth System Science, 2007, 7(6): 637-650.

[2] JABOYEDOFF M, OPPIKOFER T, ABELLN A, et al. Use of LIDAR inLandslide Investigations: AReview[J]. Natural hazards, 2012, 61(1): 5-28.

[4] GORDON S, LICHTI D, STEWART M. Application of aHigh-Resolution, Ground-Based Laser Scanner for Deformation Measurements[C]//Proceedings of 10th International FIG Symposium on Deformation Measurements, Orange, California, USA. 2001: 19-22.

[5] ALEXA M, BEHR J, COHEN -OR D, et al. Computing and rendering point set surfaces[J]. IEEE Transactions on visualization and computer graphics, 2003, 9(1): 3-15.

[6] SITHOLE G,VOSSELMAN G.ExperimentalComparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds[J].ISPRS Journal of Photogrammetry and Remote Sensing,2004(59):85 -101.

[7] 闫利,张疑.基于法向量模糊聚类的道路面点云数据滤波[J].武汉大学学报(信息科学版),2007,32(12):1119-1122.

[8] PETZOLD B,REISS P,STOSSEL W.LaserScanning Surveying and Mapping Gencies are Using a New Technique for the Derivation of Digital Terrain Models[J].ISPRS Journal of Hotogrammetry and Remote Sensing,1999(54):95 -104.

[9] PFEIFER N, KÖSTLI A, KRAUS K.Interpolation and Filtering of Laser Scanner Data-Implementation and First Results[J]. International Archives of Photogrammetry and Remote Sensing, 1998, 32: 153-159.

[10] KILIAN J,HAALA N.ENGLICH M.Capture and Evaluation of Airborne Laser Scanner Data[J].International Archives of Photogrammetry and Remote Sensing,1996,31(B3):383 -388.

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

[12] BYERS S, RAFTERY A E. Nearest-Neighbor Clutter Removal for Estimating Features in Spatial Point Processes[J]. Journal of the American Statistical Association, 1998, 93(442): 577-584.

[13] TONINI M, ABELLAN A.Rock fall Detection from Terrestrial LiDARPoint Clouds: A Clustering Approach Using R[J]. Journal of Spatial Information Science, 2014, 2014(8): 95-110.

[14] BADDELEY A J, TURNER R. Spatstat: AnR Package for Analyzing Spatial Point Pattens[J]. 2005,12(6):1-42.

[15] 李玉生, 谭开鸥, 王显华. 乌江鸡冠岭岩崩特征及成因[J].中国地质,1994(7): 26-27.

[16] 陈卫. 基于 TLS 技术的危岩体变形监测方法研究[D].西安：长安大学, 2014.

[17] 李滨,王国章,冯振,等. 地下采空诱发陡倾层状岩质斜坡失稳机制研究[J].岩石力学与工程学报, 2015, 34(6):1148-1162.

[责任编辑：张德福]

Research on TLS 3D point clouds cluster algorithm

DUAN Liyuan,ZHAO Chaoying,ZHANG Qin,NIU Yufen,ZHANG Jing

(School of Geological Engineering and Surveying, Chang’an University,Xi’an 710054, China)

The current algorithm of TLS point clouds filtering is less, and the effect is far from ideal when the terrain is complex, so a 3D point clouds cluster algorithm is proposed based on the 2D cluster algorithm to solve this problem in this paper. Taking the TLS data of Jiguanling dangerous rock mass as the example, the curvature smoothing filtering method and the 3D point clouds cluster algorithm are used to process the data respectively, thus calculating and analyzing their deformations. The experiment shows that the method proposed in this paper gets a better result, has a faster processing speed and provides a good foundation for the deformation calculation when processing data of complex terrain with dense vegetation.

TLS filtering; 3D cluster algorithm; curvature smoothing filtering; deformation

引用著录：段利媛，赵超英，张 勤，等.TLS三维点云聚类滤波算法应用研究[J].测绘工程，2017，26(4)：64-69.

10.19349/j.cnki.issn1006-7949.2017.04.012

2016-05-19

国家自然科学基金资助项目(41372375)

段利媛(1993-)，女，硕士研究生.

P221

A

1006-7949(2017)04-0064-06