凌 程，耿修瑞，楊煒暾，趙永超

(中國科學(xué)院電子學(xué)研究所 中國科學(xué)院空間信息處理與應(yīng)用系統(tǒng)技術(shù)重點(diǎn)實(shí)驗(yàn)室， 北京 100190； 中國科學(xué)院大學(xué)， 北京 100049)(2017年11月13日收稿； 2018年1月9日收修改稿)

圖像配準(zhǔn)是一項(xiàng)非常重要的圖像預(yù)處理技術(shù)。隨著數(shù)字圖像處理的發(fā)展，圖像配準(zhǔn)技術(shù)已經(jīng)被廣泛應(yīng)用于高分辨率圖像生成[1]、醫(yī)學(xué)圖像對齊[2]、高精度3D重建[3]、遙感圖像變化檢測、多源圖像融合等領(lǐng)域。這些關(guān)鍵領(lǐng)域中，需要高精度的圖像配準(zhǔn)，即亞像素級的圖像配準(zhǔn)?，F(xiàn)有的亞像素級圖像配準(zhǔn)算法主要分為3類[4-5]：基于插值的方法[6-8]、解最優(yōu)化問題法[9-10]和擴(kuò)展的相位相關(guān)法[2,11-13]。其中，基于插值的方法采用插值技術(shù)對相似性函數(shù)或者圖像進(jìn)行插值重采樣，以插值后的最大值位置作為配準(zhǔn)位置，以此達(dá)到亞像素的精度。Guizar-Sicairos等[7]對相似形函數(shù)內(nèi)插法進(jìn)行快速改進(jìn)；Yousef等[8]在此基礎(chǔ)上進(jìn)一步提高速度，在保證配準(zhǔn)精度的同時大大減少計算復(fù)雜度。本文在此將Yousef等[8]提出的快速相似性函數(shù)內(nèi)插法稱為fast-SSDFT方法。解最優(yōu)化問題法的目標(biāo)函數(shù)為衡量兩幅待配準(zhǔn)圖像之間相似程度的函數(shù)，支持各種變換模型，具有參數(shù)選擇靈活，計算精度高的特點(diǎn)，但是其計算量非常大，收斂概率和得到全局最優(yōu)解的概率也需要提高?；谙辔幌嚓P(guān)的方法利用傅里葉變換的特性，將時域上的平移量轉(zhuǎn)化為求解頻域上的相位相關(guān)，具有精度高、對噪聲魯棒和計算復(fù)雜度低等優(yōu)點(diǎn)。

目前所有擴(kuò)展的相位相關(guān)法都是基于傅里葉變換的平移特性，即兩幅待配準(zhǔn)的圖像在空間域上存在平移，轉(zhuǎn)換到頻域上則表現(xiàn)為線性相位，通過對頻域上線性相位的求解來獲得空間域上的平移量。傳統(tǒng)的相位相關(guān)法直接將歸一化互功率譜矩陣進(jìn)行傅里葉反變換，通過判斷沖激函數(shù)的峰值位置得到平移參數(shù)的值，但是這種方法只能得到整數(shù)像素值。Hoge[11]提出，通過對歸一化互功率譜矩陣進(jìn)行奇異值分解，用最大奇異值所對應(yīng)的奇異向量估算出平移量，可以達(dá)到亞像素的精度。該方法通過選取最大奇異值，拋棄小的奇異值，起到對噪聲過濾的作用，因此該方法對平移參數(shù)的估計具有較高的精度。另外為了求得亞像素精度的平移量，需要對得到的奇異向量的相位值進(jìn)行一維相位解纏，然后用最小二乘法估計出平移量。Tong等[12]利用Ransac算法[14]替代最小二乘法對解纏后的相位值進(jìn)行估計，在噪聲較大的情況下，Ransac算法能夠?yàn)V除噪聲帶來的異常值，使得估計出的平移參數(shù)精度更高。Feng等[13]針對低信噪比的圖像提出一種魯棒的亞像素配準(zhǔn)算法，該算法通過計算兩幅圖像的互相關(guān)系數(shù)[15]的質(zhì)心得到精確的平移參數(shù)，從而達(dá)到亞像素的配準(zhǔn)精度，本文在此將其稱為CC_centriod方法。

在理想無噪的情況下，纏繞相位具有一定的周期性，可以通過積分的方法進(jìn)行相位解纏。通過積分的方法求得真實(shí)相位具有一個前提，即相鄰兩個像元間的真實(shí)相位差的絕對值要小于π[16]，但是如果兩幅圖像之間的平移量超過某個閾值或者受噪聲影響時，所得到的一維纏繞相位向量相鄰的真實(shí)相位差的絕對值會超出π，因此傳統(tǒng)的積分法進(jìn)行相位解纏就會變得不再適用。本文提出一種改進(jìn)的相位解纏方法對傳統(tǒng)積分法得到的結(jié)果進(jìn)行校正，改進(jìn)的算法可以避免積分法受到此前提和噪聲的影響，使得估計出的平移參數(shù)更精確。本文把基于奇異值分解(SVD)的相位相關(guān)法結(jié)合改進(jìn)的相位解纏方法稱為SVD-IIM(SVD-improved integral method)。

1 基于SVD分解的相位相關(guān)法

1.1 方法原理和流程

令x0和y0表示兩幅大小為M×N的圖像之間在行方向和列方向上的平移量，可以表示成

g(x,y)=f(x-x0,y-y0).

(1)

式中：f(x,y)和g(x,y)為二維離散函數(shù)，其中變量x=1,2,3,…,M和變量y=1,2,3,…,N。根據(jù)傅里葉變換的性質(zhì)可以得到

G(u,v)=F(u,v)exp{-j2π(ux0/M+vy0/N)}.

(2)

式中：F(u,v)和G(u,v)分別是f(x,y)和g(x,y)的二維離散傅里葉變換函數(shù)，其中變量u=1,2,3,…,M和變量v=1,2,3,…,N。因此可以得到歸一化互功率譜函數(shù)Q(u,v)

= exp{-j2π(ux0/M+vy0/N)}.

(3)

式中:*代表共軛，傳統(tǒng)的相位相關(guān)法直接計算歸一化互功率譜函數(shù)Q(u,v)的傅里葉反變換：

在這個句子中，漢語外延文化信息“滄?！焙汀吧Ｌ铩狈g成“sea and mulberry fields”?！皽婧Ｗ兩Ｌ铩币馕吨薮蟮母淖儯⒄Z讀者很容易理解文化負(fù)載詞的文化含義，所以譯者盡可能地逐字翻譯。

Q(x,y)=δ(x-x0,y-y0).

(4)

在理想情況下，歸一化互功率譜矩陣Q在理論上是一個秩為1的矩陣，因此歸一化互功率譜函數(shù)Q(u,v)可以進(jìn)行如下分解：

Q(u,v)=exp{-j2π(ux0/M+vy0/N)}

= exp{-j2πux0/M}exp{-j2πvy0/N}

(5)

式中：H代表共軛轉(zhuǎn)置，令向量qx0和qy0為：

qx0=(e-j(2π/M)x0,e-j(2π/M)2x0,…,e-j(2π/M)Mx0)T，

(6)

qy0=(ej(2π/N)y0,ej(2π/N)2y0,…,ej(2π/N)Ny0)T.

(7)

因此歸一化互功率譜矩陣Q可以表示為

(8)

但是，在一般情況下，Q并不是一個秩為1的矩陣，因此文獻(xiàn)[11]提出可以通過對Q矩陣進(jìn)行SVD分解，用最大奇異值和相對應(yīng)的奇異向量對Q進(jìn)行秩1估計，由于奇異值最大的分量受噪聲影響較小，因而秩1估計可以起到抑制噪聲的作用。對Q進(jìn)行秩1估計：

(9)

式中：σ1為最大的奇異值，u1和v1為σ1所對應(yīng)的奇異向量，因此有u1≈qx0和v1≈qy0。對u1和v1進(jìn)行相位解纏，得到解纏后的相位px和py，并用最小二乘法擬合px和py，得到px和py的斜率kx和ky，通過得到的斜率最終求得平移參數(shù)：

(10)

圖1 基于SVD分解的相位相關(guān)法流程圖Fig.1 Flow chart of phase correlation method based on SVD

1.2 積分法相位解纏原理

在理想情況下，圖像的采樣率滿足Nyquist采樣定理，采樣頻率必須大于信號最高頻率的2倍，解纏相位中相鄰像素點(diǎn)之間的相位差值不可能超過半個周期，當(dāng)滿足此條件時必然能通過積分法，由纏繞相位解纏出正確的真實(shí)相位值。令φ(m)為周期纏繞前的真實(shí)相位值，φ(m)為相應(yīng)的纏繞相位，相位纏繞算子可以如下定義：

ω(φ(m))=φ(m)=φ(m)+2πk(m).

(11)

其中-π<ω(φ(m))<π，得到屬于(-π,π)的纏繞相位φ(m)，定義差分算子Δ，因此對于真實(shí)相位值有

Δφ(m)=φ(m+1)-φ(m).

(12)

由于真實(shí)相位值中相鄰像素點(diǎn)之間的相位差值不可能超過半個周期，所以有

-π<Δφ(m)<π.

(13)

對相鄰纏繞相位也進(jìn)行差分運(yùn)算得

Δφ(m)=φ(m+1)-φ(m)=Δφ(m)+

2πΔk(m).

(14)

對該相位差使用纏繞算子得

ω(Δφ(m))=Δφ(m)+2πΔk(m)+2πk′(m).

(15)

根據(jù)纏繞算子ω的定義，ω(Δφ(m))必然在區(qū)間(-π,π)內(nèi)，由式(13)可知，Δφ(m)也屬于區(qū)間(-π,π)，因此有

Δk(m)+k′(m)=0.

(16)

所以可將式(15)簡化為

ω(Δφ(m))=Δφ(m).

(17)

因此最終的真實(shí)相位可以表示成

ω(Δφ(n)).

(18)

從以上推導(dǎo)可以看出，使用積分法進(jìn)行相位解纏依賴于一個重要的前提，即真實(shí)相位中相鄰像素點(diǎn)之間的相位差值必須在(-π,π)范圍內(nèi)。而在實(shí)際遙感圖像配準(zhǔn)應(yīng)用中，由于有地形突變、相位混疊、平移量較大等各種因素的存在，相位信息中將會引入噪聲，在計算解纏相位的過程中，并不確定式(13)是否成立，所以通過積分法即式(18)計算所得到的相位值只是相位的估計值，其準(zhǔn)確性并不能保證[17]。

1.3 改進(jìn)的相位解纏方法

如圖2所示，如果當(dāng)真實(shí)相位差的絕對值大于π時，使用積分法進(jìn)行相位解纏會出現(xiàn)圖2(a)所示的現(xiàn)象，如果對φ(i+1)進(jìn)行加上或減去2π的校正操作后更符合趨勢斜率kx，那么執(zhí)行此校正操作，使得校正后的φ(i+1)-φ(i)更符合趨勢斜率，如圖2(b)所示。

圖2 相位校正示意圖Fig.2 Schematic of phase correlation

2 實(shí)驗(yàn)結(jié)果及分析

本文選取的實(shí)驗(yàn)數(shù)據(jù)為高分2號(GF-2)衛(wèi)星和資源3號(ZY-3)衛(wèi)星的高分辨率多光譜圖像數(shù)據(jù)(部分圖像數(shù)據(jù)如圖3 所示)，分別進(jìn)行校正效果實(shí)驗(yàn)和配準(zhǔn)精度對比實(shí)驗(yàn)。其中校正效果實(shí)驗(yàn)主要對比傳統(tǒng)積分法和改進(jìn)的相位解纏方法對纏繞相位的解纏效果以及相應(yīng)的配準(zhǔn)精度，配準(zhǔn)精度對比則實(shí)驗(yàn)對比在不同的平移量、加入不同大小的高斯白噪聲、不同混疊情況下，本文所提SVD-IIM方法與已有的SVD方法[11]、SVD-RANSAC方法[12]，CC-centriod[13]方法和fast-SSDFT[8]方法這4種方法的配準(zhǔn)精度。

2.1 校正效果實(shí)驗(yàn)

對上文所述的高分2號數(shù)據(jù)在光譜維上取均值得到其灰度圖像，在其中選取一幅大小為1 280像素×1 280像素的原始參考圖像，如圖4(a)所示。將其在原灰度圖像上分別在行方向平移541個像素和在列方向平移548個像素，得到另一幅平移后的原始待配準(zhǔn)圖像，如圖4(b)所示。那么這兩幅灰度圖像的平移關(guān)系為(541,548)，分別對兩幅圖像進(jìn)行10倍降采樣得到兩幅大小為128像素×128像素的圖像，如圖4(c)和4(d)所示，因此其平移關(guān)系變?yōu)?54.1,54.8)[18]。

按照圖1所示的流程計算圖4(c)和圖4(d)的歸一化互功率譜矩陣Q(u,v)，并濾除其中的高頻分量和小于特定閾值的頻率分量，再對Q矩陣進(jìn)行SVD分解得到最大奇異值對應(yīng)的奇異向量u和v，分別利用傳統(tǒng)積分法和改進(jìn)的相位解纏方法對u和v的相位進(jìn)行相位解纏。利用最小二乘法分別對圖5(b)、5(d)所示的解纏相位值的斜率進(jìn)行估計，其配準(zhǔn)結(jié)果如表1所示，顯然傳統(tǒng)積分法的誤差遠(yuǎn)遠(yuǎn)大于改進(jìn)的相位解纏方法的誤差。由圖5(b)、5(d)可以看出，改進(jìn)的相位解纏算法能對傳統(tǒng)積分法得到的結(jié)果進(jìn)行很好的校正，因此當(dāng)參考圖像和待配準(zhǔn)圖像之間平移量較大時，即當(dāng)相鄰像素點(diǎn)之間相位差的絕對值接近或超過p時，利用傳統(tǒng)的積分法相位解纏得到的結(jié)果并不可靠，而經(jīng)過改進(jìn)的相位解纏方法可以精確地估計出平移參數(shù)。

圖3 GF-2數(shù)據(jù)(a，b)和ZY-3數(shù)據(jù)(c，d)Fig.3 Imagery of GF-2(a，b) and ZY-3(c，d)

圖4 原始配準(zhǔn)圖和10倍降采樣后的配準(zhǔn)圖Fig.4 Original images and images after 10 times downsampling used for registration

表1 傳統(tǒng)積分法和改進(jìn)相位解纏方法配準(zhǔn)精度比較Table 1 Precision comparison between the traditional integral method and the improved phase unwrapping method

2.2 配準(zhǔn)精度對比實(shí)驗(yàn)

對上文所述的高分2號數(shù)據(jù)和資源3號數(shù)據(jù)在光譜維上取均值，利用得到的灰度圖像進(jìn)行配準(zhǔn)精度對比試驗(yàn)，選取大小為1 280像素×1 280像素的子圖為原始參考圖像，以其為基準(zhǔn)進(jìn)行不同程度的平移，得到原始待配準(zhǔn)圖像。分別對原始參考圖像和原始待配準(zhǔn)圖像進(jìn)行10倍降采樣得到實(shí)驗(yàn)所用的大小為128像素×128像素的配準(zhǔn)圖像。圖6和圖7為本文提出的SVD-IIM方法與其他4種方法分別在高分2號數(shù)據(jù)和資源3號數(shù)據(jù)上進(jìn)行的對比結(jié)果。其中圖6(a)、6(d)和圖7(a)、7(d)為對于不同平移量的絕對誤差對比示意圖，其中圖6(a)和圖7(a)平移量為57～62個像素，圖6(d)和圖7(d)平移量為55～60個像素，步長均為0.1個像素；圖6(b)、6(e)和圖7(b)、7(e)為在加入不同高斯白噪聲情況下的絕對誤差對比示意圖，將配準(zhǔn)圖像歸一化至0～256，加入標(biāo)準(zhǔn)差為6～10的高斯白噪聲，真實(shí)平移量為(59.4,7.8)；對待配準(zhǔn)圖減采樣之前的原圖進(jìn)行不同程度的高斯模糊可以模擬不同情況下的混疊，而高斯模糊等價于低通濾波器[2]，圖6(c)、6(f)和圖7(c)、7(f)為在不同混疊程度下的絕對誤差對比示意圖，高斯模糊窗大小為25×25，模糊因子為高斯函數(shù)的標(biāo)準(zhǔn)差, 真實(shí)平移量為(59.4,7.8)。圖中所有絕對誤差為100幅隨機(jī)選取的待配準(zhǔn)圖像配準(zhǔn)誤差的平均值。

圖5 相位解纏結(jié)果示意圖Fig.5 Schematic of the phase unwrapping results

圖6 GF-2數(shù)據(jù)配準(zhǔn)精度對比圖Fig.6 The contrast diagram of GF-2 data registration accuracy

圖7 ZY-3數(shù)據(jù)配準(zhǔn)精度對比圖Fig.7 The contrast diagram of ZY-3 data registration accuracy

如圖6(a)～6(c)和圖7(a)～7(c)所示，對于不同的平移量、加入不同大小的高斯白噪聲、不同混疊情況下，本文提出的SVD-IIM方法的誤差都在0附近，而原有的SVD方法和SVD-RANSAC方法的誤差都在10個像素以上，因此可以進(jìn)一步確定在平移量較大時傳統(tǒng)積分法得到的相位解纏結(jié)果是錯誤的，而改進(jìn)的相位解纏方法所得到的結(jié)果是可靠的，并且精度較高。

圖6(d)～6(f)和圖7(d)～7(f)對比在不同平移量，加入不同高斯白噪聲，不同混疊情況下的3種方法的配準(zhǔn)精度，可以看出本文提出的SVD-IIM方法的誤差相比于CC-centriod方法和fast-SSDFT方法的誤差更小，并且更穩(wěn)定。對于較大平移量的配準(zhǔn)問題，參考圖像和待配準(zhǔn)圖像間存在較大的非匹配區(qū)域，這些非匹配區(qū)域的信息在估計平移參數(shù)時如同噪聲，因此在存在大量噪聲的情況下，本文提出的SVD-IIM方法可以將配準(zhǔn)精度控制在0.1個像素范圍內(nèi)。

3 結(jié)束語

本文針對基于奇異值分解(SVD)的相位相關(guān)法中傳統(tǒng)積分法進(jìn)行相位解纏的可靠性問題，利用線性相位單調(diào)變化的特性，提出一種改進(jìn)的相位解纏方法，對傳統(tǒng)積分法得到的結(jié)果進(jìn)行校正，可以精確得到真實(shí)相位值，使得對于兩幅平移量較大的待配準(zhǔn)圖像可以穩(wěn)定并準(zhǔn)確地估計出其平移參數(shù)。本文用真實(shí)光學(xué)圖像的實(shí)驗(yàn)驗(yàn)證了本文提出的SVD-IIM方法在配準(zhǔn)問題上相比于原有的SVD方法、SVD-RANSAC方法、CC-centriod方法和fast-SSDFT方法更具魯棒性和準(zhǔn)確性，具有重要的應(yīng)用價值。