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

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

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

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

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

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

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

隨著計(jì)算機(jī)技術(shù)、光電技術(shù)與傳感器技術(shù)的發(fā)展，獲得高精度、高空間分辨率的三維地形信息不再像過去那樣遙不可及，地面三維激光掃描(Terrestrial Laser Scanner, TLS)技術(shù)憑借其受觀測環(huán)境影響小，可以實(shí)現(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è)工作強(qiáng)度的情況下，將主要的作業(yè)內(nèi)容放在內(nèi)業(yè)數(shù)據(jù)處理上，數(shù)據(jù)處理流程主要包括：點(diǎn)云濾波、DEM生成、DEM配準(zhǔn)。為了得到準(zhǔn)確的形變監(jiān)測結(jié)果，在對點(diǎn)云數(shù)據(jù)進(jìn)行后處理分析前必須先濾除不屬于地形信息的點(diǎn)( 如植被、電線、建筑物等)，所以點(diǎn)云濾波一直是影響三維激光掃描數(shù)據(jù)精度的重要因素，也是近年來學(xué)者們研究的熱點(diǎn)問題。

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

1 點(diǎn)云聚類濾波算法原理

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

1.1 建立模型

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

二維泊松過程理論認(rèn)為，在二維空間中，

(1)

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

(2)

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

(3)

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

(4)

其中，di是DK的計(jì)算值。

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

(5)

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

1.2 模型參數(shù)估計(jì)

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

(6)

M算法估計(jì)λf，λn，p為

(7)

(8)

(9)

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

1.3 三維點(diǎn)云聚類濾波算法原理

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

n維超球體的體積為

(10)

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

在三維空間中，

(11)

2 工程實(shí)例

2.1 研究區(qū)域概況

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

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

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

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

圖1 雞冠嶺危巖體

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

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

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

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

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

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

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

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

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

(12)

圖4 噪聲探測

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

3 結(jié) 論

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

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

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

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

[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] 閆利,張疑.基于法向量模糊聚類的道路面點(diǎn)云數(shù)據(jù)濾波[J].武漢大學(xué)學(xué)報(bào)(信息科學(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] 隋立春,張熠斌,柳艷,等.基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點(diǎn)云數(shù)據(jù)濾波[J].測繪學(xué)報(bào),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ī)制研究[J].巖石力學(xué)與工程學(xué)報(bào), 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三維點(diǎn)云聚類濾波算法應(yīng)用研究[J].測繪工程，2017，26(4)：64-69.

10.19349/j.cnki.issn1006-7949.2017.04.012

2016-05-19

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

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

P221

A

1006-7949(2017)04-0064-06