楊麗娟, 田錚,2, 溫金環(huán), 延偉東

(1.西北工業(yè)大學(xué) 應(yīng)用數(shù)學(xué)系, 陜西 西安 710129; 2.中國科學(xué)院 遙感科學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100101)

點(diǎn)集匹配在圖像分析、形狀匹配等許多領(lǐng)域中已獲得廣泛的應(yīng)用[1-2]，可粗略的分為剛性和非剛性，其中非剛性匹配問題備受關(guān)注。但由于點(diǎn)集中通常包含大量異常值，這嚴(yán)重影響了匹配性能，因此需要研究更穩(wěn)健的點(diǎn)集匹配算法。

近些年，基于混合高斯模型(Gaussian mixture model,GMM)的匹配方法更是受到了研究者的廣泛關(guān)注[3]。在最大似然估計(jì)框架下，利用最大化期望(expectation maximization,EM)算法求解參數(shù)的最大似然估計(jì)。Jian等[4]用混合高斯模型表示點(diǎn)集，通過最小化2個(gè)混合高斯之間的統(tǒng)計(jì)差異，進(jìn)而確定點(diǎn)集之間的對應(yīng)關(guān)系。Myronenko等[5]假設(shè)點(diǎn)集做一致性運(yùn)動，提出了基于一致性點(diǎn)漂移(coherent point drift,CPD)的點(diǎn)集匹配方法，通過增加額外的均勻分布模型化異常值。但它本質(zhì)上對應(yīng)于最小二乘方法，對異常值不穩(wěn)健。為了提高對異常值的穩(wěn)健性，一種可選的方法是選用更合適的概率模型，如具有厚尾性質(zhì)的t分布[3]。Zhou等[6]提出用混合t分布模型(student′s t mixture model,SMM)處理點(diǎn)集之間的匹配問題，通過最大化似然函數(shù)的對數(shù)期望，估計(jì)變換參數(shù)。但是它和CPD方法以及多數(shù)方法一樣，只允許各向同性協(xié)方差，而且需要手動確定正則化參數(shù)。

針對最大似然估計(jì)框架中存在的奇異性問題，可選的解決方式是貝葉斯分析[7-8]。Qu等[9]提出了在變分貝葉斯框架下的點(diǎn)集匹配方法(VBPSM)，將點(diǎn)集匹配過程分為聚類分析和回歸分析兩部分，并分別利用2個(gè)GMM模型化，通過引入轉(zhuǎn)移變量將兩部分統(tǒng)一起來。但對于非剛性變換，它沒有包含對變換函數(shù)的平滑性約束。另一個(gè)問題是估計(jì)的變換函數(shù)沒有直接反映待對準(zhǔn)點(diǎn)集之間的變換關(guān)系。在Simpson等[10]提出的非剛性配準(zhǔn)的概率配準(zhǔn)框架中，通過給變換參數(shù)指派一個(gè)合適的先驗(yàn)分布可將正則化并入概率框架中，根據(jù)數(shù)據(jù)質(zhì)量自適應(yīng)地確定正則化參數(shù)。

貝葉斯SMM已經(jīng)廣泛應(yīng)用到了密度估計(jì)、聚類和模型選擇中[11-12]。Archambeau等[12]在估計(jì)變量的后驗(yàn)分布過程中，考慮了變量之間的相關(guān)。Baldacchino等[13]將貝葉斯SMM應(yīng)用于非線性系統(tǒng)識別中的穩(wěn)健回歸。

綜上，針對存在異常值的非剛性點(diǎn)集匹配問題，本文將點(diǎn)集匹配問題表示為SMM對應(yīng)的概率密度估計(jì)問題。通過引入真實(shí)后驗(yàn)分布的變分近似，匹配問題轉(zhuǎn)化為最大化變分下界，同時(shí)允許任何形式的協(xié)方差。在參數(shù)推斷過程中，利用變分貝葉斯最大化期望(variational Bayesian EM,VBEM)算法確定變換參數(shù)。通過先驗(yàn)?zāi)Ｐ停瑢⒖臻g正則化約束并入貝葉斯SMM模型中，可自適應(yīng)地確定正則化參數(shù)。模擬點(diǎn)集和真實(shí)圖像上的實(shí)驗(yàn)結(jié)果證明所提方法可大大提高匹配性能。

1 貝葉斯混合t分布模型

(1)

式中,Θ={A,B,{Λm},{πm},{vm}},t分布的概率密度函數(shù)為:

(2)

式中，μ和Λ分別為均值向量和精度矩陣,Γ(·)為伽瑪函數(shù),當(dāng)自由度參數(shù)v→∞時(shí),t分布收斂到高斯分布。t分布可表示為無限個(gè)具有相同均值、不同精度的尺度高斯分布的混合,即:

(3)

于是,似然函數(shù)可寫為:

利用指數(shù)族分布的共軛先驗(yàn)特性,選擇模型變量的先驗(yàn):①混合比例系數(shù)π的先驗(yàn)分布服從狄利克雷分布Dir(π|κ0),其中κ0是先驗(yàn)樣本大小(稱為超參數(shù),下同);②每個(gè)混合成分的精度矩陣Λm的先驗(yàn)分布服從Wishart分布W(Λm|γ0,S0),其中γ0和S0分別為自由度和尺度矩陣;③仿射變換矩陣A和非剛性形變系數(shù)矩陣B的每列元素的先驗(yàn)分布分別服從零均值向量、精度參數(shù)υ和η的高斯分布;④精度參數(shù)向量υ和η均服從伽瑪分布,其中a0,b0和c0,d0分別是先驗(yàn)形狀和尺度參數(shù)對;⑤正則項(xiàng)以先驗(yàn)形式并入貝葉斯SMM模型中[10],即:

(8)

式中，λ>0是標(biāo)量正則化參數(shù);⑥正則化參數(shù)λ服從形狀參數(shù)為τ0和尺度參數(shù)為ζ0的伽瑪分布?？筛鶕?jù)經(jīng)驗(yàn)確定超參數(shù)的值。通過變分推斷,可獲得變量{ΘL,ΘVB}={ΘL,{A,B,Λ,π,υ,η,λ}}的后驗(yàn)分布,由于沒有自由度參數(shù)v={vm}的先驗(yàn),利用最大似然方法可確定它的值,記為ΘML={v},詳見下一節(jié)。

變量之間的聯(lián)合分布可表示為:

(9)

圖1中給出了非剛性點(diǎn)集匹配對應(yīng)的概率圖模型,即有向無環(huán)圖[7]。

圖1 非剛性點(diǎn)集匹配的概率圖模型

2 參數(shù)的變分推斷

考慮近似后驗(yàn)分布q(U,Z,ΘVB),對數(shù)似然函數(shù)可表示為:

lnp(X|Y)=L(q)+

(10)

所有變量的后驗(yàn)分布可分解為:

q(U,Z,π,Λ,A,υ,B,η)=q(U,Z)q(ΘVB)

(11)

記q(ΘVB)=q(π)q(Λ)q(A)q(υ)q(B)q(η)

q(λ)。關(guān)于q(U,Z)和q(ΘVB)分別最大化L(q)可獲得如下VBEM算法,在第k次迭代:

VBE-step:

q(U,Z)k∝exp(〈lnp(X,U,Z|Y,ΘVB)〉q(ΘVB)k-1)

(12)

VBM-step:

(13)

式中，〈·〉q(·)表示關(guān)于后驗(yàn)分布q(·)的期望,后面下標(biāo)省略。在變分貝葉斯期望步驟(VBE-step)中,固定變分變量ΘVB,計(jì)算隱變量ΘL的后驗(yàn)分布;在變分貝葉斯最大化步驟(VBM-step)中,固定隱變量ΘL,重新計(jì)算ΘVB的后驗(yàn)分布。下面列出部分變量的后驗(yàn)分布:

1) 考慮隱變量之間的相關(guān)性,其后驗(yàn)分布為:

(14)

q(znm=1)=

(15)

規(guī)范化后,示性變量元素表示為:

?n,?m

(16)

以示性變量為條件,尺度變量的后驗(yàn)分布服從伽瑪分布:q(unm|znm=1)～G(unm|αnm,βnm),其中

(17)

(18)

(19)

(20)

式中，diag(〈υ〉)表示〈υ〉中元素構(gòu)成的方陣,矩陣Ξ的元素為Ξnm:=〈znmunm〉,Λm(q)表示Λm的第q個(gè)對角元素。

(21)

(22)

4) 正則化參數(shù)λ的后驗(yàn)分布服從伽瑪分布:q(λ)～G(λ|τ,ζ),其中

(23)

(24)

根據(jù)伽瑪分布性質(zhì),有〈λ〉=τ/ζ,可根據(jù)不同的應(yīng)用自適應(yīng)地確定λ的值。

其余變量的后驗(yàn)分布見附錄。由于沒有自由度參數(shù){vm}的先驗(yàn),可直接最大化下界L(q),獲得下面非線性方程:

(25)

(26)

算法1: 基于SMM的非剛性點(diǎn)集匹配算法

2: 初始化:

先驗(yàn):a0,b0,c0,d0,κ0,γ0,S0,τ0,ζ0

變換:A,B

3: forS0上由粗到精do

4: repeat

VBE-step:

根據(jù)(15～16)推斷示性變量Z

根據(jù)(17～18)推斷尺度變量U

VBM-step:

根據(jù)(27)更新混合比例系數(shù)π的κ

根據(jù)(19～20)推斷仿射變換矩陣A

根據(jù)(30～31)推斷精度υ的a和b

根據(jù)(21～22)推斷非剛性形變系數(shù)矩陣B

根據(jù)(32～33)推斷精度η的c和d

根據(jù)(28～29)更新精度Λ的γ和S

根據(jù)(23～24)更新正則化參數(shù)λ的τ和ζ

根據(jù)(26)更新自由度v

計(jì)算第k步的L′(q)

6: end for

算法的時(shí)間復(fù)雜度主要由Z,Λ和B的后驗(yàn)更新決定。在每次迭代中,更新這些變量的時(shí)間復(fù)雜度分別近似為O(DNM2),O(D2(N+D)M)和O(DM3),一般情況下,D?N,M,故算法的總時(shí)間復(fù)雜度可近似為O(KM3)或O(KNM2),其中K為迭代次數(shù)。

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

為了驗(yàn)證本文算法的有效性,首先測試具有形狀差異的模擬數(shù)據(jù)集,然后將本文方法應(yīng)用于包含非剛性形變的圖像。并與5種算法進(jìn)行了對比:GMM-L2[4],CPD[5], SMM[6]和VBPSM (默認(rèn)各向異性協(xié)方差)[8]和TPS-RPM[15]。超參數(shù)的值設(shè)置如下:①超參數(shù)a0,b0,c0,d0,τ0和ζ0設(shè)定為0.01;②κ0的元素設(shè)定為M;③γ0設(shè)定為2×(D+1),若選擇各向異性協(xié)方差,則S0設(shè)定為單位矩陣I,否則設(shè)定為標(biāo)量1;④高斯徑向基函數(shù)的寬度參數(shù)設(shè)定為1。針對各向異性和各向同性協(xié)方差2種選擇,本文方法分別記為VBSMM(ANI)和VBSMM(ISO)。

3.1 自適應(yīng)確定正則化參數(shù)的有效性

為了檢驗(yàn)自適應(yīng)確定正則化參數(shù)的有效性,選用文獻(xiàn)[15]中公共數(shù)據(jù)集的3組存在異常值的漢字福形狀點(diǎn)集,記錄算法1在固定正則化參數(shù)時(shí)獲得的匹配召回率(recall),其中正則化參數(shù)的固定取值為0.5k,k=1,2,…,40,并與算法1在自適應(yīng)確定正則化參數(shù)時(shí)獲得的匹配召回率做對比。圖2給出了召回率關(guān)于正則化參數(shù)的變化曲線,其中空心對應(yīng)于算法1在固定正則化參數(shù)時(shí)的匹配召回率,實(shí)心對應(yīng)于算法1在自適應(yīng)確定正則化參數(shù)時(shí)的匹配召回率,其中自適應(yīng)確定的正則化參數(shù)分別為13.989 8,9.870 0和7.102 9。從圖中曲線可以看出,算法1通過自適應(yīng)方法總是能“自動”找到相對合適的正則化參數(shù)。這在批量處理點(diǎn)集匹配問題時(shí),可避免手動確定不合適的正則化參數(shù)。

圖2 召回率關(guān)于正則化參數(shù)的變化曲線(其中實(shí)心位置對應(yīng)于算法1自適應(yīng)確定的正則化參數(shù))

3.2 合成數(shù)據(jù)集實(shí)驗(yàn)

為評價(jià)對異常值的穩(wěn)健性,選用文獻(xiàn)[15]中的公共數(shù)據(jù)集,包含魚和漢字福兩種形狀點(diǎn)集,分別包含96和108個(gè)模型點(diǎn)。數(shù)據(jù)點(diǎn)樣本包含5個(gè)異常值比例水平(從0.0到2.0),對于每一個(gè)比例,均有100個(gè)數(shù)據(jù)點(diǎn)集樣本,采用平均的召回率評價(jià)匹配結(jié)果。圖3給出了6種算法應(yīng)用于2種形狀點(diǎn)集的平均召回率統(tǒng)計(jì)結(jié)果。

圖3 6種算法應(yīng)用于魚和漢字福2種形狀點(diǎn)集的召回率對比結(jié)果

從圖中可看到,隨著異常值比例的增加,所有的算法的匹配效果均迅速衰減,其中GMM-L2方法最嚴(yán)重,這主要是由于異常值破壞了點(diǎn)集的形狀結(jié)構(gòu)。但相比其他5種算法,本文提出的VBSMM算法總是能給出較好的匹配結(jié)果。圖4給出了算法1在2種形狀點(diǎn)集上的匹配結(jié)果示例。

圖4 本文算法在不同異常值比例下的匹配結(jié)果示例

3.3 CMU房子序列圖像實(shí)驗(yàn)

為驗(yàn)證本文算法對模型點(diǎn)集和數(shù)據(jù)點(diǎn)集中均存在異常值時(shí)的穩(wěn)健性,選用CMU房子序列(http:∥vasc.ri.cmu.edu/idb/html/motion/index.html),其包含不同視角的111幅房子圖像,每個(gè)圖像中標(biāo)記出了30個(gè)已知對應(yīng)關(guān)系的特征點(diǎn)。通過給每個(gè)點(diǎn)集增加不同數(shù)量服從均勻分布的異常值以生成待匹配點(diǎn)集,其中異常值比例分別為0.5和1.0?；诓蓸訁^(qū)間10∶10∶100對圖像序列進(jìn)行采樣,可獲得10組圖像對。因?yàn)?0個(gè)特征點(diǎn)的對應(yīng)關(guān)系已知,故可直接統(tǒng)計(jì)每個(gè)采樣區(qū)間內(nèi)的平均召回率。圖5比較了6種算法的平均召回率統(tǒng)計(jì)結(jié)果對比。從圖中可以看出,不同方法的性能出現(xiàn)了不同程度的震蕩,這主要是由于異常值生成的隨機(jī)性。整體來看,相比于其他5種算法,本文提出的VBSMM算法更穩(wěn)健,這與3.2部分的實(shí)驗(yàn)結(jié)果一致。圖6給出了本文算法在不同異常值比例下的匹配結(jié)果示例(正確匹配對數(shù)均為30/30)。

圖5 6種算法應(yīng)用于房子序列圖像的召回率對比結(jié)果

圖6 本文算法在不同異常值比例下得匹配結(jié)果示例

4 結(jié) 論

本文針對存在異常值的非剛性點(diǎn)集匹配問題提出了貝葉斯SMM模型方法。在變分貝葉斯框架中,該方法將點(diǎn)集匹配問題模型化為SMM對應(yīng)的概率密度估計(jì)問題,利用VBEM算法迭代更新變量的后驗(yàn)分布,而且允許各向異性和各向同性協(xié)方差2種選擇。通過先驗(yàn)?zāi)Ｐ鸵胝齽t化項(xiàng),可自適應(yīng)地確定正則化參數(shù)。最后,模擬點(diǎn)集和真實(shí)圖像上的實(shí)驗(yàn)結(jié)果驗(yàn)證了本文所提方法的有效性。

附錄:近似變分后驗(yàn)

其他變量的后驗(yàn)分布如下:

1) 對于混合比例系數(shù)π,其后驗(yàn)分布仍服從狄利克雷分布:q(π)～Dir(π|κ),其中:

(27)

2) 對于每個(gè)混合成分的精度矩陣Λm,其后驗(yàn)分布仍服從Wishart分布:q(Λm)～W(Λm|γm,Sm),其中:

(28)

(29)

當(dāng)選擇各向同性協(xié)方差時(shí),混合成分的精度服從伽瑪分布。

3) 對于精度參數(shù)υ,其后驗(yàn)分布仍服從伽瑪分布:q(υl)～G(υl|al,bl),其中:

(30)

(31)

4) 對于精度參數(shù)η,其后驗(yàn)分布仍服從伽瑪分布,q(ηl)～G(ηl|cl,dl),其中:

(32)

(33)