基于面提取的三維巖體點(diǎn)云孔洞檢測與修復(fù)方法*

馬曌月，肖俊，王穎

(中國科學(xué)院大學(xué)人工智能學(xué)院， 北京 100049)

三維重建是一種將真實(shí)物體場景刻畫為符合計(jì)算機(jī)邏輯表達(dá)的數(shù)學(xué)模型的技術(shù)，在計(jì)算機(jī)圖形學(xué)、計(jì)算機(jī)視覺以及虛擬現(xiàn)實(shí)等領(lǐng)域有廣泛的應(yīng)用。巖體三維數(shù)值模型在巖石工程領(lǐng)域可以作為有限元、離散元、非連續(xù)變形分析等數(shù)值分析方法的模型輸入，其數(shù)值分析結(jié)果對巖體穩(wěn)定性判斷、關(guān)鍵塊體識別以及在工程實(shí)施過程中做出科學(xué)決策起著重要作用。因此，利用三維重建技術(shù)建立巖體三維數(shù)值模型在巖體工程領(lǐng)域有著十分重要的地位。

在三維模型重建中，三維激光掃描技術(shù)因其采集數(shù)據(jù)速度快、精度高，逐漸成為重建數(shù)據(jù)獲取主流方式之一。三維掃描儀獲取的點(diǎn)云數(shù)據(jù)是能夠刻畫物體表面特性的空間數(shù)據(jù)，為描述三維現(xiàn)實(shí)世界提供了最直接和有效的表達(dá)方式。

然而，在使用如圖1(a)所示的三維掃描儀掃描巖體點(diǎn)云數(shù)據(jù)的過程中，由于掃描測量角度不同，掃描儀對不同距離的巖體表面進(jìn)行掃描所獲得的數(shù)據(jù)信息的準(zhǔn)確度和精度也是不同的。在目標(biāo)巖體中，由于如圖1(b)所示的障礙物的陰影和遮擋等問題都會導(dǎo)致采集的點(diǎn)云數(shù)據(jù)丟失一部分。巖體自身的復(fù)雜拓?fù)浣Y(jié)構(gòu)也會導(dǎo)致掃描數(shù)據(jù)的缺失；此外，原始點(diǎn)云的濾波處理和障礙物移除等前期工作，也會導(dǎo)致點(diǎn)云數(shù)據(jù)中出現(xiàn)無法預(yù)料的數(shù)據(jù)丟失。數(shù)據(jù)丟失會導(dǎo)致巖體表面產(chǎn)生孔洞，從而在一定程度上影響后續(xù)巖體模型重建結(jié)果的準(zhǔn)確性。

圖1 掃描巖體數(shù)據(jù)示意圖Fig.1 Scanning rock mass data

點(diǎn)云的孔洞檢測及修復(fù)對后續(xù)的重建模型具有重要意義，當(dāng)前主流的點(diǎn)云孔洞檢測與修復(fù)算法可分為基于三角網(wǎng)格和基于點(diǎn)云數(shù)據(jù)的算法。基于三角網(wǎng)格的修復(fù)算法比較成熟，典型方法包括基于體素的修復(fù)和基于表面重建的修復(fù)；而基于點(diǎn)云數(shù)據(jù)的修復(fù)算法則往往針對特定類型的模型進(jìn)行修復(fù)，如垂直建筑物和文物模型等。以上方法更適用于形狀規(guī)則、密度均勻的三維點(diǎn)云模型。雖然主流算法能夠滿足多類模型的修復(fù)需求，但是在實(shí)際工程計(jì)算中，由于掃描的巖體點(diǎn)云數(shù)據(jù)量太大、噪聲和離群點(diǎn)較多、孔洞數(shù)目多且復(fù)雜，現(xiàn)有的檢測和修復(fù)方法都不能對大型巖體點(diǎn)云數(shù)據(jù)獲得滿意的修復(fù)效果。

由于真實(shí)巖體點(diǎn)云中包含大量不規(guī)則平面，平面提取是巖體點(diǎn)云數(shù)據(jù)處理的重要環(huán)節(jié)，本文提出基于平面提取的巖體點(diǎn)云孔洞檢測與修復(fù)方法。目前巖體平面提取算法主要包括基于區(qū)域生長[1]、基于模型擬合[2]及基于特征聚類[3]等3類方法。本文首先采用一種優(yōu)化的區(qū)域生長算法進(jìn)行平面提取，然后在平面提取的結(jié)果上，將點(diǎn)云映射至對應(yīng)平面獲取鄰域夾角，通過設(shè)置點(diǎn)云鄰域夾角閾值檢測孔洞邊界。接著根據(jù)孔洞邊界點(diǎn)集所屬平面及擬合誤差劃分孔洞類型，通過擬合平面并新增采樣點(diǎn)修復(fù)孔洞。由于平面提取過程包含數(shù)據(jù)去噪、平面合并等過程，可以過濾絕大部分的噪聲點(diǎn)，對后續(xù)巖體點(diǎn)云數(shù)據(jù)孔洞修復(fù)的計(jì)算有直接的簡化作用。本文算法無需對點(diǎn)云進(jìn)行網(wǎng)格化處理，可以快速有效地解決巖體點(diǎn)云數(shù)據(jù)的修復(fù)問題。

1 相關(guān)工作

針對研究目的和適用領(lǐng)域的不同，研究者們提出了大量的孔洞檢測與修復(fù)算法。本文僅介紹應(yīng)用較廣泛、較前沿的算法，更多詳細(xì)的檢測與修復(fù)算法可參考文獻(xiàn)[4-5]。

1.1 孔洞檢測

點(diǎn)云孔洞檢測環(huán)節(jié)對后續(xù)孔洞修復(fù)流程來說十分重要。目前孔洞檢測算法主要包括基于三角網(wǎng)格與基于點(diǎn)云數(shù)據(jù)的兩類檢測算法?；谌蔷W(wǎng)格的孔洞識別通過將三維點(diǎn)云模型三角網(wǎng)格化，在三角網(wǎng)格結(jié)構(gòu)中，根據(jù)三角形面片的拓?fù)潢P(guān)系進(jìn)行孔洞識別。當(dāng)網(wǎng)格模型較大時，孔洞邊界點(diǎn)判斷計(jì)算量會隨之加大，效率隨之降低。

基于點(diǎn)云數(shù)據(jù)的孔洞檢測方法一般是通過建立點(diǎn)云空間索引，遍歷點(diǎn)云數(shù)據(jù)，根據(jù)孔洞邊界點(diǎn)的不同特征設(shè)置閾值過濾條件，從而提取孔洞邊界點(diǎn)。根據(jù)空間維度不同，基于點(diǎn)云數(shù)據(jù)的識別算法又可分為基于二維平面投影映射的識別和基于三維空間的識別?；诙S平面投影映射的邊界識別典型算法包括創(chuàng)建投影平面，計(jì)算投影點(diǎn)的角度標(biāo)準(zhǔn)閾值差，從而確定孔洞邊界點(diǎn)[6]；對投影點(diǎn)進(jìn)行參數(shù)化后，根據(jù)鄰域點(diǎn)集在采樣點(diǎn)處的場力大小之和可以表示點(diǎn)集的平均值，設(shè)置閾值，從而識別點(diǎn)云的邊界特征點(diǎn)[7]；其他還包括加權(quán)度量[8]、吊錘法[4]等判別邊界特征點(diǎn)?；谌S空間的邊界識別算法則多考慮點(diǎn)云數(shù)據(jù)本身的三維特征，基于點(diǎn)云曲率、鄰域點(diǎn)法向夾角等設(shè)置單一或混合閾值進(jìn)行檢測，該類算法對參數(shù)設(shè)置敏感度高。此外，針對特定的點(diǎn)云數(shù)據(jù)，還有基于核回歸[9]等邊界提取方法，該類算法不通過閾值參數(shù)提取孔洞邊界，適用于特定的內(nèi)部不規(guī)則尖銳孔洞的識別優(yōu)化。

1.2 孔洞修復(fù)

目前點(diǎn)云的孔洞修復(fù)算法也主要包括基于三角網(wǎng)格的孔洞修復(fù)和基于點(diǎn)云數(shù)據(jù)的孔洞修復(fù)?；谌蔷W(wǎng)格的孔洞修復(fù)主要分為基于體素的算法和基于表面的算法?；隗w素的修復(fù)算法首先將網(wǎng)格模型變換成由體素表示的模型，然后在離散空間中應(yīng)用不同的方法進(jìn)行孔洞修復(fù)，其采用的典型算法為有向體素融合，在此基礎(chǔ)上根據(jù)孔洞鄰域表面信息，使用多種偏微分方程進(jìn)行聚斂優(yōu)化等[10-12]。基于表面的典型算法則是在描述點(diǎn)云的拓?fù)浣Y(jié)構(gòu)后，根據(jù)特征線匹配、幾何特征、形狀計(jì)算和能量函數(shù)、徑向基函數(shù)等[13-16]進(jìn)行孔洞修復(fù)?？傮w來說，基于三角網(wǎng)格的修復(fù)對于一般的規(guī)則物體(如機(jī)械零件、動物點(diǎn)云模型)具有較好的修復(fù)效果，對不規(guī)則的大型點(diǎn)云模型或孔洞區(qū)域較大的網(wǎng)格模型修復(fù)效果不佳[17]。

直接對點(diǎn)云數(shù)據(jù)進(jìn)行修復(fù)的研究包括將曲線擬合的方法應(yīng)用到建筑物垂直面的點(diǎn)云孔洞修復(fù)中，對點(diǎn)云數(shù)據(jù)投影至二維平面并進(jìn)行擬合，再恢復(fù)至三維空間插入數(shù)據(jù)點(diǎn)[18]；針對由稀疏3D體積構(gòu)成的高程面中的孔洞，將點(diǎn)云數(shù)據(jù)投影到二維平面上，按內(nèi)外順序依次計(jì)算孔洞內(nèi)部填充點(diǎn)二維坐標(biāo)，再根據(jù)鄰域點(diǎn)作用關(guān)系補(bǔ)齊三維坐標(biāo)[19]；其他算法還包括利用八叉樹構(gòu)建隱式曲面[20]等?；邳c(diǎn)云數(shù)據(jù)進(jìn)行修復(fù)的算法大都依賴孔洞周圍區(qū)域信息，對較大孔洞和邊界開孔洞無法提取有效的周圍區(qū)域信息，存在提取點(diǎn)云邊界不準(zhǔn)確、錯誤識別孔洞、耗時長等問題，修復(fù)效果不佳[5]。此外，將神經(jīng)網(wǎng)絡(luò)應(yīng)用于點(diǎn)云孔洞修復(fù)中也是近年來的算法嘗試趨勢，包括在修復(fù)軟件的基礎(chǔ)上，通過改進(jìn)的BP神經(jīng)網(wǎng)絡(luò)重設(shè)置初始權(quán)值和閾值，對修復(fù)數(shù)據(jù)進(jìn)行優(yōu)化等[21]。該類方法適用范圍較小，對特定模型調(diào)整神經(jīng)網(wǎng)絡(luò)算法模型參數(shù)耗時也較多。目前已有的修復(fù)方法大多適用于一般的規(guī)則點(diǎn)云模型，本文根據(jù)巖體點(diǎn)云的平面特性，在平面提取結(jié)果的基礎(chǔ)上進(jìn)行點(diǎn)云孔洞檢測與修復(fù)。

2 本文算法

本文提出的基于巖體點(diǎn)云平面提取的孔洞檢測與修復(fù)算法流程圖如圖2所示。由于三維掃描激光掃描得到的高精度、大尺寸的巖體點(diǎn)云數(shù)據(jù)會消耗大量的存儲空間，從而影響修復(fù)的速度和效率，為了得到良好的巖體點(diǎn)云修復(fù)效果，需要解決以下問題：如何有效地對輸入的大規(guī)模且含較多噪聲點(diǎn)的巖體點(diǎn)云進(jìn)行孔洞檢測；如何快速準(zhǔn)確地修復(fù)巖體點(diǎn)云中出現(xiàn)的大尺寸孔洞。由于巖體數(shù)據(jù)量大、密度不均勻且夾雜很多離群點(diǎn)，應(yīng)用一種優(yōu)化的平面提取算法對離群點(diǎn)和噪聲點(diǎn)進(jìn)行過濾，從而減少誤差和冗余計(jì)算的時間。對巖體點(diǎn)云構(gòu)建快速空間索引，在擬合平面的基礎(chǔ)上進(jìn)行修復(fù)，保留原始巖體點(diǎn)云的幾何特征，從而達(dá)到快速準(zhǔn)確修復(fù)的目的。

因此，本文算法主要包括3個步驟：基于優(yōu)化區(qū)域生長的平面提取、點(diǎn)云孔洞檢測與孔洞分類修復(fù)。該算法通過平面提取對巖體點(diǎn)云數(shù)據(jù)進(jìn)行去噪及分類處理，提高了孔洞邊界檢測的效率；平面提取的平面擬合結(jié)果也可直接用于后續(xù)的孔洞修復(fù)中的特征面建立，使得孔洞檢測與修復(fù)算法更為高效。

2.1 平面提取

使用一種適用于巖體點(diǎn)云的優(yōu)化區(qū)域生長算法[22]進(jìn)行平面提取，研究對象為巖體點(diǎn)云模型Rock1-3。區(qū)域生長是一種穩(wěn)健高效的平面提取算法，首先算法對激光掃描儀采集到的原始點(diǎn)云進(jìn)行預(yù)處理，去除無效點(diǎn)，計(jì)算每個點(diǎn)的法向量。其次，將剩余點(diǎn)云根據(jù)分辨率劃分為體素，使用RANSAC算法初步識別每個體素中的平面。然后將所有體素中的最大平面作為生長單元，采用區(qū)域生長算法聚集相似平面的體素。最后，將未分配的點(diǎn)集進(jìn)行二次生長并映射至相鄰平面，提取出巖體點(diǎn)云所有的平面集合。

點(diǎn)云預(yù)處理通常掃描的原始點(diǎn)云中包含噪聲點(diǎn)，這些噪聲點(diǎn)會影響平面提取的準(zhǔn)確性。因此首先對輸入點(diǎn)云數(shù)據(jù)建立kd-tree空間索引，然后遍歷全部點(diǎn)云及其對應(yīng)鄰域數(shù)據(jù)，刪除孤立、異常坐標(biāo)等噪聲點(diǎn)。由于點(diǎn)的法向信息對后續(xù)不同方向的平面提取十分重要，采用主成分分析法(principal components analysis，PCA)估計(jì)法向信息。定義點(diǎn)云鄰域?yàn)橐渣c(diǎn)云中任意一點(diǎn)為中心的一定半徑尺度的鄰域球，即該點(diǎn)與其鄰域內(nèi)的其他點(diǎn)組成的局部空間。對點(diǎn)云數(shù)據(jù)中的每一個點(diǎn)P(x,y,z)，在半徑為r(考慮巖體點(diǎn)云平均密度分布，此處設(shè)置r=0.5 m)的球面上查找其k(k=40)最近鄰。構(gòu)建點(diǎn)P的協(xié)方差矩陣，然后選擇與最小特征值相對應(yīng)的特征向量,作為點(diǎn)P的法向量。計(jì)算公式如下

(1)

∑·Vj=λj·Vj,j∈{0,1,2}.

(2)

圖2 孔洞檢測與修復(fù)算法流程圖Fig.2 Flow chart of holes detection and filling algorithm

生成體素在三維空間中，一個體素及其26個相鄰體素構(gòu)建為一個三階立方體，因此，每個體素最多具有26個相鄰體素。為了從點(diǎn)云中生成體素，首先使用點(diǎn)云的三維坐標(biāo)信息計(jì)算邊界框。然后，根據(jù)分辨率設(shè)置體素大小vsize將點(diǎn)云空間劃分為體素，同時建立相鄰體素信息。最后，刪除空體素以提高體素查詢效率。此步驟中，體素大小的參數(shù)設(shè)置取決于巖體點(diǎn)云模型的尺寸規(guī)模和點(diǎn)云密度，體素設(shè)置過大，會導(dǎo)致每個體素內(nèi)部包含多個平面，影響后續(xù)合并過程；體素設(shè)置過小，可能會導(dǎo)致在擬合平面過程中容易受到噪聲干擾，且影響算法性能。通過實(shí)驗(yàn)調(diào)參，此處對實(shí)驗(yàn)?zāi)Ｐ蚏ock1、Rock2設(shè)置參數(shù)為0.8 m，對Rock3設(shè)置參數(shù)為0.65 m。

RANSAC平面檢測為了獲得初始高質(zhì)量的生長單元，在每個體素中使用健壯的RANSAC算法檢測平面，選擇體素中最大平面作為生長單元。為避免生成偽平面，將法向量作為約束條件。平面檢測算法每次從點(diǎn)云中選擇3個非共線的點(diǎn)擬合平面，然后統(tǒng)計(jì)點(diǎn)云中落在此平面上的散點(diǎn)個數(shù)，將其作為該擬合平面的分?jǐn)?shù)。若點(diǎn)到平面的距離小于給定閾值d且點(diǎn)法線與平面法線的夾角小于給定閾值angle，則認(rèn)為該點(diǎn)落在平面上，經(jīng)過T次迭代可以檢測出最優(yōu)平面。

假設(shè)點(diǎn)云數(shù)據(jù)共包含N個點(diǎn)，其中最優(yōu)平面包含n，那么在一次檢測中即可檢測出最優(yōu)平面的概率為p1，則經(jīng)過T次檢測出最優(yōu)平面的概率等于T次檢測失敗的補(bǔ)集，假設(shè)經(jīng)過T次，至少檢測出一次最優(yōu)平面的最小概率為pt，則可以計(jì)算出在此概率下的最小檢測次數(shù)Tmin。

單次掃描中可檢測到該平面的概率p1為

(3)

T次試驗(yàn)后檢測成功的概率pt為

(4)

T次試驗(yàn)至少找到一組好的檢測結(jié)果最小概率為

(5)

RANSAC平面檢測后，每個體素只包含一個最優(yōu)平面。由于后續(xù)步驟是以平面為生長單元的區(qū)域生長算法，因此只有與當(dāng)前平面相交的體素才是可生長體素。對每個體素內(nèi)檢測到的平面進(jìn)行擴(kuò)展，并檢測相鄰體素與平面是否相交，若不相交，則該體素上的點(diǎn)集位于該平面的一端，因此不屬于同一平面，將該相鄰體素在鄰域中刪除。此步驟刪除了一些無效鄰近體素，經(jīng)過一次鄰域優(yōu)化，重建了點(diǎn)云鄰域信息，可以減少后續(xù)進(jìn)行區(qū)域生長步驟時的鄰域搜索次數(shù)，提升算法性能。

區(qū)域生長本文使用區(qū)域生長算法從所有體素中檢測平面。本文首先從體素集合選擇包含最多點(diǎn)數(shù)的體素作為種子區(qū)域，將兩個體素的代表平面之間的距離(d=0.15 m)和兩平面法線夾角(angle=15°)作為生長規(guī)則；然后根據(jù)生長規(guī)則選擇合適的相鄰體素進(jìn)行擴(kuò)展，直到所有相鄰體素均已添加完畢，得到一個平面點(diǎn)集，重復(fù)這一步驟，直至沒有體素作為種子區(qū)域；最后通過不斷地添加新的種子區(qū)域，迭代計(jì)算得到所有的平面點(diǎn)集合。

剩余點(diǎn)集處理在區(qū)域生長完成后，可以粗略檢測出初始平面集U1。但是由于一個體素中可能包含多個平面，因此在區(qū)域生長過程中，一個體素中較小的平面點(diǎn)集可能未能分配至對應(yīng)平面，出現(xiàn)剩余點(diǎn)集。剩余點(diǎn)集主要存在于平面內(nèi)部以及兩個平面交界處，內(nèi)部點(diǎn)云被遺棄的主要原因是巖體表面粗糙度有變化，導(dǎo)致這些點(diǎn)難以通過共面性檢測；邊界處的剩余點(diǎn)云主要是因?yàn)榉ㄏ蛄考s束的限制以及邊界處區(qū)域的復(fù)雜性，這種現(xiàn)象導(dǎo)致了在共面性檢測過程中小的平面被舍棄，因而需要進(jìn)行二次區(qū)域生長。

若剩余點(diǎn)集與其相鄰平面滿足二次區(qū)域生長的生長規(guī)則，則將其添加到對應(yīng)平面中進(jìn)行二次區(qū)域生長，生長規(guī)則與區(qū)域生長算法規(guī)則相同：1)點(diǎn)P與平面的距離小于距離閾值；2)點(diǎn)P的法向量與平面法向夾角小于角度閾值。此處剩余點(diǎn)集未被處理的主要原因是巖體點(diǎn)云粗糙性造成的平面點(diǎn)云信息丟失，因此此處閾值參數(shù)比上文區(qū)域生長使用的閾值高出30%～100%(d2=0.25 m，angle2=25°)。

表1記錄了本文平面提取算法的參數(shù)設(shè)置。

最終，算法將巖體點(diǎn)云中的所有點(diǎn)集都分配至對應(yīng)平面，得到最終的平面集合U2。

表1 平面提取算法參數(shù)設(shè)置Table 1 Parameter setting of plane extraction algorithm

2.2 孔洞檢測

建立kd-tree點(diǎn)云空間索引后，在巖體點(diǎn)云平面提取的基礎(chǔ)上，結(jié)合基于點(diǎn)云數(shù)據(jù)的邊界提取算法，通過分析巖體點(diǎn)云鄰域點(diǎn)夾角來判斷該點(diǎn)是否為邊界點(diǎn)，并據(jù)此提取孔洞邊界。

對于結(jié)構(gòu)特征復(fù)雜的巖體點(diǎn)云數(shù)據(jù)，通過平面提取算法，刪去離群點(diǎn)等噪聲點(diǎn)，對保留的點(diǎn)云數(shù)據(jù)進(jìn)行遍歷，依據(jù)kd-tree快速檢索點(diǎn)的鄰域信息。

此步驟中，選取點(diǎn)云數(shù)據(jù)集中任意一點(diǎn)M，當(dāng)點(diǎn)M及鄰域點(diǎn)集屬于同一平面時，將點(diǎn)M及其鄰域點(diǎn)投影至對應(yīng)平面，根據(jù)下式計(jì)算其鄰域夾角α。

(6)

當(dāng)點(diǎn)M及其鄰域點(diǎn)集屬于多個平面時，將該點(diǎn)M及其對應(yīng)鄰域點(diǎn)分別映射至對應(yīng)平面并求解相應(yīng)的鄰域夾角，然后多個鄰域夾角之和作為該點(diǎn)的對應(yīng)鄰域夾角。若點(diǎn)M的鄰域夾角和小于閾值(結(jié)合巖體真實(shí)數(shù)據(jù)集特性，本文設(shè)置為220°)，則加入孔洞邊界點(diǎn)集，否則作為平面內(nèi)部點(diǎn)并舍棄。

圖3展示了針對初始的點(diǎn)云邊界提取結(jié)果，由于閾值設(shè)置問題，頂點(diǎn)M3與M6的鄰域夾角大于設(shè)置閾值，但其作為孔洞邊界點(diǎn)并未被檢測到。

圖3 未檢測邊界點(diǎn)示意圖Fig.3 Not-detected boundary point diagram

對點(diǎn)云集合進(jìn)行角度閾值過濾可以得到位于孔洞邊界上的大部分頂點(diǎn)，而初始的孔洞邊界點(diǎn)集并不連通。因此對目前的孔洞邊界點(diǎn)集進(jìn)行廣度優(yōu)先遍歷，設(shè)置鄰域半徑(r=0.5 m)及相鄰點(diǎn)數(shù)(k=30)，將不相鄰的邊界點(diǎn)集的最短路徑上的頂點(diǎn)添加為孔洞邊界點(diǎn)。然后迭代地在擴(kuò)展后的孔洞邊界點(diǎn)集中尋找最小點(diǎn)集合回路，將得到的回路點(diǎn)集合作為孔洞點(diǎn)邊界線。

2.3 孔洞修復(fù)

本文算法在檢測出的巖體點(diǎn)云孔洞邊界線點(diǎn)集上進(jìn)行孔洞修復(fù)。根據(jù)巖體點(diǎn)云平面提取結(jié)果，將點(diǎn)云孔洞分類為單一平面孔洞和復(fù)合平面孔洞。定義所有邊界點(diǎn)集位于同一平面內(nèi)的孔洞為單一平面孔洞，邊界點(diǎn)集位于多個(大于等于2)平面相交處的孔洞為復(fù)合平面孔洞。

圖4(a)展示了所有邊界點(diǎn)集位于同一平面內(nèi)的孔洞；圖4(b)展示了所有邊界點(diǎn)集位于多個平面相交處的孔洞，圖4紅色圓框框出了孔洞局部放大區(qū)域。

圖4 平面不同孔洞圖Fig.4 Different types of holes

孔洞修復(fù)步驟包括：

1)構(gòu)建孔洞坐標(biāo)系，計(jì)算孔洞鄰近域的最小二乘擬合平面，作為孔洞點(diǎn)及其鄰近域的特征面。假設(shè)孔洞點(diǎn)及其鄰近域的數(shù)據(jù)點(diǎn)可以用M0、M1、…、Mn-1來表示，這些點(diǎn)到該特征面的距離的平方和應(yīng)該最小。與前文平面提取算法中的計(jì)算方法類似，計(jì)算該擬合平面的法線相當(dāng)于求其協(xié)方差矩陣的最小特征值所對應(yīng)的特征向量n，根據(jù)空間平面方程原理即可求解。

假設(shè)在孔洞坐標(biāo)系中的3個坐標(biāo)軸分別為U軸、V軸、W軸，規(guī)定O作為該孔洞坐標(biāo)系的原點(diǎn),特征面為UOV平面。將孔洞及其鄰近域投影到該特征平面上，特征面的法向量即為孔洞坐標(biāo)系的W軸，圖5給出了孔洞坐標(biāo)系示例，其中法向量為n(xs,ys,zs)。

2)新增采樣點(diǎn)。對點(diǎn)云孔洞進(jìn)行新增點(diǎn)采樣，將孔洞點(diǎn)及其鄰近域投影至孔洞坐標(biāo)系下，進(jìn)行二維投影，在二維面內(nèi)進(jìn)行插值。為保證修復(fù)后的巖體點(diǎn)云具有密度一致性、后續(xù)重建保持完整的巖體幾何拓?fù)溥B接，此處選取原始巖體點(diǎn)云數(shù)據(jù)間的平均密度作為插值距離，在二維投影面內(nèi)進(jìn)行等距差值。計(jì)算完畢后，將孔洞坐標(biāo)系變換為原始坐標(biāo)系，得到新增點(diǎn)采樣結(jié)果。

經(jīng)上述步驟后，孔洞修復(fù)完成。針對前文不同分類孔洞，將孔洞邊界點(diǎn)集擬合平面，當(dāng)孔洞邊界點(diǎn)集合與擬合平面的法向誤差超過閾值(15°)時，判斷該孔洞為復(fù)合孔洞，否則為單一孔洞。對于單一孔洞，根據(jù)上述步驟直接進(jìn)行擬合修復(fù)。對于復(fù)合孔洞，根據(jù)孔洞邊界點(diǎn)集合選擇其所屬的多個鄰域平面作為特征平面，計(jì)算多個平面的交線，然后在對應(yīng)多個鄰域平面內(nèi)分別進(jìn)行采樣點(diǎn)集修復(fù)。最后迭代地修復(fù)提取到的所有孔洞，完成孔洞修復(fù)。

圖5 孔洞坐標(biāo)系示意圖Fig.5 Hole coordinate system diagram

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

3.1 實(shí)驗(yàn)說明

為了驗(yàn)證本文算法的有效性，本文展示了本文算法與基于點(diǎn)云數(shù)據(jù)的孟算法[18]和Hai算法[19]對多個巖體點(diǎn)云數(shù)據(jù)模型的修復(fù)效果圖及量化對比，以充分論證本文算法的良好性能。

本文算法以C++實(shí)現(xiàn)，所采用的實(shí)驗(yàn)環(huán)境為IntelCorei5-6400CPU@2.70 GHz，8 GB內(nèi)存Windows10 64位操作系統(tǒng)，開源點(diǎn)云庫PCL1.7.2(Point Cloud Library，PCL)。實(shí)驗(yàn)數(shù)據(jù)為真實(shí)巖體點(diǎn)云模型Rock1-3,點(diǎn)云數(shù)據(jù)中包含100 000～600 000個點(diǎn)。

3.2 實(shí)驗(yàn)結(jié)果與分析

為驗(yàn)證算法的有效性，分別在點(diǎn)云規(guī)模和孔洞復(fù)雜性不同的巖體點(diǎn)云上進(jìn)行修復(fù)實(shí)驗(yàn)。

圖6分別展示了本文算法與其他算法對Rock1-3模型的修復(fù)對比。其中，圖6(a)為各算法的原始輸入模型Rock1-3的點(diǎn)云數(shù)據(jù)圖，6(b)為孟慶年等[18]的修復(fù)算法效果，6(c)為Hai[19]的修復(fù)算法效果，6(d)為本文修復(fù)算法效果。圖6(b)、6(c)、6(d)中，白色點(diǎn)云區(qū)域?yàn)樵键c(diǎn)云數(shù)據(jù)，彩色區(qū)域?yàn)楦魉惴▽?yīng)的實(shí)際修復(fù)點(diǎn)云，所有實(shí)驗(yàn)結(jié)果在Cloud Compare中進(jìn)行可視化顯示。

如圖6所示，在整體點(diǎn)云修復(fù)效果對比中，本文算法修復(fù)效果最佳。在圖6(b)中，孟算法[18]孔洞修復(fù)結(jié)果存在部分小孔洞未識別，黃色圓框框出部分展示了在巖體點(diǎn)云密度不均勻處，有多處錯誤修復(fù)，在點(diǎn)云外部區(qū)域添加了多余的修復(fù)點(diǎn)集。圖6(c)中，Hai[19]也存在少量孔洞未識別的問題，并且往往將幾個孔洞合并為同一大孔洞修復(fù)，破壞了原有的巖體點(diǎn)云結(jié)構(gòu)，導(dǎo)致修復(fù)結(jié)果錯誤。圖6(d)展示了本文算法的修復(fù)效果，在巖體點(diǎn)云平面提取的結(jié)果上，可以檢測并修復(fù)絕大部分孔洞，修復(fù)局部細(xì)節(jié)更為平滑，與周圍區(qū)域融合效果更佳，并且保留了原有巖體點(diǎn)云特征和鄰域信息，邊界特征清晰。

表2展示了本文與其他算法的定量統(tǒng)計(jì)對比，Rock1-3原始模型原有孔洞初始值分別為11、26、37。對比項(xiàng)包括修復(fù)后的巖體點(diǎn)云數(shù)據(jù)頂點(diǎn)數(shù)Q(原始模型初始點(diǎn)云數(shù)記為Original)，算法對孔洞的檢測數(shù)N、正確修復(fù)孔洞數(shù)M、修復(fù)錯誤數(shù)E，其中E為未識別孔洞數(shù)與修復(fù)錯誤數(shù)的總和，修復(fù)算法運(yùn)行時間t/s。

圖6 Rock1-3修復(fù)效果對比圖Fig.6 Rock1-3 filling effect comparison

表2 Rock1-3模型修復(fù)定量對比Table 2 Rock1-3 filling effect statistical comparison

由表2知，從修復(fù)正確率上分析，Rock1模型中初始孔洞數(shù)11個，本文算法全部正確識別并正確修復(fù)；Rock2模型中初始孔洞數(shù)26個，本文算法正確識別25個并正確修復(fù)；Rock3模型中初始孔洞數(shù)37個，本文算法正確識別35個并正確修復(fù)34個，對比其他算法正確率更高。針對位于巖體外邊界的孔洞，本文算法存在部分未能正確修復(fù)的情況，其原因?yàn)椴糠治挥谕膺吔缟系目锥礋o法有效檢測。修復(fù)后的點(diǎn)云頂點(diǎn)數(shù)Hai[19]最多，體現(xiàn)其修復(fù)質(zhì)量不佳，這是由于其算法更適用于垂直物高程面，修復(fù)孔洞區(qū)域過大導(dǎo)致的；孟慶年等[18]修復(fù)后的頂點(diǎn)數(shù)也比本文算法更多，這是由于其部分孔洞未識別，且修復(fù)點(diǎn)云錯誤導(dǎo)致的，這在圖6中均有展示。

從算法運(yùn)行效率上分析，本文算法運(yùn)行時間最短，且隨點(diǎn)云數(shù)據(jù)量的增加，耗時增加較少，這是由于算法考慮了巖體點(diǎn)云平面特性，在平面提取的基礎(chǔ)上不需要網(wǎng)格化和再次調(diào)參來提取點(diǎn)云特征，因此在運(yùn)行效率上有較大優(yōu)勢。孟慶年等[18]和Hai[19]由于需多次擬合誤差，在計(jì)算過程中耗時較長。由運(yùn)行時間定量對比可知，本文修復(fù)算法運(yùn)行效率更高。

總體來說，對于Rock1-3巖體點(diǎn)云模型，本文算法對比其他基于點(diǎn)云數(shù)據(jù)的修復(fù)方法(孟慶年等[18]、Hai[19])，能在最大程度上識別巖體點(diǎn)云原有的孔洞，并對識別到的孔洞進(jìn)行正確修復(fù)，孔洞識別正確率最高，同時修復(fù)錯誤率最小。

3.3 算法分析

本文提出的算法針對巖體點(diǎn)云特性進(jìn)行檢測修復(fù)，提高了巖體點(diǎn)云數(shù)據(jù)修復(fù)的準(zhǔn)確性，同時算法效率較高，魯棒性也較強(qiáng)。相較于其他基于點(diǎn)云數(shù)據(jù)修復(fù)的算法，利用了巖體平面提取的特性，平面提取時過濾了噪聲點(diǎn)，簡化了后續(xù)邊界提取步驟；平面擬合過程簡化了后續(xù)孔洞修復(fù)過程中的擬合孔洞平面步驟，因此效率更高，且針對大型不規(guī)則巖體點(diǎn)云數(shù)據(jù)修復(fù)適用性更好。其他算法則存在需要多次誤差擬合導(dǎo)致計(jì)算變慢、對結(jié)構(gòu)復(fù)雜的巖體點(diǎn)云模型不適用、對密度不均勻處的點(diǎn)云計(jì)算錯誤導(dǎo)致巖體孔洞誤識別誤修復(fù)等問題。

因此，本文算法能較好較快地解決巖體點(diǎn)云數(shù)據(jù)孔洞檢測與修復(fù)的問題，對比已有修復(fù)方法，修復(fù)準(zhǔn)確性、運(yùn)行效率更高，對巖體點(diǎn)云模型適用性更佳，優(yōu)勢明顯。

4 總結(jié)

提出一種基于巖體平面提取的點(diǎn)云數(shù)據(jù)檢測修復(fù)方案，算法的重要貢獻(xiàn)之一在于首次將點(diǎn)云平面提取的結(jié)果應(yīng)用于孔洞檢測與修復(fù)，理論上平面提取的效果越好，算法最終的效果越好。從實(shí)驗(yàn)結(jié)果來看，本方法針對于巖體點(diǎn)云模型，孔洞修復(fù)識別率和正確率更高，運(yùn)行時間更快，整體修復(fù)效果更佳。由于考慮巖體點(diǎn)云特性進(jìn)行孔洞檢測修復(fù)，本文算法目前針對大型巖體點(diǎn)云模型效果良好，對其他特性的點(diǎn)云數(shù)據(jù)集的適用性未做深入討論，這也是后續(xù)工作的一項(xiàng)重點(diǎn)，另外，由于在點(diǎn)云外邊界處缺失數(shù)據(jù)未能有效地檢測孔洞，這部分工作也需要后續(xù)進(jìn)一步研究。