融合kmeans 聚類與Hausdorff 距離的點(diǎn)云精簡(jiǎn)算法改進(jìn)

彭海駒，嚴(yán)科文*，林 松，賴浩源，張澤鑫

（1. 惠州市自然資源規(guī)劃勘測(cè)院，廣東 惠州 516000）

隨著三維激光掃描技術(shù)的不斷發(fā)展，點(diǎn)云精簡(jiǎn)成為高效處理海量的點(diǎn)云數(shù)據(jù)重要技術(shù)手段。常用的點(diǎn)云精簡(jiǎn)方法有體素下采樣法[1]、曲率采樣法[2]、基于聚類的精簡(jiǎn)方法、隨機(jī)采樣法。相關(guān)研究人員對(duì)幾種方法分別進(jìn)行探討，歸納出其面臨的問(wèn)題有容易丟失特征點(diǎn)數(shù)據(jù)、精度較差、曲率估算精度較差問(wèn)題等[3-4]。針對(duì)上述問(wèn)題，本文對(duì)融合kmeans 聚類與Hausdorff距離的點(diǎn)云精簡(jiǎn)算法進(jìn)行改進(jìn)，基于其對(duì)特征區(qū)域和非特征區(qū)域采用不同特征點(diǎn)提取策略的原理，在非特征點(diǎn)區(qū)域仍然采用kmeans 聚類算法提取特征點(diǎn)，但采用手肘法確定聚類數(shù)，保證聚類精度。在特征區(qū)域采用維數(shù)特征Hausdorff 距離代替主曲率Hausdorff 距離提取特征點(diǎn)，避免了曲率的估算和在曲率值過(guò)小時(shí)設(shè)定Hausdorff 距離閾值，最后融合kmeans 聚類簇心與采用維數(shù)特征Hausdorff 距離提取的特征點(diǎn)實(shí)現(xiàn)數(shù)據(jù)精簡(jiǎn)。

1 融合kmeans與Hausdorff距離的點(diǎn)云精簡(jiǎn)算法

1.1 kmeans聚類算法

kmeans 聚類算法中的k代表聚類個(gè)數(shù)，means 代表取每一個(gè)聚類中數(shù)值的平均值作為該簇的中心，即每一個(gè)簇都用這個(gè)類的中心描述[5]。這種聚類算法容易實(shí)現(xiàn)，其具體實(shí)現(xiàn)流程如下：

1）確定簇的個(gè)數(shù)k。

2）計(jì)算各個(gè)采樣點(diǎn)到簇中心的距離，一般是歐氏距離。

3）根據(jù)新劃分的簇，更新“簇中心”。

4）重復(fù)步驟（2）和步驟（3），直到“簇中心”不再移動(dòng)。

1.2 Hausdorff距離

Hausdorff距離是一種描述兩組點(diǎn)集之間相似程度的量度[6-7]，若存在兩組數(shù)據(jù)點(diǎn)A={a1,a2,…,an} ，B={b1,b2,…,bn} ，則這兩組數(shù)據(jù)的單向Hausdorff 距離為：

1.3 算法描述

融合kmeans 與Hausdorff 距離的點(diǎn)云精簡(jiǎn)算法首先建立空間索引，然后基于局部點(diǎn)云數(shù)據(jù)擬合局部曲面，再通過(guò)擬合曲面估計(jì)采樣點(diǎn)的2 個(gè)主曲率值。對(duì)于某一采樣點(diǎn)pi和其鄰域內(nèi)某點(diǎn)qj的主曲率分別為k1、k2和k1'、k2'，則采樣點(diǎn)與其鄰域內(nèi)點(diǎn)曲率的Hausdorff距離計(jì)算公式如下：

因此，采樣點(diǎn)pi的Hausdorff距離可以表示為：

該方法通過(guò)設(shè)定Hausdorff距離閾值實(shí)現(xiàn)特征點(diǎn)云的提取，然后對(duì)于非特征區(qū)域采用kmean 算法聚類，將每個(gè)聚類中心作為非特征區(qū)域的特征點(diǎn)與通過(guò)Hausdorff距離閾值提取的特征點(diǎn)融合完成精簡(jiǎn)。

2 改進(jìn)的融合kmeans與Hausdorff距離的點(diǎn)云精簡(jiǎn)算法

2.1 手肘法

手肘法是一種確定數(shù)據(jù)聚類最優(yōu)聚類數(shù)的方法[8]，李健[4]等采用平均密度計(jì)算kmeans聚類數(shù)，并沒(méi)有考慮聚類精度，對(duì)于聚類結(jié)果通常采用SSE作為精度檢驗(yàn)指標(biāo)，其計(jì)算公式如式（8）所示。而手肘法的核心思想就是使得SSE 最小。對(duì)于kmeans 來(lái)說(shuō)，k值取得越大，則樣本的劃分更加精細(xì)，每個(gè)簇的聚合程度更高，則SSE越小。并且，當(dāng)k的取值小于真實(shí)聚類數(shù)時(shí)，SSE 會(huì)隨著k值的增大迅速下降；當(dāng)k達(dá)到真實(shí)聚類數(shù)時(shí)，SSE 的下降幅度會(huì)驟減。隨后繼續(xù)增大k值，SEE 會(huì)慢慢趨于平緩。SSE 和k值的關(guān)系圖與一個(gè)手肘相似，而這個(gè)肘部對(duì)應(yīng)的k值就是真實(shí)聚類數(shù)。

式中，Ci為第i個(gè)簇；p為Ci中的采樣點(diǎn)；mi為Ci的質(zhì)心。

以斯坦福大學(xué)兔子點(diǎn)云為例，其原始點(diǎn)云數(shù)據(jù)如圖1所示，共8 171個(gè)點(diǎn)。采用手肘法，以100個(gè)聚類數(shù)作為步距，計(jì)算聚類個(gè)數(shù)在100～8 000 的SSE 值，結(jié)果如圖2 所示。圖中隨著k值的不斷增大，SSE 不斷減小，且整個(gè)圖形符合手肘法的趨勢(shì)。分析圖2 可知，當(dāng)k小于1 000時(shí)，SSE的下降幅度較大；當(dāng)k大于1 000時(shí)，SSE開(kāi)始緩慢下降，故對(duì)于該兔子點(diǎn)云數(shù)據(jù)的聚類個(gè)數(shù)可取值為1 000，圖3 為k取1 000 時(shí)兔子點(diǎn)云聚類效果圖。

圖1 原始兔子點(diǎn)云

圖2 聚類個(gè)數(shù)與SSE間關(guān)系

圖3 k=1 000 聚類結(jié)果

2.2 維數(shù)特征

主成分分析法（principal component analysis，PCA）是一種數(shù)學(xué)降維方法[9]，能夠?qū)⒃瓉?lái)變量重新組合成一組新的互不相關(guān)的幾個(gè)綜合變量，從綜合變量中選取部分綜合變量能較完整的反映原來(lái)變量的統(tǒng)計(jì)信息。對(duì)于采樣點(diǎn)P，其鄰域點(diǎn)集為Pi=[xi yi zi]，則可以通過(guò)式（8）得到P點(diǎn)的協(xié)方差矩陣M。

式中，k為P點(diǎn)鄰域點(diǎn)個(gè)數(shù)；Pˉ為P點(diǎn)鄰域點(diǎn)集的質(zhì)心點(diǎn)坐標(biāo)；M是一個(gè)對(duì)稱半正定矩陣，故其矩陣的3 個(gè)特征值都為非負(fù)值，對(duì)于點(diǎn)云數(shù)據(jù)，則這3 個(gè)特征值所對(duì)應(yīng)的特征向量為兩兩正交的向量。

對(duì)采樣點(diǎn)P鄰域內(nèi)點(diǎn)集進(jìn)行主成分分析可以得到3 個(gè)特征值λ1、λ2、λ3，其中λ1≥λ2≥λ3，這3 個(gè)特征值之間有以下關(guān)系：

1）λ1≥λ2≥λ3時(shí)，局部區(qū)域表現(xiàn)為線型，如圖4a所示。

2）λ1≈λ2≥λ3時(shí)，局部區(qū)域表現(xiàn)為面型，如圖4b所示。

3）λ1≈λ2≈λ3時(shí)，局部區(qū)域表現(xiàn)為球形，如圖4c所示。

圖4 特征分析示意圖

依據(jù)3個(gè)特征值按式（9）定義了維數(shù)特征[10]，其中，a1D+a2D+a3D=1。特征維數(shù)描述了局部區(qū)域點(diǎn)云數(shù)據(jù)的空間分布特征。

2.3 改進(jìn)后算法

在李健[4]等提出的融合kmeans 與Hausdorff 距離點(diǎn)云精簡(jiǎn)算法的基礎(chǔ)上，先采用維數(shù)特征的Hausdorff 距離提取特征點(diǎn)，再對(duì)剩余非特征點(diǎn)利用手肘法確定最優(yōu)聚類數(shù)，以每個(gè)聚類中心點(diǎn)作為非特征區(qū)域內(nèi)的特征點(diǎn)，保證數(shù)據(jù)均勻，不造成空洞。改進(jìn)后的算法首先避免了在平面區(qū)域曲率值過(guò)小而需要設(shè)定閾值的問(wèn)題，同時(shí)在確定聚類數(shù)中也保證了聚類精度。

對(duì)于某一采樣點(diǎn)pi和其鄰域內(nèi)某點(diǎn)qj的特征維數(shù)分別為a1D、a2D、a3D和則采樣點(diǎn)與其鄰域內(nèi)點(diǎn)維數(shù)特征的Hausdorff距離計(jì)算公式如下：

因此，采樣點(diǎn)pi的維數(shù)特征Hausdorff 距離可以表示為：

3 實(shí)例驗(yàn)證與分析

采用Rigel VZ-1000掃描儀對(duì)某古亭建筑進(jìn)行高密度多站掃描，依據(jù)標(biāo)靶球粗配準(zhǔn)后，再采用最近點(diǎn)迭代算法（iterative closest point，ICP）精配準(zhǔn)，然后利用SOR算法去除少量體外孤點(diǎn)，利用拉普拉斯濾波算法平滑表面毛刺點(diǎn)，配準(zhǔn)去噪后點(diǎn)云數(shù)據(jù)如圖5a所示。

本文通過(guò)將改進(jìn)的精簡(jiǎn)算法與李健提出的精簡(jiǎn)算法、基于曲率的精簡(jiǎn)算法、隨機(jī)采樣精簡(jiǎn)算法和體素下采樣精簡(jiǎn)算法的實(shí)驗(yàn)結(jié)果進(jìn)行比較，來(lái)評(píng)價(jià)改進(jìn)后的方法，5 種精簡(jiǎn)方法實(shí)驗(yàn)結(jié)果如圖5b、c、d、e、f所示。

圖5 不同精簡(jiǎn)方法實(shí)驗(yàn)結(jié)果

原始點(diǎn)云個(gè)數(shù)為726 417個(gè)，采用5種精簡(jiǎn)方法將古亭點(diǎn)云都精簡(jiǎn)至20%左右?；谇什捎煤碗S機(jī)采樣2 種方法在檐頂剔除了大量數(shù)據(jù)點(diǎn)，特別是基于曲率的精簡(jiǎn)方法造成了數(shù)據(jù)整體分布不均勻，本文方法、文獻(xiàn)[5]方法和體素下采樣方法都能夠獲得較為均勻的數(shù)據(jù)。圖6為本文方法和文獻(xiàn)[5]方法提取特征點(diǎn)和非特征區(qū)域內(nèi)特征點(diǎn)的結(jié)果圖，其中對(duì)于文獻(xiàn)[5]中基于曲率Hausdorff距離提取68 415個(gè)特征點(diǎn)，如圖6a所示；采用本文方法中基于維數(shù)特征Hausdorff距離提取70 791個(gè)特征點(diǎn)，如圖6d所示；依據(jù)手肘法確定最優(yōu)聚類個(gè)數(shù)約為70 000個(gè)，為保證相同精簡(jiǎn)率，對(duì)于文獻(xiàn)[5]中剩余非特征區(qū)域依據(jù)剩余非特征點(diǎn)密度提取聚類中心點(diǎn)73 057個(gè)點(diǎn)，如圖6b所示；本文方法對(duì)于剩余非特征區(qū)提取聚類中心7 000 個(gè)點(diǎn)，如圖6e 所示；文獻(xiàn)[5]方法和本文方法融合特征點(diǎn)和非特征區(qū)域特征點(diǎn)結(jié)果如圖6c、f 所示?？梢钥闯鲭m然提取的特征個(gè)數(shù)相同，但是本文方法的特征點(diǎn)相較于文獻(xiàn)[5]方法在檐頂?shù)奶卣鼽c(diǎn)更多。

圖6 2種Hausdorff距離提取特征點(diǎn)及非特征區(qū)域特征點(diǎn)圖

為了定量分析5種精簡(jiǎn)算法的精簡(jiǎn)精度，對(duì)精簡(jiǎn)后的點(diǎn)云數(shù)據(jù)建立三角網(wǎng)，通過(guò)統(tǒng)計(jì)三角形個(gè)數(shù)和三角形總面積來(lái)評(píng)定精簡(jiǎn)精度[11-12]。統(tǒng)計(jì)結(jié)果如表1 所示。以原始點(diǎn)云三角形總面積近似為真值，采用其他5 種精簡(jiǎn)方法精簡(jiǎn)后的點(diǎn)云構(gòu)建的三角網(wǎng)總面積都小于真值，符合精簡(jiǎn)后的數(shù)據(jù)精度有損失的規(guī)律，其中采用本文方法精簡(jiǎn)的結(jié)果最接近真值，其表面精度損失較少，即精簡(jiǎn)精度相對(duì)于其他4種方法精簡(jiǎn)精度較高。

表1 不同精簡(jiǎn)方法的結(jié)果比較

4 結(jié) 語(yǔ)

本文改進(jìn)了融合kmeans 聚類和Hausdorff 距離的點(diǎn)云精簡(jiǎn)方法，通過(guò)對(duì)聚類數(shù)的確定和Hausdorff距離計(jì)算參數(shù)加以改進(jìn)，提出了先通過(guò)計(jì)算特征維數(shù)的Hausdorff距離提取特征點(diǎn)，然后采用手肘法確定剩余非特征點(diǎn)的最優(yōu)聚類數(shù)，將特征點(diǎn)和聚類中心數(shù)據(jù)融合作為點(diǎn)云精簡(jiǎn)結(jié)果的技術(shù)路線；最后通過(guò)實(shí)際掃描的古亭建筑的點(diǎn)云數(shù)據(jù)，分別采用本文方法，文獻(xiàn)[5]方法、基于曲率采樣、隨機(jī)采樣和體素下采樣5 種方法進(jìn)行精簡(jiǎn)實(shí)驗(yàn)，結(jié)果表明改進(jìn)后的方法在相同的精簡(jiǎn)率下精度更高。