張韜，肖俊，王穎

(中國科學(xué)院大學(xué)人工智能學(xué)院， 北京 100049) (2020年2月20日收稿； 2020年4月17日收修改稿)

隨著信息化程度的不斷提升，在巖體工程中利用數(shù)值模擬技術(shù)對巖體進(jìn)行力學(xué)計(jì)算與分析，對預(yù)防地質(zhì)災(zāi)害的意義和價值越來越大。然而對巖體進(jìn)行精確的三維重建是對巖體進(jìn)行后續(xù)分析的基礎(chǔ)，因此如何獲取完整的巖體地表信息對巖體重建至關(guān)重要。

三維激光掃描是獲取巖體三維信息的重要設(shè)備，但由于掃描環(huán)境以及測量范圍的限制，對于大型巖體，往往不能通過一次掃描而獲得全部的場景，需要在不同的視角對巖體進(jìn)行多站掃描，而將多個視角獲得的巖體點(diǎn)云數(shù)據(jù)融合成一個點(diǎn)云數(shù)據(jù)是巖體點(diǎn)云配準(zhǔn)的主要工作。因此，點(diǎn)云配準(zhǔn)是三維數(shù)據(jù)處理中關(guān)鍵的步驟之一，也是后續(xù)進(jìn)行點(diǎn)云數(shù)據(jù)應(yīng)用的基礎(chǔ)。

點(diǎn)云配準(zhǔn)作為計(jì)算機(jī)視覺的基本問題，近年來，隨著研究人員深入的研究，取得了很多成果，其中最為常用的方法是迭代最近點(diǎn)算法[1](iterative closest point，ICP)。ICP算法是通過迭代的思想，每次迭代使目標(biāo)點(diǎn)云與源點(diǎn)云距離最小，不用尋找特征描述子，只需要優(yōu)化一個目標(biāo)函數(shù)，并使算法收斂。雖然這種算法簡單，但是初始位置以及噪聲點(diǎn)對該算法影響很大，容易陷入局部最優(yōu)，最終導(dǎo)致配準(zhǔn)失敗。王飛鵬等[2]提出基于高斯曲率的改進(jìn)算法，通過高斯曲率過濾掉部分點(diǎn)云，將剩余的點(diǎn)云應(yīng)用ICP算法進(jìn)行配準(zhǔn)，改善了ICP算法的效率。

此外，還有許多基于對應(yīng)特征點(diǎn)的算法被提出，這類算法是尋找源點(diǎn)云與目標(biāo)點(diǎn)云中的對應(yīng)點(diǎn)對，通過這些對應(yīng)的點(diǎn)對來計(jì)算點(diǎn)云變換參數(shù)。例如，Aiger等[3]提出的四點(diǎn)法(4-points congruent sets，4PCS)是一種基于剛體變換比例不變量的算法，利用剛體變換不會改變兩點(diǎn)的距離以及交叉點(diǎn)的對應(yīng)比例的性質(zhì)；Yang和Zang[4]提出基于曲線的配準(zhǔn)方法，通過提取點(diǎn)云中的曲線進(jìn)行配準(zhǔn)；Xian等[5]提出一個通過將三維無序點(diǎn)云轉(zhuǎn)換為二維數(shù)字圖像的方法，在二維圖像中通過尺度不變特征變換(scale invariant feature transform，SIFT)算法提取特征點(diǎn)進(jìn)行配準(zhǔn)；Wang等[6]提出基于高斯曲率的N階完全圖算法，通過選擇高斯曲率較大的點(diǎn)構(gòu)造一個N階完全圖的特征描述子進(jìn)行匹配；Hu等[7]利用巖體點(diǎn)云表面大部分為平面的特點(diǎn)，通過提取巖體點(diǎn)云的平面進(jìn)行配準(zhǔn)；Biber和StraBer[8]提出一種三維正態(tài)分布變換算法(3D-normal distribution transform，3 D-NDT)，利用點(diǎn)云表面點(diǎn)分布的統(tǒng)計(jì)信息以及迭代更新的策略實(shí)現(xiàn)配準(zhǔn)。

近幾年，由于機(jī)器學(xué)習(xí)方法的興起，一些基于機(jī)器學(xué)習(xí)的點(diǎn)云配準(zhǔn)算法被提出。Vongkulbhisal 等[9]提出一種判別優(yōu)化方法(discriminative optimization，DO)，該方法與傳統(tǒng)構(gòu)建損失函數(shù)和求解損失函數(shù)不同，它通過學(xué)習(xí)一個矩陣的更新序列求解損失函數(shù)從而得到最終的變換；Elbaz 等[10]提出一種用自編碼器進(jìn)行局部定位的算法(localization by registration using a deep auto-encoder reduced cover set，LORAX)，該算法通過將分割的點(diǎn)云投影成為深度圖，再通過自編碼器構(gòu)造一個特征描述子進(jìn)行配準(zhǔn)；Lu 等[11]提出一種通過端到端的神經(jīng)網(wǎng)絡(luò)尋找匹配點(diǎn)進(jìn)行配準(zhǔn)的方法；舒程珣等[12]提出通過基于圖像的卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行配準(zhǔn)的算法；劉鳴等[13]提出一種基于獨(dú)立成分分析的配準(zhǔn)方法,該方法通過提取源點(diǎn)云與目標(biāo)點(diǎn)云的獨(dú)立成分進(jìn)行點(diǎn)云配準(zhǔn)。

然而巖體點(diǎn)云不同于一般的點(diǎn)云，具有密度大、數(shù)據(jù)量大，且大部分區(qū)域?yàn)槠矫娴忍攸c(diǎn)，使得一些常用的點(diǎn)云配準(zhǔn)算法不能滿足巖體點(diǎn)云配準(zhǔn)的精度要求。本文基于巖體點(diǎn)云表面大部分區(qū)域?yàn)槠矫娴奶攸c(diǎn)，并根據(jù)巖體工程對于精度的較高要求，提出一種基于幾何特征逐層過濾的巖體點(diǎn)云配準(zhǔn)算法。該方法通過高斯曲率過濾掉絕大部分點(diǎn)云，再通過主曲率、主方向、法向量、鄰域的協(xié)方差矩陣的特征值和特征向量等進(jìn)行匹配點(diǎn)對的層層篩選，進(jìn)而準(zhǔn)確地找到對應(yīng)的匹配點(diǎn)對計(jì)算變換參數(shù)，最后利用ICP算法進(jìn)行精配準(zhǔn)。算法具有較高的精度。

1 背景知識

點(diǎn)云配準(zhǔn)即是將源點(diǎn)云與目標(biāo)點(diǎn)云融合在一起，設(shè)源點(diǎn)云為P，目標(biāo)點(diǎn)云為Q，則可將點(diǎn)云配準(zhǔn)抽象為：P經(jīng)過剛體變換，變換到Q的坐標(biāo)系下，即

Q=RP+T.

(1)

其中:R為一個正交矩陣，T為平移向量。點(diǎn)云配準(zhǔn)的主要工作就是計(jì)算R和T，如果求出2個點(diǎn)云中的匹配點(diǎn)對，則可以通過解一個方程組求出R和T，其中方程組為：

(2)

(3)

對于一個曲面，經(jīng)過剛體變換之后，部分幾何特征不會發(fā)生變化，而部分幾何特征會發(fā)生有規(guī)律的變化?？梢酝ㄟ^這些不變的幾何特征以及有規(guī)律變化的幾何特征進(jìn)行匹配點(diǎn)的篩選，從而進(jìn)行配準(zhǔn)。其中曲率作為剛體變換的不變量，可以用來篩選初始匹配點(diǎn)對。對于一個點(diǎn)云中的點(diǎn)變換到其對應(yīng)點(diǎn)，該點(diǎn)與其對應(yīng)點(diǎn)的鄰域的協(xié)方差矩陣相似，從而可以通過判斷鄰域協(xié)方差矩陣是否相似進(jìn)一步篩選對應(yīng)的匹配點(diǎn)。如果2個點(diǎn)互為對應(yīng)點(diǎn)，則可以通過計(jì)算各自的協(xié)方差矩陣的特征向量矩陣來計(jì)算旋轉(zhuǎn)矩陣。而曲面的法向量和主方向經(jīng)過剛體變換之后，也會發(fā)生相應(yīng)的旋轉(zhuǎn)，則可以利用法向量和主方向進(jìn)一步篩選。

2 基于匹配點(diǎn)逐層過濾的巖體點(diǎn)云配準(zhǔn)算法

基于對于巖體點(diǎn)云密度大、數(shù)量大以及巖體表面大部分為區(qū)域平面的特點(diǎn)，提出基于幾何特征逐層過濾的巖體點(diǎn)云配準(zhǔn)算法，該算法通過幾何特征來篩選匹配點(diǎn)對以完成粗配準(zhǔn)，最后用ICP算法進(jìn)行精配準(zhǔn)，算法框圖如圖1所示。

圖1 算法流程圖Fig.1 Algorithm flowchart

該算法包括4個主要步驟：1)計(jì)算高斯曲率等幾何特征，通過高斯曲率過濾源點(diǎn)云與目標(biāo)點(diǎn)云；2)逐層過濾匹配點(diǎn)對；3)矩陣聚類；4)ICP精配準(zhǔn)。下面介紹算法的關(guān)鍵步驟。

2.1 幾何特征的計(jì)算

點(diǎn)云的高斯曲率通過該點(diǎn)的鄰域點(diǎn)進(jìn)行計(jì)算。本文通過Zhang等[14]提出的算法計(jì)算高斯曲率。該算法計(jì)算點(diǎn)x的主曲率是通過該點(diǎn)的鄰域{x1,x2,…,xN}來計(jì)算。在計(jì)算出主曲率的同時，還計(jì)算了對應(yīng)的主方向以及法向量。最終通過k=k1×k2求出高斯曲率。該算法計(jì)算高斯曲率大小與鄰域大小有關(guān)，幾何特征計(jì)算受點(diǎn)云表面噪聲點(diǎn)影響較大，選取鄰域較大，則高斯曲率計(jì)算通常偏小，且耗時較多，但是魯棒性較強(qiáng)；選取鄰域較小，高斯曲率計(jì)算偏大，耗時較小，魯棒性較弱。對于配準(zhǔn)來說，如果噪聲較大，則應(yīng)該選擇較大的鄰域進(jìn)行計(jì)算，這樣魯棒性較強(qiáng)，噪聲較小，選擇較小的鄰域進(jìn)行計(jì)算，鄰域大小的選擇應(yīng)該根據(jù)點(diǎn)云的具體情況而定。一般而言，設(shè)置鄰域大小為100個點(diǎn)有較強(qiáng)的魯棒性，能夠得到較好的結(jié)果。

2.2 幾何特征逐層過濾

(4)

(5)

2.3 矩陣聚類以及精配準(zhǔn)階段

前面的步驟可以得到很多匹配的對應(yīng)點(diǎn)對，但是有的點(diǎn)對不是真正的對應(yīng)點(diǎn)對，需要進(jìn)行篩選。前面步驟求出的大部分矩陣應(yīng)該符合源點(diǎn)云到目標(biāo)點(diǎn)云的變換，少部分矩陣不符合該剛體變換，則可以對符合條件的矩陣進(jìn)行聚類，將矩陣數(shù)量最多的類別的對應(yīng)點(diǎn)對通過最小二乘法求出剛體變換的參數(shù)。目前，主流的精配準(zhǔn)算法為ICP算法，將算法前面階段計(jì)算的旋轉(zhuǎn)矩陣和平移向量作為ICP算法的初始值進(jìn)行迭代優(yōu)化，最終得到更精確的結(jié)果。

3 算法分析與比較

為了測試算法的魯棒性，分別在不同強(qiáng)度高斯噪聲、不同初始位置、不同重合度的條件下進(jìn)行實(shí)驗(yàn)。為了比較算法性能，用ICP算法和Wang等[6]提出的基于N階完全圖的算法(N-point complete graphs，NCG)進(jìn)行對比實(shí)驗(yàn)，其中ICP算法是在點(diǎn)云庫[15](point cloud library，PCL)中實(shí)現(xiàn)。本文在公共數(shù)據(jù)集RockBench[16]和在北京香山公園采集的數(shù)據(jù)上進(jìn)行實(shí)驗(yàn)。其中利用RockBench中的3個點(diǎn)云數(shù)據(jù)，為數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3，在香山公園中采集的數(shù)據(jù)為數(shù)據(jù)4、數(shù)據(jù)5、數(shù)據(jù)6，分別為不同角度采集的數(shù)據(jù)。實(shí)驗(yàn)中通過計(jì)算根均方差(root mean square，RMS)衡量算法精度。其中RMS公式為

(6)

3.1 魯棒性檢測實(shí)驗(yàn)

魯棒性檢測分為對噪聲、初始位置和重合度的魯棒性檢測。對于噪聲的魯棒性檢測是對源點(diǎn)云中的每個點(diǎn)替換成不同強(qiáng)度的高斯噪聲點(diǎn)，再將該點(diǎn)云變換到另一位置，成為該點(diǎn)云對應(yīng)的目標(biāo)點(diǎn)云，其中高斯噪聲的強(qiáng)度分別為(0，0.012m)，(0，0.042m)，(0，0.12m)，(0，0.152m)，源點(diǎn)云和目標(biāo)點(diǎn)云的初始位置以及匹配結(jié)果如圖2所示。

圖2 噪聲魯棒性實(shí)驗(yàn)結(jié)果Fig.2 Result of noise robustness test

對于初始位置魯棒性檢測，是將點(diǎn)云變換到不同的位置進(jìn)行匹配，源點(diǎn)云和目標(biāo)點(diǎn)云的初始位置以及匹配結(jié)果如圖3所示。對于不同重合度的魯棒性檢測，是在兩組部分重疊的點(diǎn)云上進(jìn)行實(shí)驗(yàn)，第1組點(diǎn)云重疊的點(diǎn)數(shù)為63 785，其中源點(diǎn)云中有118 230個點(diǎn)，目標(biāo)點(diǎn)云中有201 738個點(diǎn)。第2組點(diǎn)云重疊的點(diǎn)的個數(shù)為42 151，其中源點(diǎn)云中有185 308個點(diǎn)，目標(biāo)點(diǎn)云有164 143個點(diǎn)，源點(diǎn)云和目標(biāo)點(diǎn)云的初始位置以及匹配結(jié)果如圖4，配準(zhǔn)RMS如表1所示。

圖3 不同初始位置的配準(zhǔn)結(jié)果Fig.3 Registration results for different initial positions

圖4 不同重合度配準(zhǔn)結(jié)果Fig.4 Registration results with different degrees of coincidence

表1 魯棒性測試RMSTable 1 RMS of robustness test

從圖2、圖3、圖4及表1可以看出，初始位置對算法無影響，不同的重疊率對算法無影響，在存在高斯噪聲的情況下，算法比較穩(wěn)定，在高斯噪聲較小的情況下，算法有較高的精度。

3.2 算法對比實(shí)驗(yàn)

本文對比實(shí)驗(yàn)包括ICP算法、NCG算法，以及本文提出的經(jīng)過ICP精配準(zhǔn)的算法和未經(jīng)過ICP精配準(zhǔn)的算法，并通過計(jì)算RMS比較算法精度。分別在RockBench數(shù)據(jù)集中的數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3和在北京香山公園采集的數(shù)據(jù)4、數(shù)據(jù)5、數(shù)據(jù)6中進(jìn)行了對比實(shí)驗(yàn)。實(shí)驗(yàn)數(shù)據(jù)的基本信息如表2所示。數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3的初始位置如圖5，目標(biāo)點(diǎn)云是通過對源點(diǎn)云進(jìn)行變換得到，數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3的配準(zhǔn)結(jié)果如圖6所示，RMS如表3所示。數(shù)據(jù)4、數(shù)據(jù)5、數(shù)據(jù)6的初始位置如圖7所示，配準(zhǔn)RMS如表4所示，配準(zhǔn)結(jié)果如圖8所示。

表2 點(diǎn)云數(shù)據(jù)基本信息Table 2 Information of point cloud data

由表3可知，ICP算法和NCG算法在理想的情形下精度很高，而本文提出的算法精度能夠與ICP算法和NCG算法媲美。由表4可以看出，本

表3 數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3的RMSTable 3 RMS of Data1, Data2， and Data3

文的算法在香山的3個巖石數(shù)據(jù)上表現(xiàn)都很穩(wěn)定，本文的粗配準(zhǔn)算法和經(jīng)過精細(xì)配準(zhǔn)的算法精度相差很小，而且本文粗配準(zhǔn)算法在精度上比其他2種算法要高。ICP算法由于存在局部最優(yōu)的問題，導(dǎo)致ICP在數(shù)據(jù)4和數(shù)據(jù)6上未匹配成功。而NCG算法在3個點(diǎn)云數(shù)據(jù)上雖然都匹配成功，但是精度都沒有本文未經(jīng)過ICP精配準(zhǔn)的算法精度高。根據(jù)實(shí)驗(yàn)可知，本算法能夠有效避免初始位置對算法的影響。在篩選匹配點(diǎn)對的過程中，本算法是通過匹配鄰域來匹配對應(yīng)點(diǎn)，這使得算法對噪聲比較魯棒，而且算法是通過幾何特征逐層篩選匹配點(diǎn)對，最終能夠剔除誤匹配點(diǎn)對，從而得到精確的匹配點(diǎn)對。

圖5 數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3初始位置Fig.5 Initial position of Data1， Data2, and Data3

圖6 數(shù)據(jù)1、數(shù)據(jù)2、數(shù)據(jù)3匹配結(jié)果Fig.6 Registration result of Data1, Data2, and Data3

圖7 數(shù)據(jù)4、數(shù)據(jù)5、數(shù)據(jù)6的初始位置Fig.7 Initial position of Data4, Data5, and Data6

圖8 算法在數(shù)據(jù)4、數(shù)據(jù)5、數(shù)據(jù)6的匹配結(jié)果Fig.8 Registration results of Data4, Data5, and Data6

3.3 運(yùn)行時間分析

本文算法總體時間主要分為4部分：Tg表示計(jì)算高斯曲率等幾何特征的時間復(fù)雜度，Tr表示逐層過濾篩選匹配點(diǎn)對的時間復(fù)雜度，Tc表示矩陣聚類的時間復(fù)雜度，Ticp表示精配準(zhǔn)的時間復(fù)雜度。這里，

Tg=O(nlogn+mlogm)，

(7)

T=Tg+Tr+Tc+Ticp.

(8)

若選取合適的參數(shù)，可以使K1，K2，K3較小，最終使得

T≈O(NlogN).

(9)

其中，N=max(m,n)。

對于NCG算法，若選取合適的參數(shù)，可使NCG算法時間復(fù)雜度為：TNCG≈O(NlogN)，其中，N=max(m,n)。

ICP算法時間復(fù)雜度為：Ticp=O(mlogn)。因此當(dāng)源點(diǎn)云與目標(biāo)點(diǎn)云的數(shù)量相差不大時，3個算法的時間復(fù)雜度相近。

4 總結(jié)

本方法針對巖體點(diǎn)云的數(shù)據(jù)量大、密度大，大部分區(qū)域?yàn)槠矫娴奶攸c(diǎn)，提出一種基于幾何特征逐層過濾的巖體點(diǎn)云配準(zhǔn)算法，該算法通過高斯曲率過濾掉絕大部分點(diǎn)云，再通過主曲率、主方向、法向量、鄰域的協(xié)方差矩陣的特征值和特征向量逐層篩選匹配點(diǎn)對，最終精確地找到匹配點(diǎn)對。實(shí)驗(yàn)結(jié)果表明，本文算法能夠很好地匹配巖體點(diǎn)云，對于初始位置和高斯噪聲都具有較好的魯棒性。本文方法目前只針對大型巖體點(diǎn)云，在包含其他特性的點(diǎn)云上的適用性未做深入討論，這也是后續(xù)工作的重點(diǎn)。