王愷 陳世平 曾

（1 北京空間機(jī)電研究所，北京 100094）

（2 中國(guó)空間技術(shù)研究院，北京 100081）

（3 中國(guó)資源衛(wèi)星應(yīng)用中心，北京 100094）

0 引言

攝影測(cè)量中控制點(diǎn)的提取是攝影測(cè)量系統(tǒng)檢校的基礎(chǔ)，也是提取物體三維信息的基礎(chǔ)。在數(shù)字?jǐn)z影測(cè)量中以影像匹配代替?zhèn)鹘y(tǒng)的人工觀測(cè)，可以達(dá)到自動(dòng)確立同名像點(diǎn)的目的[1]。影像匹配分為基于模板和基于特征的兩種方法?；谀０宓挠跋衿ヅ浞椒ㄓ址譃榛叶认嚓P(guān)法、相位相關(guān)法和互信息法等[2]。由于相位相關(guān)法計(jì)算速度快，對(duì)變化成像條件或頻域噪聲有很高的魯棒性，適合于航天遙感攝影測(cè)量，是本文將要研究的方法。

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

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

1 優(yōu)化相位匹配的原理

則

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

對(duì)輸入影像加窗及傅里葉變換后，按照式（5）計(jì)算相位相關(guān)函數(shù)，進(jìn)而得到相位差。如果輸入影像中存在混疊，不等于，即幅值并不等于 1[6]?？梢詫⒌姆底鳛楹饬炕殳B大小的標(biāo)志，幅值越接近1，受混疊影響越小。據(jù)此得權(quán)函數(shù)為

式中 參數(shù)α通過實(shí)驗(yàn)確定，具體方法見參考文獻(xiàn)[6]。將權(quán)函數(shù)與相乘，即可求得受混疊影響較小的屬于相位差平面的點(diǎn)集。

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

式中 r、c為相位差矩陣的行列標(biāo)號(hào)；M、N為行列數(shù)；J為任意整數(shù)。為了得到相位差平面，需要進(jìn)行相位差解纏繞。根據(jù)式（7）、（8），相位差矩陣可以按行或列進(jìn)行解纏繞，從而將二維解纏繞簡(jiǎn)化為一維解纏繞。解纏繞前需要對(duì)相位差濾波去噪，去噪效果直接決定解纏繞的品質(zhì)[13]。如果要實(shí)現(xiàn)優(yōu)良的去噪效果，則需針對(duì)具體影像的相位差設(shè)計(jì)濾波器。

2 算法實(shí)現(xiàn)

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

整個(gè)算法如下：

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

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

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

3 實(shí)驗(yàn)驗(yàn)證

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

表1、2分別列出不同景物影像在不同混疊時(shí)的平均匹配誤差和最大匹配誤差，混疊越小、景物細(xì)節(jié)越豐富，則匹配誤差越小。紋理不豐富的場(chǎng)景，如農(nóng)田等，隨著濾波器帶寬變窄，降采樣后圖像中的混疊變少，匹配誤差隨之變小。當(dāng)高斯濾波器帶寬較小、降采樣后圖像紋理減少較多時(shí)，匹配誤差隨著帶寬的減小反而由減小變?yōu)樵龃螅姳?、2中σ=3時(shí)的情況。對(duì)基本無紋理的場(chǎng)景，如湖水等，自身高頻分量較少，匹配誤差隨著濾波器帶寬的減小沒有顯著下降。當(dāng)混疊較小時(shí)（σ=3），基本無紋理的場(chǎng)景匹配誤差優(yōu)于1/10個(gè)像元；紋理不豐富的場(chǎng)景匹配誤差優(yōu)于1/20像元；紋理豐富的場(chǎng)景（如村莊等）匹配誤差優(yōu)于1/100像元。

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

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

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

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

表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 結(jié)束語

經(jīng)以上實(shí)驗(yàn)驗(yàn)證，本文提出的匹配方法對(duì)混疊有一定的抑制能力，對(duì)噪聲有較強(qiáng)的魯棒性。對(duì)于紋理豐富的、存在少量混疊的影像匹配誤差優(yōu)于1/100個(gè)像元。本方法是一種高精度的、可靠的影像匹配方法，可以將其應(yīng)用到數(shù)字?jǐn)z影測(cè)量軟件中實(shí)現(xiàn)高精度、高可靠性地確定同名像點(diǎn)。

References)

[1]張祖勛, 張劍清. 數(shù)字?jǐn)z影測(cè)量學(xué)[M]. 武漢: 武漢大學(xué)出版社, 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.