李 燦 孫 未 張力云

(成都理工大學(xué)核技術(shù)與自動化工程學(xué)院 四川 成都 610059)

0 引 言

點云配準(zhǔn)技術(shù)是一種重要的數(shù)字檢測技術(shù)，現(xiàn)已蓬勃發(fā)展并廣泛運用于各個行業(yè)，如計算機視覺、高精度加工、無損檢測、虛擬現(xiàn)實等[1-4]。通過三維掃描儀對物體表面進行多角度掃描，得到不同視角下的三維點云數(shù)據(jù)，再將所得的點云數(shù)據(jù)進行配準(zhǔn)，即可得到完整的模型。Besl等[5]提出的最近點迭代(Iterative Closest Point,ICP)算法廣泛運用于點云配準(zhǔn)中，作為經(jīng)典的配準(zhǔn)算法，其主要手段是通過尋找兩個點集的對應(yīng)點，并計算其變換矩陣，但該算法容易陷入局部最小值并且效率不高[6-7]。

由于ICP算法存在的部分缺陷，眾多學(xué)者都提出了各種各樣的改進算法來減小配準(zhǔn)誤差和提高配準(zhǔn)效率。Bae[8]提出基于曲率和法向變化率的幾何基本最近點迭代(GP-ICP)算法，該算法使得ICP算法對迭代初始值的要求降低。Ying等[9]提出基于七維空間迭代的Scale-ICP算法，雖然其收斂速度較快，但對迭代過程非常依賴且僅能應(yīng)用于點云整體存在相似變換的情況。Sharp等[10]提出的最近點迭代(ICPIF)算法主要利用不變特征配準(zhǔn)，該算法可實現(xiàn)局部重疊的點云配準(zhǔn)，但是在理論上對點云數(shù)據(jù)缺失情況下的配準(zhǔn)很難實現(xiàn)。

本文提出一種基于高斯似然估計因子分析的點云配準(zhǔn)算法。將點云配準(zhǔn)的數(shù)學(xué)公式模型擴展為因子分析模型，這樣便將點云配準(zhǔn)轉(zhuǎn)化為對因子模型參數(shù)的求解；采用多元高斯函數(shù)對點云數(shù)據(jù)進行逼近，并通過EM算法求解出因子載荷矩陣，再利用所求得的因子載荷矩陣完成對點云的配準(zhǔn)。實驗表明，與其他幾種算法相比，即使在點云無序、數(shù)據(jù)存在遮擋、缺失不完整，以及噪聲環(huán)境下，本文算法仍可實現(xiàn)快速精確配準(zhǔn)，收斂速度快，具有良好的魯棒性。

1 基本原理

1.1 因子分析

一組D維的觀測數(shù)據(jù)X={x1,x2,…,xN}的正交因子模型為：

(1)

(2)

通過上述正交因子數(shù)學(xué)模型可知，觀測數(shù)據(jù)X的協(xié)方差矩陣為：

∑=AAT+φ

(3)

式中：φ為特殊因子矩陣ε的協(xié)方差。

1.2 點云配準(zhǔn)原理

將源點云P={p1,p2,…,pn}和目標(biāo)點云Q={q1,q2,…,qm}進行剛性幾何變換：

pi=Rqj+ti=1,2,…,n;j=1,2,…,m

(4)

式中：R∈R3×3為旋轉(zhuǎn)矩陣,且RTR=I；t∈R3×1為平移向量。由式(3)可知，兩點云的協(xié)方差矩陣可表示為:

由式(4)可得到兩點云協(xié)方差矩陣的關(guān)系為：

∑P=R∑QRT

(5)

點云數(shù)據(jù)在實際運用中可能會有噪聲干擾，可以通過正交因子模型對數(shù)據(jù)處理，計算出兩點云的因子載荷矩陣，可以將特殊因子協(xié)方差矩陣看作是噪聲項。所以降噪后的點云協(xié)方差矩陣為：

(6)

ΛP=ΛQ

(7)

綜上所述，可得：

(8)

為了計算出平移向量t需要先計算出點云P和點云Q的均值點，即：

(9)

進一步得到：

(10)

2 參數(shù)求解

2.1 極大似然估計

對樣本X={x1,x2,…,xN}進行極大似然估計：

(11)

2.2 EM迭代

為了求解2.1節(jié)的最大值，這里采用最大期望算法(Expectation Maximization Algorithm，EM)對上述問題進行求解，EM算法分為兩個部分。E步：計算完整數(shù)據(jù)的對數(shù)似然函數(shù)的期望；M步：通過將中間量最大化的值來獲得新的參數(shù)值。

由2.1節(jié)可知為了估算旋轉(zhuǎn)矩陣，需要求解的參數(shù)為AP、AQ、φP、φQ，在EM算法下需要計算每個觀測值的概率屬性γi，同樣為了計算γi需要假設(shè)一個中間隱藏變量fi，計算結(jié)束之后可消掉。則代價函數(shù)為：

(12)

式中：θ={μ,A,φ}表示待求解參數(shù)；fi是第i個中間隱藏變量且f～N(O,I)；γi(fi)是第i個中間隱藏變量的觀測值。E步，計算γi(fi)，由高斯邊緣分布條件可得X和z的聯(lián)合分布為：

(13)

因此，有：

(14)

根據(jù)多元高斯分布公式，得到：

(15)

M步,將代價函數(shù)進一步表示為：

(16)

對A求偏導(dǎo)：

(17)

令式(17)等于0，可得：

(18)

由式(13)可知：

(19)

綜上所述：

(20)

求解過程中μ值不變。同理，用此方法對φ求偏導(dǎo)，省略推導(dǎo)過程，得到推導(dǎo)結(jié)果：

A(μfi|xi(μfi|xi)T+∑fi|xi)AT]

(21)

本文算法流程如圖1所示。

圖1 本文算法流程

本文算法是在MATLAB 2017a版本下實現(xiàn)的，具體步驟如下：

Step1通過式(1)將點云P和點云Q表示為正交因子模型。

%E步：

for i=1:N

muf(:,i)=A′/(A*A′+phi)*(X(:,i)-mu);

%更新每一個中間均值變量muf

Ef(:,:,i)=eye(K)-A′/(A*A′+phi)*A;

%更新每一個中間協(xié)方差變量Ef

end

Step3由式(20)分別更新點云P和點云Q的正交因子載荷矩陣AP、AQ和殘差矩陣φP、φQ,具體代碼如下：

%M步：

%更新A

A1=zeros(D,K);

%初始化中間變量矩陣為D×K的矩陣

A2=zeros(K,K);

%初始化中間變量矩陣為K×K的矩陣

for i=1:N

A1=A1+(X(:,i)-mu)*muf(:,i)′;

%計算(X(:,i)-mu)*muf(:,i)′的累加和A1

A2=A2+muf(:,i)*muf(:,i)′+Ef(:,:,i);

%計算muf(:,i)*muf(:,i)′+Ef(:,:,i)的累加和A2

end

A=A1/A2;

%利用累加和A1和A2更新因子載荷矩陣A更新phi

P=zeros(D,D);

%中間殘差矩陣初始化為D×D的矩陣

for i=1:N

P=P+X(:,i)*X(:,i)′-X(:,i)*muf(:,i)′*

A′-A*muf(:,i)*X(:,i)′+A*(muf(:,i)*muf(:,i)′+

Ef(:,:,i))*A′;

%采用累加和的方式更新中間殘差矩陣P

end

phi=diag(P/N);%取1/N倍中間殘差矩陣的對角線元素

%得到1×D的向量

phi=diag(phi);%將1×D的向量轉(zhuǎn)化為對角線殘差矩陣，

%完成對對角線殘差矩陣phi的更新

Step4重復(fù)Step2和Step3直至收斂，由式(6)-式(8)利用AP、AQ估算旋轉(zhuǎn)矩陣R。

Step6由式(10)估算出平移向量t。

3 實驗與仿真分析

仿真主要是通過對來自斯坦福大學(xué)的公開點云數(shù)據(jù)Bunny(分辨率：31 607個點)和Horse(分辨率：48 485個點)三維點云數(shù)據(jù)進行仿真配準(zhǔn)，和對來自三維激光掃描儀HandySCAN 700分別從正面與側(cè)面對現(xiàn)場實際物體進行掃描的實驗數(shù)據(jù)進行配準(zhǔn)，噪聲為高斯白噪聲。三維激光掃描儀HandySCAN 700配置參數(shù)如下：

產(chǎn)品用途：行業(yè)掃描

掃描元件：CCD

產(chǎn)品類型為3D掃描儀

掃描范圍：自由掃描為0.1～4.0 m

掃描速度：480 000次測量/秒

產(chǎn)品尺寸：122×77×294 mm

產(chǎn)品重量：0.85 kg

使用條件如下：

電源類型：100～220 V

工作溫度：5℃～40℃

工作濕度：10%～90%

使用方式：手持

激光類別：2級

精度：最高0.03 mm

實驗預(yù)仿真環(huán)境為：MATLAB 2017a版本、i7-6700HQ四核處理器和GTX965M。為了確保配準(zhǔn)精度，實驗與仿真是在不能受到有其他軟件的影響的條件下進行的。

3.1 經(jīng)典配準(zhǔn)

通過對Bunny和Horse點云的隨機旋轉(zhuǎn)和平移得到目標(biāo)點云，然后對Bunny目標(biāo)點云做隨機丟失25%處理，對Horse的目標(biāo)點云添加50 dB的高斯白噪聲。點云初始狀態(tài)如圖2所示。

圖2 Bunny和Horse點云的初始狀態(tài)

分別利用本文算法、CPD算法、Go-ICP算法對三維點云數(shù)據(jù)配準(zhǔn)，比較三種算法的配準(zhǔn)精度和時間，配準(zhǔn)效果如圖3所示。

圖3 三種算法對Bunny和Horse點云的配準(zhǔn)效果

三種算法配準(zhǔn)精度和配準(zhǔn)時間如表1所示。

表1 三種算法的配準(zhǔn)時間和配準(zhǔn)誤差

由圖3和表1可知，CPD算法對Bunny和Horse點云配準(zhǔn)在時間上明顯多于本文算法和Go-ICP算法。在配準(zhǔn)精度方面，本文算法對Bunny的配準(zhǔn)誤差要略大于Go-ICP算法，但對Horse的配準(zhǔn)精度更優(yōu)于其他兩種算法。

3.2 實物配準(zhǔn)

為證明本文算法具有實用性，采用三維激光掃描儀HandySCAN 700分別從正面與側(cè)面對兩組物體進行掃描，從正面與側(cè)面分別掃描的目的是為了得到不同視覺下的點云數(shù)據(jù)集。得到Bottle和Banana點云的初始狀態(tài)，如圖4所示。其中：Bottle源點云(右)分辨率為21 450個點，目標(biāo)點云(左)分辨率為21 430 個點；Banana源點云(左)分辨率為45 271個點，目標(biāo)點云(右)分辨率為46 370 個點。

(a) Bottle (b) Banana

將本文算法與CPD算法和Go-ICP算法的配準(zhǔn)效果進行對比，結(jié)果如圖5所示。

(a) (b) (c)

三種算法對Bottle和Banana點云的配準(zhǔn)時間和配準(zhǔn)誤差如表2所示。

表2 三種算法對Bottle和Banana點云的配準(zhǔn)誤差和時間

實驗結(jié)果表明，對Bottle和Banana進行點云配準(zhǔn)時，本文算法所需的時間明顯少于CPD和Go-ICP兩種算法。三種算法在對Bottle進行點云配準(zhǔn)時的配準(zhǔn)誤差相差甚微，但本文算法對Banana配準(zhǔn)時的誤差要明顯小于其他兩種算法的精度。

3.3 現(xiàn)場掃描機械點云配準(zhǔn)

為充分驗證本文算法對實物掃描數(shù)據(jù)進行配準(zhǔn)的有效性,用便攜式激光掃描儀 (HandySCAN 700TM) 對現(xiàn)場機械器件從正面與側(cè)面分別進行掃描，從正面與側(cè)面分別掃描的目的是為了得到不同視覺下的機械點云數(shù)據(jù)集，進行實物數(shù)據(jù)采集，然后對掃描數(shù)據(jù)做相應(yīng)處理后，再進行配準(zhǔn)。機械器件目標(biāo)點云(右)分辨率28 437個點，和待配準(zhǔn)點云(左)分辨率21 328個點，初始狀態(tài)如圖6所示。

圖6 機械器件的初始狀態(tài)

將本文算法與CPD算法和Go-ICP算法分別配準(zhǔn)，配準(zhǔn)效果如圖7所示。

圖7 配準(zhǔn)效果圖

三種算法的配準(zhǔn)時間和精度如表3所示。

表3 三種算法的配準(zhǔn)時間和配準(zhǔn)誤差

由圖6和表3可知，對機械器件進行點云配準(zhǔn)時，本文算法在配準(zhǔn)時間上與Go-ICP算法相差甚微且遠(yuǎn)優(yōu)于CPD算法，但在配準(zhǔn)精度方面，本文算法明顯優(yōu)于其他兩種算法。

綜合分析發(fā)現(xiàn)，本文算法具有較好的實用性。該算法的應(yīng)用前景主要是在歷史文物保護、工程測量、醫(yī)學(xué)圖像重建等方面。

4 結(jié) 語

本文針對數(shù)據(jù)存在白噪聲、數(shù)據(jù)缺失、無序，以及仿射情況下的三維點云提出一種基于主成分因子法的點云配準(zhǔn)算法。將帶白噪聲的點云數(shù)學(xué)模型擴展為正交因子模型，從而將點云配準(zhǔn)問題轉(zhuǎn)化為對模型參數(shù)的求解；通過運用極大似然法得到因子載荷矩陣，進一步精確化模型參數(shù)，從而實現(xiàn)點云配準(zhǔn)。實驗表明，本文算法有效提升了配準(zhǔn)精度和效率；與Go-ICP和CPD兩種算法對比，本文算法的穩(wěn)定性較高。