俞 浩, 高 飛

(合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥 230009)

0 引 言

三維激光掃描技術(shù)是20世紀(jì)90年代中葉出現(xiàn)的新興技術(shù),它可以不接觸被測物體得到被測物體表面大量的密度高的三維點信息以及反射率等信息,由于三維激光掃描的快速性、準(zhǔn)確性,因此該技術(shù)又被稱為“實景復(fù)制技術(shù)”。該技術(shù)被廣泛運用在逆向工程、虛擬現(xiàn)實、古建筑保護(hù)與修復(fù)等領(lǐng)域。

通常在外業(yè)采集數(shù)據(jù)時,往往對某一物體不能通過一次掃描得到完整的三維點云數(shù)據(jù),因此在實際作業(yè)時,需要從不同角度對物體進(jìn)行數(shù)據(jù)采集,因為在不同角度下進(jìn)行測量時坐標(biāo)系不同,所以必須要通過一定的方法對將不同角度下的點云數(shù)據(jù)進(jìn)行點云配準(zhǔn),即將不同坐標(biāo)系下的點云數(shù)據(jù)轉(zhuǎn)換到同一個坐標(biāo)系之下,從而得到物體完整的三維點云數(shù)據(jù)。常見的粗配準(zhǔn)方法有七參數(shù)轉(zhuǎn)換法和四元數(shù)法。

七參數(shù)轉(zhuǎn)換法就是要求得兩點云所在坐標(biāo)系之間的3個平移量、3個旋轉(zhuǎn)量以及1個尺度變換參數(shù)。轉(zhuǎn)換公式如下:

X=ΔX+(1+k)RP

(1)

R=R(α)R(β)R(γ)

(2)

(3)

(4)

(5)

由(1)式、(5)式可得:

(6)

其中,X、Y、Z為目標(biāo)點云坐標(biāo);Δx、Δy、Δz為3個平移參數(shù);k為尺度變換參數(shù);α、β、γ為3個旋轉(zhuǎn)參數(shù);px、py、pz為源點云坐標(biāo)。

由(6)式可知，得到7個轉(zhuǎn)換參數(shù)就可以實現(xiàn)源點云坐標(biāo)到目標(biāo)點云的轉(zhuǎn)化。

四元數(shù)是Hamilton在1843年找到的一種擴(kuò)展的復(fù)數(shù),于1986年由Faugeras和Hebert用四元數(shù)進(jìn)行配準(zhǔn)研究[1]。一個四元數(shù)h擁有1個實部和3個虛部,即

h=q0+q1i+q2j+q3k

(7)

其中,i、j、k為四元數(shù)的3個虛部單位。四元數(shù)是一個包含4個矢量的陣列:

(8)

(9)

(10)

定義源點云與目標(biāo)點云的質(zhì)心為:

(11)

兩點云之間的協(xié)方差矩陣為:

(12)

由協(xié)方差矩陣構(gòu)成的4×4對稱矩陣為:

(13)

計算Q(R)的特征向量與特征值,其最大特征值對應(yīng)的特征向量為qR=[q0q1q2q3]T。

qT=μx-R(qR)μp

(14)

本文通過引入特征向量的歐氏距離來解決兩片點云間的特征匹配問題,通過七參數(shù)轉(zhuǎn)換法對點云坐標(biāo)進(jìn)行轉(zhuǎn)換,確保在進(jìn)行最近點迭代法前兩點云有著比較好的初始位置。本文通過應(yīng)用實例驗證該方法的可行性。

1 特征點匹配

要實現(xiàn)不同坐標(biāo)系下的點云配準(zhǔn),就要求解出不同坐標(biāo)系之間的轉(zhuǎn)換參數(shù)(旋轉(zhuǎn)矩陣R和平移矩陣T),那么就必須要對不同坐標(biāo)系下的相同點進(jìn)行匹配。

目前特征點匹配主要分為2類:① 通過控制點來進(jìn)行坐標(biāo)系之間的轉(zhuǎn)換[2];② 根據(jù)點云重疊部分的點以及點與點之間的固有關(guān)系來進(jìn)行配準(zhǔn)。這種配準(zhǔn)方法主要有求解點與點之間的配準(zhǔn)算法[3-4]、求解點與平面之間距離的配準(zhǔn)算法、計算豪斯多夫距離來衡量點云相似度以及基于幾何特征的配準(zhǔn)算法[5-6]。著名的最近點迭代法就是通過求解點與點之間的歐式距離的最小值來進(jìn)行點對匹配,文獻(xiàn)[7]通過計算點鄰域的切平面與點之間的距離來進(jìn)行匹配,但是這種方法的魯棒性不是太好。文獻(xiàn)[8]進(jìn)行了一些豪斯多夫距離在點云配準(zhǔn)的理論研究,但未進(jìn)行實際配準(zhǔn)應(yīng)用。

總的來說,第2種配準(zhǔn)步驟主要分為尋找對應(yīng)點和求解對應(yīng)關(guān)系2步。

尺度不變特征變換(scale-invariant feature transform,SIFT)算法可以將不同坐標(biāo)系下的點云里的特征點提取出來。SIFT算法是David Lowe在1999年提出,2004年整理完善的[9],最初是一種根據(jù)電腦視覺的算法來偵測與描述影像的局部特征點,這些局部特征點可以幫助識別物體。該算法同樣可以運用到點云數(shù)據(jù)處理中來,其實質(zhì)是計算出特征點的方向。它所找出來的關(guān)鍵點其實是一些在點云數(shù)據(jù)中比較突出的點,這些點不會因視角、光照等因素而改變。

在將特征點提取出來后,就可以用KDTree對這些特征點進(jìn)行其法向量的求解。通過對其鄰域點集求法向量來近似代替其法向量[10]。對鄰域點的三維坐標(biāo)進(jìn)行中心化,求得協(xié)方差矩陣并對角化,求得3個特征值、最小特征值對應(yīng)的特征向量即為法向量。

(15)

C·nj=λj·nj

(16)

通過對不同坐標(biāo)系下的點云數(shù)據(jù)P的所有特征向量ni(i=1,2,3,…,n)和點云數(shù)據(jù)Q的特征向量nj(j=1,2,3,…,n)求解歐式距離,若最小歐式距離小于一個閾值ε,則這2個特征向量所對應(yīng)的特征點即被認(rèn)為是在不同坐標(biāo)系下的相同點；若最小歐式距離大于這個閾值ε,則不進(jìn)行匹配,進(jìn)行下一個特征向量的最小歐式距離的求解。

運用求解歐式距離對2個不同坐標(biāo)系下的特征點進(jìn)行匹配的具體流程如圖1所示。

圖1 計算歐式距離匹配點對流程

2 點云粗配準(zhǔn)

用匹配好的兩組不同坐標(biāo)系下的點云數(shù)據(jù)進(jìn)行七參數(shù)求解[11]，即

A=(1+k)RB+T

(17)

加上最小二乘的約束條件,得出一個最符合的七參數(shù)。其中

(17)式中,A為目標(biāo)點云;B為源點云;R為旋轉(zhuǎn)矩陣;T為平移矩陣;k為縮放比例。由于點云坐標(biāo)式基于掃描儀的極坐標(biāo)系求出,因此縮放比例k可以忽略。

(17)式可以變?yōu)?

A=RB+T

(18)

結(jié)合最小二乘的約束條件：

VTPV=min

(19)

可以得到一個最為準(zhǔn)確的旋轉(zhuǎn)矩陣R和平移矩陣T。

3 點云精配準(zhǔn)

運用粗配準(zhǔn)求得的旋轉(zhuǎn)矩陣R和平移矩陣T可以將源點云數(shù)據(jù)轉(zhuǎn)換成目標(biāo)點云附近,實現(xiàn)2片點云的粗配準(zhǔn)。但是不同點云之間會有分叉的情況出現(xiàn),基于這種情況,引入最近點迭代法對點云位置進(jìn)行精配準(zhǔn)。

最近點迭代法(ICP算法)最先由Besl和Mckay提出[12-13],其本質(zhì)就是在最小二乘的約束條件下,求出2個點集之間的變換關(guān)系,假設(shè){Pi,i=1,2,3,…,n}為源點集,{Qi,i=1,2,3,…,m}為目標(biāo)點集,則構(gòu)建函數(shù)為：

(20)

通過迭代確定這2個點集的最優(yōu)匹配,直到達(dá)到設(shè)定的某個收斂準(zhǔn)則為止,也就是通過不停地求解這2個點集的變換關(guān)系，并且根據(jù)這個變換更新2個點集的位置,直到變換關(guān)系滿足這個收斂準(zhǔn)則,這時2個點集之間的變換關(guān)系就認(rèn)為是最優(yōu)解，即旋轉(zhuǎn)矩陣R與平移矩陣T最優(yōu),其具體步驟為:

(1) 首先從源點云及目標(biāo)點云中分別提取點集Pi∈P,Qi∈Q。

(2) 計算源點集中的點與目標(biāo)點集中的點的對應(yīng)關(guān)系,使得‖Qi-Pi‖2=min。

(3) 計算旋轉(zhuǎn)矩陣R與平移矩陣T,滿足關(guān)系如下：

(21)

(4) 運用旋轉(zhuǎn)矩陣R與平移矩陣T變換P點云的位置,即Pi′={Pi′=RPi+T}。

(5) 計算均方差δ,具體公式如下：

(22)

(6) 比較δ與最初設(shè)定的閾值ε。若δ<ε,則停止計算；若δ>ε,則重復(fù)步驟(2)～步驟(6)。

這樣經(jīng)過每一次迭代,Pi就會更靠近點云Qi,基于這種思想,Besl證明了這種方法的收斂性。

4 應(yīng)用實例分析

實驗數(shù)據(jù)為2個不同坐標(biāo)系下的某古文物鼎的點云數(shù)據(jù)，初始點云如圖2所示,左側(cè)點云為源點云,右側(cè)點云為目標(biāo)點云,目的是要將左側(cè)點云向右側(cè)點云轉(zhuǎn)化。先用SIFT算法提取出特征點后再對特征點進(jìn)行法向量求解，并進(jìn)行點對配對,得到的粗配準(zhǔn)的旋轉(zhuǎn)矩陣R和平移矩陣T分別為:

圖2 兩點云初始位置

粗配準(zhǔn)的結(jié)果如圖3所示。用最近點迭代法進(jìn)行迭代50次實現(xiàn)精配準(zhǔn),再迭代的過程中每一次迭代都會生成一個旋轉(zhuǎn)矩陣及平移矩陣,相對應(yīng)的源點云也就一次一次地更加逼近目標(biāo)點云。其中最后一次的旋轉(zhuǎn)矩陣和平移矩陣分別為:

圖3 粗配準(zhǔn)后的點云位置

將初始點云直接進(jìn)行ICP算法陷入的局部最優(yōu)的情況如圖4所示。

圖4 直接運用ICP算法陷入的局部最優(yōu)情況

直接進(jìn)行ICP算法每一次迭代的具體誤差值以及運行時間如圖5所示,經(jīng)過粗配準(zhǔn)之后再運用ICP算法進(jìn)行精配準(zhǔn)后的效果如圖6所示。

圖6 經(jīng)過ICP算法迭代后的兩片點云

從圖4可以看出,直接進(jìn)行ICP算法進(jìn)行點云配準(zhǔn)會導(dǎo)致點云并不能很好地拼接在一起,甚至?xí)霈F(xiàn)局部最優(yōu)的情況,如圖5所示,雖然迭代后的誤差值逐漸趨于收斂,但是實際位置并沒有很好地配準(zhǔn)。因此本文的粗配準(zhǔn)算法有一定的必要性,通過粗配準(zhǔn)將兩點云調(diào)整至大致重疊的位置,它們之間僅僅相差較小的旋轉(zhuǎn)和平移。

圖5 直接ICP算法運行時間及誤差值

從圖6可以看出,經(jīng)過粗配準(zhǔn)之后再進(jìn)行點云數(shù)據(jù)精配準(zhǔn),兩片點云幾乎完全重疊,迭代50次后均方誤差達(dá)到0.000 111 194 m,在運行時間方面,文中程序所用的計算機(jī)配置為Intel(R)Core(TM)i5-7500處理器,16 GB內(nèi)存,直接ICP算法花費時間為63 289 ms,經(jīng)過粗配準(zhǔn)之后再進(jìn)行ICP算法的花費時間為48 011 ms。具體差異見表1所列。

表1 2種方法差異

5 結(jié) 論

由于多站數(shù)據(jù)的坐標(biāo)系不同,初始點云的位置差異比較大,直接進(jìn)行ICP算法無法準(zhǔn)確配準(zhǔn)拼接,甚至?xí)萑刖植孔顑?yōu)的問題。而僅僅用四元數(shù)法或七參數(shù)轉(zhuǎn)換坐標(biāo)進(jìn)行粗配準(zhǔn)并不能達(dá)到所要的精度。本文針對這個問題提出新的解決方法,即先利用特征點對應(yīng)的特征向量尋找點對并求解七參數(shù)的方法進(jìn)行粗配準(zhǔn),將點云位置之間的剛性變換變小,再通過ICP算法對點云數(shù)據(jù)進(jìn)行精確配準(zhǔn)。這樣既可以解決直接進(jìn)行ICP算法陷入的局部最優(yōu)問題,也可以達(dá)到所需要的精度。