段利媛，趙超英，張 勤，牛玉芬,張靜

(長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安 710054)

TLS三維點云聚類濾波算法應(yīng)用研究

段利媛，趙超英，張 勤，牛玉芬,張靜

(長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安 710054)

現(xiàn)有地面三維激光掃描點云數(shù)據(jù)濾波算法較少，針對地形復(fù)雜區(qū)域的點云濾波效果更是不甚理想，因此對二維聚類算法進行改進，提出三維點云聚類濾波算法，并對其在地形復(fù)雜區(qū)域的TLS數(shù)據(jù)濾波中的應(yīng)用進行研究。以重慶雞冠嶺危巖體的TLS數(shù)據(jù)為例，分別采用曲率平滑濾波方法和文中提出的點云聚類濾波方法處理，并對兩種方法處理過的數(shù)據(jù)進行形變量計算和分析。實驗證明，針對植被覆蓋茂密、地形復(fù)雜的山體，該方法的點云濾波效果較好，且處理速度有較大提升，能為點云后期形變量計算提供較好的基礎(chǔ)。

TLS；濾波；三維聚類算法；曲率平滑濾波；形變量

隨著計算機技術(shù)、光電技術(shù)與傳感器技術(shù)的發(fā)展，獲得高精度、高空間分辨率的三維地形信息不再像過去那樣遙不可及，地面三維激光掃描(Terrestrial Laser Scanner, TLS)技術(shù)憑借其受觀測環(huán)境影響小，可以實現(xiàn)非接觸、快速、高精度獲取研究區(qū)域逼真三維形態(tài)的優(yōu)勢迅速崛起，被廣泛應(yīng)用在地形測繪[1]、滑坡調(diào)查與變形監(jiān)測[2]、危巖體位移監(jiān)測[3]、地面沉降監(jiān)測[4]等領(lǐng)域。三維激光掃描技術(shù)在降低了外業(yè)工作強度的情況下，將主要的作業(yè)內(nèi)容放在內(nèi)業(yè)數(shù)據(jù)處理上，數(shù)據(jù)處理流程主要包括：點云濾波、DEM生成、DEM配準。為了得到準確的形變監(jiān)測結(jié)果，在對點云數(shù)據(jù)進行后處理分析前必須先濾除不屬于地形信息的點( 如植被、電線、建筑物等)，所以點云濾波一直是影響三維激光掃描數(shù)據(jù)精度的重要因素，也是近年來學(xué)者們研究的熱點問題。

目前針對TLS數(shù)據(jù)的點云濾波方法較少，主要有平滑濾波法、移動最小二乘曲面法[5]、雙邊濾波法等，這些方法原理簡單、較易實現(xiàn)，計算速度快，卻容易出現(xiàn)過度光順的現(xiàn)象，造成細部特征損失，模型扭曲變形。針對機載、車載三維激光掃描(也稱LiDAR，Light Detection And Ranging)的點云濾波方法較多，也相對成熟，TLS數(shù)據(jù)濾波也可以借鑒使用，主要有：坡度或法向濾波算法[6-7]、移動窗口法[8]、迭代線性最小二乘濾波[9]、形態(tài)學(xué)及其改進方法[10-11]等。這些算法依賴于掃描數(shù)據(jù)的空間位置信息，所處理點云的尺度較大，并且要求其空間坐標系為平面—高程系統(tǒng)，對自由坐標系下的地面三維激光掃描數(shù)據(jù)處理并不完全適用，還存在計算量大、計算速度慢，對于低矮草叢、灌木的處理較困難等問題，特別是在研究區(qū)域地形復(fù)雜的情況下?；诖?，本文提出一種基于R語言的以點云距離為判別特征的三維點云聚類濾波算法，該算法在原有的二維聚類算法[12-14]上進行改進，有效避免二維聚類算法使用中進行點云投影時不可避免的點云空間結(jié)構(gòu)信息的損失，充分利用TLS點云數(shù)據(jù)的高維度和高密度特點。為此，本文選取重慶武隆雞冠嶺危巖體為研究對象，進行點云濾波處理，與常用曲率平滑濾波方法進行對比分析，并對濾波處理后的數(shù)據(jù)進行形變量計算與分析。該危巖體與地面呈近90°分布，巖體表面風(fēng)化嚴重，凹凸不平，并常年覆有茂密植被，大大增加該危巖體形變監(jiān)測的難度，這些情況能夠更好地驗證本文所提出濾波算法對植被茂密的復(fù)雜地形點云數(shù)據(jù)濾波處理的效果。

1 點云聚類濾波算法原理

該算法主要包括建立模型和模型參數(shù)估計兩部分。

1.1 建立模型

將研究區(qū)域的特征點與噪聲點看作疊加在一起且分別服從齊次泊松分布[12]，首先計算每個數(shù)據(jù)點到其周圍K個最近點的歐氏距離DK，由于高密度區(qū)點的DK值較低密度區(qū)小，所以特征點的DK值比噪聲點的DK值小，由此做基于距離值的空間聚類及特征點提取。

二維泊松過程理論認為，在二維空間中，

(1)

只考慮一個齊次泊松過程時，可以得出距離DK的分布，對于x∈[0,)，

(2)

式(2)指假設(shè)每個數(shù)據(jù)點都有一個以該點為中心的半徑為x的圓。如果DK比x大，這個圓里肯定有0,1,…,K-1個點，因此可以得到密度算式為

(3)

這是一個轉(zhuǎn)換的伽馬隨機變量γ～Γ(K,λπ)，其中γ=(DK)2，所以DK～Γ1/2(K,λπ)，即一個廣義伽馬分布。在二維齊次泊松過程中，強度(過程速率)λ的極大似然估計量為

(4)

其中，di是DK的計算值。

當(dāng)DK為混合分布時：將特征點和噪聲點的分布模型做簡單疊加，假設(shè)DK的分布(大約)是兩個分布組合的泊松過程，DK的模型為

(5)

假設(shè)每個數(shù)據(jù)點都有一個δi∈{0,1}，如果第i個數(shù)據(jù)點在特征聚類群中，δi=1，反之(第i個點在噪聲點聚類群中)則δi=0。因此每個數(shù)據(jù)點都有一個DK的計算值di和一個未知量δi，以此將數(shù)據(jù)點聚類到特征點群和噪聲點群中。

1.2 模型參數(shù)估計

用期望最大化(Expectation Maximization，EM)算法估計DK模型參數(shù)λf，λn，p。EM算法可以更詳細地分為E(Expectation)算法和M(Maximization)算法。E算法估計缺失值δi為

(6)

M算法估計λf，λn，p為

(7)

(8)

(9)

根據(jù)在混合數(shù)據(jù)點中每個點的觀測值DK是否有更高的密度這一標準，再用上述估計值，就可以對點云數(shù)據(jù)進行聚類分析。

1.3 三維點云聚類濾波算法原理

在n維過程，可以用n維超球體的分布結(jié)果代替圓的分布。

n維超球體的體積為

(10)

然后用Vnxn代替πx2，重復(fù)上面的過程，得出結(jié)果：γ～Γ(k,λVn) ,其中γ=(DK)n。

在三維空間中，

(11)

2 工程實例

2.1 研究區(qū)域概況

雞冠嶺危巖體位于四川省涪陵區(qū)武隆縣興順鄉(xiāng)核桃村烏江邊灘峽左岸。本文掃描區(qū)域為緊鄰1994年雞冠嶺巖崩段的新生危巖體[15-16]，分布標高為540～740 m，長約300 m，寬約70 m，分布面積約 2.1×104m2，方量約為 45.06×104m3，如圖1所示。烏江在此段切割強烈，地形陡峭，且該巖體地處桐麻灣背斜核部地層中[15,17]，由灰?guī)r和頁巖夾煤層組成， 具有上硬下軟的地質(zhì)結(jié)構(gòu)特征。由于背斜核部受到強烈的褶皺作用，應(yīng)力集中，導(dǎo)致上覆層狀巖體呈現(xiàn)出彎曲變形的特征。且此段危巖下部巖溶十分發(fā)育，對上部巖體的穩(wěn)定性非常不利，極易形成大型崩塌。

項目組使用徠卡HDS8800掃描儀于2012年10月、2013年8月對該區(qū)域進行掃描，掃描范圍如圖1所示，本文選取圖1中兩塊矩形區(qū)域(下文稱面積較小區(qū)域為區(qū)域一，面積較大區(qū)域為區(qū)域二)驗證所設(shè)計濾波算法對植被茂密的復(fù)雜地形點云數(shù)據(jù)濾波處理效果。所選研究區(qū)域植被茂密(高大樹木與低矮植被并存)，且山體表面凹凸不平，能有效驗證濾波算法針對地形復(fù)雜多植被區(qū)域的濾波效果。

2.2 數(shù)據(jù)處理

對區(qū)域一和區(qū)域二的兩期點云數(shù)據(jù)進行處理，主要步驟包括：濾波、DEM建模、模型配準及形變量計算。本文主要工作為濾波及形變量計算兩步，其中濾波方法分為兩種：曲率平滑濾波方法、三維點云聚類濾波方法。后期形變量使用Geomagic軟件獲取。

圖1 雞冠嶺危巖體

2.2.1 濾波結(jié)果與分析

區(qū)域一和區(qū)域二的濾波結(jié)果分別如圖2、圖3所示。圖中矩形區(qū)域為主要噪聲點所在處，處理過程中主要比較矩形內(nèi)的點云濾波效果。圖2(a)為區(qū)域一原始點云，點云數(shù)量為200,278，圖2(b)為經(jīng)人工和曲率平滑濾波處理的點云，點云數(shù)量為184,338，從圖2(b)可以看出，人工和曲率平滑濾波只是濾除了明顯的噪聲點(高大樹木的樹冠部分)，對于矮小植被和高大樹木剩余的樹干部分則無法進一步處理，在地形復(fù)雜的part1和part4兩處尤其明顯。圖2(c)為經(jīng)三維點云聚類算法處理的點云，K=8時濾波效果最好，此時點云數(shù)量為148,923。用EM算法估計區(qū)域一DK模型的擬合，結(jié)果如表1所示，得出模型參數(shù)的估計值λf=308.15，λn=71.41，p=0.94。圖3(a)為區(qū)域二原始點云，點云數(shù)量為581,601，圖3(b)為經(jīng)人工和曲率平滑濾波處理的點云，點云數(shù)量為547,003，圖中part3和part4兩部分殘余噪聲點明顯，其中part4部分皆是矮小的草木，且地形精細復(fù)雜，常用曲率平滑濾波方法對此幾乎失效，圖3(c)為經(jīng)三維點云聚類算法處理的點云，點云數(shù)量為436,767，part4部分噪聲點全部被濾除，且完整保存該處的山體特征。估計區(qū)域二DK模型的擬合，結(jié)果如表2所示，模型參數(shù)估計值λf=580.13，λn=76.47，p=0.93。

注：a：原始點云；b:人工、曲率平滑濾波方法處理點云;c:三維點云聚類算法處理點云(K=8)圖2 點云濾波結(jié)果對比(區(qū)域一)

注：a:原始點云;b:人工、曲率平滑濾波處理點云;c:三維點云聚類算法處理點云(K=10)圖3 點云濾波結(jié)果對比(區(qū)域二)

表1 EM估計模型擬合的迭代計算(區(qū)域一)

表2 EM估計模型擬合的迭代計算(區(qū)域二)

圖4為兩區(qū)域的噪聲探測結(jié)果，經(jīng)試驗，K=8時區(qū)域一的濾波效果最好，K=10時區(qū)域二的濾波效果最好。從圖4中可以看出，探測出的噪聲點主要分布于矩形中，與噪聲的實際分布一致，零星分布的其余紅色點則是由于其掃描點云質(zhì)量不好，太過稀疏引起的噪聲點誤判，所以外業(yè)掃描時得到質(zhì)量較好的點云數(shù)據(jù)是非常關(guān)鍵的。

2.2.2 形變量結(jié)果與分析

對濾波后的點云進行DEM建模、模型配準后，對兩期DEM模型進行3D誤差分析，計算山體形變量，結(jié)果見圖4，圖中白色部分為兩期數(shù)據(jù)都缺失部分，灰色部分為因一期數(shù)據(jù)缺失而無法進行比較的部分(2013年8月數(shù)據(jù)局部點云過于稀疏，質(zhì)量不好)。從圖可以看出：區(qū)域一使用曲率平滑濾波方法處理后的最大形變量為±0.070 m，平均形變量為-0.015～0.015 m，標準偏差(見式(12))為0.019 m；使用三維點云聚類算法處理后的最大形變量降為±0.040 m，平均形變量為-0.014～0.014 m，標準偏差為0.017 m。區(qū)域二使用曲率平滑濾波方法處理后的最大形變量為±0.2 m，平均形變量為-0.025～0.024 m，標準偏差為0.037 m，使用三維點云聚類算法處理后的最大形變量為±0.06 m，平均形變量為-0.014～0.014 m，標準偏差為0.018 m。兩種處理方法的形變量結(jié)果中，最大形變量相差較大，平均形變量、標準偏差值相差較小，這是因為曲率平滑濾波無法將植被茂密區(qū)域的噪聲點全部濾除，其中仍然存在較明顯噪聲點(見圖2(b)、圖3(b))，故圖4(a)中該區(qū)域的形變量值較大，同時影響整個研究區(qū)域的平均形變量和標準偏差值，反映到圖5中則表現(xiàn)為：兩種方法處理后的形變量分布大體相同，只在矩形部分差異較大。從形變量結(jié)果中可以看出，三維點云聚類濾波算法針對地形復(fù)雜區(qū)域的TLS數(shù)據(jù)濾波效果較好，可以為后期形變量計算提供更好的基礎(chǔ)，在TLS點云數(shù)據(jù)濾波處理中具有明顯優(yōu)勢。模型配準時，由于2013年8月數(shù)據(jù)點云過于稀疏，導(dǎo)致該無標靶數(shù)據(jù)配準誤差較大，為1.4 cm，對形變量結(jié)果影響較大，故針對復(fù)雜地形的無標靶點云數(shù)據(jù)的配準亦是今后努力的方向。

(12)

圖4 噪聲探測

注：(a-1)：曲率平滑濾波處理后的區(qū)域一形變量；(b-1)：三維點云聚類算法處理后的區(qū)域一形變量；(a-2)：曲率平滑濾波處理后的區(qū)域二形變量；(b-2)：三維點云聚類算法處理后的區(qū)域二形變量圖5 研究區(qū)域形變量

3 結(jié) 論

針對現(xiàn)有濾波算法在實現(xiàn)地形復(fù)雜、植被茂密區(qū)域的TLS數(shù)據(jù)濾波處理方面的不足提出改進的三維點云聚類濾波算法，更充分地利用TLS數(shù)據(jù)的高維度、高密度特點。以雞冠嶺危巖體數(shù)據(jù)為例，通過對其進行濾波處理，并與常用曲率平滑濾波方法進行對比分析，證明三維點云聚類算法對多植被的復(fù)雜地形點云數(shù)據(jù)進行濾波處理的優(yōu)越性，通過計算與分析研究區(qū)域的形變量，驗證該算法對植被茂密的復(fù)雜地形點云數(shù)據(jù)濾波的適用性。在此基礎(chǔ)上對該算法的優(yōu)點及不足之處做出以下結(jié)論：

1)三維點云聚類算法能夠?qū)Φ匦螐?fù)雜、多植被區(qū)域的TLS數(shù)據(jù)進行快速濾波，且能較好保存區(qū)域的細部特征，濾波效果較好。

2)該算法中K值的選擇沒有可依據(jù)的準則，需按地形的復(fù)雜程度試驗確定，但一般情況下K≥3。

3)該算法對外業(yè)掃描數(shù)據(jù)質(zhì)量要求較高，點云過于稀疏時可能會引起部分誤判。

[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] 閆利,張疑.基于法向量模糊聚類的道路面點云數(shù)據(jù)濾波[J].武漢大學(xué)學(xué)報(信息科學(xué)版),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] 隋立春,張熠斌,柳艷,等.基于改進的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點云數(shù)據(jù)濾波[J].測繪學(xué)報,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].中國地質(zhì),1994(7): 26-27.

[16] 陳衛(wèi). 基于 TLS 技術(shù)的危巖體變形監(jiān)測方法研究[D].西安：長安大學(xué), 2014.

[17] 李濱,王國章,馮振,等. 地下采空誘發(fā)陡傾層狀巖質(zhì)斜坡失穩(wěn)機制研究[J].巖石力學(xué)與工程學(xué)報, 2015, 34(6):1148-1162.

[責(zé)任編輯：張德福]

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三維點云聚類濾波算法應(yīng)用研究[J].測繪工程，2017，26(4)：64-69.

10.19349/j.cnki.issn1006-7949.2017.04.012

2016-05-19

國家自然科學(xué)基金資助項目(41372375)

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

P221

A

1006-7949(2017)04-0064-06