一种快速的In SAR图像精配准方法

薛海伟,冯大政

(西安电子科技大学雷达信号处理国家重点实验室,陕西西安 710071)

一种快速的In SAR图像精配准方法

薛海伟,冯大政

(西安电子科技大学雷达信号处理国家重点实验室,陕西西安 710071)

针对传统的配准方法通过对图像过采样插值来实现亚像素精度配准,但是搜索精度受限于插值单元尺寸,且计算量会随精度要求的提高而大幅增加的问题,提出了一种快速解析搜索亚像素偏移量的图像精配准方法.首先,该方法以干涉频谱信噪比为匹配测度,通过整合过采样插值和匹配测度搜索过程,将亚像素偏移量搜索过程转化为连续函数优化问题,实现了对亚像素偏移量的解析搜索,提高了配准精度;然后,采用双迭代算法对代价函数进行优化,能够快速准确地获得代价函数最大值对应的偏移量.实验结果表明,该方法不仅具有较高的配准精度,且在计算复杂度方面也有较大改善.

干涉合成孔径雷达;亚像素配准;信噪比;图像插值;双迭代算法

干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)是一种利用两幅或多幅合成孔径雷达(Synthetic Aperture Radar,SAR)图像来获取地面数字高程模型(Digital Elevation Model,DEM)的遥感技术[1-2].图像配准是将含有相同场景或目标的SAR图像进行空间几何对准的过程,是InSAR数据处理中的一个关键步骤[3-9].由于SAR图像存在相干斑,基于特征点的配准方法不再适用于亚像素精度的干涉SAR配准[10].实际中主要利用小块图像对的相似性测度来确定足够数量的像素偏移点对,并采用低阶多项式来近似干涉SAR图像间的配准函数[11].现有的基于图像数据的匹配测度主要有互相关函数[1,3,9]、相位平均波动函数[4]和干涉频谱信噪比[5]等.其中,互相关函数是计算图像间偏移量的最大似然估计方法,但是在计算互相关函数之前需要对局部条纹频率进行补偿,具有一定的局限性[12].相位平均波动函数[4]是基于缠绕相位进行的,在一定程度上降低了偏移量的估计精度.干涉频谱信噪比由复干涉图频谱中的最大值与剩余频谱模值之和相比得到,潜在地补偿了局部条纹频率.为保证干涉相位测量的精度,干涉SAR图像配准的精度要达到亚像素级[1,13].传统方法一般采用插值技术来达到亚像素级精度[14],但是在由插值获得的离散网格上进行离散搜索获得的亚像素偏移量是亚最优的.一方面偏移量的计算精度取决于插值单元尺寸,插值单元尺寸越小,偏移量精度就越高,但是不能完全消除配准误差,理论上偏移量计算精度只能达到插值单元尺寸的一半;另一方面,逐点插值运算及后续的匹配测度计算会导致很大的计算量,特别是在配准精度要求较高时[8];另外匹配测度最大点的位置也会随插值位置的不同而改变[3].

近年来,人们提出了许多方法来提高配准精度.其中,文献[3]基于互相关函数,采用二次B样条函数逼近截断Sinc函数来构造代价函数,并利用双迭代算法交替地在水平方向和垂直方向搜索代价函数的最大值,实现亚像素级精度配准.文献[14]分析研究了干涉SAR图像配准中几种常用的插值函数及其对配准精度的影响.文献[9]从原始图像计算互相关系数,利用插值得到互相关系数曲面进而确定与峰值对应的亚像素偏移量.文献[8]利用分式布朗运动模型来拟合SAR图像,并利用拟合模型的统计特性来提取两幅图像间的相对偏移量.此外,通过联合运动平台轨迹的精确信息及辅助地形信息也可达到提高配准精度的目的[11].笔者提出了一种快速解析搜索亚像素偏移量的方法.该方法以干涉频谱信噪比作为匹配测度,通过整合插值操作和匹配测度搜索过程构造一个可解析的代价函数,并采用双迭代算法在连续域内搜索代价函数的极值,得到对应的亚像素偏移量.由于该代价函数对噪声具有良好的稳健性,双迭代算法能够快速收敛,大大降低了运算量,同时在连续域搜索提高了偏移量的搜索精度,因而能够实现更高的配准精度.

1　干涉频谱经验信噪比

当两幅干涉SAR图像之间准确配准时,其相位差图像即干涉相位图中会出现干涉条纹.干涉SAR通过对干涉相位图进行处理来提取地面高程信息,因此配准后干涉相位图的质量可作为配准精度的准则[4-5],其中干涉相位图的质量由干涉图频谱的经验信噪比(Signal-to-Noise Ratio,SNR)来衡量.

假设两幅需要配准的SAR复图像分别为I1(x,y)和I2(x,y),当两幅图像之间存在整数像素偏移(u,v)时,其干涉图表示为

其中,上标*表示共轭.对R(x,y;u,v)进行傅里叶变换,得到二维干涉图频谱的表达式为

其中,M和N分别表示干涉图频谱的行数和列数.干涉图频谱峰值处的频率代表了干涉图中的主要条纹频率,由复干涉图频谱中主要条纹频率的模值与其余频率模值之和的比值即可得到SNR,其表达式为

当两幅图像精确配准时,经验信噪比达到最大值;反之,则较小.这样通过计算在所有可能偏移位置上的SNR,并根据最大SNR值对应的偏移位置,就能够确定图像对之间的真实偏移量.

由于SAR成像将目标以散射点模型表示,由大量散射点随机组成的分辨单元受到相干斑干扰,具有较大的能量和较高的空间频率,给干涉SAR图像配准带来困难.与互相关函数相比,干涉频谱信噪比具有更高的准确性及对相干斑噪声的稳健性.图1给出了对存在偏移的干涉SAR图像进行配准时的两种匹配测度比较,采用一组Etna火山口的仿真数据[15].从图1可看出,与对局部条纹频率补偿后的互相关函数相比,信噪比具有更加尖锐的峰值和更低的旁瓣,更有利于后续的快速优化搜索.

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

2.1解析代价函数的构造

图1　两种匹配测度比较图

其中,ϕ(x,y)是插值函数,lϕ表示插值函数长度的一半.由于理想的插值Sinc函数需要无穷采样点,无法实现.通常采用对称和可分离的近似插值函数,这里采用4×4三次卷积插值函数[14].

在干涉SAR中,当地形变化平缓或基线较短时,散射体是局部平稳的[4],图像窗口内其他像素也可用式(4)的形式表示.和之间存在亚像素偏移(dx,dy)时的干涉图可表示为

将R(x,y;dx,dy)进行二维离散傅里叶变换,得到复干涉图频谱为

值得指出的是,与式(5)相比,式(9)的突出特点是在插值过程之前先计算图像间的复干涉图频谱,其亚像素偏移量dx和dy只涉及后面的插值过程.同时随dx和dy连续变化.也就是说,可预先计算图像间存在整数像素偏移时的复干涉图频谱,再利用少量插值操作就可得到图像间存在任意亚像素偏移量时的频谱,大幅度降低代价函数的计算量.另一方面,由于插值过程是后续进行的,可根据精度要求对插值函数进行选择,具有较高的灵活性.

根据式(3),图像间存在亚像素偏移量(dx,dy)时的SNR可计算如下:

基于以上分析,这里建立了一个包含亚像素偏移量dx和dy的代价函数,其最大值位置对应最优的亚像素偏移量.

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

如上所述,亚像素偏移量搜索过程可转化为如下连续函数优化问题:

图1(a)所示为采用一组仿真数据[15]进行亚像素偏移量搜索时的代价函数示意图,代价函数的最大值位置对应最优的亚像素偏移量,可采用双迭代算法来解决这个最大值寻优问题.双迭代算法的思想来源于循环最大算法[16],分别在行方向与列方向上交替搜索代价函数最大值,使得代价函数值交替上升.双迭代算法采用下面的迭代形式:

其中,第1个不等式由式(12a)得到,后面的不等式由式(12b)得到.这样,与代价函数最大值对应的偏移量dx和dy就被交替地搜索,直到收敛到两个固定值,即为最优的亚像素偏移量.

3　实验结果与分析

在像素级配准的基础上,利用文中方法确定图像窗口对之间的亚像素偏移量,并从配准精度和计算量两方面,将文中方法与几种传统方法进行了对比.实验结果表明,该方法具有更好的性能.

图2　Etna数据和实验结果

首先,采用一组Etna火山口的仿真数据[15](数据基于SIR-C/X-SAR在X波段产生)来验证所提出亚像素偏移量估计方法的性能.为便于实验和分析,分别截取了512×512像素的图像,如图2所示.互相关法、平均波动函数法和最大谱估计法都是在对副图像进行亚像素插值的基础上,分别计算在插值网格点上的匹配测度,进而确定亚像素偏移量.在上述方法中,亚像素插值单元尺寸均为0.1像素.将文中方法的性能与上述方法进行比较,首先将图像对划分为小块图像窗口对,窗口大小为51×51像素,采用文中方法计算窗口对之间的亚像素偏移量,利用图2(a)中所示窗口对的实验结果来比较分析各方法的性能.实验中计算机的配置如下:双核2.50 GHz Pentium(R)处理器,2 GB内存,使用Matlab7.13编程.图2(b)所示为双迭代优化算法的收敛性能.表1列出了分别由几种方法计算得到的窗口间的偏移量、对应的相干系数均值、RSNR和运算时间.

表1　不同方法得到的偏移量、相干系数均值、RSNR和运算时间比较

从表1可以看出,文中方法得到的偏移量对应的相干系数均值和RSNR均优于其他几种方法的,表明文中方法具有更高的偏移量搜索精度,同时运算时间大大降低.在计算得到所有窗口对间的偏移量并去除坏点后,利用最小二乘拟合就能获得全部像素点的亚像素偏移量[11],进而将副图像进行重采样与主图像实现配准.图2(b)所示为采用不同方法进行配准后的相干直方图比较,采用最大谱估计法、平均波动函数法、互相关法和文中方法配准后得到的相干系数均值分别为0.493 5、0.488 7、0.490 4和0.504 6,表明文中方法实现了对亚像素偏移量的解析搜索,具有更高的配准精度,更重要的是,运算速度也得到大幅度提升.

进一步采用文中方法进行实测干涉SAR图像亚像素精度配准.图像对由TerraSAR-X分别在2009年2月12和23日获得,成像区域是澳大利亚的艾尔斯岩石,如图3所示.笔者取出对应不同场景特征区域的子图像对来进行实验,分别是:岩石一角(区域A)、岩石一侧的陡坡(区域B)和一块平坦区域(区域C),如图3(a)所示.子图像对大小为900×900像素.在像素级配准的基础上,亚像素精度配准过程如下:首先,分别将主副图像划分为大小51×5 1像素的窗口对,采用文中方法确定这些窗口对之间的亚像素偏移量(由于与表1中窗口大小一样,运算时间可参考表1);然后,对这些偏移量进行最小均方拟合来确定每一个像素间的偏移关系[11].最后根据偏移关系对副图像重采样,并进一步计算干涉相位.区域A、B和C的实验结果分别如图4、图5和图6所示,配准后的相干系数均值和RSNR如表2所示.从表2可看出,上述3个区域采用文中方法均能够获得更高的相干系数均值和RSNR,表明文中方法能获得更高的配准精度和更加清晰的干涉图;同时文中方法具有更快的运算速度.

图3　TerraSAR-X实测干涉SAR图像对

图4　区域A实验结果

图5　区域B实验结果

图6　区域C实验结果

表2　文中方法与其他方法的相干系数均值、RSNR比较

4　结束语

笔者提出了一种快速解析搜索亚像素偏移量的图像精配准方法.该方法以干涉频谱信噪比作为匹配测度,将亚像素偏移量搜索过程转化为连续函数优化问题,实现了对亚像素偏移量的解析搜索;并采用双迭代算法对代价函数进行优化,能够快速准确地获得代价函数最大值对应的亚像素偏移量.实验结果表明,该方法不仅具有较高的配准精度,在计算复杂度方面也有较大改善.

[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.

(编辑:齐淑娟)

简 讯

我校先进材料与纳米科技学院举办材料科学发展论坛.2015年12月11日～12日,先进材料与纳米科技学院第二届材料科学发展论坛在我校召开.期间,中国科学院金属研究所成会明院士、西安交通大学汪宏教授、湘潭大学周益春教授、吉林大学孙洪波教授和清华大学任天令教授等分别作了题为《高质量石墨烯及其他二维材料的CVD生长》、《新型无机电介质材料研究》、《耦合下铁电电畴演变规律:实验观测及理论分析》、《飞秒激光微纳制造》和《基于石墨烯材料的新型器件》的报告.

摘自《西电科大报》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.