一種頻域亞像元影像匹配方法

王愷 陳世平 曾

（1 北京空間機電研究所，北京 100094）

（2 中國空間技術研究院，北京 100081）

（3 中國資源衛星應用中心，北京 100094）

0 引言

攝影測量中控制點的提取是攝影測量系統檢校的基礎，也是提取物體三維信息的基礎。在數字攝影測量中以影像匹配代替傳統的人工觀測，可以達到自動確立同名像點的目的[1]。影像匹配分為基于模板和基于特征的兩種方法。基于模板的影像匹配方法又分為灰度相關法、相位相關法和互信息法等[2]。由于相位相關法計算速度快，對變化成像條件或頻域噪聲有很高的魯棒性，適合于航天遙感攝影測量，是本文將要研究的方法。

文獻[3]于 1975年提出了相位相關匹配方法，但只能求解出兩函數間的整數級平移；文獻[4-5]通過對相位相關函數進行多相位分解，將相位相關法擴展到亞像元級，但如果待匹配影像存在嚴重混疊時，匹配誤差很大；文獻[6]分析了混疊對相位匹配的影響，并對相位相關函數首先進行掩模處理，將相位差平面中受混疊影響較大的點去除，然后使用最小二乘法在頻域中直接求解平移分量；文獻[7]指出無噪聲相位相關函數是一個秩1矩陣，可以使用奇異值分解得出有噪聲相位相關函數的一個秩1子空間，進而對其在頻域中使用最小二乘法求解平移分量，這種方法對噪聲和混疊都有較高的魯棒性和精度，但計算量較大；文獻[8]在相位差求解過程中使用快速最大概率密度功率法來估計平移參數，在計算速度和可靠性間取得了一定平衡；文獻[9]引入兩點梯度法來求解平移參數，可靠性較高，運算速度較快。

從相位相關法的發展可見，其大致可分為兩類：空域Dirac Delta函數法和頻域直接求解法。前者計算過程簡單快速，對噪聲抑制性好，但只能求出整數像元的平移[3]，后經改進，雖然可以實現亞像元匹配，但受混疊影響大[9]；后者計算過程復雜，匹配結果的可靠性較前者低，但可以在一定噪聲、混疊存在的情況下實現高精度匹配。為了實現影像精確、可靠地匹配，本文提出一種將上述兩種方法有機結合的頻域匹配新方法：優化相位匹配法，先使用空域Dirac Delta函數法進行粗匹配，然后使用頻域直接求解法進行精匹配，并且在精匹配時使用優化法，進一步提高了算法的精度和可靠性。

1 優化相位匹配的原理

則

為實現亞像元精度的匹配，采用在頻域直接求解相位差平面斜率的方法。使用這種方法的關鍵有兩點：1）要正確地找出屬于相位差平面的點集；2）通過此點集估計出相位差平面斜率參數。

對輸入影像加窗及傅里葉變換后，按照式（5）計算相位相關函數，進而得到相位差。如果輸入影像中存在混疊，不等于，即幅值并不等于 1[6]。可以將的幅值作為衡量混疊大小的標志，幅值越接近1，受混疊影響越小。據此得權函數為

式中 參數α通過實驗確定，具體方法見參考文獻[6]。將權函數與相乘，即可求得受混疊影響較小的屬于相位差平面的點集。

以上得到的相位差位于[- π ,π]之間，相位差面并非平面，相位差是纏繞的。文獻[11-12]給出此時的相位差為

式中 r、c為相位差矩陣的行列標號；M、N為行列數；J為任意整數。為了得到相位差平面，需要進行相位差解纏繞。根據式（7）、（8），相位差矩陣可以按行或列進行解纏繞，從而將二維解纏繞簡化為一維解纏繞。解纏繞前需要對相位差濾波去噪，去噪效果直接決定解纏繞的品質[13]。如果要實現優良的去噪效果，則需針對具體影像的相位差設計濾波器。

2 算法實現

根據上述算法原理，先使用空域 Dirac Delta函數法對兩幅影像進行粗配準，得出整數級別的平移，真實平移必定位于中，再通過以上的優化法求取最佳亞像元平移估計值，算法流程如圖1所示。

整個算法如下：

3）在參考影像上截取子影像 p1，p1左上角在原影像中的位置為，大小為×；在待匹配影像上截取子影像p2，p2左上角在原影像中的位置為(0, 0)，大小為；

4）對p1和p2使用Blackman-Harris窗進行加窗處理，對處理后的影像進行傅里葉變換得和；

圖1 算法流程圖Fig.1 Flow chart of the algorithm

3 實驗驗證

為了驗證算法的匹配精度，本文采用0.5m分辨率的航空影像表征地面真實情況，通過濾波和降采樣生成待匹配影像對。如果待匹配影像對在原高分辨率影像中的平移差為，降采樣率為 l，則降采樣后待匹配影像對理論平移差為。選擇42幅航空影像，影像內容為村莊、農田、湖、山等，每種影像7幅，每幅影像大小為256×256像元，每像元的量化位數為8bit。降采樣率l可取任意值，本文取為4。或取整數值-2～1像元。為了測試不同混疊程度下的算法精度，采用高斯低通濾波器對原圖像進行濾波，空域高斯核函數的“標準差”σ取值為{1.0, 1.5, 2.0, 2.5, 3.0}，模板大小為9×9像元，對應頻域帶寬為{0.477 5,0.318 3,0.238 7,0.191 0,0.159 2}。標準差越大，降采樣后影像中的混疊越小，降采樣后共生成3 360幅影像。

表1、2分別列出不同景物影像在不同混疊時的平均匹配誤差和最大匹配誤差，混疊越小、景物細節越豐富，則匹配誤差越小。紋理不豐富的場景，如農田等，隨著濾波器帶寬變窄，降采樣后圖像中的混疊變少，匹配誤差隨之變小。當高斯濾波器帶寬較小、降采樣后圖像紋理減少較多時，匹配誤差隨著帶寬的減小反而由減小變為增大，見表1、2中σ=3時的情況。對基本無紋理的場景，如湖水等，自身高頻分量較少，匹配誤差隨著濾波器帶寬的減小沒有顯著下降。當混疊較小時（σ=3），基本無紋理的場景匹配誤差優于1/10個像元；紋理不豐富的場景匹配誤差優于1/20像元；紋理豐富的場景（如村莊等）匹配誤差優于1/100像元。

表1 平均匹配誤差Tab.1 The average matching errors 像元

表2 最大匹配誤差Tab.2 The maximum matching errors 像元

不同信噪比下算法的匹配誤差見表3。場景選擇為村莊，共7幅高分辨航空影像，每幅影像大小為256×256像元，降采樣率l=4，δx或δy取整數值-2～1，高斯濾波器σ=3，降采樣后圖像加入高斯白噪聲，共生成 560幅待匹配影像。隨著信噪比（采用峰值信噪比，PSNR）增大，匹配誤差隨之減小。當PSNR＞30dB，匹配誤差受PSNR影響較小，算法對存在噪聲的影像的匹配性能較好。

其他研究者的相位匹配算法的匹配誤差見表4。場景選擇為村莊，共7幅高分辨航空影像，每幅影像大小為 256×256像元，降采樣率 l=4，δx或 δy取整數值-2～1，高斯濾波器 σ=3，降采樣后圖像加入高斯白噪聲（PSNR=30dB），共生成112幅待匹配影像。除了COSI-Corr算法匹配精度和本文相近外，其他算法的精度都相對較低。在實驗中發現，COSI-Corr算法有時會出現不收斂的情況，可靠性比本文提出的算法低。

表3 不同峰值信噪比下的匹配誤差（σ=3）Tab.3 The matching errors in the different PSNR（σ=3）像元

表4 不同相位匹配算法的匹配誤差（σ=3，PSNR=30dB）Tab.4 The matching errors about the different phase matching methods（σ=3，PSNR=30dB）像元

4 結束語

經以上實驗驗證，本文提出的匹配方法對混疊有一定的抑制能力，對噪聲有較強的魯棒性。對于紋理豐富的、存在少量混疊的影像匹配誤差優于1/100個像元。本方法是一種高精度的、可靠的影像匹配方法，可以將其應用到數字攝影測量軟件中實現高精度、高可靠性地確定同名像點。

References)

[1]張祖勛, 張劍清. 數字攝影測量學[M]. 武漢: 武漢大學出版社, 2012.ZHANG Zuxun, ZHANG Jieqing. Digital Photogrammetry[M]. Wuhan: Wuhan University Press, 2012. (in Chinese)

[2]Barbara Z, Jan F. Image Registration Methods: A Survey[J]. Image and Vision Computing, 2003, 21(11): 977-1000.

[3]Kuglin C D, Hines D C. The Phase Correlation Image Alignment Method[C]. In Proc. Int. Conf. on Cybernetics and Society.New York: IEEE, 1975:163-165.

[4]Shekarforoush H, Berthod M, Zerubia J. Subpixel Image Registration by Estimating the Polyphone Decomposition of Cross Power Spectrum[C]. IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 1996: 532-537.

[5]Foroosh H, Zerubia J, Berthod M. Extension of Phase Correlation to Subpixel Registration[J]. IEEE Trans. Image Process,2002, 22(2): 188-200.

[6]Stone H S, Orchard M, Change E C, etal. A Fast Direct Fourier-based Algorithm for Subpixel Registration of Images[J].IEEE Trans. Geosci. Remote Sensing, 2001, 39(10): 2235-2243.

[7]Hoge W S. A Subspace Identification Extension to the Phase Correlation Method[J]. IEEE Trans. Medical Imaging, 2003,22(2): 277-280.

[8]Liu J G, Yan H S. Robust Phase Correlation Methods for Sub-pixel Feature Matching[C]. lst SEAS DTC Technical Conference. Edinburgh, 2006.

[9]Leprince S, Member S, Ayoub F, etal. Automatic and Precise Ortho Rectification, Coregistration, and Subpixel Correlation of Satellite Images, Application to Ground Deformation Measurements[J]. IEEE Trans. Geosci. Remote Sensing, 2007,45(6): 1529-1558.

[10]Fredric J H. On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform[J]. Proceedings of the IEEE, 1978, 66(1): 51-83.

[11]Murat B, Hassan F. Inferring Motion From the Bank Constraint of the Phase Matrix[J]. IEEE ICASSP 2005 Proc., 2005(2):925-928.

[12]Murat B, Hassan F. Estimating Subpixel Shifts Directly from the Phase Difference[C]. ICIP, 2005: 1057 -1060.

[13]Dennis C G, Mark D P. Two-dimensional Phase Unwrapping Theory, Algorithms and Software[M]. New York: John Wiley amp; Sons Inc., 1998.