符艷軍，孫開鋒

(1 西安外事學(xué)院，西安 710077； 2 西安精密機械研究所，西安 710075)

0 引言

在基于視覺的飛行器導(dǎo)航與制導(dǎo)、遙感衛(wèi)星災(zāi)害監(jiān)控與環(huán)境監(jiān)測、醫(yī)學(xué)圖像分析等應(yīng)用中， 常常需要對不同成像傳感器獲取的異源圖像進行匹配。

基于互信息測度的匹配方法是一種完全基于圖像灰度統(tǒng)計概率的匹配方法，不需要對原圖像間的灰度關(guān)系作任何假設(shè)，非常適合于異源圖像的匹配，但互信息測度最大的缺陷是計算量大， 匹配耗時較長[1]。為此，研究者們就圍繞減少互信息測度計算量以及改進搜索策略等方面提出了相應(yīng)的加速方法。文獻(xiàn)[2-5]分別采用模擬退火、粒子群等優(yōu)化搜索算法提高匹配速度；文獻(xiàn)[6]采用灰度壓縮和粒子群優(yōu)化算法相結(jié)合加快匹配速度；文獻(xiàn)[7]通過提取多層次特征點減少了互信息測度的計算量；文獻(xiàn)[8-10]采用由粗到精的多層匹配方法對算法進行加速；文獻(xiàn)[11-15]將分層匹配與各種優(yōu)化算法相結(jié)合提高匹配速度。

文中利用匹配過程中相鄰基準(zhǔn)子圖間灰度統(tǒng)計特征之間的相關(guān)性，通過差量法減少每一個匹配位置互信息的計算量來加快匹配速度。與已有的以灰度壓縮、特征提取或采取優(yōu)化搜索策略等快速互信息匹配方法相比，該方法既沒有減少參與匹配的像素數(shù)，也無需對圖像灰度等級進行處理，不會對匹配精度造成影響。

1 圖像熵及互信息匹配

圖像熵描述了圖像信源的平均信息量，反映了圖像灰度的統(tǒng)計信息，相似的兩幅圖像其圖像熵也相近。利用互信息進行圖像匹配的實質(zhì)是：當(dāng)兩幅圖像在空間位置配準(zhǔn)時，其重疊部分所對應(yīng)像素對的灰度互信息達(dá)到最大值。 兩幅圖像U和V的互信息I(U,V)可以用下式表示：

I(U,V)=H(U)+H(V)-H(U,V)

(1)

其中:H(U)、H(V) 分別表示圖像U及V的熵；H(U,V)表示兩幅圖像的聯(lián)合熵。熵通常由變量的概率密度來表示。針對灰度圖像，可以分別用兩幅圖像的灰度直方圖和聯(lián)合直方圖進行估計。 計算公式分別為：

(2)

(3)

(4)

式(2)～式(4)中，hU(i)表示圖像U中灰度值為i的像素出現(xiàn)的次數(shù)；hV(i)表示圖像V中灰度值為i的像素出現(xiàn)的次數(shù)；hUV(i,j)表示同一像素位置上圖像U中灰度值為i、圖像V中灰度值為j的像素對出現(xiàn)的次數(shù)；T為圖像的總像素數(shù)。

研究發(fā)現(xiàn)[16]，互信息本身的大小與待配準(zhǔn)兩圖像間的重疊度有一定的關(guān)聯(lián)性，為了消除這種關(guān)聯(lián)關(guān)系，文獻(xiàn)[17]提出了利用歸一化互信息(normalized mutual information， NMI)作為相似性測度。NMI能夠減少對圖像重疊部分的敏感性，實驗表明它比標(biāo)準(zhǔn)的互信息方法更具魯棒性。NMI表達(dá)式為：

NMI(U,V)=(H(U)+H(V))/H(U,V)

(5)

在飛行器導(dǎo)航中，匹配的目的是確定飛行器的當(dāng)前位置，通過計算實測圖與基準(zhǔn)圖中每一個匹配位置對應(yīng)基準(zhǔn)子圖的互信息并找出最大互信息值對應(yīng)位置，即可確定基準(zhǔn)圖與實測圖的相對位置, 實現(xiàn)對飛行器的定位。

2 快速互信息匹配算法

與基于特征的圖像匹配算法相比，基于最大歸一化互信息的圖像匹配最大缺點是計算速度慢。 每計算一次互信息測度，需要遍歷圖像中的所有像素，這使得計算量大大增加。由式(5)可以看出，歸一化互信息的計算量主要集中在兩幅圖像各自的熵及其聯(lián)合熵的計算。從式(2)～式(4)可以看出，圖像熵及聯(lián)合熵均為求和運算，其計算量主要集中于每一個求和項的計算，因此，如果能減少每次匹配過程中需要計算的求和項的項數(shù)， 則可以加快匹配速度。因為圖像分布具有塊狀結(jié)構(gòu)，在遍歷式搜索匹配過程中，各相鄰基準(zhǔn)子圖的灰度統(tǒng)計量具有很強的相關(guān)性，基準(zhǔn)圖像(i,j)位置處對應(yīng)子圖的灰度直方圖與其四鄰域位置處對應(yīng)子圖的灰度直方圖非常相近。

圖1 相鄰基準(zhǔn)子圖及其灰度直方圖

如圖1所示，圖1(a)為原始基準(zhǔn)圖像，圖1(b)為圖1(a)的一個子圖，圖1(c)為圖1(b)在原基準(zhǔn)圖上右移一個像素對應(yīng)的子圖，圖1(d)和圖1(e)分別為圖1(b)和圖1(c)的灰度直方圖。從圖中可以看出，相鄰兩個子圖的直方圖非常相似，也即兩相鄰子圖存在大量的灰度值其出現(xiàn)的頻率沒有發(fā)生變化，相應(yīng)的，在式(2)～式(4)的求和項中分別有很多項數(shù)沒有發(fā)生變化，因此， 當(dāng)前匹配位置處基準(zhǔn)子圖的熵以及與實測圖的聯(lián)合熵的計算可以前一匹配位置處基準(zhǔn)子圖的熵及與實測圖的聯(lián)合熵的計算為基準(zhǔn)，通過比較相鄰兩個匹配位置處基準(zhǔn)子圖對應(yīng)的直方圖及聯(lián)合直方圖的變化情況，用差量法進行計算。即在整個匹配過程中，除了要完整計算基準(zhǔn)圖上第一匹配位置(1,1)處對應(yīng)基準(zhǔn)子圖的熵以及與實測圖的聯(lián)合熵外，其它所有匹配位置的基準(zhǔn)子圖的熵及與實測圖的聯(lián)合熵的計算主要集中在與前一匹配位置相比灰度直方圖矩陣和聯(lián)合直方圖矩陣中發(fā)生變化的那些元素對應(yīng)的求和項上。據(jù)此， 文中互信息匹配算法的流程圖如圖2所示。具體步驟如下:

1)獲取實測圖像A和基準(zhǔn)圖像R，將兩幅圖像的灰度值調(diào)整到同一灰度區(qū)間。

圖2 文中匹配算法流程圖

2)按照式(2)或式(3)計算A的信息熵H(A)并存儲。

3)從基準(zhǔn)圖像R的左上角(1,1)點截取與實測圖像A大小相等的第1幅基準(zhǔn)子圖S1,統(tǒng)計S1各灰度值出現(xiàn)次數(shù)并存入初始灰度直方圖矩陣hS1，按照式(2)或式(3)計算S1的信息熵H(S1)并存儲，最后將hS1和H(S1)分別賦值給hS0和H(S0)。

4)統(tǒng)計實測圖A與第1幅基準(zhǔn)子圖S1對應(yīng)像素位置的灰度對出現(xiàn)次數(shù)并存入初始聯(lián)合直方圖矩陣hAS1，按照式(4)計算初始聯(lián)合熵H(AS1)并保存，將hAS1和H(AS1)分別賦值給hAS0和H(AS0)；再按照式(5)計算實測圖與基準(zhǔn)子圖的歸一化互信息值并存儲。

5)從基準(zhǔn)圖上(1,2)點位置開始，逐行逐像素按照“Z”字形進行迭代匹配。匹配過程中， 當(dāng)前匹配位置記為(i,j)，當(dāng)j>1時，當(dāng)前匹配位置的前一匹配位置應(yīng)為(i,j-1)；當(dāng)j=1時，則當(dāng)前匹配位置的前一匹配位置應(yīng)為(i-1,j)，則當(dāng)前位置(i,j)處的歸一化互信息可按如下步驟計算:

①通過比較當(dāng)前匹配位置(i,j)處基準(zhǔn)子圖和前一匹配位置(i,j-1)或(i-1,j)處基準(zhǔn)子圖對應(yīng)列或行像素變化情況，按照差量法將前一基準(zhǔn)子圖的直方圖矩陣hS0更新為當(dāng)前子圖的直方圖矩陣hS。具體做法為: 在行方向匹配時，當(dāng)前匹配位置的前一匹配位置為(i,j-1)，則以hS0為基準(zhǔn)，減去前一基準(zhǔn)子圖第一列對應(yīng)像素灰度值出現(xiàn)次數(shù)，加上當(dāng)前基準(zhǔn)子圖最后一列對應(yīng)像素灰度值出現(xiàn)次數(shù)，即得到hS；在列方向匹配時，當(dāng)前匹配位置的前一匹配位置為(i-1,j)，則以hS0為基準(zhǔn)，減去前一基準(zhǔn)子圖第一行對應(yīng)像素灰度值出現(xiàn)次數(shù)，加上當(dāng)前基準(zhǔn)子圖最后一行對應(yīng)像素灰度值出現(xiàn)次數(shù)，即得到hS。

②找出hS與hS0的不同元素，分別組成矩陣hΔS和hΔS0，按照式(6)、式(7)分別計算hΔS和hΔS0對應(yīng)的信息熵ΔS及ΔS0；其中，normS0(S)為與hΔS0中元素對應(yīng)的前一基準(zhǔn)子圖熵值計算過程中的求和項，因為這些求和項在前一匹配位置已經(jīng)計算過，此次無需再重復(fù)計算，而只需計算這些項的和即可。

(6)

(7)

③以前一基準(zhǔn)子圖的信息熵H(S0)為基準(zhǔn)，按照式(8)通過差量法將H(S0)更新為當(dāng)前基準(zhǔn)子圖的信息熵H(S)。

H(S)=H(S0)-ΔS0+ΔS

(8)

④計算實測圖與當(dāng)前基準(zhǔn)子圖的聯(lián)合直方圖矩陣hAS，找出hAS和前一基準(zhǔn)子圖聯(lián)合直方圖矩陣hAS0的不同元素，分別組成矩陣hΔAS和hΔAS0，按照式(8)、式(9)分別計算hΔAS和hΔAS0對應(yīng)的聯(lián)合熵ΔAS及ΔAS0；其中，normAS0為與hΔAS0中元素對應(yīng)的前一基準(zhǔn)子圖聯(lián)合熵計算過程中的求和項，因為這些求和項在前一匹配位置已經(jīng)計算過，此次無需再重復(fù)計算，而只需計算這些項的和即可。

(9)

(10)

⑤以前一基準(zhǔn)子圖的聯(lián)合熵H(AS0)為基準(zhǔn)，按照式(11)通過差量法將H(AS0)更新為當(dāng)前基準(zhǔn)子圖的聯(lián)合熵H(AS)。

H(AS)=H(AS0)-ΔAS0+ΔAS

(11)

⑥按照式(12)計算實測圖與當(dāng)前基準(zhǔn)子圖的歸一化互信息。

NMI(A,S)=(H(A)+H(S))/H(A,S)

(12)

6) 遍歷整個基準(zhǔn)圖，選取歸一化互信息最大值對應(yīng)的位置作為最終匹配點。

3 算法時間復(fù)雜性分析

文中算法與傳統(tǒng)方法相比，計算量的差別主要體現(xiàn)在每一個匹配位置處基準(zhǔn)子圖各統(tǒng)計量的計算上，為此只分析基準(zhǔn)子圖各統(tǒng)計量計算的時間復(fù)雜度。

設(shè)基準(zhǔn)圖尺寸為N×N，實測圖尺寸為n×n，則基準(zhǔn)子圖大小與實測圖大小一致均為n×n。

1)基準(zhǔn)子圖直方圖計算

每一個匹配位置處，傳統(tǒng)算法需要遍歷整個基準(zhǔn)子圖來計算該子圖的直方圖，而文中方法只須遍歷基準(zhǔn)子圖的第一行(或列)和最后一行(或列)即可獲得該子圖直方圖。對于整個匹配過程來說，傳統(tǒng)方法其時間復(fù)雜度為O((N-n)2n2)，文中方法時間復(fù)雜度為O((N-n)2n)。

2)基準(zhǔn)子圖熵及聯(lián)合熵計算

由于相鄰基準(zhǔn)子圖其直方圖發(fā)生變化情況跟圖像內(nèi)容有關(guān)，很難從理論上推導(dǎo)兩子圖間有多少個灰度值其出現(xiàn)頻率會發(fā)生變化，也很難推導(dǎo)出其聯(lián)合直方圖中有多少個灰度值對的頻率會發(fā)生變化。以直方圖為例，最壞情況時有2n(n為基準(zhǔn)子圖各邊長所含像素個數(shù))個灰度值頻率發(fā)生變化，且2n不超過基準(zhǔn)圖灰度等級數(shù)。

為此，文中在標(biāo)準(zhǔn)圖像集上進行了大量的仿真實驗，通過對仿真結(jié)果進行統(tǒng)計得出：相鄰子圖間灰度直方圖中發(fā)生變化的元素平均為50%左右，聯(lián)合灰度直方圖中發(fā)生變化的元素平均為75%左右，這就意味著文中方法相比傳統(tǒng)方法，其基準(zhǔn)子圖熵值的計算量減少了近50%，聯(lián)合熵計算量減少了約25%，時間復(fù)雜度或時間頻度有所降低。部分實驗結(jié)果列于表1～表4。

另一方面，因文中方法需要存儲前一匹配位置基準(zhǔn)子圖的多個統(tǒng)計量，算法的空間復(fù)雜度有所增加，但由于硬件技術(shù)的發(fā)展大大提高了計算機的存儲容量，使得存儲容量的局限性對于算法的影響大大降低。

4 仿真實驗

為了驗證算法的有效性，將文中方法應(yīng)用于SAR與可見光圖像匹配，以Matlab R2015為仿真平臺，進行兩類實驗：一類是通過比較每個匹配位置處文中方法與傳統(tǒng)遍歷式互信息匹配方法在計算量方面的差別以及總耗時來說明文中方法的有效性；另一類是將文中方法與多分辨率分層匹配相結(jié)合，與單純的多分辨率方法進行匹配耗時及匹配精度的比較，其中的多分辨率分層匹配采用d3小波基，小波分解級數(shù)為1。限于篇幅，文中只列出部分結(jié)果。

匹配過程中設(shè): 基準(zhǔn)子圖直方圖計算中需要統(tǒng)計的像素灰度值的次數(shù)為M，基準(zhǔn)子圖熵值計算中需要計算的求和項的項數(shù)為M1，基準(zhǔn)子圖與實測圖聯(lián)合熵計算中需要計算的求和項項數(shù)為M2。

選取可見光與SAR圖像進行匹配，如圖3所示，其中圖3(a)為可見光基準(zhǔn)圖，大小為422×358；圖3(b)為SAR實測圖像1，大小為150×150，圖3(c)為SAR實測圖像2，大小為100×90，圖3(d)為基于文中差量法匹配時實測圖在基準(zhǔn)圖上的匹配定位結(jié)果；表1為差量法與傳統(tǒng)標(biāo)準(zhǔn)算法列方向匹配時部分匹配位置計算量統(tǒng)計結(jié)果的比較，表2為差量法與傳統(tǒng)標(biāo)準(zhǔn)算法行方向匹配時部分匹配位置計算量統(tǒng)計結(jié)果的比較，表3為文中差量法與傳統(tǒng)標(biāo)準(zhǔn)方法匹配結(jié)果和匹配耗時的比較，表4為差量法結(jié)合多分辨率分層匹配法與單純多分辨率分層匹配法的匹配結(jié)果和匹配耗時的比較。

圖3 SAR與可見光圖像匹配

由圖3可以看出，采用歸一化互信息作為匹配測度，可以實現(xiàn)異源圖像之間的匹配；由表1～表3可以看出，相比傳統(tǒng)方法，文中差量法能夠有效減少除第一個匹配位置外其它所有匹配位置互信息測度的計算量，總匹配耗時平均減少20%以上；由表4可以看出，文中差量法可以與多分辨率分層匹配技術(shù)相結(jié)合，其匹配精度與單純的多分辨率分層匹配算法精度相同(在有些情況下匹配精度有所下降)，但其匹配耗時平均下降20%左右。

5 結(jié)語

與已有的以灰度壓縮、特征提取或采取優(yōu)化搜索策略等快速互信息匹配方法不同,文中基于圖像的塊狀結(jié)構(gòu)特性以及匹配過程中相鄰子圖統(tǒng)計特性的相關(guān)性, 通過差量法減少互信息測度的計算量來加快匹配速度, 其互信息測度計算精度及匹配精度與原始的基于標(biāo)準(zhǔn)互信息計算方法的計算精度及匹配精度相同, 但計算量大為減少。 另外,如果參與匹配的兩幅圖像為高分辨率清晰圖像，則可以把文中方法與灰度壓縮、多分辨率分層匹配等方法相結(jié)合，即先對待匹配的圖像進行灰度壓縮或多分辨率分解，然后采用遍歷式方法進行搜索，而在每一個搜索點處采用文中的差量法計算互信息測度，這樣可以進一步減少匹配耗時。

表1 SAR與可見光圖像在列方向匹配時部分位置計算量統(tǒng)計結(jié)果

表2 SAR與可見光圖像在行方向匹配時部分位置計算量統(tǒng)計結(jié)果

表3 差量法與傳統(tǒng)標(biāo)準(zhǔn)法在SAR與可見光圖像上的匹配結(jié)果

表4 差量法結(jié)合多分辨率分層算法與單純多分辨率分層算法在SAR與可見光圖像上的匹配結(jié)果