基于核密度估計(jì)的點(diǎn)云魯棒配準(zhǔn)算法

林洪彬 劉 彬 張玉存

1.燕山大學(xué)，秦皇島，066004

2.河北省測(cè)試計(jì)量技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室，秦皇島，066004

0 引言

點(diǎn)云配準(zhǔn)是計(jì)算機(jī)視覺(jué)領(lǐng)域的基本問(wèn)題之一，廣泛應(yīng)用于距離數(shù)據(jù)融合、醫(yī)學(xué)圖像對(duì)準(zhǔn)、目標(biāo)定位、跟蹤及識(shí)別［1］。ICP算法是經(jīng)典的點(diǎn)云配準(zhǔn)算法，然而，ICP算法存在如下不足：①代價(jià)函數(shù)的不可微性要求點(diǎn)云之間有足夠的重疊區(qū)域，否則算法容易陷入局部極?。虎贗CP算法需要一個(gè)緊的初始條件。上述兩點(diǎn)要求在應(yīng)用ICP算法之前必須通過(guò)一定的預(yù)處理方法估計(jì)點(diǎn)云之間的重疊區(qū)域，并進(jìn)行粗配準(zhǔn)，以獲得ICP算法的初始值。點(diǎn)云配準(zhǔn)的難點(diǎn)在于異常值、背景和噪聲等對(duì)配準(zhǔn)精度的影響。為了消除上述因素對(duì)配準(zhǔn)算法精度和魯棒性的影響，文獻(xiàn)［2－3］相繼對(duì)ICP算法進(jìn)行了改進(jìn)。除此之外，人們還提出了兩類(lèi)方法對(duì)配準(zhǔn)算法的魯棒性進(jìn)行改善。第一類(lèi)方法通過(guò)對(duì)異常值、背景和噪聲進(jìn)行統(tǒng)計(jì)學(xué)建模，進(jìn)而采用 ML/MAP估計(jì)的方法進(jìn)行配準(zhǔn)。如Hasler等［4］研究了圖像配準(zhǔn)中的異常值建模問(wèn)題；Ying等［5］研究了基于統(tǒng)計(jì)模型的局部重疊對(duì)象識(shí)別問(wèn)題。然而，這類(lèi)方法的問(wèn)題在于，在許多實(shí)際應(yīng)用場(chǎng)合，對(duì)干擾因素的精確統(tǒng)計(jì)學(xué)建模是不可行的。第二類(lèi)方法采用變換不變的魯棒特征結(jié)合代價(jià)函數(shù)進(jìn)行配準(zhǔn)，以獲得針對(duì)異常值、背景和噪聲等干擾魯棒性更好的配準(zhǔn)算法。如Viola等［6］提出了一種基于融合信息的配準(zhǔn)算法；Tsin等［7］提出了一種基于統(tǒng)計(jì)相關(guān)的魯棒點(diǎn)云配準(zhǔn)算法（KC算法）；Huttenlocher等［8］提出了一種基于局部Hausdorff距離的配準(zhǔn)算法；Chui等［9］提出了一種基于混合模型的特征配準(zhǔn)算法；Jian等［10］提出了一種基于高斯混合模型的魯棒配準(zhǔn)算法（GMM（gaussian mixture model）算法）等。楊靜等［11］提出了一種基于簡(jiǎn)單Schur凹函數(shù)的圖像配準(zhǔn)測(cè)度；李暉等［12］提出了塔式分解和模糊梯度場(chǎng)的圖像配準(zhǔn)算法。這類(lèi)算法的問(wèn)題在于不同的代價(jià)函數(shù)和優(yōu)化方法對(duì)算法性能影響較大。因此，設(shè)計(jì)合適的代價(jià)函數(shù)并根據(jù)代價(jià)函數(shù)的特點(diǎn)探尋全局合適的優(yōu)化方法成為此類(lèi)算法的關(guān)鍵，本文即屬于此類(lèi)方法。

本文基于核密度估計(jì)理論提出了一種新的測(cè)度函數(shù)，對(duì)不同參數(shù)下測(cè)度的性能進(jìn)行分析，論證了利用該測(cè)度進(jìn)行尋優(yōu)配準(zhǔn)時(shí)存在的局部極值現(xiàn)象和極值漂移現(xiàn)象。在此基礎(chǔ)上提出了一種基于BFGS擬牛頓法的變尺度優(yōu)化配準(zhǔn)算法，解決了算法魯棒性與配準(zhǔn)精度之間的矛盾。

1 點(diǎn)云的核密度分布模型建立

設(shè)M＝｛mi，i＝1，2，…，Nm｝表示模型點(diǎn)云，S＝｛si，i＝1，2，…，Ns｝表示場(chǎng)景點(diǎn)云，兩個(gè)點(diǎn)云可以視為對(duì)目標(biāo)所對(duì)應(yīng)的實(shí)向量空間Rd進(jìn)行觀測(cè)所得到的樣本。

其中，H為d×d階正定帶寬矩陣，為簡(jiǎn)化計(jì)算，取H＝hId，Id為d×d階單位矩陣；K：Rd→R為滿(mǎn)足如下條件的非負(fù)核函數(shù)：

對(duì)于三維點(diǎn)云，d＝3。出于一般性考慮，選取高斯函數(shù)作為概率估計(jì)的核函數(shù)，即

由此，可以通過(guò)上述過(guò)程分別建立離散點(diǎn)云M和S的連續(xù)核密度分布模型。

2 相似性測(cè)度的構(gòu)造

核密度函數(shù)相似性測(cè)度對(duì)配準(zhǔn)算法性能的提高具有重要作用，通常用于描述兩個(gè)核密度分布相似性的測(cè)度有積分絕對(duì)值誤差測(cè)度、Kullback－Liebler測(cè)度、MI（mutual information）測(cè)度、歐氏測(cè)度等［10］。不同的測(cè)度函數(shù)優(yōu)化過(guò)程具有不同的魯棒性和計(jì)算效率，如歐氏距離計(jì)算方便、意義明確，在相似性評(píng)價(jià)方面得到了廣泛應(yīng)用，但相比于Kullback－Liebler測(cè)度，其魯棒性稍差。在文獻(xiàn)［10］提出的GMM測(cè)度函數(shù)的基礎(chǔ)上，本文提出了如下測(cè)度函數(shù)：

當(dāng)α→0時(shí)，有

當(dāng)α→1時(shí)，有

此時(shí)，d1（f，g）為與之間的歐氏距離。因此，利用本文提出的測(cè)度函數(shù)可以通過(guò)單一參數(shù)α可以調(diào)整算法的魯棒性和漸進(jìn)效率。

當(dāng)配準(zhǔn)對(duì)象為剛體時(shí)，空間變換矩陣T可以表示為：T（M）＝Rmi＋t，（i＝ 1，2，…，Nm），其中R為3×3階旋轉(zhuǎn)矩陣，t為3×1階平移向量。此時(shí)，

令x＝Rx′＋t，則

由式（3）可得

由于R為剛體變換矩陣，故，RRT＝I3，detR＝1，則

將式（10）代入式（8）可得

當(dāng)α＝1時(shí)，測(cè)度函數(shù)F（T）的解析表達(dá)式為

當(dāng)α＝2時(shí)，測(cè)度函數(shù)F（T）的解析表達(dá)式為

3 相似性測(cè)度的性能分析

首先，隨機(jī)生成一個(gè)3×100的服從正態(tài)分布的隨機(jī)矩陣作為場(chǎng)景點(diǎn)云M。其次，對(duì)場(chǎng)景點(diǎn)云M沿x軸方向平移3個(gè)單位，并隨機(jī)選取其中50列作為模型點(diǎn)云S。h分別取0.1，10，20，50，100時(shí)，測(cè)度函數(shù)隨偏移量tx的變化關(guān)系如圖1所示（說(shuō)明：本文除角度的量綱?。ā悖┗騬ad，其他均為量綱一量）。再采用同樣的方法測(cè)試測(cè)度函數(shù)對(duì)旋轉(zhuǎn)參數(shù)的性能，設(shè)定模型點(diǎn)云和場(chǎng)景點(diǎn)云之間的旋轉(zhuǎn)關(guān)系為θx＝π/3。此時(shí)，測(cè)度函數(shù)隨旋轉(zhuǎn)量θx的變化關(guān)系如圖2所示。

圖1 不同參數(shù)下測(cè)度函數(shù)與偏移量tx關(guān)系曲線(xiàn)

圖2 不同參數(shù)下測(cè)度函數(shù)與繞x軸旋轉(zhuǎn)角θx關(guān)系曲線(xiàn)

由式（14）可知，α＝1的情況下本文測(cè)度等價(jià)于文獻(xiàn)［7］和文獻(xiàn)［14］所采用的測(cè)度。圖1表明，針對(duì)平移變換，通過(guò)搜索測(cè)度函數(shù)的極大值可以確定平移參量，實(shí)現(xiàn)點(diǎn)云的平移配準(zhǔn)。圖1a表明：α＝1時(shí)，尺度參數(shù)越大測(cè)度函數(shù)曲線(xiàn)越寬，曲線(xiàn)越光滑，峰值位置偏離真實(shí)值，采用大尺度參數(shù)的測(cè)度函數(shù)進(jìn)行配準(zhǔn)時(shí)將存在殘余誤差，我們把這種現(xiàn)象稱(chēng)為過(guò)尺度參數(shù)導(dǎo)致的極值漂移現(xiàn)象；尺度參數(shù)越小測(cè)度函數(shù)曲線(xiàn)越窄，曲線(xiàn)可能出現(xiàn)局部極值點(diǎn)，采用小尺度參數(shù)的測(cè)度函數(shù)進(jìn)行配準(zhǔn)時(shí)算法可能因陷入局部極小而導(dǎo)致失配，我們把這種現(xiàn)象稱(chēng)為欠尺度參數(shù)導(dǎo)致的局部極值現(xiàn)象。圖1b表明，當(dāng)選取α＝2時(shí)，測(cè)度函數(shù)曲線(xiàn)在極值漂移和局部極值方面性能更好，采用此測(cè)度函數(shù)的配準(zhǔn)算法具有更好的魯棒性和配準(zhǔn)精度。圖2表明，在進(jìn)行旋轉(zhuǎn)參數(shù)的配準(zhǔn)時(shí)，選取參數(shù)α＝1也會(huì)導(dǎo)致極大值漂移和局部極值的問(wèn)題，而選取參數(shù)α＝2可以改善這兩種現(xiàn)象，提高算法的魯棒性和準(zhǔn)確性。

為了檢驗(yàn)測(cè)度函數(shù)極值點(diǎn)的統(tǒng)計(jì)穩(wěn)定性，將實(shí)驗(yàn)重復(fù)進(jìn)行100次，統(tǒng)計(jì)各種測(cè)度函數(shù)全局極大值所在位置的均值μ和標(biāo)準(zhǔn)差σ，實(shí)驗(yàn)結(jié)果如表1和表2所示。

表1 測(cè)度函數(shù)對(duì)平移參數(shù)估計(jì)結(jié)果統(tǒng)計(jì)表

表2 測(cè)度函數(shù)對(duì)旋轉(zhuǎn)參數(shù)估計(jì)結(jié)果統(tǒng)計(jì)

對(duì)比三種測(cè)度可以發(fā)現(xiàn)，與MI測(cè)度相比，本文提出的測(cè)度和GMM測(cè)度對(duì)于平移參數(shù)和旋轉(zhuǎn)參數(shù)的估計(jì)均具有更好的準(zhǔn)確性和有效性。通過(guò)對(duì)比GMM測(cè)度實(shí)驗(yàn)結(jié)果和本文測(cè)度實(shí)驗(yàn)結(jié)果可以發(fā)現(xiàn)，當(dāng)α＝1時(shí)，本文測(cè)度的估計(jì)性能與GMM測(cè)度估計(jì)性能相同，而當(dāng)α＝2時(shí)，測(cè)度對(duì)參數(shù)估計(jì)的性能優(yōu)于GMM測(cè)度的估計(jì)性能。

4 基于BFGS優(yōu)化的變尺度配準(zhǔn)算法

通過(guò)表1和表2可以發(fā)現(xiàn)，應(yīng)用本文提出的測(cè)度可以改善配準(zhǔn)參數(shù)估計(jì)的準(zhǔn)確性和有效性。然而，對(duì)比不同尺度參數(shù)下GMM測(cè)度和本文測(cè)度的性能可以發(fā)現(xiàn)，隨著尺度的增大，對(duì)平移參數(shù)和旋轉(zhuǎn)參數(shù)估計(jì)的偏倚現(xiàn)象明顯，導(dǎo)致參數(shù)估計(jì)準(zhǔn)確度降低。圖1和圖2也表明，尺度參數(shù)過(guò)小可能會(huì)導(dǎo)致局部極值現(xiàn)象的產(chǎn)生，從而導(dǎo)致點(diǎn)云失配。為此，本文提出了一種基于BFGS優(yōu)化的變尺度配準(zhǔn)算法，其優(yōu)化過(guò)程如圖3所示。

首先，隨機(jī)選取初始參數(shù)T0作為極值點(diǎn)搜索的初始條件，同時(shí)，為了保證算法的收斂域和平滑性，選取一個(gè)較大的初始參數(shù)h0。再利用BFGS擬牛頓法搜索該尺度下測(cè)度函數(shù)極大值所在位置T1，在通過(guò)某種機(jī)制減小尺度參數(shù)h并以T1為搜索初始條件，獲得T2依次類(lèi)推，直至達(dá)到預(yù)定的配準(zhǔn)精度。

圖3 變尺度優(yōu)化過(guò)程示意圖

基于BFGS優(yōu)化方法的變尺度配準(zhǔn)算法流程如下：

（1）初始化。h0＝3λmax和T0＝［00000 0］T，并令k＝0，其中，λmax為兩個(gè)點(diǎn)云自相關(guān)矩陣的最大特征值，I4表示4×4階單位矩陣。

（2）更新變換矩陣T。以Tk為初始條件，以hk為尺度，利用BFGS擬牛頓法搜索測(cè)度函數(shù)F（Tk，hk）的極大值。并將此極大值所在位置記為T(mén)k＋1。

（3）更新尺度參數(shù)h。hk＋1＝μhk，選取μ∈（0.9，0.96），k←k＋1。

（4）循環(huán)步驟（2）、（3），直至 ‖Tk＋1－Tk‖ ＜ε。

采用這種衰減尺度的BFGS擬牛頓搜索法進(jìn)行配準(zhǔn)主要有三個(gè)優(yōu)點(diǎn)：① 拓展了算法的收斂域；②克服了欠尺度參數(shù)導(dǎo)致的局部極值現(xiàn)象，減小了算法失配的可能；③克服了過(guò)尺度參數(shù)導(dǎo)致的極大值漂移現(xiàn)象，提高了配準(zhǔn)精度。

5 實(shí)驗(yàn)研究

為了驗(yàn)證算法的性能，本節(jié)分別通過(guò)仿真信號(hào)和實(shí)測(cè)信號(hào)的配準(zhǔn)實(shí)驗(yàn)對(duì)算法的可靠性和魯棒性進(jìn)行分析和比較。

5.1 受白噪聲干擾點(diǎn)云配準(zhǔn)實(shí)驗(yàn)（實(shí)驗(yàn)1）

（1）模型點(diǎn)云M的構(gòu)造。通過(guò)如下公式構(gòu)造一條三維空間曲線(xiàn)，并在該曲線(xiàn)均勻選取200個(gè)點(diǎn)作為模型點(diǎn)云M：

（2）場(chǎng)景點(diǎn)云S的構(gòu)造。對(duì)模型點(diǎn)云M進(jìn)行旋轉(zhuǎn)和平移變換，變換參數(shù)

再增加信噪比為20dB的白噪聲作為干擾，由此構(gòu)造場(chǎng)景點(diǎn)云S。原始點(diǎn)云以及通過(guò)ICP、GMM和本文方法的配準(zhǔn)結(jié)果如圖4所示。

圖4表明，在給定實(shí)驗(yàn)條件下ICP算法、GMM算法和本文算法均可實(shí)現(xiàn)點(diǎn)云的精確配準(zhǔn)，說(shuō)明這幾種方法相對(duì)與白噪聲干擾都具有魯棒性。

圖4 帶噪聲點(diǎn)云配準(zhǔn)實(shí)驗(yàn)結(jié)果

5.2 配準(zhǔn)算法收斂性能實(shí)驗(yàn)（實(shí)驗(yàn)2）

為了檢驗(yàn)算法的收斂性能，設(shè)計(jì)了實(shí)驗(yàn)2，實(shí)驗(yàn)中模型點(diǎn)云M的生成方法和實(shí)驗(yàn)1相同。在此基礎(chǔ)上分別改變繞x軸方向轉(zhuǎn)角θx∈（－π，π）和沿x方向平移tx∈（－20，20），并附加20dB白噪聲作為場(chǎng)景點(diǎn)云S。分別用ICP算法，GMM算法和本文算法（α＝1，μ＝0.9）進(jìn)行配準(zhǔn)，求取各種算法的配準(zhǔn)誤差。配準(zhǔn)誤差隨變換參數(shù)變化關(guān)系曲線(xiàn)如圖5所示。

圖5 配準(zhǔn)誤差隨變換參數(shù)變化關(guān)系曲線(xiàn)

觀察圖5a可以發(fā)現(xiàn)，在針對(duì)單純平移參數(shù)的配準(zhǔn)中，ICP算法和本文算法在（－20，20）范圍內(nèi)可獲得全局收斂的配準(zhǔn)效果，而GMM算法的收斂范圍為（－4.3，4.3）。這說(shuō)明在處理單一平移參數(shù)的配準(zhǔn)時(shí)，GMM算法要求兩個(gè)點(diǎn)云之間的相對(duì)平移量較小，限制了算法的實(shí)際應(yīng)用。由圖5b可知，針對(duì)單純旋轉(zhuǎn)變換，ICP算法的收斂區(qū) 間 為（－85°，85°），GMM 算 法 的 收 斂 區(qū) 間為（－90°，105°），而 本 文 算 法 的 收 斂 區(qū) 間為（－180°，170°）。收斂區(qū)間越寬，說(shuō)明使用該算法時(shí)兩個(gè)點(diǎn)云之間的相對(duì)變換參數(shù)的允許范圍就越大，從而算法的適用性就越好。綜合分析圖5可以發(fā)現(xiàn)，不論是針對(duì)旋轉(zhuǎn)參數(shù)的配準(zhǔn)還是針對(duì)平移參數(shù)的配準(zhǔn)，相對(duì)于ICP算法和GMM算法，本文提出的算法都具有更好的收斂性能。

5.3 局部重疊點(diǎn)云配準(zhǔn)實(shí)驗(yàn)（實(shí)驗(yàn)3）

為了驗(yàn)證點(diǎn)云局部重疊情況下本文算法與傳統(tǒng)算法的性能，設(shè)計(jì)了實(shí)驗(yàn)3。實(shí)驗(yàn)原始數(shù)據(jù)為一組人臉激光掃描散亂點(diǎn)云，記為PtsFace，掃描點(diǎn)數(shù)NPtsFace＝392，點(diǎn)云PtsFace的坐標(biāo)分布范圍如下：（x×y×z）∈（［－1.18，1.11］×［－0.58，1.15］×［－1.82，2.56］）。取散亂點(diǎn)PtsFace的前80%的點(diǎn)作為模型點(diǎn)云M。再對(duì)PtsFace進(jìn)行變換，參數(shù)［θxθyθztxtytz］T＝［π/5 0 0 0.2 0 0］T，取變換后點(diǎn)云的后70% 作為場(chǎng)景點(diǎn)云S。實(shí)驗(yàn)結(jié)果如圖6所示。

圖6 局部重疊點(diǎn)云配準(zhǔn)實(shí)驗(yàn)結(jié)果

圖6表明，在給定實(shí)驗(yàn)條件下，ICP算法陷入局部極小，不能實(shí)現(xiàn)給定點(diǎn)云的配準(zhǔn)。圖6c表明，當(dāng)參數(shù)選取為σ＝0.56時(shí)，GMM算法存在殘余誤差。由前面的理論分析可知，該殘余誤差是由過(guò)尺度參數(shù)引起的極值點(diǎn)漂移引起的，不能通過(guò)增加迭代次數(shù)和降低門(mén)限值的方法加以改善。

圖6d表明，當(dāng)參數(shù)選取為σ＝0.02時(shí)，GMM算法配準(zhǔn)失敗，由上文分析可知，這是由于欠尺度參數(shù)導(dǎo)致的局部極值現(xiàn)象所引起的，而這種現(xiàn)象主要存在于點(diǎn)云之間存在非重疊區(qū)域的情況。然而，在實(shí)際應(yīng)用中點(diǎn)云之間往往存在明顯的非重疊區(qū)域。因此，GMM算法在實(shí)際應(yīng)用中也受到限制。圖6e和圖6f表明，本文方法可以實(shí)現(xiàn)具有非重疊區(qū)域點(diǎn)云的精確配準(zhǔn)。

6 結(jié)束語(yǔ)

針對(duì)傳統(tǒng)點(diǎn)云配準(zhǔn)算法魯棒性差、收斂區(qū)間窄的缺點(diǎn)，基于核密度估計(jì)的思想，提出了一種新的魯棒的點(diǎn)云配準(zhǔn)算法。本文從點(diǎn)云的概率分布模型建立入手，將點(diǎn)云配準(zhǔn)問(wèn)題表述為核密度函數(shù)相似性尋優(yōu)問(wèn)題，提出了用于描述核密度分布相似性的測(cè)度函數(shù)，研究了該測(cè)度函數(shù)的性能及其與其它測(cè)度函數(shù)之間的關(guān)系。針對(duì)算法實(shí)際應(yīng)用中過(guò)尺度參數(shù)和欠尺度參數(shù)兩種情況造成的配準(zhǔn)殘余誤差和配準(zhǔn)失敗的情況，提出了一種變尺度BFGS尋優(yōu)方法，實(shí)現(xiàn)了點(diǎn)云的魯棒、精確配準(zhǔn)。

［1］Salvi J，Matabosch C，F(xiàn)ofi D，et al.A Review of Recent Range Image Registration Methods with Accuracy Evaluation［J］.Image and Vision Computing，2007，25（5）：578－596.

［2］Fitsgibbon A W.Robust Registration of 2Dand 3D Point Sets［J］.Image and Vision Computing，2001，21（13/14）：1145－1153.

［3］Granger S，Pennec X.Multi－scale EM－ICP：A Fast and Robust Approach for Surface Registration［C］//In ECCV’02.London：Springer－ Verlag，2002：418－432.

［4］Hasler D，Sviaz L，Süsstrunk S.Outlier Modeling in Image Matching［J］.PAMI，IEEE Trans.，2003，25（3）：301－315.

［5］Ying Z，Castanon D.Partially Occluded Object Recognition Using Statistical Models［J］.International Journal of Computer Vision，2002，49（1）：57－78.

［6］Viola P，WellsⅢ W M.Alignement by Maximization of Mututal Information［J］.International Journal of Computer Vision，1997，24（2）：137－154.

［7］Tsin Y，Kanade T.A Correlation－based Approach to Robust Point Set Registration［C］//ECCV’04.London：Springer－Verlag，2004：558－569.

［8］Huttenlocher D P，Klanderman G A，Gregory A K，et al.Comparing Images Using the Hausdorff Distance［J］.IEEE Trans.Pattern Anal.Mach.Intell.，1993，15（9）：850－863.

［9］Chui H，Rangarajan A.A Feature Registration Framework Using Mixture Models［C］//IEEE Workshop on Mathematical Methods in Biomedical Image Analysis，Washington：IEEE Computer Society，2000：190－197.

［10］Jian B，Vemuri B C.A Robust Algorithm for Point Set Registration Using Mixture of Gaussians［C］//ICCV’05.Beijing，2005：1246－1251.

［11］楊靜，胡順波，劉常春，等.基于簡(jiǎn)單Schur凹函數(shù)的圖像配準(zhǔn)測(cè)度研究［J］.電子學(xué)報(bào)，2008，36（12）：2328－2332.

［12］李暉，彭玉華，尹勇.基于平移旋轉(zhuǎn)不變的塔式分解和模糊梯度場(chǎng)的醫(yī)學(xué)圖像配準(zhǔn)［J］.電子學(xué)報(bào)，2009，37（4）：854－859.

［13］Silverman B W.Density Estimation for Statistics and Data Analysis［M］.London：Chapman and Hall，1986：170－180.

［14］Boughorbel F，Mercimek M.A New Method for the Registration of Three－dimensional Pointsets：The Gaussian Fields Framework［J］.Image and Vision Computing，2010，28（1）：124－137.