薛海偉,馮大政

(西安電子科技大學(xué)雷達(dá)信號處理國家重點(diǎn)實(shí)驗(yàn)室,陜西西安 710071)

一種快速的In SAR圖像精配準(zhǔn)方法

薛海偉,馮大政

(西安電子科技大學(xué)雷達(dá)信號處理國家重點(diǎn)實(shí)驗(yàn)室,陜西西安 710071)

針對傳統(tǒng)的配準(zhǔn)方法通過對圖像過采樣插值來實(shí)現(xiàn)亞像素精度配準(zhǔn),但是搜索精度受限于插值單元尺寸,且計(jì)算量會(huì)隨精度要求的提高而大幅增加的問題,提出了一種快速解析搜索亞像素偏移量的圖像精配準(zhǔn)方法.首先,該方法以干涉頻譜信噪比為匹配測度,通過整合過采樣插值和匹配測度搜索過程,將亞像素偏移量搜索過程轉(zhuǎn)化為連續(xù)函數(shù)優(yōu)化問題,實(shí)現(xiàn)了對亞像素偏移量的解析搜索,提高了配準(zhǔn)精度;然后,采用雙迭代算法對代價(jià)函數(shù)進(jìn)行優(yōu)化,能夠快速準(zhǔn)確地獲得代價(jià)函數(shù)最大值對應(yīng)的偏移量.實(shí)驗(yàn)結(jié)果表明,該方法不僅具有較高的配準(zhǔn)精度,且在計(jì)算復(fù)雜度方面也有較大改善.

干涉合成孔徑雷達(dá);亞像素配準(zhǔn);信噪比;圖像插值;雙迭代算法

干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar,InSAR)是一種利用兩幅或多幅合成孔徑雷達(dá)(Synthetic Aperture Radar,SAR)圖像來獲取地面數(shù)字高程模型(Digital Elevation Model,DEM)的遙感技術(shù)[1-2].圖像配準(zhǔn)是將含有相同場景或目標(biāo)的SAR圖像進(jìn)行空間幾何對準(zhǔn)的過程,是InSAR數(shù)據(jù)處理中的一個(gè)關(guān)鍵步驟[3-9].由于SAR圖像存在相干斑,基于特征點(diǎn)的配準(zhǔn)方法不再適用于亞像素精度的干涉SAR配準(zhǔn)[10].實(shí)際中主要利用小塊圖像對的相似性測度來確定足夠數(shù)量的像素偏移點(diǎn)對,并采用低階多項(xiàng)式來近似干涉SAR圖像間的配準(zhǔn)函數(shù)[11].現(xiàn)有的基于圖像數(shù)據(jù)的匹配測度主要有互相關(guān)函數(shù)[1,3,9]、相位平均波動(dòng)函數(shù)[4]和干涉頻譜信噪比[5]等.其中,互相關(guān)函數(shù)是計(jì)算圖像間偏移量的最大似然估計(jì)方法,但是在計(jì)算互相關(guān)函數(shù)之前需要對局部條紋頻率進(jìn)行補(bǔ)償,具有一定的局限性[12].相位平均波動(dòng)函數(shù)[4]是基于纏繞相位進(jìn)行的,在一定程度上降低了偏移量的估計(jì)精度.干涉頻譜信噪比由復(fù)干涉圖頻譜中的最大值與剩余頻譜模值之和相比得到,潛在地補(bǔ)償了局部條紋頻率.為保證干涉相位測量的精度,干涉SAR圖像配準(zhǔn)的精度要達(dá)到亞像素級[1,13].傳統(tǒng)方法一般采用插值技術(shù)來達(dá)到亞像素級精度[14],但是在由插值獲得的離散網(wǎng)格上進(jìn)行離散搜索獲得的亞像素偏移量是亞最優(yōu)的.一方面偏移量的計(jì)算精度取決于插值單元尺寸,插值單元尺寸越小,偏移量精度就越高,但是不能完全消除配準(zhǔn)誤差,理論上偏移量計(jì)算精度只能達(dá)到插值單元尺寸的一半;另一方面,逐點(diǎn)插值運(yùn)算及后續(xù)的匹配測度計(jì)算會(huì)導(dǎo)致很大的計(jì)算量,特別是在配準(zhǔn)精度要求較高時(shí)[8];另外匹配測度最大點(diǎn)的位置也會(huì)隨插值位置的不同而改變[3].

近年來,人們提出了許多方法來提高配準(zhǔn)精度.其中,文獻(xiàn)[3]基于互相關(guān)函數(shù),采用二次B樣條函數(shù)逼近截?cái)郤inc函數(shù)來構(gòu)造代價(jià)函數(shù),并利用雙迭代算法交替地在水平方向和垂直方向搜索代價(jià)函數(shù)的最大值,實(shí)現(xiàn)亞像素級精度配準(zhǔn).文獻(xiàn)[14]分析研究了干涉SAR圖像配準(zhǔn)中幾種常用的插值函數(shù)及其對配準(zhǔn)精度的影響.文獻(xiàn)[9]從原始圖像計(jì)算互相關(guān)系數(shù),利用插值得到互相關(guān)系數(shù)曲面進(jìn)而確定與峰值對應(yīng)的亞像素偏移量.文獻(xiàn)[8]利用分式布朗運(yùn)動(dòng)模型來擬合SAR圖像,并利用擬合模型的統(tǒng)計(jì)特性來提取兩幅圖像間的相對偏移量.此外,通過聯(lián)合運(yùn)動(dòng)平臺軌跡的精確信息及輔助地形信息也可達(dá)到提高配準(zhǔn)精度的目的[11].筆者提出了一種快速解析搜索亞像素偏移量的方法.該方法以干涉頻譜信噪比作為匹配測度,通過整合插值操作和匹配測度搜索過程構(gòu)造一個(gè)可解析的代價(jià)函數(shù),并采用雙迭代算法在連續(xù)域內(nèi)搜索代價(jià)函數(shù)的極值,得到對應(yīng)的亞像素偏移量.由于該代價(jià)函數(shù)對噪聲具有良好的穩(wěn)健性,雙迭代算法能夠快速收斂,大大降低了運(yùn)算量,同時(shí)在連續(xù)域搜索提高了偏移量的搜索精度,因而能夠?qū)崿F(xiàn)更高的配準(zhǔn)精度.

1　干涉頻譜經(jīng)驗(yàn)信噪比

當(dāng)兩幅干涉SAR圖像之間準(zhǔn)確配準(zhǔn)時(shí),其相位差圖像即干涉相位圖中會(huì)出現(xiàn)干涉條紋.干涉SAR通過對干涉相位圖進(jìn)行處理來提取地面高程信息,因此配準(zhǔn)后干涉相位圖的質(zhì)量可作為配準(zhǔn)精度的準(zhǔn)則[4-5],其中干涉相位圖的質(zhì)量由干涉圖頻譜的經(jīng)驗(yàn)信噪比(Signal-to-Noise Ratio,SNR)來衡量.

假設(shè)兩幅需要配準(zhǔn)的SAR復(fù)圖像分別為I1(x,y)和I2(x,y),當(dāng)兩幅圖像之間存在整數(shù)像素偏移(u,v)時(shí),其干涉圖表示為

其中,上標(biāo)*表示共軛.對R(x,y;u,v)進(jìn)行傅里葉變換,得到二維干涉圖頻譜的表達(dá)式為

其中,M和N分別表示干涉圖頻譜的行數(shù)和列數(shù).干涉圖頻譜峰值處的頻率代表了干涉圖中的主要條紋頻率,由復(fù)干涉圖頻譜中主要條紋頻率的模值與其余頻率模值之和的比值即可得到SNR,其表達(dá)式為

當(dāng)兩幅圖像精確配準(zhǔn)時(shí),經(jīng)驗(yàn)信噪比達(dá)到最大值;反之,則較小.這樣通過計(jì)算在所有可能偏移位置上的SNR,并根據(jù)最大SNR值對應(yīng)的偏移位置,就能夠確定圖像對之間的真實(shí)偏移量.

由于SAR成像將目標(biāo)以散射點(diǎn)模型表示,由大量散射點(diǎn)隨機(jī)組成的分辨單元受到相干斑干擾,具有較大的能量和較高的空間頻率,給干涉SAR圖像配準(zhǔn)帶來困難.與互相關(guān)函數(shù)相比,干涉頻譜信噪比具有更高的準(zhǔn)確性及對相干斑噪聲的穩(wěn)健性.圖1給出了對存在偏移的干涉SAR圖像進(jìn)行配準(zhǔn)時(shí)的兩種匹配測度比較,采用一組Etna火山口的仿真數(shù)據(jù)[15].從圖1可看出,與對局部條紋頻率補(bǔ)償后的互相關(guān)函數(shù)相比,信噪比具有更加尖銳的峰值和更低的旁瓣,更有利于后續(xù)的快速優(yōu)化搜索.

2　亞像素偏移量的快速解析搜索算法

2.1解析代價(jià)函數(shù)的構(gòu)造

圖1　兩種匹配測度比較圖

其中,?(x,y)是插值函數(shù),l?表示插值函數(shù)長度的一半.由于理想的插值Sinc函數(shù)需要無窮采樣點(diǎn),無法實(shí)現(xiàn).通常采用對稱和可分離的近似插值函數(shù),這里采用4×4三次卷積插值函數(shù)[14].

在干涉SAR中,當(dāng)?shù)匦巫兓骄徎蚧€較短時(shí),散射體是局部平穩(wěn)的[4],圖像窗口內(nèi)其他像素也可用式(4)的形式表示.和之間存在亞像素偏移(dx,dy)時(shí)的干涉圖可表示為

將R(x,y;dx,dy)進(jìn)行二維離散傅里葉變換,得到復(fù)干涉圖頻譜為

值得指出的是,與式(5)相比,式(9)的突出特點(diǎn)是在插值過程之前先計(jì)算圖像間的復(fù)干涉圖頻譜,其亞像素偏移量dx和dy只涉及后面的插值過程.同時(shí)隨dx和dy連續(xù)變化.也就是說,可預(yù)先計(jì)算圖像間存在整數(shù)像素偏移時(shí)的復(fù)干涉圖頻譜,再利用少量插值操作就可得到圖像間存在任意亞像素偏移量時(shí)的頻譜,大幅度降低代價(jià)函數(shù)的計(jì)算量.另一方面,由于插值過程是后續(xù)進(jìn)行的,可根據(jù)精度要求對插值函數(shù)進(jìn)行選擇,具有較高的靈活性.

根據(jù)式(3),圖像間存在亞像素偏移量(dx,dy)時(shí)的SNR可計(jì)算如下:

基于以上分析,這里建立了一個(gè)包含亞像素偏移量dx和dy的代價(jià)函數(shù),其最大值位置對應(yīng)最優(yōu)的亞像素偏移量.

2.2亞像素偏移量的解析搜索

如上所述,亞像素偏移量搜索過程可轉(zhuǎn)化為如下連續(xù)函數(shù)優(yōu)化問題:

圖1(a)所示為采用一組仿真數(shù)據(jù)[15]進(jìn)行亞像素偏移量搜索時(shí)的代價(jià)函數(shù)示意圖,代價(jià)函數(shù)的最大值位置對應(yīng)最優(yōu)的亞像素偏移量,可采用雙迭代算法來解決這個(gè)最大值尋優(yōu)問題.雙迭代算法的思想來源于循環(huán)最大算法[16],分別在行方向與列方向上交替搜索代價(jià)函數(shù)最大值,使得代價(jià)函數(shù)值交替上升.雙迭代算法采用下面的迭代形式:

其中,第1個(gè)不等式由式(12a)得到,后面的不等式由式(12b)得到.這樣,與代價(jià)函數(shù)最大值對應(yīng)的偏移量dx和dy就被交替地搜索,直到收斂到兩個(gè)固定值,即為最優(yōu)的亞像素偏移量.

3　實(shí)驗(yàn)結(jié)果與分析

在像素級配準(zhǔn)的基礎(chǔ)上,利用文中方法確定圖像窗口對之間的亞像素偏移量,并從配準(zhǔn)精度和計(jì)算量兩方面,將文中方法與幾種傳統(tǒng)方法進(jìn)行了對比.實(shí)驗(yàn)結(jié)果表明,該方法具有更好的性能.

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

首先,采用一組Etna火山口的仿真數(shù)據(jù)[15](數(shù)據(jù)基于SIR-C/X-SAR在X波段產(chǎn)生)來驗(yàn)證所提出亞像素偏移量估計(jì)方法的性能.為便于實(shí)驗(yàn)和分析,分別截取了512×512像素的圖像,如圖2所示.互相關(guān)法、平均波動(dòng)函數(shù)法和最大譜估計(jì)法都是在對副圖像進(jìn)行亞像素插值的基礎(chǔ)上,分別計(jì)算在插值網(wǎng)格點(diǎn)上的匹配測度,進(jìn)而確定亞像素偏移量.在上述方法中,亞像素插值單元尺寸均為0.1像素.將文中方法的性能與上述方法進(jìn)行比較,首先將圖像對劃分為小塊圖像窗口對,窗口大小為51×51像素,采用文中方法計(jì)算窗口對之間的亞像素偏移量,利用圖2(a)中所示窗口對的實(shí)驗(yàn)結(jié)果來比較分析各方法的性能.實(shí)驗(yàn)中計(jì)算機(jī)的配置如下:雙核2.50 GHz Pentium(R)處理器,2 GB內(nèi)存,使用Matlab7.13編程.圖2(b)所示為雙迭代優(yōu)化算法的收斂性能.表1列出了分別由幾種方法計(jì)算得到的窗口間的偏移量、對應(yīng)的相干系數(shù)均值、RSNR和運(yùn)算時(shí)間.

表1　不同方法得到的偏移量、相干系數(shù)均值、RSNR和運(yùn)算時(shí)間比較

從表1可以看出,文中方法得到的偏移量對應(yīng)的相干系數(shù)均值和RSNR均優(yōu)于其他幾種方法的,表明文中方法具有更高的偏移量搜索精度,同時(shí)運(yùn)算時(shí)間大大降低.在計(jì)算得到所有窗口對間的偏移量并去除壞點(diǎn)后,利用最小二乘擬合就能獲得全部像素點(diǎn)的亞像素偏移量[11],進(jìn)而將副圖像進(jìn)行重采樣與主圖像實(shí)現(xiàn)配準(zhǔn).圖2(b)所示為采用不同方法進(jìn)行配準(zhǔn)后的相干直方圖比較,采用最大譜估計(jì)法、平均波動(dòng)函數(shù)法、互相關(guān)法和文中方法配準(zhǔn)后得到的相干系數(shù)均值分別為0.493 5、0.488 7、0.490 4和0.504 6,表明文中方法實(shí)現(xiàn)了對亞像素偏移量的解析搜索,具有更高的配準(zhǔn)精度,更重要的是,運(yùn)算速度也得到大幅度提升.

進(jìn)一步采用文中方法進(jìn)行實(shí)測干涉SAR圖像亞像素精度配準(zhǔn).圖像對由TerraSAR-X分別在2009年2月12和23日獲得,成像區(qū)域是澳大利亞的艾爾斯巖石,如圖3所示.筆者取出對應(yīng)不同場景特征區(qū)域的子圖像對來進(jìn)行實(shí)驗(yàn),分別是:巖石一角(區(qū)域A)、巖石一側(cè)的陡坡(區(qū)域B)和一塊平坦區(qū)域(區(qū)域C),如圖3(a)所示.子圖像對大小為900×900像素.在像素級配準(zhǔn)的基礎(chǔ)上,亞像素精度配準(zhǔn)過程如下:首先,分別將主副圖像劃分為大小51×5 1像素的窗口對,采用文中方法確定這些窗口對之間的亞像素偏移量(由于與表1中窗口大小一樣,運(yùn)算時(shí)間可參考表1);然后,對這些偏移量進(jìn)行最小均方擬合來確定每一個(gè)像素間的偏移關(guān)系[11].最后根據(jù)偏移關(guān)系對副圖像重采樣,并進(jìn)一步計(jì)算干涉相位.區(qū)域A、B和C的實(shí)驗(yàn)結(jié)果分別如圖4、圖5和圖6所示,配準(zhǔn)后的相干系數(shù)均值和RSNR如表2所示.從表2可看出,上述3個(gè)區(qū)域采用文中方法均能夠獲得更高的相干系數(shù)均值和RSNR,表明文中方法能獲得更高的配準(zhǔn)精度和更加清晰的干涉圖;同時(shí)文中方法具有更快的運(yùn)算速度.

圖3　TerraSAR-X實(shí)測干涉SAR圖像對

圖4　區(qū)域A實(shí)驗(yàn)結(jié)果

圖5　區(qū)域B實(shí)驗(yàn)結(jié)果

圖6　區(qū)域C實(shí)驗(yàn)結(jié)果

表2　文中方法與其他方法的相干系數(shù)均值、RSNR比較

4　結(jié)束語

筆者提出了一種快速解析搜索亞像素偏移量的圖像精配準(zhǔn)方法.該方法以干涉頻譜信噪比作為匹配測度,將亞像素偏移量搜索過程轉(zhuǎn)化為連續(xù)函數(shù)優(yōu)化問題,實(shí)現(xiàn)了對亞像素偏移量的解析搜索;并采用雙迭代算法對代價(jià)函數(shù)進(jìn)行優(yōu)化,能夠快速準(zhǔn)確地獲得代價(jià)函數(shù)最大值對應(yīng)的亞像素偏移量.實(shí)驗(yàn)結(jié)果表明,該方法不僅具有較高的配準(zhǔn)精度,在計(jì)算復(fù)雜度方面也有較大改善.

[1]LI F,GOLDSTEIN R.Studies of Multibaseline Spaceborne Interferometric Synthetic Aperture Radars[J].IEEE Transactions on Geoscience and Remote Sensing,1990,28(1):88-97.

[2]KUGLER F,SCHULZE D,HAJNSEK I,et al.TanDEM-X Pol-InSAR Performance for Forest Height Estimation[J]. IEEE Transactions on Geoscience and Remote Sensing,2014,52(10):6404-6422.

[3]LIU B Q,FENG D Z,SHUI P L,et al.Analytic Search Method for Interferometric SAR Image Registration[J].IEEE Geoscience and Remote Sensing Letters,2008,5(2):294-298.

[4]LIN Q,VESECKY J F,ZEBKER H A,et al.New Approaches in Interferometric SAR Data Processing[J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(3):560-567.

[5]GABRIEL A K,GOLDSTEIN R M.Crossed Orbit Interferometry:Theory and Experimental Results from SIR-B[J]. International Journal of Remote Sensing,1988,9(5):857-872.

[6]王國力,周偉,柴勇,等.基于單演信號的SAR圖像配準(zhǔn)[J].電子與信息學(xué)報(bào),2013,35(8):1779-1785. WANG Guoli,ZHOU Wei,CHAI Yong,et al.SAR Image Registration Based on Monogenic Signal Theory[J]. Journal of Electronics&Information Technology,2013,35(8):1779-1785.

[7]WANG T,JONSSON S,HANSSEN R F.Improved SAR Image Coregistration Using Pixel-offset Series[J].IEEE Geoscience and Remote Sensing Letters,2014,11(9):1465-1469.

[8]DANUDIRDJO D,HIROSE A.Local Subpixel Coregistration of Interferometric Synthetic Aperture Radar Images Based on Fractal Models[J].IEEE Transactions on Geoscience and Remote Sensing,2013,51(7):4292-4301.

[9]YAGUE-MARTINEZ N,EINEDER M,CONG X Y,et al.Ground Displacement Measurement by TerraSAR-X Image Correlation:The 2011 Tohoku-Oki Earthquake[J].IEEE Geoscience and Remote Sensing Letters,2012,9(4): 539-543.

[10]LAPINI A,BIANCHI T,ARGENTI F,et al.Blind Speckle Decorrelation for SAR Image Despeckling[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(2):1044-1058.

[11]NITTI D O,HANSSEN R F,REFICE A,et al.Impact of DEM-assisted Coregistration on High-resolution SAR Interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(3):1127-1143.

[12]BALMLER R,EINEDER M.Accuracy of Differential Shift Estimation by Correlation and Split-bandwidth Interferometry for Wideband and Delta-k SAR Systems[J].IEEE Geoscience and Remote Sensing Letters,2005,2(2): 151-155.

[13]CHEN A C,ZEBKER H A.Reducing Ionospheric Effects in InSAR Data Using Accurate Coregistration[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(1):60-70.

[14]HANSSEN R,BAMLER R.Evaluation of Interpolation Kernels for SAR Interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,1999,37(1):318-321.

[15]EPSILON NOUGHT.Radar Remote Sensing:Sample Data[DB/OL].[2014-10-21].http://epsilon.nought.de/.

[16]STOICA P,SELEN Y.Cyclic Minimizers,Majorization Techniques,and the Expectation Maximization Algorithm:a Refresher[J].IEEE Signal Processing Magazine,2004,21(1):112-114.

(編輯:齊淑娟)

簡 訊

我校先進(jìn)材料與納米科技學(xué)院舉辦材料科學(xué)發(fā)展論壇.2015年12月11日～12日,先進(jìn)材料與納米科技學(xué)院第二屆材料科學(xué)發(fā)展論壇在我校召開.期間,中國科學(xué)院金屬研究所成會(huì)明院士、西安交通大學(xué)汪宏教授、湘潭大學(xué)周益春教授、吉林大學(xué)孫洪波教授和清華大學(xué)任天令教授等分別作了題為《高質(zhì)量石墨烯及其他二維材料的CVD生長》、《新型無機(jī)電介質(zhì)材料研究》、《耦合下鐵電電疇演變規(guī)律:實(shí)驗(yàn)觀測及理論分析》、《飛秒激光微納制造》和《基于石墨烯材料的新型器件》的報(bào)告.

摘自《西電科大報(bào)》2015.12.28

Fast subpixel registration method for InSAR images

XUE Haiwei,FENG Dazheng
(National Key Lab.of Radar Signal Processing,Xidian Univ.,Xi’an 710071,China)

In conventional registration methods,image interpolation is used to extract subpixel offsets which lead to the limitation that the registration accuracy is restricted by the interpolation unit.Moreover,the computational burden is heavy when high accuracy is demanded.In order to reduce the computational load,we propose an efficient offset estimation method for interferometric synthetic aperture radar(InSAR) image subpixel registration.Firstly,a novel cost function continuously varying with offsets is established by integrating the empirical signal-to-noise ratio(SNR)and the interpolation operation.This suggests a more accurate registration since the accuracy of the estimated offsets does not depend on the interpolation unit. Secondly,an efficient bi-iterative algorithm is employed to solve the cost function in the continuous domain. The subpixel offsets associated with the maximum of the cost function can be exactly obtained with low computational complexity.Simulated and real data are tested to illustrate the good performance and computational efficiency of the proposed method.

interferometric synthetic aperture radar(InSAR);subpixel image registration;signal-tonoise ratio(SNR);image interpolation;bi-iterative algorithm

TN957

A

1001-2400(2016)03-0172-07

10.3969/j.issn.1001-2400.2016.03.030

2015-08-15

國家自然科學(xué)基金資助項(xiàng)目(61271293)

薛海偉(1985-),男,西安電子科技大學(xué)博士研究生,E-mail:xuehw001@163.com.