解析搜索偏移量的InSAR图像精配准方法

薛海伟, 冯大政

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



解析搜索偏移量的InSAR图像精配准方法

薛海伟, 冯大政

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

基于插值的干涉合成孔径雷达(interferometric synthetic aperture radar,InSAR)图像配准方法能够达到亚像素精度,但配准精度取决于插值单元尺度,且计算量会随精度要求提高大幅增加。提出了一种将图像配准过程转化为连续函数优化问题的精配准方法。首先,该方法通过整合插值和互相关函数搜索过程构造了一个解析代价函数,在连续域内搜索的偏移量精度不依赖于插值单元,提高了配准精度;然后利用该代价函数对于偏移量参数的梯度信息,采用拟牛顿法来优化代价函数,快速收敛并得到与最优值对应的亚像素偏移量,具有较低的运算量。实验结果表明该方法与传统方法相比有更高的配准精度,且在计算复杂度方面有较大改善。

干涉合成孔径雷达; 亚像素配准; 解析搜索; 互相关; 图像插值; 拟牛顿法

0　引　言

干涉合成孔径雷达(interferometric synthetic aperture radar,InSAR)测量技术是一种利用两幅或多幅合成孔径雷达(synthetic aperture radar,SAR)图像来获取地面高程信息的遥感技术[1-6]。图像配准是将含有相同场景或目标的图像进行空间几何对准的过程[7-11],是InSAR数据处理中的一个关键步骤。当图像间的配准误差超过一个像素时,构成相位差的SAR信号间的相干性大大降低,不能获得准确的干涉相位信息,因此InSAR图像配准精度必须达到亚像素级[7,12]。

当地形变化缓慢或干涉基线不是很大时,通常用低阶多项式来近似InSAR图像间的配准函数[6]。其中,多项式参数可通过对一系列像素偏移点对进行最小二乘拟合得到。在SAR图像中,由于相干斑的存在,基于地面特征的控制点变得不可靠[13],因此一般利用小块图像窗口的匹配测度来确定足够数量的像素偏移点对[14]。现有的基于图像数据的匹配测度主要有互相关函数[5,9]、干涉频谱信噪比[10]和相位平均波动函数[11]等。亚像素级的配准精度要求图像窗口间的偏移量估计精度也要达到亚像素级。传统的配准方法一般采用插值技术来达到亚像素级精度[4],即首先利用插值技术将图像的像素分解为亚像素,再计算以插值单元尺寸为单位的离散网格点上的匹配测度值,最后根据测度最优值的位置确定亚像素偏移量。但是,在离散域搜索获得的亚像素偏移量并不是最优的:一方面偏移量的计算精度取决于插值单元尺寸,插值单元尺寸越小,偏移量精度就越高,但是不能完全消除配准误差,理论上偏移量计算精度只能达到插值单元尺寸的一半;另一方面,逐点插值运算及后续的匹配测度计算会导致很大的计算量,特别是在配准精度要求较高时[5-6];另外匹配测度最优值的位置也会随插值位置的不同而改变。近年来,人们提出了许多方法来实现亚像素配准精度。其中,文献[1]对图像配准方法进行了综述性介绍。文献[4]分析研究了干涉SAR图像中插值方法的优化,及其对配准性能的影响。文献[5]采用二次B样条函数逼近截断Sinc函数来构造代价函数,并利用双迭代算法交替地在水平方向和垂直方向搜索代价函数的最大值,实现亚像素级精度配准。文献[6]利用分式布朗运动模型来拟合SAR图像,并利用拟合模型的统计特性来提取两幅图像间的相对偏移量。此外,通过联合运动平台轨迹的精确信息及辅助地形信息也可以达到提高配准精度的目的[14]。

本文在互相关函数的基础上构造了随亚像素偏移量连续变化的代价函数,并采用拟牛顿方法来优化代价函数,实现了对亚像素偏移量的快速解析搜索,提高了配准精度。仿真和实测数据的实验结果表明该方法能够有效的实现图像配准,并相对于基于插值的配准方法具有更高的配准精度和更低的运算量。

1　亚像素配准解析推导

1.1互相关函数(cross-correlation,CC)

利用InSAR测量技术可获取两幅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是计算图像窗口间偏移量的最优测度[15]。I1(x,y)和I2(x,y)间存在偏移(u,v)时的CC为

(2)

式中，*表示共轭;M和N分别为图像行方向和列方向的像素点数;e-jΦ(x,y)为补偿局部地形的相位因子。由Cauchy-Schwarz不等式可知,互相关函数在I1与I2精确配准时达到最大值[16]。CC达到的精度由标准差σCC给出[15]

(3)

式中，γ表示相干系数。式(3)可作为选取配准窗口大小的依据。式(2)中分母为两幅图像的能量函数,分别用P1和P2(u,v)来表示。对于局部地形补偿,采用如下的方法:将局部地形假设为一个斜面,对应斜面的局部条纹频率由干涉图频谱最大值的位置来确定,并在互相关之前将局部相位的复指数与主图像相乘。则InSAR图像配准表示为如下的优化问题

(4)

通过计算在所有可能偏移位置上的图像间互相关函数值,并找到最大值所对应的位置,就可以确定两幅图像之间的真实偏移量。在式(4)的基础上,本文通过联合插值过程构造了代价函数在连续域的解析表达式进而实现对亚像素偏移量的解析搜索。

1.2解析代价函数推导

(5)

式中，φ(x,y)是插值函数,其值由插值函数的形式和相对距离决定;l表示插值函数长度的一半。由式(5)可以看出,m和n的取值范围由(x,y)点的位置和l共同决定。由于理想的Sinc插值函数需要无穷采样点,无法实现。通常采用对称和可分离的近似插值函数以降低运算量[4]。由于插值函数是可分离的,可先后进行两个方向的一维插值,示意图见图1。

图1　二维连续插值的一维分解示意图

(6)

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

(7)

令x-m=p,y-n=q,对式(7)进行推导,则目标函数变为

(8)

(9)

通过上述推导,InSAR图像亚像素精度配准可表示为如下的优化问题

(10)

这样,通过将插值操作和亚像素搜索过程进行联合就构造了一个新的解析代价函数,同时亚像素偏移量搜索问题转换为连续函数优化问题。值得指出,经过变量代换后的式(10)的一个显著特点是在插值步骤之前先进行整数像素偏移时的互相关操作;同时注意到亚像素偏移参数只涉及插值这一个过程。也就是说,与式(6)不同,可以先计算图像间存在整数像素偏移时的互相关测度,进而通过插值得到图像间存在任意亚像素偏移时的代价函数。利用优化方法寻找代价函数最大值的同时实现了对亚像素偏移量的连续搜索,提高了搜索偏移量的精确度,也意味着更高的配准精度。式(10)的另一个重要优势是可以直接获得代价函数|Rdx,dy|2分别对亚像素偏移参数dx和dy的一阶偏导信息,从而使拟牛顿法这一高效的优化方法变得可行,实际上方法经过数步迭代就可实现收敛,大大降低了运算量。

1.3代价函数的梯度表达式推导

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

(11)

(12)

为了便于推导,令

(13)

则代价函数|Rdx,dy|2对变量dx的一阶偏导为

(14)

代价函数|Rdx,dy|2对变量dy的一阶偏导为

(15)

2　拟牛顿法优化代价函数

(16)

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

基于以上分析,应用BFGS拟牛顿方法优化|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出发,沿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方法中,由当前迭代点的一阶导数及可用信息构造出来的S具有正定性,因此保证产生的搜索方向均为上升方向,每次迭代使代价函数有所上升,且有较高的收敛速度。

3　实验结果与计算复杂度分析

3.1实验结果

本节采用一组仿真和实测数据来验证所提出方法的有效性。首先采用一组Etna火山口的仿真数据[18](数据是基于SIR-C/X-SAR在X波段产生)来验证所提出亚像素偏移量估计方法的性能。为了便于实验和分析,分别截取了512×512像素的子图像,如图2所示。采用互相关系数和信噪比(signal to noise ratio,SNR)两个指标来比较配准方法的性能。SNR由干涉相位频谱中的最大项与剩余项之和相比得到,用来衡量干涉条纹的清晰度。平均波动函数法、互相关法和最大谱估计法都是在对副图像进行亚像素插值的基础上,分别计算平均波动函数、互相关系数和SNR以确定亚像素偏移量。在上述3种方法中,亚像素插值单元尺寸均为0.1像素。将图像对划分为小块图像窗口对,窗口大小为51像素×51像素,并利用图2(a)中所示窗口对的实验结果来比较分析各方法的性能。实验中计算机的配置如下:双核2.50 GHz Pentium (R)处理器,2GB内存,使用Matlab7.13编程。图2(b)所示为本文方法的收敛性能。表1列出了分别由几种方法计算得到的窗口间的偏移量、对应的互相关系数、SNR和运算时间。从表1可以看出,本文方法在得到的互相关系数和SNR上优于其他方法,表明本文方法具有更高的偏移量搜索精度,且具有更快的运算速度。在利用本文方法得到全部可用图像窗口间的亚像素偏移量后,通过最小二乘拟合确定所有像素点的偏移量[13],最后将副图像进行重采样实现与主图像精确配准。图2(d)所示为采用不同方法进行配准后的相干直方图比较,采用平均波动函数法、互相关法、最大谱估计法和本文方法配准后得到的相干系数均值分别为0.486 2、0.466 6、0.482 9和0.491 7,表明本文方法具有较高的配准精度。

图2　Etna数据和实验结果

窗口对方法参数最优偏移量互相关系数SNR/dB运算时间/s1平均波动函数法(-0.4,0.3)0.5149-30.22316.1535互相关法(-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平均波动函数法(-0.4,0.4)0.6396-29.22525.9356互相关法(-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图像数据进行验证。图像对由TerraSAR-X分别在2009年2月12和23日获得,成像区域是澳大利亚的艾尔斯岩石。我们取出对应不同场景特征区域的子图像对来进行实验,分别是:岩石一角(区域A)、岩石一侧的陡坡(区域B)和一块平坦区域(区域C),子图像对大小为800像素×800像素,如图3中所示。在像素级配准的基础上,亚像素精度配准过程如下:首先,分别将主副图像分解为大小为51像素×51像素的窗口对。采用本文方法确定这些窗口对之间的亚像素偏移量。然后通过最小二乘拟合确定所有像素点的偏移量[13]。最后将副图像进行重采样,并进一步计算干涉相位。区域A、B和C的实验结果如图4所示,配准后的相干系数均值和SNR如表2所示。从表2中可以看出,采用本文方法均能够获得更高的相干系数均值和SNR值,表明本文方法能获得更高的配准精度和更加清晰的干涉图。

3.2计算复杂度分析

配准方法的计算量主要来自于亚像素插值和匹配测度计算。当图像对大小为m×n像素,插值因子为r(如插值因子为10时对应0.1亚像素插值单元)时,互相关法计算复杂度为O(r2m2n2)。采用FFT计算互相关函数可以将计算复杂度降为O(r2mnlog2(mn))。本文提出的方法的计算复杂度为O(mnlog2(mn))+O(l2K),其中K为所采用拟牛顿法的迭代次数。从表1中可以看出,本文方法的运算时间要低于其他方法。

表2　本文方法与其他方法的相干系数、SNR比较结果

图4　TerraSAR-X实测数据实验结果

4　结　论

本文提出了一种干涉SAR复图像亚像素精度配准方法。该方法构造了可解析搜索的代价函数并获得其梯度信息的表达式,利用拟牛顿法优化代价函数,加快了运算速度,同时实现了在连续域搜索偏移量,使其达到插值意义上的最优,提高了配准精度。实验结果表明该方法有更高的配准精度,且在计算复杂度方面有较大改善。

[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；网络优先出版日期：2016-06-07。

国家自然科学基金(61271293)资助课题

TN 957

A

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

薛海伟(1985-),男,博士研究生,主要研究方向为干涉SAR数据处理、SAR成像。

E-mail:xuehw01@163.com

冯大政(1959-),男,教授,博士,主要研究方向为雷达成像、阵列信号处理、盲信号处理、神经网络。

E-mail:dzfeng@xidian.edu.cn

网络优先出版地址：http://www.cnki.net/kcms/detail/11.2422.TN.20160607.1141.008.html