周玫君，趙遼英，厲小潤

(1.杭州電子科技大學計算機學院，浙江 杭州 310018；2.浙江大學信息與電子工程系，浙江 杭州 310027)

0 引 言

圖像配準是影響圖像拼接、圖像變化檢測和圖像融合等圖像處理的關(guān)鍵技術(shù)[1]，廣泛應用于醫(yī)學圖像處理、高分辨率圖像生成以及圖像重建等實際應用中，要求配準精度達到亞像素級[2]?，F(xiàn)有的亞像素級圖像配準算法主要有基于灰度插值法[3]、解最優(yōu)化問題法[4]和變換域求解法[5-7]。基于灰度插值法首先通過再采樣使2幅圖像放大并且具有相同的分辨率，然后進行配準，但在再采樣時必須進行插值，而插值算法的質(zhì)量決定了插值方法的精確度。解最優(yōu)化問題法可歸結(jié)為求配準參數(shù)最小化的代價函數(shù)問題，適用于多光譜圖像及亮度變化的配準，但算法復雜度較高。變換域求解法對噪聲、光照等魯棒性強，對圖像灰度信息依賴較低，因此受到廣泛關(guān)注。早期的變換域求解法主要采用擴展相位相關(guān)法，如文獻[5]利用曲線擬合的思想將相位相關(guān)算法推廣到亞像素精度，但達到的精度有限并且只能解決平移問題。為求解旋轉(zhuǎn)和縮放參數(shù)，常用的方法是相位相關(guān)法結(jié)合擴展傅里葉梅林變換(Fourier-Mellin，F(xiàn)M)法，如文獻[6]提出FM結(jié)合sinc函數(shù)來近似表示圖像間的相位相關(guān)函數(shù)，在一定程度上提高了運動參數(shù)的估計精度，但該方法受限于sinc函數(shù)和相位相關(guān)函數(shù)之間的誤差，并且沒有考慮縮放問題。為進一步提高算法的精度，袁恒等[7]將矩陣乘法離散傅里葉變換與FM結(jié)合提出上采樣快速FM變換法，但該方法不能實現(xiàn)任意精度的亞像素參數(shù)估計。針對高精度旋轉(zhuǎn)縮放參數(shù)估計問題，本文綜合考慮配準精度和運算速度，提出一種改進的曲面擬合修正上采樣傅里葉梅林算法，在不降低執(zhí)行效率的情況下，可獲得任意高精度的亞像素級運動參數(shù)。

1 基礎(chǔ)理論

基于傅里葉變換求圖像間的旋轉(zhuǎn)、縮放和平移參數(shù)主要使用2個技術(shù)：相位相關(guān)匹配[8]和傅里葉梅林變換[9]。

1.1 相位相關(guān)匹配

設(shè)大小均為M×N圖像f1(x,y)和f2(x,y)存在如下的平移關(guān)系：

f2(x,y)=f1(x-x0,y-y0)

(1)

則f1(x,y)和f2(x,y)之間的歸一化互功率譜為：

(2)

對式(2)進行傅里葉反變換可得f1(x,y)和f2(x,y)之間的脈沖函數(shù)δ：

q(x,y)=δ(x-x0,y-y0)

(3)

搜索其峰值位置，即可確定圖像f2相對圖像f1的像素級平移參數(shù)(x0,y0)。

1.2 Fourier-Mellin變換

假設(shè)2幅圖像f1(x,y)，f2(x,y)存在平移、旋轉(zhuǎn)和縮放變換：

(4)

式中，φ，s，Δx，Δy分別為旋轉(zhuǎn)、縮放和平移參數(shù)，則f1(x,y)，f2(x,y)的傅里葉變換的幅度譜M1，M2滿足如下關(guān)系：

(5)

將幅度譜轉(zhuǎn)換到對數(shù)-極坐標中，可以得到：

M2(ε,θ)=M1(ε-d,θ-φ)

(6)

2 曲面擬合修正上采樣傅里葉梅林法

在實際的應用場景中，圖像間的運動都是連續(xù)的，如果采用第1節(jié)中的相位相關(guān)法和傅里葉梅林法估計圖像運動，只能得到對數(shù)-極坐標下的整像素級運動參數(shù)，不能解決圖像的實際問題。為了得到任意精度的亞像素級運動參數(shù)，本文提出改進的曲面擬合修正上采樣傅里葉梅林算法，首先，對預處理后的2幅圖像進行傅里葉變換，并對其歸一化互功率譜Q(u,v)執(zhí)行2倍零填充，反變換檢測其峰值坐標得到旋轉(zhuǎn)縮放矢量的1/2像素級估計值(dε2,dθ2)；然后，應用矩陣乘法離散傅里葉變換快速計算n倍上采樣相位相關(guān)峰作為旋轉(zhuǎn)縮放矢量的亞像素初始估計值；最后，采用曲面擬合法對亞像素初始估計值進行修正，進而得到任意精度的運動估計結(jié)果。

2.1 亞像素級旋轉(zhuǎn)縮放初始估計值的快速計算

(7)

2.2 曲面擬合計算旋轉(zhuǎn)縮放參數(shù)調(diào)整值

在n倍上采樣相位相關(guān)曲面上，采用二元二次多項式曲面擬合函數(shù)在以初始估計值(ε0,θ0)為中心的3×3領(lǐng)域進行擬合：

q(ε,θ)=a0+a1ε+a2θ+a3ε2+a4εθ+a5θ2

(8)

式中，q(ε,θ)表示坐標為(ε,θ)的相關(guān)函數(shù)值。將點(ε0,θ0)與其周圍8點的相關(guān)函數(shù)值代入式(8)，并利用最小二乘法得到該多項式的系數(shù)a0,a1,a2,a3,a4,a5，即可求出擬合曲面的極值點位置：

(9)

用擬合曲面的極值點坐標(Δε,Δθ)除以n的結(jié)果(dεn,dθn)對1/2像素級的估計值(dε2,dθ2)進行更新，得到對數(shù)-極坐標下的任意精度的亞像素級旋轉(zhuǎn)、縮放估計值，ε=dε2+dεn,θ=dθ2+dθn。再對ε,θ進行坐標變換公式計算得到旋轉(zhuǎn)參數(shù)φ和縮放參數(shù)s。

圖1 算法總體流程

2.3 算法總體流程

按照上述分析，算法總體流程如圖1所示，主要步驟如下：

(1)對圖像進行預處理。為解決頻譜混疊和傅里葉變換的邊緣效應對配準精度的影響，首先對圖像進行梯度[10]和加窗處理[11]。

(2)計算初始旋轉(zhuǎn)縮放像素級參數(shù)。對預處理后的2幅圖像進行傅里葉變換并對其歸一化互功率譜執(zhí)行2倍零填充，反變換檢測其峰值坐標即可得到旋轉(zhuǎn)縮放矢量的1/2像素級估計值(dε2,dθ2)。

(3)求n倍上采樣相位相關(guān)峰。利用矩陣乘法離散傅里葉變換[9]快速計算n倍上采樣相位相關(guān)曲面，求得的峰值坐標(ε0,θ0)即為坐標反變換前的旋轉(zhuǎn)縮放亞像素級初始估計參數(shù)。

(4)曲面擬合計算旋轉(zhuǎn)縮放參數(shù)調(diào)整值。在上采樣相位相關(guān)曲面上，對以峰值坐標為中心的3×3區(qū)域進行曲面擬合，修正峰值以得到任意精度的亞像素級調(diào)整值(dεn,dθn)。

(5)使用(dεn,dθn)對初始估計值(dε2,dθ2)進行修正，從而實現(xiàn)任意精度級別的旋轉(zhuǎn)縮放運動估計。

最后根據(jù)對數(shù)-極坐標和笛卡爾坐標間的映射，計算笛卡爾坐標下的旋轉(zhuǎn)縮放參數(shù)，對待配準圖進行反變換，再利用基于擬合的擴展相位相關(guān)法求亞像素平移參數(shù)。

3 實驗結(jié)果與分析

分別采用仿真圖像配準實驗和真實圖像配準實驗來驗證本文算法的性能。選用FM結(jié)合sinc曲線擬合的方法[6](簡稱FM_sinc)、快速上采樣FM變換法[7](簡稱Fast_FM)和本文提出的曲面擬合修正上采樣FM算法進行實驗比較。由于本文算法的重點在于旋轉(zhuǎn)和縮放參數(shù)的精確估計，因此，實驗中所有方法采用相同的算法估計平移參數(shù)，預處理步驟一致。實驗環(huán)境配置為：Inter CPU(i7-7700,3.60 GHz)、8.0 GB內(nèi)存、64位的Windows10操作系統(tǒng)。

3.1 仿真序列圖像實驗

圖2(a)為某地區(qū)上空航拍遙感影像中截取的大小為330×330的子圖。以圖2(a)為參考圖，取隨機值(3.6，2.4)平移變換后，依次對其進行步長為5°的旋轉(zhuǎn)操作，縮放系數(shù)每隔2幅圖增加0.1，共得到6幅序列圖像，圖2(b)和(c)分別為第3和第6幅序列圖。

圖2 仿真實驗的部分圖像

用配準后的均方根誤差(Root Mean Square Error，RMSE)評價配準精度，其表達式為：

(10)

圖3 3種算法的配準均方根誤差

分別采用3種算法求得的配準后圖像的RRMSE如圖3所示。經(jīng)計算，F(xiàn)M_sinc，F(xiàn)ast_FM和本文算法的平均RRMSE分別為0.862 0，0.427 3和0.154 1，本文算法的配準精度分別提高了0.707 9像素和0.273 2像素。從圖3可以看出：FM_sinc算法的RRMSE最大，并隨著圖像旋轉(zhuǎn)縮放的增大而增大。Fast_FM算法對旋轉(zhuǎn)縮放參數(shù)的估計精度有明顯的提升，但其精度受限于n的取值，只能得到1/n精度估計結(jié)果。而本文算法在1/n精度估計的基礎(chǔ)上，利用擬合修正估計值，實現(xiàn)了任意精度估計結(jié)果。

實驗中，通過計算可得，F(xiàn)M_sinc，F(xiàn)ast_FM和本文算法的平均運行時間分別為0.763 7 s，0.790 1 s和0.805 3 s，算法運行時間分別增加約0.041 6 s和0.015 2 s。說明本文算法能在增加較少運行時間的情況下獲得更高的配準精度。

仿真圖像在平移固定，旋轉(zhuǎn)縮放(15°，1.2)下，采用3種算法配準后的圖像與參考圖像之間的聯(lián)合概率直方圖如圖4所示。

圖4 旋轉(zhuǎn)縮放(15°，1.2)時，3種算法配準結(jié)果的聯(lián)合直方圖

從圖4可以看出：圖4(a)中有大量的點分散在對角線周圍，相比圖4(a)，圖4(b)中有部分點更靠近對角線，圖4(c)中點的分布最靠近直角線。聯(lián)合直方圖中的點分布越靠近對角線，誤差越小，圖4結(jié)果表明本文算法的配準誤差最小。

3.2 真實圖像實驗

圖5 真實圖像

選取某地區(qū)2001年拍攝的SPOT4圖像和2004年拍攝的Landsat TM圖像作為真實圖像進行實驗。分別從中心截取大小為330×330的子圖作為參考圖和待配準圖，如圖5所示。

由于真實遙感圖像間的旋轉(zhuǎn)、縮放和平移量未知，無法對結(jié)果客觀評價，所以根據(jù)配準后的拼接效果進行主觀評價。圖6給出3種算法配準后的拼接結(jié)果圖，其中每幅圖右邊的小圖依次為左圖中從左到右的圓圈中拼接區(qū)域的局部放大圖。從圖6可以看出，3種算法都得到了很好的配準結(jié)果，但仔細分析各圖的第一個局部放大區(qū)域可以看出，F(xiàn)M_sinc求得的配準后圖像在拼接處有完全錯開的斷層，F(xiàn)ast_FM算法求得的圖像拼接處斷層錯開近1/4，本文算法比Fast_FM的拼接效果稍微好一些，斷層處錯開了近1/2。由此可見，本文算法的配準效果最佳。

圖6 3種算法配準拼接結(jié)果

4 結(jié)束語

針對亞像素級的旋轉(zhuǎn)、縮放運動的高精度估計問題，本文提出一種基于曲面擬合修正矩陣乘法的傅里葉梅林變換算法，計算結(jié)果表明本文算法的計算精度誤差為0.15像素級別，是一種計算精度高的配準方法，具有廣泛的應用前景。如何進一步提高計算精度，增強其魯棒性，減少計算時間是下一步研究的重點。