CN102955157A - 一种用于干涉合成孔径雷达图像精配准的快速相干系数法 - Google Patents

Abstract

本发明涉及一种用于干涉合成孔径雷达图像精配准的快速相干系数法，其包括以下步骤：步骤1)读入经过粗配准达到像素级精度的干涉合成孔径雷达图像对；步骤2)进行图像预处理，将读入的图像对划分为若干子图像对，并对每组子图像对进行去均值化，使得每幅子图像的均值为零；步骤3)解优化问题估计亚像素级图像偏移量并建立点对应；步骤4)反演图像间的几何扭曲函数，对辅图像进行坐标变换。本发明将传统相干系数法中的偏移量搜索问题转化为一个多变量非线性优化问题，可以在不对图像进行任何过采样处理的情况下精确估计出图像间的亚像素级偏移量，具有速度快，精度高等优点。

Description

一种用于干涉合成孔径雷达图像精配准的快速相干系数法

技术领域

本发明涉及雷达信号处理领域，特别涉及一种用于干涉合成孔径雷达图像精配准的快速相干系数法。

背景技术

干涉合成孔径雷达(Interferometric Synthetic Aperture Radar，InSAR)是一种有效的雷达成像体制，其已被广泛用于数字高程反演(P.A.Rosen，S.Hensley，I.R.Joughin，F.K.Li，S.N.Madsen，E.Rodriquez and R.M.Goldstein，“Synthetic ApertureRadar Interferometry，”Proc.IEEE，vol.88，no.3，pp.333-382，Mar.2000)和地表变形测量(D.Massonnet，K.Feigl，M.Rossi and F.Adragna，“Radar Interferometric Mapping ofDeformation in the Year after the Landers Earthquake，”Nature，vol.369，pp.227-230，May 1994)等领域。为了获得干涉相位，InSAR通过同一平台上的两个天线或者同一个天线在不同的几何位置对目标的散射回波进行测量。由于两次测量几何不同，导致获得的两幅图像(即图像对)之间存在一定的几何扭曲。图像配准就是求取图像间的几何扭曲函数，使得两幅图像在几何上精确对准。

对于相对平坦区域的短基线干涉图像对，其几何扭曲函数常用低阶多项式来描述(D.O.Nitti，R.F.Hanssen，A.Refice，F.Bovenga and R.Nutricato，“Impact ofDEM-Assisted Coregistration on High-Resolution SAR Interferometry，”IEEE Trans.Geosci.Remote Sens.，vol.49，no.3，pp.1127-1143，Mar.2011)。四参数多项式模型(例如相似变换)，六参数多项式模型(仿射变换)和十二参数多项式模型(二阶多项式)等为经常使用的多项式模型(Z.Li and J.Bethel，“Image Coregistration in SARInterferometry，”in Proc.Int.Archives of Photogrammetry，Remote Sensing and SpatialInformation Sciences，Beijing，2008，pp.433-438)。

对于大起伏区域和大基线干涉情况，目标高程引起的几何扭曲将不能被忽视(E.Sansosti，P.Berardino，M.Manunta，F.Serafino and G.Fornaro，“Geometrical SAR ImageRegistration，”IEEE Trans.Geosci.Remote Sens.，vol.44，no.10，pp.2861-2870，Oct.2006)，低阶多项式关系将无法精确描述，在本申请中也不予以考虑。

为了反演几何扭曲函数，一种常用的方法是在图像间构建一系列点对应(tie pointcorrespondences)，通过对点对应的位置进行最小二乘拟合估计出多项式参数。点对应可以通过将原始图像对划分为若干组子图像对，估计出每组子图像对间的精确偏移量来构建。为了实现高精度配准以降低反演误差，InSAR图像配准精度必须优于十分之一个像素(G.Krieger，K.P.Papathanassiou and S.R.Cloude，“SpacebornePolarimetric SAR Interferometry：Performance Analysis and Mission Concepts，”EURASIP J.App.Sig.Pr，vol.2005，no.20，pp.3272-3292，2005)，这意味着估计出的偏移量精度至少优于十分之一个像素。

传统的相干系数法(F.K.Li and R.M.Goldstein，“Studies of MultibaselineSpaceborne Interferometric Synthetic Aperture Radars，”IEEE Trans.Geosci.RemoteSens.，vol.28，no.1，pp.88-97，Jan.1990)通过以下两种方法来估计图像间的亚像素级偏移量(Z.Li and J.Bethel，“Image Coregistration in SAR Interferometry，”in Proc.Int.Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，Beijing，2008，pp.433-438)：

●过采样干涉图像对(G.Rufino，A.Moccia and S.Esposito，“DEM Generation byMeans of ERS Tandem Data，”IEEE Trans.Geosci.Remote Sens.，vol.36，no.6，pp.1905-1912，Nov.1998)。即对原始图像对进行过采样，搜索使过采样后的图像对具有最大相干性的偏移量。

●过采样相干系数(F.K.Li and R.M.Goldstein，“Studies of Multibaseline SpaceborneInterferometric Synthetic Aperture Radars，”IEEE Trans.Geosci.Remote Sens.，vol.28，no.1，pp.88-97，Jan.1990)。即插值粗配准得到的相干系数以获得亚像素偏移量。

对于上述两种方法，偏移量估计精度完全取决于过采样率。为了获得高精度配准，过采样率就必须增大，由此导致的计算复杂度将增加。因此，设计出一种快速的、无需过采样的高精度图像配准算法是InSAR信号处理及应用领域一个亟需解决的问题。

发明内容

本发明的目的在于，提供一种用于干涉合成孔径雷达图像精配准的快速相干系数法，以解决现有的配准算法所存在的不足。

为了估计InSAR图像间的几何扭曲函数，通常需要建立一系列的点对应。这可以通过将原始图像对划分为若干子图像对，估计出每组子图像对间的亚像素级偏移量实现。传统的相干系数法在实现精配准时精度受限，且计算量大。为了克服这些缺点，本发明的技术方案提出一种用于InSAR图像精配准的快速相干系数法，该方法只考虑相对平坦区域的小基线干涉情况，基于双线性插值算法，其为图像处理领域最常用的插值算法，可以实现精度和计算复杂度间的最佳折中(B.Zitova and J.Flusser，“Image Registration Methods：A Survey，”Image Vision Comput.，vol.21，no.11，pp.977-1000，Oct.2003)。

如果I_ΔxΔy(x₀+Δx，y₀+Δy)为图像I中一个待插值的非整数位置像素点，I₀₀(x₀，y₀)、I₀₁(x₀，y₀+1)、I₁₀(x₀+1，y₀)和I₁₁(x₀+1，y₀+1)为与其相邻的四个整数像素位置，则由双线性插值算法，I_ΔxΔy点的值可以表示为：

I_ΔxΔy＝I₀₀(1-Δx)(1-Δy)+I₁₀Δx(1-Δy)+I₀₁(1-Δx)Δy+I₁₁ΔxΔy    (1)

其中Δx和Δy为I_ΔxΔy在x轴和y轴上偏离于I₀₀的亚像素级偏移量。

式(1)可以进一步整理为：

I_ΔxΔy＝m₁+m₂Δx+m₃Δy+m₄ΔxΔy   (2)

其中：

$\left\{\begin{array}{c}{m}_1=I_{00},m_2=I_{10}-I_{00}\\ m_3=I_{01}-I_{00},m_4=I_{00}-I_{10}-I_{01}+I_{11}\end{array}\right.---\left(3\right)$

若I_m和I_s为任意一组InSAR子图像对，其中I_m为主图像(master image)，I_s为辅图像(slave image)。精配准的目的就是寻求最优的Δx和Δy使得主图像I_m与(Δx，Δy)偏移后的辅图像I_sΔ间的相干系数达到最大值。由上式(2)，I_sΔ可以表示为：

I_sΔ＝A₀+A₁Δx+Δ₂Δy+A₃ΔxΔy    (4)

由上式(3)，参数A0、A1、A2和A3可以由二维卷积快速得以计算：

$\left\{\begin{array}{cc}{A}_0=I_s& A_1=I_s\⊗\left[\begin{array}{cc}1& -1\end{array}\right]\\ A_2=I_s\⊗\left[\begin{array}{c}1\\ -1\end{array}\right]& A_3=I_s\⊗\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}\right.---\left(5\right)$

上式中

表示二维卷积，卷积操作使得到的矩阵大小发生了变化，必须选取有效部分使得到的A矩阵大小相等并且一一对应。

对于平坦区域的短基线干涉情况，目标高程起伏引起的相位非平稳性很小，可以认为得到的图像对具有平稳性。对于长度为L的两个均值为零的平稳复高斯信号z₁和z₂，在各态历经性假设下，其相干系数定义为(R.Touzi，A.Lopes，J.Bruniquel andP.W.Vachon，“Coherence Estimation for SAR Imagery，”IEEE Trans.Geosci.RemoteSens.，vol.37，no.1，pp.135-149，Jan.1999)：

$\mathrm{Corr}(s_1,s_2)=\left|\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{s}_{1i}{s}_{2i}^{*}\right|/\sqrt{\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|s_{1i}\right|}^2\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|s_{2i}\right|}^2}---\left(6\right)$

式中“*”表示复共轭。然而，实际的SAR图像其均值不为零。因此，为了求取图像对间的相干系数，我们须对图像进行去均值处理，否则得到的相干系数将不具有平移不变性，即Corr(X，Y)＝Corr(Y，X)＝Corr(X，bY)≠Corr(X，a+bY)，a≠0，b≠0。由式(4)和(5)可以看到，双线性插值不影响图像的均值，所以去均值处理可以直接针对子图像对I_m和I_s，此时，I_m与I_sΔ间的相干系数可表示为：

$\mathrm{Corr}(I_m,I_{\mathrm{s\Δ}})=\left|\mathrm{\Σ}(I_m.\×I_{\mathrm{s\Δ}}^{*})\right|/\sqrt{\mathrm{\Σ}{\left|I_m\right|}^2\mathrm{\Σ}{\left|I_{\mathrm{s\Δ}}\right|}^2}$

$=\frac{|b_0+b_1\mathrm{\Δx}+b_2\mathrm{\Δy}+b_3\mathrm{\Δx\Δy}|}{\sqrt{c(a_0+a_1\mathrm{\Δx}+a_2\mathrm{\Δy}+a_3\mathrm{\Δ}{x}^2+a_4\mathrm{\Δ}{y}^2+a_5\mathrm{\Δx\Δy}+a_6\mathrm{\Δ}{x}^2\mathrm{\Δy}+a_7\mathrm{\Δx\Δ}{y}^2+a_8\mathrm{\Δ}{x}^2\mathrm{\Δ}{y}^2)}}---\left(7\right)$

其中“.×”表示标量乘积，式中相关参数为：

$\left\{\begin{array}{c}{b}_0=\mathrm{\Σ}(I_m.\×A_0^{*}),b_1=\mathrm{\Σ}(I_m.\×A_1^{*}),b_2=\mathrm{\Σ}(I_m.\×A_2^{*}),b_3=\mathrm{\Σ}(I_m.\×A_3^{*})\\ a_0=\mathrm{\Σ}(A_0.\×A_0^{*}),a_1=\mathrm{\Σ}(A_0.\×A_1^{*}+A_0^{*}.\×A_1),a_2=\mathrm{\Σ}(A_0.\×A_2^{*}+A_0^{*}.\×A_2)\\ a_3=\mathrm{\Σ}(A_1.\×A_1^{*}),a_4=\mathrm{\Σ}(A_2.\×A_2^{*}),a_5=\mathrm{\Σ}(A_0.\×A_3^{*}+A_0^{*}.\×A_3+A_1.\×A_2^{*}+A_1^{*}.\×A_2)\\ a_6=\mathrm{\Σ}(A_1.\×A_3^{*}+A_1^{*}.\×A_3),a_7=\mathrm{\Σ}(A_2.\×A_3^{*}+A_2^{*}.\×A_3),a_8=\mathrm{\Σ}(A_3.\×A_3^{*}),c=\mathrm{\Σ}{\left|I_m\right|}^2\end{array}\right.---\left(8\right)$

由此，精配准过程可转化为求解以下优化问题：

$\left\{\begin{array}{c}(\mathrm{\Δx},\mathrm{\Δy})=\arg \max \frac{|b_0+b_1\mathrm{\Δx}+b_2\mathrm{\Δy}+b_3\mathrm{\Δx\Δy}|}{\sqrt{c(a_0+a_1\mathrm{\Δx}+a_2\mathrm{\Δy}+a_3\mathrm{\Δ}{x}^2+a_4\mathrm{\Δ}{y}^2+a_5\mathrm{\Δx\Δy}+a_6\mathrm{\Δ}{x}^2\mathrm{\Δy}+a_7\mathrm{\Δx\Δ}{y}^2+a_8\mathrm{\Δ}{x}^2\mathrm{\Δ}{y}^2)}}\\ s.t.0\≤\mathrm{\Δx}\≤\mathrm{1,0}\≤\mathrm{\Δy}\≤1\end{array}\right.---\left(9\right)$

上式是一个带约束多变量非线性优化问题，可以通过序列二次规划算法(P.T.Boggs and J.W.Tolle，“Sequential Quadratic Programming，”in Acta Numerica，vol.4，A.Iserles，Ed.Cambridge：Cambridge University Press，1995，pp.1-51)快速得以求解。通过求解式(9)，我们可以精确得到InSAR图像对间的最优亚像素级偏移量。此时，根据子图像对中心位置间坐标偏移关系，我们可建立一对点对应。

利用上述算法对每一组子图像对进行处理，可以得到一系列的偏移量以及点对应，通过最小二乘算法对这些点对应进行多项式拟合，我们可以反演出原始InSAR图像对间的几何扭曲函数。根据得到的扭曲函数，我们对辅图像进行坐标变换，即可实现其与主图像间的几何对准。

本发明的一种用于干涉合成孔径雷达图像精配准的快速相干系数法，其包括以下步骤：

步骤1)：读入经过粗配准达到像素级精度的干涉合成孔径雷达图像对；本发明在粗配准基础上实现快速、高精度的图像亚像素级配准，即精配准。

步骤2)：进行图像预处理，将读入的图像对划分为若干子图像对，并对每组子图像对进行去均值化，使得每幅子图像的均值为零；

步骤3)：解优化问题估计亚像素级图像偏移量并建立点对应；

步骤4)：反演图像间的几何扭曲函数，对辅图像进行坐标变换。

为了得到高精度配准，所述的步骤2)将步骤1)读入粗配准后图像对划分为若干组子图像对，并对每组子图像对进行去均值化处理。本专利的重点是精配准，实现粗配准的算法很多，只要读入的图像是粗配准过了就行。

作为上述技术方案的一种改进，所述的步骤3)，通过求解一个以相干系数最大化为目标函数、以图像在两个方向上亚像素偏移为变量的多元非线性优化方程来快速、精确地估计出步骤2)所得到的预处理后的子图像对间的偏移量，并基于此偏移量构建点对应。

作为上述技术方案的一种改进，所述的步骤4)中是对步骤3)得到的点对应进行多项式拟合，反演出图像间的几何扭曲函数，并基于此几何扭曲函数，对辅图像进行坐标变换，使之与主图像在几何上对准。

作为上述技术方案的一种改进，设主图像I_m和辅图像I_s为任意一组去均值后的InSAR子图像对，所述的步骤3)通过求解优化问题估计I_m和I_s间的亚像素级偏移量的步骤包括：

31)先计算A参数：

$\left\{\begin{array}{cc}{A}_0=I_s& A_1=I_s\⊗\left[\begin{array}{cc}1& -1\end{array}\right]\\ A_2=I_s\⊗\left[\begin{array}{c}1\\ -1\end{array}\right]& A_3=I_s\⊗\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}\right.---\left(2\right)$

其中

表示二维卷积；

32)基于A参数，计算以下优化问题中涉及到的参数：

$\left\{\begin{array}{c}{b}_0=\mathrm{\Σ}(I_m.\×A_0^{*}),b_1=\mathrm{\Σ}(I_m.\×A_1^{*}),b_2=\mathrm{\Σ}(I_m.\×A_2^{*}),b_3=\mathrm{\Σ}(I_m.\×A_3^{*})\\ a_0=\mathrm{\Σ}(A_0.\×A_0^{*}),a_1=\mathrm{\Σ}(A_0.\×A_1^{*}+A_0^{*}.\×A_1),a_2=\mathrm{\Σ}(A_0.\×A_2^{*}+A_0^{*}.\×A_2)\\ a_3=\mathrm{\Σ}(A_1.\×A_1^{*}),a_4=\mathrm{\Σ}(A_2.\×A_2^{*}),a_5=\mathrm{\Σ}(A_0.\×A_3^{*}+A_0^{*}.\×A_3+A_1.\×A_2^{*}+A_1^{*}.\×A_2)\\ a_6=\mathrm{\Σ}(A_1.\×A_3^{*}+A_1^{*}.\×A_3),a_7=\mathrm{\Σ}(A_2.\×A_3^{*}+A_2^{*}.\×A_3),a_8=\mathrm{\Σ}(A_3.\×A_3^{*}),c=\mathrm{\Σ}{\left|I_m\right|}^2\end{array}\right.---\left(3\right)$

33)基于以上得到的各参数，我们利用序列二次规划算法求解以下优化问题：

$\left\{\begin{array}{c}(\mathrm{\Δx},\mathrm{\Δy})=\arg \max \frac{|b_0+b_1\mathrm{\Δx}+b_2\mathrm{\Δy}+b_3\mathrm{\Δx\Δy}|}{\sqrt{c(a_0+a_1\mathrm{\Δx}+a_2\mathrm{\Δy}+a_3\mathrm{\Δ}{x}^2+a_4\mathrm{\Δ}{y}^2+a_5\mathrm{\Δx\Δy}+a_6\mathrm{\Δ}{x}^2\mathrm{\Δy}+a_7\mathrm{\Δx\Δ}{y}^2+a_8\mathrm{\Δ}{x}^2\mathrm{\Δ}{y}^2)}}\\ s.t.0\≤\mathrm{\Δx}\≤\mathrm{1,0}\≤\mathrm{\Δy}\≤1\end{array}\right.---\left(4\right)$

式中Δx和Δy即为所求的亚像素级偏移量；

34)基于得到的偏移量，得到每组子图像对中心位置的坐标对应关系，并以此作为点对应。

本发明的优点在于，本发明提出一种用于干涉合成孔径雷达图像精配准的快速相干系数法，仅考虑到相对平坦区域的短基线干涉情况，将传统相干系数法中的偏移量搜索问题转化为一个多变量非线性优化问题，以解决传统干涉合成孔径雷达图像配准算法精配准时的大时间消耗。

InSAR高程反演或变形测量对图像的配准精度要求很高，传统的相干系数法通过过采样粗配准相干系数或者过采样图像对的方式实现图像精配准，配准精度受限于过采样率。过采样率越高，配准精度越高，而高过采样率使得处理的数据量增大，影响到配准的速度。本发明将传统相干系数法的相干搜索转化成一个以相干系数最大化为目标函数、以图像在两个方向上亚像素偏移为变量的多元非线性优化问题，其可通过序列二次规划算法快速精确求解，估计出的偏移量不再受采样率的限制。因此，该方法可在不对图像进行任何过采样处理的情况下得到两幅图像间的亚像素级偏移量，具有快速、高精度等优点。

附图说明

图1为本发明的实现方案；

图2为InSAR图像对(左：主图像，右：辅图像)；

图3为配准后的干涉条纹图；

图4为配准后的相干系数图。

具体实施步骤

为了更好地理解本发明的技术方案，以下结合实际InSAR图像对本发明的实施步骤做进一步的描述。

算法的具体实施步骤如图1所示。我们使用该方法对一组大小为1000×1000的Radarsat-2干涉图像对进行精配准，如图2所示，其中左图为主图像，右图为辅图像。两幅图像数据分别采集于2008年5月4日和28日，成像地点位于美国亚利桑那州的凤凰城(South Phoenix，Arizona，USA)。为了视觉上更好的理解所选的区域，参考谷歌地球(Google Earth)光学图像，我们可将成像区域划分为五部分：裸露的土地、建筑物、居民区、停车场和植被覆盖的土地，分别在主图像中标记为1到5。由图2可以看到，这组图像对间的几何扭曲很小，因此，我们采用四参数的相似变换模型来描述其几何扭曲。如果

为主图像与复图像间的一组点对应，那么此时的几何扭曲函数可表示为：

$\left(\begin{array}{c}{x}_s\\ y_s\end{array}\right)=\left(\begin{array}{cc}s\cos \mathrm{\θ}& -s\sin \mathrm{\θ}\\ s\sin \mathrm{\θ}& s\cos \mathrm{\θ}\end{array}\right)\left(\begin{array}{c}{x}_m\\ y_m\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)---\left(1\right)$

其中s和θ分别表示图像间的缩放和旋转因子，t_x和t_y分别表示图像在x和y方向上的平移量。式(1)具有四个自由度，因此，为了反演参数，我们至少得建立两组点对应。

本例中，我们用相干系数法对原始图像对进行粗配准，根据步骤1，读入粗配准后的图像对。根据步骤2，对读入的图像对进行预处理。为了精确反演出几何扭曲函数，我们将图像对划分为四组500×500的子图像对，并对每个子图像进行去均值处理，即减去其均值。设I_m(主图像)和I_s(辅图像)为任意一组去均值后的InSAR子图像对，按照步骤3，我们通过求解优化问题估计I_m和I_s间的亚像素级偏移量。我们先计算A参数：

$\left\{\begin{array}{cc}{A}_0=I_s& A_1=I_s\⊗\left[\begin{array}{cc}1& -1\end{array}\right]\\ A_2=I_s\⊗\left[\begin{array}{c}1\\ -1\end{array}\right]& A_3=I_s\⊗\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}\right.---\left(2\right)$

其中

表示二维卷积。基于A参数，我们计算以下优化问题中涉及到的参数：

$\left\{\begin{array}{c}{b}_0=\mathrm{\Σ}(I_m.\×A_0^{*}),b_1=\mathrm{\Σ}(I_m.\×A_1^{*}),b_2=\mathrm{\Σ}(I_m.\×A_2^{*}),b_3=\mathrm{\Σ}(I_m.\×A_3^{*})\\ a_0=\mathrm{\Σ}(A_0.\×A_0^{*}),a_1=\mathrm{\Σ}(A_0.\×A_1^{*}+A_0^{*}.\×A_1),a_2=\mathrm{\Σ}(A_0.\×A_2^{*}+A_0^{*}.\×A_2)\\ a_3=\mathrm{\Σ}(A_1.\×A_1^{*}),a_4=\mathrm{\Σ}(A_2.\×A_2^{*}),a_5=\mathrm{\Σ}(A_0.\×A_3^{*}+A_0^{*}.\×A_3+A_1.\×A_2^{*}+A_1^{*}.\×A_2)\\ a_6=\mathrm{\Σ}(A_1.\×A_3^{*}+A_1^{*}.\×A_3),a_7=\mathrm{\Σ}(A_2.\×A_3^{*}+A_2^{*}.\×A_3),a_8=\mathrm{\Σ}(A_3.\×A_3^{*}),c=\mathrm{\Σ}{\left|I_m\right|}^2\end{array}\right.---\left(3\right)$

基于以上得到的各参数，我们利用序列二次规划算法求解以下优化问题：

$\left\{\begin{array}{c}(\mathrm{\Δx},\mathrm{\Δy})=\arg \max \frac{|b_0+b_1\mathrm{\Δx}+b_2\mathrm{\Δy}+b_3\mathrm{\Δx\Δy}|}{\sqrt{c(a_0+a_1\mathrm{\Δx}+a_2\mathrm{\Δy}+a_3\mathrm{\Δ}{x}^2+a_4\mathrm{\Δ}{y}^2+a_5\mathrm{\Δx\Δy}+a_6\mathrm{\Δ}{x}^2\mathrm{\Δy}+a_7\mathrm{\Δx\Δ}{y}^2+a_8\mathrm{\Δ}{x}^2\mathrm{\Δ}{y}^2)}}\\ s.t.0\≤\mathrm{\Δx}\≤\mathrm{1,0}\≤\mathrm{\Δy}\≤1\end{array}\right.---\left(4\right)$

式中Δx和Δy即为所求的亚像素级偏移量。基于该方法我们可快速计算出每组子图像对间的最佳亚像素偏移量及其所对应的最大相干系数，结果如表1所示。基于得到的偏移量，我们可得到每组子图像对中心位置的坐标对应关系，并以此作为点对应，这样我们可得到四组点对应。基于得到的点对应，我们利用最小二乘法反演式(1)中的四个参数。然后根据得到的几何扭曲函数对原始图像对中的辅图像(图2右)进行坐标变换，使得其与主图像(图2左)在几何上对准。基于配准后的图像对，即可计算出干涉相位差和相干系数。为了降低噪声的影响，计算的干涉相位差和相干系数采用了三视估计，即利用每个像素点及其周围3×3邻域内的所有像素点的信息来估计该像素点的相位差和相干系数，得到的干涉条纹图和相干系数图分别如图3和图4所示。

表1

表1提出的快速算法和相干系数法对各组子图像对配准所得到的最佳亚像素级偏移量、最佳偏移量所对应的最大相干系数以及处理时间。

注：我们将原始的InSAR图像对分成四组子图像对，按照每组图像对在原始图像中的位置，我们分别在第一列中用左上、右上、左下和右下来标记。第三列括号内的第一个和第二个元素分别表示两幅图像在x方向和y方向上的偏移量，图像坐标系定义如图2所示。

表2提出的快速算法和相干系数法配准性能比较

表2

算法	全局相干系数	平均三视相干系数	信噪比(dB)
				相干系数法	0.0109	0.5420	-37.4674
提出的快速算法	0.0109	0.5448	-37.3571

从图3和图4可以看到，对于稳定的区域，如裸露的土地和建筑物，得到的干涉条纹很清晰，相干系数很大；对于居民区，干涉条纹变得模糊，相干系数较小，这可能由于在干涉时间基线(24天)内，居民区发生了一些小的变化而引起了时间去相干；对于停车场和植被覆盖的土地，由于时间基线内车辆随时的变动，植被的生长及晃动等引起的时间去相干非常大，以至于在这些区域我们几乎看不到干涉条纹，相干系数也很低。

为了说明提出的快速算法的优势，我们也用相干系数法对各组子图像对进行精配准：利用双线性插值算法分别对每组子图像对进行十倍过采样，搜索使得每组采样后的子图像对相干性最大时的偏移量，得到的最佳偏移量及其对应的最大相干系数如表1所示。基于得到的偏移量，构建点对应，反演几何扭曲函数，对辅图像进行坐标变换。

为了比较两种算法的配准性能，我们分别考察时间消耗、配准后图像的相干系数、配准后图像的信噪比等评价因子。时间消耗描述算法的运行速度，本实验中算法的时间消耗是通过在一台具有4.0GB内存，2.4GHz主频的电脑上，运行各个算法的Matlab代码处理各组子图像对，记录运行时间而得。两种算法处理每组子图像对需要的时间消耗如表1所示。相干性描述两幅图像间的一致性和相似性，除了表1中所列的最大相干系数，我们还计算出配准后的两幅图像间的全局相干系数和平均三视相干系数，以衡量每种算法的全局相干性和局部相干性。若I_m和I_sr为配准后的图像对，那么全局相干系数可表示为：

${\mathrm{Corr}}_{\mathrm{Global}}=\frac{\left|\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\right[I_m(m,n)-{\mathrm{\μ}}_m\left]{[I_{\mathrm{sr}}(m,n)-{\mathrm{\μ}}_{\mathrm{sr}}]}^{*}\right|}{\sqrt{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{|I_m(m,n)-{\mathrm{\μ}}_m|}^2\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{|I_{\mathrm{sr}}(m,n)-{\mathrm{\μ}}_{\mathrm{sr}}|}^2}}---\left(5\right)$

平均三视相干系数定义为：

${\stackrel{\‾}{\mathrm{Corr}}}_{3-\mathrm{Look}}=\frac{\underset{m=2}{\overset{M-1}{\mathrm{\Σ}}}\underset{n=2}{\overset{N-1}{\mathrm{\Σ}}}\frac{\left|\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{j=-1}{\overset{1}{\mathrm{\Σ}}}\right[I_m(m+i,n+j)-{\mathrm{\μ}}_{\mathrm{mi},j}\left]{[I_{\mathrm{sr}}(m+i,n+j)-{\mathrm{\μ}}_{\mathrm{sri},j}]}^{*}\right|}{\sqrt{\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{h=-1}{\overset{1}{\mathrm{\Σ}}}{|I_m(m+i,n+j)-{\mathrm{\μ}}_{\mathrm{mi},j}|}^2\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{j=-1}{\overset{1}{\mathrm{\Σ}}}{|I_{\mathrm{sr}}(m+i,n+j)-{\mathrm{\μ}}_{\mathrm{sri},j}|}^2}}}{(M-2)\×(N-2)}---\left(6\right)$

这里M和N分别为配准后的图像对在x和y方向上的大小；μ_m和μ_sr分别为I_m和I_sr的均值；μ_mi，j和μ_sri，j分别为图像I_m和I_sr中以像素点(i，j)为中心的3×3邻域内所有元素的均值。

信噪比(SNR)描述干涉条纹的清晰度，其定义为干涉条纹二维频谱中最大模值分量与其它频率分量模值之和的比值(A.K.Gabriel and R.M.Goldstein，“Crossed OrbitInterferometry：Theory and Experimental Results from SIR-B，”Int.J.Remote Sens.，vol.9，no.5，pp.857-872，May 1988)：

$\mathrm{SNR}=\frac{f_{\max }}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}f(m,n)-f_{\max }}---\left(7\right)$

两种算法配准后得到的全局相干系数、平均三视相干系数及信噪比如表2所示。由表1和表2所列的各项评价结果我们可以看到，对于子图像对的配准，提出的快速算法可以得到比相干系数法更好的相干性。整体配准效果上，提出的算法可以实现和相干系数法同样的全局相干性，但具有更好的局部相干性。在干涉条纹清晰度方面，提出的算法也优于相干系数法。最重要的是，提出的算法在运行速度上具有绝对的优势：其比相干系数法平均快了8525倍。因此我们可以得出结论：提出的快速算法可以获得更快的配准速度、更好的相干性以及更清晰的干涉条纹。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (4)

1.一种用于干涉合成孔径雷达图像精配准的快速相干系数法，其包括以下步骤：

步骤1)：读入经过粗配准达到像素级精度的干涉合成孔径雷达图像对；

步骤2)：进行图像预处理，将读入的图像对划分为若干子图像对，并对每组子图像对进行去均值化，使得每幅子图像的均值为零；

步骤3)：解优化问题估计亚像素级图像偏移量并建立点对应；

步骤4)：反演图像间的几何扭曲函数，对辅图像进行坐标变换。

2.根据权利要求1所述的用于干涉合成孔径雷达图像精配准的快速相干系数法，其特征在于，所述的步骤3)，通过求解一个以相干系数最大化为目标函数、以图像在两个方向上亚像素偏移为变量的多元非线性优化方程来快速、精确地估计出步骤2)所得到的预处理后的子图像对间的偏移量，并基于此偏移量构建点对应。

3.根据权利要求2所述的用于干涉合成孔径雷达图像精配准的快速相干系数法，其特征在于，所述的步骤4)中是对步骤3)得到的点对应进行多项式拟合，反演出图像间的几何扭曲函数，并基于此几何扭曲函数，对辅图像进行坐标变换，使之与主图像在几何上对准。

4.根据权利要求2所述的用于干涉合成孔径雷达图像精配准的快速相干系数法，其特征在于，

设主图像I_m和辅图像I_s为任意一组去均值后的InSAR子图像对，所述的步骤3)通过求解优化问题估计I_m和I_s间的亚像素级偏移量的步骤包括：

31)先计算A参数：

$\left\{\begin{array}{cc}{A}_0=I_s& A_1=I_s\⊗\left[\begin{array}{cc}1& -1\end{array}\right]\\ A_2=I_s\⊗\left[\begin{array}{c}1\\ -1\end{array}\right]& A_3=I_s\⊗\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}\right.---\left(2\right)$

其中

表示二维卷积；

32)基于A参数，计算以下优化问题中涉及到的参数：

$\left\{\begin{array}{c}{b}_0=\mathrm{\Σ}(I_m.\×A_0^{*}),b_1=\mathrm{\Σ}(I_m.\×A_1^{*}),b_2=\mathrm{\Σ}(I_m.\×A_2^{*}),b_3=\mathrm{\Σ}(I_m.\×A_3^{*})\\ a_0=\mathrm{\Σ}(A_0.\×A_0^{*}),a_1=\mathrm{\Σ}(A_0.\×A_1^{*}+A_0^{*}.\×A_1),a_2=\mathrm{\Σ}(A_0.\×A_2^{*}+A_0^{*}.\×A_2)\\ a_3=\mathrm{\Σ}(A_1.\×A_1^{*}),a_4=\mathrm{\Σ}(A_2.\×A_2^{*}),a_5=\mathrm{\Σ}(A_0.\×A_3^{*}+A_0^{*}.\×A_3+A_1.\×A_2^{*}+A_1^{*}.\×A_2)\\ a_6=\mathrm{\Σ}(A_1.\×A_3^{*}+A_1^{*}.\×A_3),a_7=\mathrm{\Σ}(A_2.\×A_3^{*}+A_2^{*}.\×A_3),a_8=\mathrm{\Σ}(A_3.\×A_3^{*}),c=\mathrm{\Σ}{\left|I_m\right|}^2\end{array}\right.---\left(3\right)$

33)基于以上得到的各参数，我们利用序列二次规划算法求解以下优化问题：

$\left\{\begin{array}{c}(\mathrm{\Δx},\mathrm{\Δy})=\arg \max \frac{|b_0+b_1\mathrm{\Δx}+b_2\mathrm{\Δy}+b_3\mathrm{\Δx\Δy}|}{\sqrt{c(a_0+a_1\mathrm{\Δx}+a_2\mathrm{\Δy}+a_3\mathrm{\Δ}{x}^2+a_4\mathrm{\Δ}{y}^2+a_5\mathrm{\Δx\Δy}+a_6\mathrm{\Δ}{x}^2\mathrm{\Δy}+a_7\mathrm{\Δx\Δ}{y}^2+a_8\mathrm{\Δ}{x}^2\mathrm{\Δ}{y}^2)}}\\ s.t.0\≤\mathrm{\Δx}\≤\mathrm{1,0}\≤\mathrm{\Δy}\≤1\end{array}\right.---\left(4\right)$

式中Δx和Δy即为所求的亚像素级偏移量；

34)基于得到的偏移量，得到每组子图像对中心位置的坐标对应关系，并以此作为点对应。

Images