楊麗娟, 田錚,2, 溫金環(huán), 延偉東
(1.西北工業(yè)大學(xué) 應(yīng)用數(shù)學(xué)系, 陜西 西安 710129; 2.中國(guó)科學(xué)院 遙感科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100101)
點(diǎn)集匹配在圖像分析、形狀匹配等許多領(lǐng)域中已獲得廣泛的應(yīng)用[1-2],可粗略的分為剛性和非剛性,其中非剛性匹配問(wèn)題備受關(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)集,通過(guò)最小化2個(gè)混合高斯之間的統(tǒng)計(jì)差異,進(jìn)而確定點(diǎn)集之間的對(duì)應(yīng)關(guān)系。Myronenko等[5]假設(shè)點(diǎn)集做一致性運(yùn)動(dòng),提出了基于一致性點(diǎn)漂移(coherent point drift,CPD)的點(diǎn)集匹配方法,通過(guò)增加額外的均勻分布模型化異常值。但它本質(zhì)上對(duì)應(yīng)于最小二乘方法,對(duì)異常值不穩(wěn)健。為了提高對(duì)異常值的穩(wěn)健性,一種可選的方法是選用更合適的概率模型,如具有厚尾性質(zhì)的t分布[3]。Zhou等[6]提出用混合t分布模型(student′s t mixture model,SMM)處理點(diǎn)集之間的匹配問(wèn)題,通過(guò)最大化似然函數(shù)的對(duì)數(shù)期望,估計(jì)變換參數(shù)。但是它和CPD方法以及多數(shù)方法一樣,只允許各向同性協(xié)方差,而且需要手動(dòng)確定正則化參數(shù)。
針對(duì)最大似然估計(jì)框架中存在的奇異性問(wèn)題,可選的解決方式是貝葉斯分析[7-8]。Qu等[9]提出了在變分貝葉斯框架下的點(diǎn)集匹配方法(VBPSM),將點(diǎn)集匹配過(guò)程分為聚類分析和回歸分析兩部分,并分別利用2個(gè)GMM模型化,通過(guò)引入轉(zhuǎn)移變量將兩部分統(tǒng)一起來(lái)。但對(duì)于非剛性變換,它沒(méi)有包含對(duì)變換函數(shù)的平滑性約束。另一個(gè)問(wèn)題是估計(jì)的變換函數(shù)沒(méi)有直接反映待對(duì)準(zhǔn)點(diǎn)集之間的變換關(guān)系。在Simpson等[10]提出的非剛性配準(zhǔn)的概率配準(zhǔn)框架中,通過(guò)給變換參數(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ò)程中,考慮了變量之間的相關(guān)。Baldacchino等[13]將貝葉斯SMM應(yīng)用于非線性系統(tǒng)識(shí)別中的穩(wěn)健回歸。
綜上,針對(duì)存在異常值的非剛性點(diǎn)集匹配問(wèn)題,本文將點(diǎn)集匹配問(wèn)題表示為SMM對(duì)應(yīng)的概率密度估計(jì)問(wèn)題。通過(guò)引入真實(shí)后驗(yàn)分布的變分近似,匹配問(wèn)題轉(zhuǎn)化為最大化變分下界,同時(shí)允許任何形式的協(xié)方差。在參數(shù)推斷過(guò)程中,利用變分貝葉斯最大化期望(variational Bayesian EM,VBEM)算法確定變換參數(shù)。通過(guò)先驗(yàn)?zāi)P?,將空間正則化約束并入貝葉斯SMM模型中,可自適應(yīng)地確定正則化參數(shù)。模擬點(diǎn)集和真實(shí)圖像上的實(shí)驗(yàn)結(jié)果證明所提方法可大大提高匹配性能。
(1)
式中,Θ={A,B,{Λm},{πm},{vm}},t分布的概率密度函數(shù)為:
(2)
式中,μ和Λ分別為均值向量和精度矩陣,Γ(·)為伽瑪函數(shù),當(dāng)自由度參數(shù)v→∞時(shí),t分布收斂到高斯分布。t分布可表示為無(wú)限個(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ù)對(duì);⑤正則項(xiàng)以先驗(yàn)形式并入貝葉斯SMM模型中[10],即:
(8)
式中,λ>0是標(biāo)量正則化參數(shù);⑥正則化參數(shù)λ服從形狀參數(shù)為τ0和尺度參數(shù)為ζ0的伽瑪分布??筛鶕?jù)經(jīng)驗(yàn)確定超參數(shù)的值。通過(guò)變分推斷,可獲得變量{ΘL,ΘVB}={ΘL,{A,B,Λ,π,υ,η,λ}}的后驗(yàn)分布,由于沒(méi)有自由度參數(shù)v={vm}的先驗(yàn),利用最大似然方法可確定它的值,記為ΘML={v},詳見下一節(jié)。
變量之間的聯(lián)合分布可表示為:
(9)
圖1中給出了非剛性點(diǎn)集匹配對(duì)應(yīng)的概率圖模型,即有向無(wú)環(huán)圖[7]。
圖1 非剛性點(diǎn)集匹配的概率圖模型
考慮近似后驗(yàn)分布q(U,Z,ΘVB),對(duì)數(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è)對(duì)角元素。
(21)
(22)
4) 正則化參數(shù)λ的后驗(yàn)分布服從伽瑪分布:q(λ)~G(λ|τ,ζ),其中
(23)
(24)
根據(jù)伽瑪分布性質(zhì),有〈λ〉=τ/ζ,可根據(jù)不同的應(yīng)用自適應(yīng)地確定λ的值。
其余變量的后驗(yàn)分布見附錄。由于沒(méi)有自由度參數(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ù)。
為了驗(yàn)證本文算法的有效性,首先測(cè)試具有形狀差異的模擬數(shù)據(jù)集,然后將本文方法應(yīng)用于包含非剛性形變的圖像。并與5種算法進(jìn)行了對(duì)比: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。針對(duì)各向異性和各向同性協(xié)方差2種選擇,本文方法分別記為VBSMM(ANI)和VBSMM(ISO)。
為了檢驗(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í)獲得的匹配召回率做對(duì)比。圖2給出了召回率關(guān)于正則化參數(shù)的變化曲線,其中空心對(duì)應(yīng)于算法1在固定正則化參數(shù)時(shí)的匹配召回率,實(shí)心對(duì)應(yīng)于算法1在自適應(yīng)確定正則化參數(shù)時(shí)的匹配召回率,其中自適應(yīng)確定的正則化參數(shù)分別為13.989 8,9.870 0和7.102 9。從圖中曲線可以看出,算法1通過(guò)自適應(yīng)方法總是能“自動(dòng)”找到相對(duì)合適的正則化參數(shù)。這在批量處理點(diǎn)集匹配問(wèn)題時(shí),可避免手動(dòng)確定不合適的正則化參數(shù)。
圖2 召回率關(guān)于正則化參數(shù)的變化曲線(其中實(shí)心位置對(duì)應(yīng)于算法1自適應(yīng)確定的正則化參數(shù))
為評(píng)價(jià)對(duì)異常值的穩(wěn)健性,選用文獻(xiàn)[15]中的公共數(shù)據(jù)集,包含魚和漢字福兩種形狀點(diǎn)集,分別包含96和108個(gè)模型點(diǎn)。數(shù)據(jù)點(diǎn)樣本包含5個(gè)異常值比例水平(從0.0到2.0),對(duì)于每一個(gè)比例,均有100個(gè)數(shù)據(jù)點(diǎn)集樣本,采用平均的召回率評(píng)價(jià)匹配結(jié)果。圖3給出了6種算法應(yīng)用于2種形狀點(diǎn)集的平均召回率統(tǒng)計(jì)結(jié)果。
圖3 6種算法應(yīng)用于魚和漢字福2種形狀點(diǎn)集的召回率對(duì)比結(jié)果
從圖中可看到,隨著異常值比例的增加,所有的算法的匹配效果均迅速衰減,其中GMM-L2方法最嚴(yán)重,這主要是由于異常值破壞了點(diǎn)集的形狀結(jié)構(gòu)。但相比其他5種算法,本文提出的VBSMM算法總是能給出較好的匹配結(jié)果。圖4給出了算法1在2種形狀點(diǎn)集上的匹配結(jié)果示例。
圖4 本文算法在不同異常值比例下的匹配結(jié)果示例
為驗(yàn)證本文算法對(duì)模型點(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è)已知對(duì)應(yīng)關(guān)系的特征點(diǎn)。通過(guò)給每個(gè)點(diǎn)集增加不同數(shù)量服從均勻分布的異常值以生成待匹配點(diǎn)集,其中異常值比例分別為0.5和1.0。基于采樣區(qū)間10∶10∶100對(duì)圖像序列進(jìn)行采樣,可獲得10組圖像對(duì)。因?yàn)?0個(gè)特征點(diǎn)的對(duì)應(yīng)關(guān)系已知,故可直接統(tǒng)計(jì)每個(gè)采樣區(qū)間內(nèi)的平均召回率。圖5比較了6種算法的平均召回率統(tǒng)計(jì)結(jié)果對(duì)比。從圖中可以看出,不同方法的性能出現(xiàn)了不同程度的震蕩,這主要是由于異常值生成的隨機(jī)性。整體來(lái)看,相比于其他5種算法,本文提出的VBSMM算法更穩(wěn)健,這與3.2部分的實(shí)驗(yàn)結(jié)果一致。圖6給出了本文算法在不同異常值比例下的匹配結(jié)果示例(正確匹配對(duì)數(shù)均為30/30)。
圖5 6種算法應(yīng)用于房子序列圖像的召回率對(duì)比結(jié)果
圖6 本文算法在不同異常值比例下得匹配結(jié)果示例
本文針對(duì)存在異常值的非剛性點(diǎn)集匹配問(wèn)題提出了貝葉斯SMM模型方法。在變分貝葉斯框架中,該方法將點(diǎn)集匹配問(wèn)題模型化為SMM對(duì)應(yīng)的概率密度估計(jì)問(wèn)題,利用VBEM算法迭代更新變量的后驗(yàn)分布,而且允許各向異性和各向同性協(xié)方差2種選擇。通過(guò)先驗(yàn)?zāi)P鸵胝齽t化項(xiàng),可自適應(yīng)地確定正則化參數(shù)。最后,模擬點(diǎn)集和真實(shí)圖像上的實(shí)驗(yàn)結(jié)果驗(yàn)證了本文所提方法的有效性。
其他變量的后驗(yàn)分布如下:
1) 對(duì)于混合比例系數(shù)π,其后驗(yàn)分布仍服從狄利克雷分布:q(π)~Dir(π|κ),其中:
(27)
2) 對(duì)于每個(gè)混合成分的精度矩陣Λm,其后驗(yàn)分布仍服從Wishart分布:q(Λm)~W(Λm|γm,Sm),其中:
(28)
(29)
當(dāng)選擇各向同性協(xié)方差時(shí),混合成分的精度服從伽瑪分布。
3) 對(duì)于精度參數(shù)υ,其后驗(yàn)分布仍服從伽瑪分布:q(υl)~G(υl|al,bl),其中:
(30)
(31)
4) 對(duì)于精度參數(shù)η,其后驗(yàn)分布仍服從伽瑪分布,q(ηl)~G(ηl|cl,dl),其中:
(32)
(33)