薛海偉, 馮大政

(西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071)

?

解析搜索偏移量的InSAR圖像精配準方法

薛海偉, 馮大政

(西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071)

基于插值的干涉合成孔徑雷達(interferometric synthetic aperture radar,InSAR)圖像配準方法能夠達到亞像素精度,但配準精度取決于插值單元尺度,且計算量會隨精度要求提高大幅增加。提出了一種將圖像配準過程轉(zhuǎn)化為連續(xù)函數(shù)優(yōu)化問題的精配準方法。首先,該方法通過整合插值和互相關(guān)函數(shù)搜索過程構(gòu)造了一個解析代價函數(shù),在連續(xù)域內(nèi)搜索的偏移量精度不依賴于插值單元,提高了配準精度;然后利用該代價函數(shù)對于偏移量參數(shù)的梯度信息,采用擬牛頓法來優(yōu)化代價函數(shù),快速收斂并得到與最優(yōu)值對應的亞像素偏移量,具有較低的運算量。實驗結(jié)果表明該方法與傳統(tǒng)方法相比有更高的配準精度,且在計算復雜度方面有較大改善。

干涉合成孔徑雷達; 亞像素配準; 解析搜索; 互相關(guān); 圖像插值; 擬牛頓法

0　引　言

干涉合成孔徑雷達(interferometric synthetic aperture radar,InSAR)測量技術(shù)是一種利用兩幅或多幅合成孔徑雷達(synthetic aperture radar,SAR)圖像來獲取地面高程信息的遙感技術(shù)[1-6]。圖像配準是將含有相同場景或目標的圖像進行空間幾何對準的過程[7-11],是InSAR數(shù)據(jù)處理中的一個關(guān)鍵步驟。當圖像間的配準誤差超過一個像素時,構(gòu)成相位差的SAR信號間的相干性大大降低,不能獲得準確的干涉相位信息,因此InSAR圖像配準精度必須達到亞像素級[7,12]。

當?shù)匦巫兓徛蚋缮婊€不是很大時,通常用低階多項式來近似InSAR圖像間的配準函數(shù)[6]。其中,多項式參數(shù)可通過對一系列像素偏移點對進行最小二乘擬合得到。在SAR圖像中,由于相干斑的存在,基于地面特征的控制點變得不可靠[13],因此一般利用小塊圖像窗口的匹配測度來確定足夠數(shù)量的像素偏移點對[14]?，F(xiàn)有的基于圖像數(shù)據(jù)的匹配測度主要有互相關(guān)函數(shù)[5,9]、干涉頻譜信噪比[10]和相位平均波動函數(shù)[11]等。亞像素級的配準精度要求圖像窗口間的偏移量估計精度也要達到亞像素級。傳統(tǒng)的配準方法一般采用插值技術(shù)來達到亞像素級精度[4],即首先利用插值技術(shù)將圖像的像素分解為亞像素,再計算以插值單元尺寸為單位的離散網(wǎng)格點上的匹配測度值,最后根據(jù)測度最優(yōu)值的位置確定亞像素偏移量。但是,在離散域搜索獲得的亞像素偏移量并不是最優(yōu)的:一方面偏移量的計算精度取決于插值單元尺寸,插值單元尺寸越小,偏移量精度就越高,但是不能完全消除配準誤差,理論上偏移量計算精度只能達到插值單元尺寸的一半;另一方面,逐點插值運算及后續(xù)的匹配測度計算會導致很大的計算量,特別是在配準精度要求較高時[5-6];另外匹配測度最優(yōu)值的位置也會隨插值位置的不同而改變。近年來,人們提出了許多方法來實現(xiàn)亞像素配準精度。其中,文獻[1]對圖像配準方法進行了綜述性介紹。文獻[4]分析研究了干涉SAR圖像中插值方法的優(yōu)化,及其對配準性能的影響。文獻[5]采用二次B樣條函數(shù)逼近截斷Sinc函數(shù)來構(gòu)造代價函數(shù),并利用雙迭代算法交替地在水平方向和垂直方向搜索代價函數(shù)的最大值,實現(xiàn)亞像素級精度配準。文獻[6]利用分式布朗運動模型來擬合SAR圖像,并利用擬合模型的統(tǒng)計特性來提取兩幅圖像間的相對偏移量。此外,通過聯(lián)合運動平臺軌跡的精確信息及輔助地形信息也可以達到提高配準精度的目的[14]。

本文在互相關(guān)函數(shù)的基礎(chǔ)上構(gòu)造了隨亞像素偏移量連續(xù)變化的代價函數(shù),并采用擬牛頓方法來優(yōu)化代價函數(shù),實現(xiàn)了對亞像素偏移量的快速解析搜索,提高了配準精度。仿真和實測數(shù)據(jù)的實驗結(jié)果表明該方法能夠有效的實現(xiàn)圖像配準,并相對于基于插值的配準方法具有更高的配準精度和更低的運算量。

1　亞像素配準解析推導

1.1互相關(guān)函數(shù)(cross-correlation,CC)

利用InSAR測量技術(shù)可獲取兩幅SAR復圖像I1和I2,分別稱為主圖像和副圖像,配準后其干涉相位提供地面高度信息。I1和I2在點(x,y)處的復值可表示為

(1)

式中,|·|表示取模;|I1(x,y)|和|I2(x,y)|分別為主副圖像在點(x,y)的幅度;φ1(x,y)和φ2(x,y)分別為主副圖像在點(x,y)的相位。

在最大似然估計準則下,CC是計算圖像窗口間偏移量的最優(yōu)測度[15]。I1(x,y)和I2(x,y)間存在偏移(u,v)時的CC為

(2)

式中，*表示共軛;M和N分別為圖像行方向和列方向的像素點數(shù);e-jΦ(x,y)為補償局部地形的相位因子。由Cauchy-Schwarz不等式可知,互相關(guān)函數(shù)在I1與I2精確配準時達到最大值[16]。CC達到的精度由標準差σCC給出[15]

(3)

式中，γ表示相干系數(shù)。式(3)可作為選取配準窗口大小的依據(jù)。式(2)中分母為兩幅圖像的能量函數(shù),分別用P1和P2(u,v)來表示。對于局部地形補償,采用如下的方法:將局部地形假設(shè)為一個斜面,對應斜面的局部條紋頻率由干涉圖頻譜最大值的位置來確定,并在互相關(guān)之前將局部相位的復指數(shù)與主圖像相乘。則InSAR圖像配準表示為如下的優(yōu)化問題

(4)

通過計算在所有可能偏移位置上的圖像間互相關(guān)函數(shù)值,并找到最大值所對應的位置,就可以確定兩幅圖像之間的真實偏移量。在式(4)的基礎(chǔ)上,本文通過聯(lián)合插值過程構(gòu)造了代價函數(shù)在連續(xù)域的解析表達式進而實現(xiàn)對亞像素偏移量的解析搜索。

1.2解析代價函數(shù)推導

(5)

式中，φ(x,y)是插值函數(shù),其值由插值函數(shù)的形式和相對距離決定;l表示插值函數(shù)長度的一半。由式(5)可以看出,m和n的取值范圍由(x,y)點的位置和l共同決定。由于理想的Sinc插值函數(shù)需要無窮采樣點,無法實現(xiàn)。通常采用對稱和可分離的近似插值函數(shù)以降低運算量[4]。由于插值函數(shù)是可分離的,可先后進行兩個方向的一維插值,示意圖見圖1。

圖1　二維連續(xù)插值的一維分解示意圖

(6)

將式(5)代入式(6)中,可得

(7)

令x-m=p,y-n=q,對式(7)進行推導,則目標函數(shù)變?yōu)?/p>

(8)

(9)

通過上述推導,InSAR圖像亞像素精度配準可表示為如下的優(yōu)化問題

(10)

這樣,通過將插值操作和亞像素搜索過程進行聯(lián)合就構(gòu)造了一個新的解析代價函數(shù),同時亞像素偏移量搜索問題轉(zhuǎn)換為連續(xù)函數(shù)優(yōu)化問題。值得指出,經(jīng)過變量代換后的式(10)的一個顯著特點是在插值步驟之前先進行整數(shù)像素偏移時的互相關(guān)操作;同時注意到亞像素偏移參數(shù)只涉及插值這一個過程。也就是說,與式(6)不同,可以先計算圖像間存在整數(shù)像素偏移時的互相關(guān)測度,進而通過插值得到圖像間存在任意亞像素偏移時的代價函數(shù)。利用優(yōu)化方法尋找代價函數(shù)最大值的同時實現(xiàn)了對亞像素偏移量的連續(xù)搜索,提高了搜索偏移量的精確度,也意味著更高的配準精度。式(10)的另一個重要優(yōu)勢是可以直接獲得代價函數(shù)|Rdx,dy|2分別對亞像素偏移參數(shù)dx和dy的一階偏導信息,從而使擬牛頓法這一高效的優(yōu)化方法變得可行,實際上方法經(jīng)過數(shù)步迭代就可實現(xiàn)收斂,大大降低了運算量。

1.3代價函數(shù)的梯度表達式推導

式(10)中Cp,q是復值,可令

(11)

(12)

為了便于推導,令

(13)

則代價函數(shù)|Rdx,dy|2對變量dx的一階偏導為

(14)

代價函數(shù)|Rdx,dy|2對變量dy的一階偏導為

(15)

2　擬牛頓法優(yōu)化代價函數(shù)

(16)

式中，δk=xk+1-xk;γk=gk+1-gk。

基于以上分析,應用BFGS擬牛頓方法優(yōu)化|Rx|2的具體計算步驟如下。

步驟 1給定初始點x0=[0,0]T,S0=I2(2×2單位矩陣),允許誤差0<ε?1。

步驟 2置k=0,利用式(14)和式(15)計算出|Rx|2在xk處的梯度gk。

步驟 3令dk=-Skgk。從xk出發(fā),沿dk方向搜索,求步長αk,使它滿足

(17)

令δk=αkdk,xk+1=xk+δk。

步驟 4計算|Rx|2在xk+1處的梯度gk+1,若‖gk+1‖≤ε,則停止迭代,得到點xopt=xk+1;否則,進行步驟5。

步驟 5令γk=gk+1-gk。利用式(16)計算Sk+1,置k=k+1,返回步驟3。

在BFGS方法中,由當前迭代點的一階導數(shù)及可用信息構(gòu)造出來的S具有正定性,因此保證產(chǎn)生的搜索方向均為上升方向,每次迭代使代價函數(shù)有所上升,且有較高的收斂速度。

3　實驗結(jié)果與計算復雜度分析

3.1實驗結(jié)果

本節(jié)采用一組仿真和實測數(shù)據(jù)來驗證所提出方法的有效性。首先采用一組Etna火山口的仿真數(shù)據(jù)[18](數(shù)據(jù)是基于SIR-C/X-SAR在X波段產(chǎn)生)來驗證所提出亞像素偏移量估計方法的性能。為了便于實驗和分析,分別截取了512×512像素的子圖像,如圖2所示。采用互相關(guān)系數(shù)和信噪比(signal to noise ratio,SNR)兩個指標來比較配準方法的性能。SNR由干涉相位頻譜中的最大項與剩余項之和相比得到,用來衡量干涉條紋的清晰度。平均波動函數(shù)法、互相關(guān)法和最大譜估計法都是在對副圖像進行亞像素插值的基礎(chǔ)上,分別計算平均波動函數(shù)、互相關(guān)系數(shù)和SNR以確定亞像素偏移量。在上述3種方法中,亞像素插值單元尺寸均為0.1像素。將圖像對劃分為小塊圖像窗口對,窗口大小為51像素×51像素,并利用圖2(a)中所示窗口對的實驗結(jié)果來比較分析各方法的性能。實驗中計算機的配置如下:雙核2.50 GHz Pentium (R)處理器,2GB內(nèi)存,使用Matlab7.13編程。圖2(b)所示為本文方法的收斂性能。表1列出了分別由幾種方法計算得到的窗口間的偏移量、對應的互相關(guān)系數(shù)、SNR和運算時間。從表1可以看出,本文方法在得到的互相關(guān)系數(shù)和SNR上優(yōu)于其他方法,表明本文方法具有更高的偏移量搜索精度,且具有更快的運算速度。在利用本文方法得到全部可用圖像窗口間的亞像素偏移量后,通過最小二乘擬合確定所有像素點的偏移量[13],最后將副圖像進行重采樣實現(xiàn)與主圖像精確配準。圖2(d)所示為采用不同方法進行配準后的相干直方圖比較,采用平均波動函數(shù)法、互相關(guān)法、最大譜估計法和本文方法配準后得到的相干系數(shù)均值分別為0.486 2、0.466 6、0.482 9和0.491 7,表明本文方法具有較高的配準精度。

圖2　Etna數(shù)據(jù)和實驗結(jié)果

窗口對方法參數(shù)最優(yōu)偏移量互相關(guān)系數(shù)SNR/dB運算時間/s1平均波動函數(shù)法(-0.4,0.3)0.5149-30.22316.1535互相關(guān)法(-0.3,0.3)0.5202-30.15976.0830最大譜估計法(-0.3,0.3)0.5202-30.15978.0623本文方法(-0.3270,0.3334)0.5232-30.15630.89482平均波動函數(shù)法(-0.4,0.4)0.6396-29.22525.9356互相關(guān)法(-0.4,0.3)0.6375-29.21305.8070最大譜估計法(-0.3,0.3)0.6388-29.15927.9458本文方法(-0.3751,0.3201)0.6414-29.18570.7672

圖3　TerraSAR-X實測圖像

進一步地,本文方法通過實測InSAR圖像數(shù)據(jù)進行驗證。圖像對由TerraSAR-X分別在2009年2月12和23日獲得,成像區(qū)域是澳大利亞的艾爾斯巖石。我們?nèi)〕鰧煌瑘鼍疤卣鲄^(qū)域的子圖像對來進行實驗,分別是:巖石一角(區(qū)域A)、巖石一側(cè)的陡坡(區(qū)域B)和一塊平坦區(qū)域(區(qū)域C),子圖像對大小為800像素×800像素,如圖3中所示。在像素級配準的基礎(chǔ)上,亞像素精度配準過程如下:首先,分別將主副圖像分解為大小為51像素×51像素的窗口對。采用本文方法確定這些窗口對之間的亞像素偏移量。然后通過最小二乘擬合確定所有像素點的偏移量[13]。最后將副圖像進行重采樣,并進一步計算干涉相位。區(qū)域A、B和C的實驗結(jié)果如圖4所示,配準后的相干系數(shù)均值和SNR如表2所示。從表2中可以看出,采用本文方法均能夠獲得更高的相干系數(shù)均值和SNR值,表明本文方法能獲得更高的配準精度和更加清晰的干涉圖。

3.2計算復雜度分析

配準方法的計算量主要來自于亞像素插值和匹配測度計算。當圖像對大小為m×n像素,插值因子為r(如插值因子為10時對應0.1亞像素插值單元)時,互相關(guān)法計算復雜度為O(r2m2n2)。采用FFT計算互相關(guān)函數(shù)可以將計算復雜度降為O(r2mnlog2(mn))。本文提出的方法的計算復雜度為O(mnlog2(mn))+O(l2K),其中K為所采用擬牛頓法的迭代次數(shù)。從表1中可以看出,本文方法的運算時間要低于其他方法。

表2　本文方法與其他方法的相干系數(shù)、SNR比較結(jié)果

圖4　TerraSAR-X實測數(shù)據(jù)實驗結(jié)果

4　結(jié)　論

本文提出了一種干涉SAR復圖像亞像素精度配準方法。該方法構(gòu)造了可解析搜索的代價函數(shù)并獲得其梯度信息的表達式,利用擬牛頓法優(yōu)化代價函數(shù),加快了運算速度,同時實現(xiàn)了在連續(xù)域搜索偏移量,使其達到插值意義上的最優(yōu),提高了配準精度。實驗結(jié)果表明該方法有更高的配準精度,且在計算復雜度方面有較大改善。

[1] Zitová B, Flusser J. Image registration methods: a survey[J].ImageandVisionComputing, 2003, 21(3): 977-1000.

[2] Wang G L, Zhou W, Chai Y, et al. SAR image registration based on monogenic signal theory[J].JournalofElectronics&InformationTechnology, 2013, 35 (8): 1779-1785. (王國力,周偉,柴勇,等.基于單演信號的SAR圖像配準[J].電子與信息學報,2013,35(8):1779-1785.)

[3] Wang T, Jonsson S, Hanssen R F. Improved SAR image coregistration using pixel-offset series[J].IEEEGeoscienceandRemoteSensingLetters, 2014, 11(9): 1465-1469.

[4] Hanssen R, Bamler R. Evaluation of interpolation kernels for SAR interferometry[J].IEEETrans.onGeoscienceandRemoteSensing, 1999, 37(1): 318-321.

[5] Liu B Q, Feng D Z, Shui P L, et al. Analytic search method for interferometric SAR image registration[J].IEEEGeoscienceandRemoteSensingLetters, 2008, 5(2): 294-298.

[6] Danudirdjo D, Hirose A. Local subpixel coregistration of interferometric synthetic aperture radar images based on fractal models[J].IEEETrans.onGeoscienceandRemoteSensing, 2013,51(7): 4292-4301.

[7] Rosen P A, Hensley S, Joughin I R, et al. Synthetic aperture radar interferometry[J].ProceedingsofIEEE, 2000, 88(3): 333-382.

[8] Kugler F, Schulze D, Hajnsek I, et al. TanDEM-X Pol-InSAR performance for forest height estimation[J].IEEETrans.onGeoscienceandRemoteSensing, 2014, 52(10):6404-6422.

[9] Li F, Goldstein R. Studies of multibaseline spaceborne interferometric synthetic aperture radars[J].IEEETrans.onGeoscienceandRemoteSensing, 1990, 28(1): 88-97.

[10] Gabriel A K, Goldstein R M. Crossed orbit interferometry: theory and experimental results from SIR-B[J].InternationalJournalofRemoteSensing, 1988, 9(5): 857-872.

[11] Lin Q, Vesecky J F, Zebker H A, et al. New approaches in interferometric SAR data processing[J].IEEETrans.onGeoscienceandRemoteSensing, 1992, 30(3): 560-567.

[12] Chen A C, Zebker H A. Reducing ionospheric effects in InSAR data using accurate coregistration[J].IEEETrans.onGeoscienceandRemoteSensing, 2014, 52(1): 60-70.

[13] Lapini A, Bianchi T, Argenti F, et al. Blind speckle decorrelation for SAR image despeckling[J].IEEETrans.onGeoscienceandRemoteSensing, 2014, 52(2): 1044-1058.

[14] Nitti D O, Hanssen R F, Refice A, et al. Impact of DEM-assisted coregistration on high-resolution SAR interferometry[J].IEEETrans.onGeoscienceandRemoteSensing,2011,49(3):1127-1143.

[15] Zan F D. Accuracy of incoherent speckle tracking for circular Gaussian signals[J].IEEEGeoscienceandRemoteSensingLetters, 2014, 11(1): 264-267.

[16] Rosenfeld A, Kak A C.Digitalpictureprocessing[M]. New York: Academic Press, 1976: 297-302.

[17] Rohde G K, Aldroubi A, Healy D M, et al. Interpolation artifacts in sub-pixel image registration[J].IEEETrans.onImageProcessing, 2009, 18(2): 333-345.

[18] Epsilon Nought. Radar Remote Sensing.[EB/OL]. http:∥epsilon.nought.de/

Analytic offset search method for InSAR image precise registration

XUE Hai-wei, FENG Da-zheng

(National Lab of Radar Signal Processing, Xidian University, Xi’an 710071, China)

In the interpolation-based methods for subpixel registration of interferometric synthetic aperture radar (InSAR) images, but the registration accuracy is restricted by the interpolation unit, and the computational burden is heavy when high accuracy is demanded. The proposed method implements the analytical optimization of cross correlation to overcome the limitation aforementioned. Firstly, a novel analytic cost function is derived by combining the interpolation operation with the cross-correlation searching process. The offset is searched in the continuous domains, which suggests a more accurate registration since the estimation accuracy does not depend on the interpolation unit. Then, by incorporating the gradient information of the cost function, the quasi-Newton method is employed to optimize the cost function. The subpixel offset associated with the maximum of the cost function can be quickly obtained with low computational complexity. The experimental results show that the proposed method has higher registration accuracy and superior performance in terms of registration accuracy and computational complexity when compared with several conventional methods.

interferometric synthetic aperture radar (InSAR); subpixel image registration; analytic search; cross-correlation; image interpolation; quasi-Newton method

2015-09-01；

2016-02-17；網(wǎng)絡(luò)優(yōu)先出版日期：2016-06-07。

國家自然科學基金(61271293)資助課題

TN 957

A

10.3969/j.issn.1001-506X.2016.08.12

薛海偉(1985-),男,博士研究生,主要研究方向為干涉SAR數(shù)據(jù)處理、SAR成像。

E-mail:xuehw01@163.com

馮大政(1959-),男,教授,博士,主要研究方向為雷達成像、陣列信號處理、盲信號處理、神經(jīng)網(wǎng)絡(luò)。

E-mail:dzfeng@xidian.edu.cn

網(wǎng)絡(luò)優(yōu)先出版地址：http://www.cnki.net/kcms/detail/11.2422.TN.20160607.1141.008.html