薛海偉,馮大政

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

一種快速的In SAR圖像精配準方法

薛海偉,馮大政

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

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

干涉合成孔徑雷達;亞像素配準;信噪比;圖像插值;雙迭代算法

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

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

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

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

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

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

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

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

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

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

2.1解析代價函數(shù)的構造

圖1　兩種匹配測度比較圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3　TerraSAR-X實測干涉SAR圖像對

圖4　區(qū)域A實驗結(jié)果

圖5　區(qū)域B實驗結(jié)果

圖6　區(qū)域C實驗結(jié)果

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

4　結(jié)束語

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

[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圖像配準[J].電子與信息學報,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.

(編輯:齊淑娟)

簡 訊

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

摘自《西電科大報》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

國家自然科學基金資助項目(61271293)

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