CN102270341B - 一种自适应的高精度干涉sar相位估计方法 - Google Patents

Abstract

本发明公开了一种自适应的高精度干涉SAR相位估计方法，它结合维纳滤波器理论构造最佳加权矢量，通过对其组成的最佳协方差矩阵特征分解得到信号子空间和噪声子空间，充分利用相应像素对及其相邻像素的相干信息，并根据MUSIC算法中信号子空间和噪声子空间的正交性构造空间谱函数，通过谱峰搜索准确地估计相应像素间的干涉相位。它不需要确定配准误差大小和方向来得到最优权值，仅仅利用维纳滤波器获得最优权值，从而解决了传统InSAR干涉相位估计运算量大的问题。本发明适用于InSAR复杂场景地表参数精确反演等领域。

Description

一种自适应的高精度干涉SAR相位估计方法

技术领域：

本技术发明属于雷达技术领域，它特别涉及了干涉合成孔径雷达(InSAR)图像配准技术领域。

背景技术：

干涉合成孔径雷达(InSAR)技术是目前SAR遥感技术发展的重要方向，能够全天时、全天候获取大面积数字高程图像。它不仅具有对散射体的空间分布和高度的敏感性，而且具有成本低、连续性和远程遥感探测的能力。InSAR可以提取观测对象的空间三维结构特征信息，为微波定量遥感、高精度数字高程信息和观察对象细微形变信息的提取提供了可能性。

干涉合成孔径雷达技术中的三大关键处理步骤是SAR图像配准、干涉相位估计和相位展开。SAR图像配准的目地是使配准后两幅SAR图像中的一对像素对应于地面上同一分辨单元。当图像配准的精度较差时，两幅图像同一位置的像素可能对应于地面上不同的散射体，其相位差不能反映地面高度起伏的情况，将使后面的相位展开难以获得令人满意的结果。配准精度的高低直接影响干涉相位图的质量，不精确的配准会导致干涉相位估计及展开不准确。根据本人了解以及已发表的文献，目前InSAR干涉相位估计精度很高的方法有：Hai L,Guisheng L”An estimation method for InSAR interferometric phase combined with imageauto-coregistration”Science in China,Series F,2006.Zhenfang L,Zheng B”Imageauto-Coregistration and InSAR interferogram estimation using joint subspace Projection”IEEETrans On GRS,2006，是一种基于加权联合单像素模型的空间投影干涉相位估计方法，能够根据配准误差大小和方向来构造最优加权观测矢量，同时利用相邻像素的相关信息，具有自适应图像配准和降低相位噪声功能，并且可以在SAR图像配准精度很差(可以允许达到一个分辨单元)的条件下准确的估计相应像素间的干涉相位。

然而该方法在估计干涉相位时虽然在估计干涉相位不需要确定噪声子空间维数，但需要首先确定配准误差方向和大小，并且确定最优权值时需要搜索，因此计算量较大。

发明内容：

本发明为了克服InSAR干涉相位估计精度不高及运算量大的问题，提供了一种自适应的高精度干涉SAR相位估计方法。该方法结合维纳滤波器理论构造最佳加权矢量，通过对其组成的最佳协方差矩阵特征分解得到信号子空间和噪声子空间，充分利用相应像素对及其相邻像素的相干信息，并根据MUSIC算法中信号子空间和噪声子空间的正交性构造空间谱函数，通过谱峰搜索准确地估计相应像素间的干涉相位。该方法不需要首先确定配准误差大小和方向来得到最优权值，仅仅利用维纳滤波器获得最优权值，从而很好地解决了InSAR干涉相位估计运算量大的问题。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、干涉合成孔径雷达(InSAR)

干涉SAR(Interferometric SAR)是近年来得到迅速发展的地表形变探测新技术，具有大面积生成地表数字地形高度模型的能力。InSAR主要用于获取地物的空间垂直结构信息，通过该技术可以获取两个重要的参数分别为干涉相位和相干系数。InSAR信号处理是在经典SAR成像的基础上，再经过二副图像的配准、干涉相位估计、相位展开等关键步骤，得到地表三维信息。详见L.C.Graham.“Synthetic interferometer radar for topographic mapping”,Proceedings of the IEEE,1974

定义2、干涉相位

干涉相位是两幅InSAR复图像的相位差值，通过将配准后的两幅复图像共轭相乘得到干涉相位。

定义3、相干系数

相干系数是用来衡量主副图像之间的相似程度和干涉图质量最基本、最直接的衡量标准，既包含幅度相似性的信息，也包含相位相似性的信息。M₁和M₂表示零均值的复高斯随机变量，则相干系数γ定义为：

$\mathrm{\γ}=\frac{\left|E\right[M_1{M_2}^{*}\left]\right|}{\sqrt{E\left[{\left|M_1\right|}^2\right]E\left[{\left|M_2\right|}^2\right]}}$

式中，*表示取共轭运算，E[]表示统计平均值，||表示求模值。γ为相干系数，表示两个复变量之间的互相关系数，其值越大说明配准精度越高。

定义4、图像配准

所谓图像配准，就是用来计算干涉相位的两幅复图像的一对像素点必须对应地面上同一分辨单元。在图像没有配准的情况下，两幅SAR图像中同一位置的像素可能对应于地面上不同的散射体，其相位差不能反映地面高度起伏的真实情况。详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社。

定义5、粗配准

粗配准是指在主图像和副图像上选取相同的区域进行像素级的配准，图像配准的精度达到一个分辨单元。传统的图像粗配准的方法有很多种，如相大相干系数法，其基本思想是根据图像能量互相关函数的统计特性，通过寻找图像间复相关函数的最大值来进行配准，首先以其中的一幅图像作为模板，称为主图像，另外一幅作为副图像，然后通过主图像在副图像中上下左右的搜索，获得一系列的相干系数值，其最大值所对应的位置即为两幅复图像正确匹配的位置。详见Rosen PA,Hensley S,Joughin IR,et al.“Synthetic aperture radarinterferometry”Proceedings of the IEEE,2000.

定义6、精配准

精配准是指粗配准后对主图像和副图像进行亚像素级的配准，其中亚像素级的配准精度要达到1/8个分辨单元。本发明结合维纳滤波理论构造最佳加权矢量，实现SAR图像的精配准过程，并具有自适应图像配准功能。

定义7、配准误差

配准误差是指两幅SAR图像经过粗配准和精配准后同一像素点的偏移量大小。当配准误差优于1/8像素时，所造成的去相干很小，符合SAR干涉处理的精度需求。

定义8、相位噪声

相位噪声是指天线接收的回波数据中的噪声，一般认为是加性高斯噪声。

定义9、MUSIC算法

MUSIC(Multiple Signal Classification)算法——多重信号分类算法是一种基于相关矩阵分解的信号频率估计方法，其基本思想是直接对估计的随机过程相关矩阵进行特征分解得到信号子空间和噪声子空间，利用两者的正交性构造空间谱函数，通过谱峰搜索，估计信号频率。该算法最早于1979年由R.O.Schmid提出的。

定义10、信号子空间

根据MUSIC算法的相关定义，信号子空间是指对自相关矩阵R＝E[x(n)x^H(n)]特征分解，得到与信号有关的特征值λ₁…λ_L对应的特征向量u₁…u_L生成的子空间。其中E[]表示统计平均值，x(n)表示随机信号，H表示矩阵的共轭转置，L表示x(n)中的信号频率的个数。

定义11、噪声子空间

根据MUSIC算法的相关定义，噪声子空间是指对自相关矩阵R特征分解，得到与噪声有关的特征值λ_L+1…λ_M对应的特征向量u_L+1…u_M生成的子空间，其中L表示信号频率的个数，M表示矩阵R的维数。

定义12、空间谱函数

根据信号频率向量(j＝1,2,…,L)与噪声子空间的正交性，构造空间谱的函数P_MUSIC，这里ω_j表示信号频率，L表示信号频率的个数，T表示矩阵的转置；

$P_{\mathrm{MUSIC}}=\frac{1}{\underset{j=L+1}{\overset{M}{\mathrm{\Σ}}}{\left|{\mathrm{\α}}^H{u}_j\right|}^2}$

其中α＝[1,e^-jω,…,e^-j(M-1)ω]^T(ω∈[0,2π))表示信号频率估计向量，ω表示信号频率估计值，u_j(j＝L+1,…,M)是噪声子空间中对应的特征向量，H表示矩阵的共轭转置，||表示求模值，M表示矩阵R的维数，L表示信号频率的个数。在[0,2π)内改变ω，利用空间谱的函数P_MUSIC的表达式求得谱函数值，峰值点位置对应的ω'即为信号频率的估计值。详见文献“现代数字信号处理及其应用”，何子述等编著，清华大学出版社。

定义13、维纳滤波器

维纳滤波器是Rorbert Wiener于1949年首先根据最小均方误差(MMSE)准则提出的一种线性滤波器，具有自适应滤波功能。

定义14、最小均方误差准则

最小均方误差准则是一种认为滤波器输出与需要信号之差的均方值最小为最佳的准则。

ξ＝E[|f(n)|²]， $f\left(n\right)=d\left(n\right)-\hat{d}\left(n\right)=d\left(n\right)-{\stackrel{\→}{\mathrm{\ω}}}^H\stackrel{\→}{x}\left(n\right)$

其中ξ为最小均方误差准则的性能函数，f(n)为实际信号d(n)与滤波器输出信号的估计误差，为滤波器的输入信号矢量，为滤波器参量，E[]表示统计平均值，||表示求模值，H表示矩阵的共轭转置。

定义15、最佳滤波器参量

根据维纳滤波的相关定义，最佳滤波器参量是指维纳滤波器输出与需要信号之差的均方值最小时的滤波器参量。具体的数学公式表达见说明书附图1，详见文献“自适应滤波”，龚耀寰编著，电子工业出版社。

定义16、Hadamard积

Hadamard积是指两个相同维数的矩阵同一位置对应元素相乘。

定义17、统计平均

统计平均是指对一随机过程的所有样本求平均值，记为E[]。

定义18、样本平均

样本平均是指对随机过程的一个样本求平均值，由于无法描述平稳的随机过程的所有样本，通常情况下利用样本平均来代替统计平均。

定义19、自适应的高精度干涉SAR相位估计方法

该方法首先利用维纳滤波器获得最佳滤波器参量利用维纳滤波器获得最佳滤波器参量，再利用该参量构造联合最佳加权矢量，由该矢量获得协方差矩阵和相干系数矩阵，然后结合MUSIC算法，特征分解协方差矩阵得到信号子空间和噪声子空间，根据两者的正交性构造空间谱函数，通过谱峰搜索准确地估计主副图像的干涉相位。该方法流程图见说明书附图2。

本发明提供了一种自适应的高精度干涉SAR相位估计方法，它包含以下几个步骤：

步骤1、采用传统的最大相干系数法进行InSAR的图像粗配准：

采用传统的最大相干系数法对InSAR副图像进行图像粗配准，得到粗配准后副图像S₂。具体方法是：在主图像S₁中选定任一像素点作为参考点，以参考点为中心取一个P×P大小的窗口；在粗配准前的副图像S'₂中对应位置处选取一个Q×Q大小的搜索窗口，这里Q>>P保证在副图像中搜索到与主图像相匹配的P×P大小的窗口，P表示主图像S₁中参考窗口的大小，Q表示粗配准前的副图像S'₂中搜索窗口的大小；该方法见说明书附图3。

在配准前副图像S'₂中的搜索窗口中，选取一个P×P大小的窗口，在整个Q×Q大小的搜索窗口进行滑窗移动，计算每次滑窗后主副图像中P×P大小的窗口各像素点的复相干系数γ：

$\mathrm{\γ}=\frac{\underset{m=1}{\overset{P}{\mathrm{\Σ}}}\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{S}_1^{\′}(m,n){S_2^{\′}}^{*}(m+u,n+v)}{\sqrt{\underset{m=1}{\overset{P}{\mathrm{\Σ}}}\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{\left|S_1^{\′}(m,n)\right|}^2\underset{m=1}{\overset{P}{\mathrm{\Σ}}}\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{\left|S_2^{\′}(m+u,n+v)\right|}^2}}$

这里m，n为整数，其中S′₁(m,n)表示在主图像S₁中选取P×P大小的窗口,在窗口中位于第m行、第n列的像素值，这里m，n为整数，S'₂(m+u,n+v)表示副图像S'₂中滑窗移动u行、v列的搜索窗口的像素值，这里m，n，u，v为整数，γ表示每次滑窗后主副图像中P×P大小的窗口中各像素点的复相干系数，*表示共轭，||表示求模值；最后通过每次滑窗后主副图像中P×P大小的窗口中各像素点的复相干系数γ表达式计算参考窗口与搜索得到的各窗口的相干系数，每次滑窗后主副图像中P×P大小的窗口中各像素点的复相干系数γ取最大值所对应的位置即为粗配准的误差偏移位置(u,v)；根据得到的误差偏移位置(u,v)，对配准前的副图像S'₂中所有像素点移动u行、v列，得到配准前的副图像S'₂的粗配准图像S₂。

步骤2、建立维纳滤波器信号模型：

由步骤1得到的粗配准副图像S₂，建立维纳滤波器的信号模型。具体过程如下：取粗配准副图像S₂中对应第i个像素为中心的一个3×3的矩阵作为维纳滤波器的输入信号，建立如下滤波器模型：

构造回波信号矢量

${\stackrel{\‾}{S}}_2\left(i\right)={[S_2(i-4),S_2(i-3),...,S_2(i+4)]}^T$

其中S₂(i-4),S₂(i-3),…,S₂(i+4)表示粗配准后副图像S₂中3×3大小矩阵中的像素点，T表示矩阵的转置，见说明书附图4。

构造维纳滤波器参量ω(i)

ω(i)＝[ω₁,ω₂,…,ω₉]^H

其中ω₁,ω₂,…,ω₉分别表示维纳滤波器9个抽头的权系数，H表示矩阵的共轭转置。

将维纳滤波器参量ω(i)与回波信号矢量相乘得到s₂(i),将s₂(i)作为主图像S₁中第i个像素的估计值，即将主图像S₁中第i个像素的估计值s₂(i)与主图像S₁中对应的像素值s₁(i)进行相减，得到主副图像第i个像素点的估计误差值e(i)

$e\left(i\right)=s_2\left(i\right)-s_1\left(i\right)=\mathrm{\ω}{\left(i\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)-s_1\left(i\right)$

其中H表示矩阵的共轭转置。e(i)的均方值最小时对应着最佳滤波器参量ω_opt(i)。见说明书附图1。

步骤3、计算维纳滤波器的最佳滤波器参量ω_opt(i)：

求取粗配准后副图像回波信号矢量的自相关矩阵 求取粗配准后副图像回波信号矢量与主图像像素点s₁(i)的互相关矢量 上式中，E[]表示统计平均，H表示矩阵的共轭转置，*表示矩阵的共轭。

对粗配准后副图像回波信号矢量的自相关矩阵求逆，并与互相关矢量相乘，得到最佳滤波器参量ω_opt(i)，即其中-1表示求矩阵的逆。

步骤4、实现InSAR图像精配准

将步骤3得到的最佳滤波器参量ω_opt(i)与由步骤2构造的回波信号矢量相乘，得到副图像S₂中第i个像素点的精配准后的值s_opt2(i)

$s_{\mathrm{opt}2}\left(i\right)={\left({\mathrm{\ω}}_{\mathrm{opt}}\left(i\right)\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)$

其中ω_opt(i)表示像素点i的最佳滤波器参量，是步骤2中粗配准后副图像中对应像素i的3×3的矩阵构成的回波信号矢量，H表示矩阵的共轭转置。

步骤5、构造接收信号的协方差矩阵

利用步骤4精配准后副图像的相应像素点s_opt2(i)，构造联合最佳加权矢量s(i)，s(i)＝[s₁(i),s_opt2(i)]^T，这里T表示矩阵的转置。

利用最佳加权矢量s(i)，构造接收信号协方差矩阵C_s(i)为：

其中E[]表示统计平均，表示干涉相位导向矢量，表示主副图像中对应第i个像素点之间的干涉相位，x(i)表示加权矢量s(i)中的信号矢量，v(i)表示加权矢量s(i)中加性高斯噪声矢量，R(i)＝E[x(i)x(i)^H]表示信号矢量x(i)的自相关矩阵，σ²表示噪声矢量v(i)的方差，I₂表示二阶的单位矩阵，□表示Hadamard积，H表示矩阵的共轭转置。

步骤6、利用MUSIC算法分离信号和噪声子空间

将步骤5得到的接收信号的协方差矩阵C_s(i)，利用MUSIC算法分解得到信号和噪声子空间。具体步骤如下：对协方差矩阵C_s(i)进行特征值分解，求取它的最小特征值λ_min，与最小特征值λ_min对应的特征向量为噪声子空间记为u₂。

步骤7、构造谱函数估计干涉相位

将噪声子空间记为E_N＝u₂，利用干涉相位导向矢量与噪声子空间相乘，构造谱函数这里表示干涉相位估计值，T表示矩阵的转置，在[0,2π)内改变相位的值，P_MUSIC峰值位置对应的相位即为主副图像中对应第i个像素点的干涉相位值

本发明的原理：结合维纳滤波理论构造最佳加权矢量，通过对其组成的最佳协方差矩阵特征分解得到信号子空间和噪声子空间，充分利用相应像素对及其相邻像素的相干信息，并根据MUSIC算法中信号子空间和噪声子空间的正交性构造空间谱函数，通过谱峰搜索准确地估计出主副图像中相应像素间的干涉相位。

本发明的创新点在于结合维纳滤波理论构造最佳加权矢量，利用最佳加权矢量构造干涉信号模型，通过对其组成的最佳协方差矩阵特征分解得到信号子空间和噪声子空间，并根据MUSIC算法中信号子空间和噪声子空间的正交性来构造谱函数，准确的估计出干涉相位，并具有很强的相位噪声抑制功能，。本发明利用维纳滤波实现自适应图像精配准，不需要首先确定配准误差方向和大小，并且避免了最佳滤波器参量的搜索过程，从而很好地解决了InSAR干涉相位估计运算量大的问题。

本发明的优点在于与传统算法相比，本文算法整个流程简单，较小的运算量精确的估计出干涉相位。与传统InSAR干涉相位估计相比，本发明特别适用于InSAR复杂场景地表参数精确反演等领域。

附图说明：

图1是本发明维纳滤波器模型。

其中s₁(i)为主图像中像素i，S₂(i)是由粗配准后副图像中对应像素i的3×3的矩阵所构成的回波信号矢量，ω(i)为对应的维纳滤波器参量，s₂(i)为对主图像中像素i的估计，e(i)为主副图像对应像素点i的估计误差，表示乘积，表示相加，T表示矩阵的转置，H表示矩阵的共轭转置。

图2是本发明所提供方法的流程框图。

图3是最大相干系数法粗配准示意图。

其中P表示主图像中参考窗口的大小，Q表示副图像中搜索窗口的大小。

图4是粗配准后主副图像对应像素点选取图。

其中i表示主副图像对应位置的像素点，i-4,…,i-1,i+1,…,i+4表示粗配准后副图像中3×3大小矩阵中的像素点。

具体实施方式

本发明主要采用实测数据进行验证，所有步骤、结论都在VC++、MATLAB7.0上验证正确。具体实施步骤如下：

步骤1、利用相干系数法进行InSAR的图像粗配准：

读取主副复图像，主图像记为S₁；副图像记为S'₂。利用最大相干系数法对副图像S'₂进行粗配准，得到副图像的误差偏移位置为(30,18)，对副图像所有像素点向下移动30行、向右移动18列，得到副图像的粗配准图像S₂。配准后主副图像大小为417×430。

步骤2、建立维纳滤波器信号模型：

取主图像S₁的像素点i作为滤波器的需要信号，粗配准后副图像S₂中对应像素i的一个3×3的矩阵作为维纳滤波器的输入信号，建立如下滤波器模型：

${\stackrel{\‾}{S}}_2\left(i\right)={[s_2(i-4),s_2(i-3),...,s_2(i+4)]}^T$

ω(i)＝[ω₁,ω₂,…，ω₉]^H

$s_2\left(i\right)=\mathrm{\ω}{\left(i\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)$

$e\left(i\right)=s_2\left(i\right)-s_1\left(i\right)={\mathrm{\ω}\left(i\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)-s_1\left(i\right)$

具体实施方法：主图像为一个417×430大小的矩阵，取第m行，第n列元素s₁(m,n)，对应的粗配准后副图像也取第m行，第n列元素s₂(m-1,n)、s₂(m-1,n+1)、s₂(m,n-1)、s₂(m,n)、s₂(m,n+1)、s₂(m+1,n-1)、s₂(m+1,n)和s₂(m+1,n+1)作为滤波器的输入信号。这里主图像第4行，第4列元素为例：s₁(4,4)＝-0.0105-0.0740i， ${\stackrel{\‾}{S}}_2\left(\mathrm{4,4}\right)=[-0.0774+0.0025i,-0.0844-0.0701i,$ $-0.0364-0.0405i,-0.0931-0.0078i,-0.0713-0.0554i,-0.0859-0.0157i,-0.0340+0.0079i,$ ${-0.0070-0.0277i,-0.0289-0.0126i]}^T,$ 依次取m＝5,…,414，n＝5,…,427。

步骤3、计算维纳滤波器的最佳滤波器参量ω_opt(i)：

根据前面步骤3的推导，维纳滤波器中最佳维纳滤波器参量ω_opt(i)的求解方法如下所示：

$C_{S_2}\left(i\right){\mathrm{\ω}}_{\mathrm{opt}}=r_{S_{2}d}\⇒{\mathrm{\ω}}_{\mathrm{opt}}={C_{S_2}}^{-1}\left(i\right)r_{S_{2}d}$

其中 $C_{S_2}\left(i\right)=E\left[{\stackrel{\‾}{S}}_2\left(i\right){{\stackrel{\‾}{S}}_2}^H\left(i\right)\right]$ 表示输入矢量的自相关矩阵， $r_{S_{2}d}=E\left[{\stackrel{\‾}{S}}_2\left(i\right)s_1^{*}\left(i\right)\right]$

表示输入矢量与需要信号s₁(i)的互相关矢量，E[]表示统计平均。

在实际中一般用样本平均来估计统计平均，从相邻的像素中获得独立同分布的样本，对应的和由下式进行估计：

$C_{S_2}\left(i\right)=\frac{1}{2l+1}\underset{j=-l}{\overset{l}{\mathrm{\Σ}}}{\stackrel{\‾}{S}}_2(i+j){\stackrel{\‾}{S}}_2{(i+j)}^H$

$r_{S_{2}d}=\frac{1}{2l+1}\underset{j=-l}{\overset{l}{\mathrm{\Σ}}}{\stackrel{\‾}{S}}_2(i+j)s_1{(i+j)}^{*}$

这里选取像素i周围5×5矩阵大小的像素点作为样本点，对应的s₁(i+j)＝s₁(m+kk,n+ll)，是由以s₂(m+kk,n+ll)为中心点的一个3×3的矩阵构成。由像素点s₁(4,4)和求得的样本互相关矩阵 $r_{S_{2}d}=[0.0007-0.0006i,$ $0.0010+0.0005i,0.0004+0.0007i,0.0009-0.0010i,0.0015-0.0002i,0.0004-0.0005i,0.0004$ $-0.0006i,0.0018-0.0001i,0.0005-0.0010i{]}^T,$ 对应的样本协方差矩阵 最后得到ω(4,4)＝[0.3964，0.3370，0.2931，0.3732，0.8134，0.3237，0.5178，0.8472，0.5293]^T。

步骤4、实现InSAR图像精配准

SAR图像每一像素i经过步骤3处理后，可以得到每个像素点i的最佳滤波器参量ω_opt(i)。然后利用最佳滤波器参量，得到精配准后副图像的相应像素点s₂(i)

$s_2\left(i\right)={\left({\mathrm{\ω}}_{\mathrm{opt}}\left(i\right)\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)$

根据步骤3得到的最佳滤波器参量ω(4,4)，，计算出精配准后相应像素点s₂(4,4)＝-0.0190+0.0554i。

步骤5、构造接收信号的协方差矩阵

由最佳加权矢量估计最佳协方差矩阵C_s(i)，利用样本平均来估计统计平均，即

$C_s\left(i\right)=\frac{1}{2M+1}\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}s(i+j)s{(i+j)}^H$

这里选取像素i周围5×5矩阵大小的像素点作为样本点，得到的C_s(4,4)＝[0.0030，0.0044+0.0015i；0.0044-0.0015i，0.0319]。

步骤6、利用MUSIC算法分离信号和噪声子空间

根据MUSIC算法中的相关定义，C_s(4,4)经过特征分解，得到的两个特征值为0.0023和0.0115，对应的特征向量分别为[0.8623-0.1648i，-0.4789]^T和[0.4704-0.0899i，0.8779]^T。小特征值对应的特征向量为噪声子空间，即E_N＝[0.8623-0.1648i,-0.4789]^T。

步骤7、构造谱函数估计干涉相位

构造的谱函数为：

在[0,2π)内改变使达到峰值对应的相位即为干涉相位值。将步骤7中E_N代入谱函数，主副图像第4行，第4列像素的干涉相位为5.8905。其他像素点干涉相位依照每行每列递增来获得。

通过本发明具体实施方式可以看出，本发明所提供的一种自适应的高精度干涉SAR相位估计方法能够准确的估计出干涉相位，具有自适应图像配准功能和很强的噪声抑制能力，且避免了最佳滤波器参量的搜索过程，大大降低了运算量。

Claims (1)

1.一种自适应的高精度干涉SAR相位估计方法，其特征是它包含以下几个步骤：

步骤1、采用传统的最大相干系数法进行InSAR的图像粗配准：

采用传统的最大相干系数法对InSAR副图像进行图像粗配准，得到粗配准后副图像S₂；具体方法是：在主图像S₁中选定任一像素点作为参考点，以参考点为中心取一个P×P大小的窗口；在配准前的副图像S'₂中对应位置处选取一个Q×Q大小的搜索窗口，这里Q>>P保证在副图像中搜索到与主图像相匹配的P×P大小的窗口，P表示主图像S₁中参考窗口的大小，Q表示配准前的副图像S'₂中搜索窗口的大小；

在配准前副图像S'₂中选定任一像素点作为参考点，以参考点为中心取一个P×P大小的窗口，在整个Q×Q大小的搜索窗口进行滑窗移动，计算每次滑窗后主副图像中P×P大小的窗口各像素点的复相干系数γ：

$\mathrm{\γ}=\frac{\underset{m=1}{\overset{P}{\mathrm{\Σ}}}\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{S}_1^{\′}(m,n){S_2^{\′}}^{*}(m+u,n+v)}{\sqrt{\underset{m=1}{\overset{P}{\mathrm{\Σ}}}\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{\left|S_1^{\′}(m,n)\right|}^2\underset{m=1}{\overset{P}{\mathrm{\Σ}}}\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{\left|S_2^{\′}(m+u,n+v)\right|}^2}}$

这里m，n为整数，其中S′₁(m,n)表示在主图像S₁中选取P×P大小的窗口,在窗口中位于第m行、第n列的像素值，这里m，n为整数，S'₂(m+u,n+v)表示副图像S'₂中滑窗移动u行、v列的搜索窗口的像素值，这里m，n，u，v为整数，γ表示每次滑窗后主副图像中P×P大小的窗口中各像素点的复相干系数，*表示共轭，||表示求模值；最后通过每次滑窗后主副图像中P×P大小的窗口中各像素点的复相干系数γ表达式计算参考窗口与搜索得到的各窗口的相干系数，每次滑窗后主副图像中P×P大小的窗口中各像素点的复相干系数γ取最大值所对应的位置即为粗配准的误差偏移位置(u,v)；根据得到的误差偏移位置(u,v)，对配准前的副图像S'₂中所有像素点移动u行、v列，得到配准前的副图像S'₂的粗配准图像S₂；

步骤2、建立维纳滤波器信号模型：

由步骤1得到的粗配准副图像S₂，建立维纳滤波器的信号模型；具体过程如下：取粗配准副图像S₂中对应第i个像素为中心的一个3×3的矩阵作为维纳滤波器的输入信号，建立如下滤波器模型：

构造回波信号矢量

${\stackrel{\‾}{S}}_2\left(i\right)={[S_2(i-4),S_2(i-3),...,S_2(i+4)]}^T$

其中S₂(i-4),S₂(i-3),…,S₂(i+4)表示粗配准图像S₂中3×3大小矩阵中的像素点，T表示矩阵的转置；

构造维纳滤波器参量ω(i)

ω(i)＝[ω₁,ω₂,…,ω₉]^H

其中ω₁,ω₂,…,ω₉分别表示维纳滤波器9个抽头的权系数，H表示矩阵的共轭转置；

将维纳滤波器参量ω(i)与回波信号矢量相乘得到s₂(i),将s₂(i)作为主图像S₁中第i个像素的估计值，即将主图像S₁中第i个像素的估计值s₂(i)与主图像S₁中对应的像素值s₁(i)进行相减，得到主副图像第i个像素点的估计误差值e(i)

$e\left(i\right)=s_2\left(i\right)-s_1\left(i\right)=\mathrm{\ω}{\left(i\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)-s_1\left(i\right)$

其中H表示矩阵的共轭转置；e(i)的均方值最小时对应着最佳滤波器参量ω_opt(i)；

步骤3、计算维纳滤波器的最佳滤波器参量ω_opt(i)：

求取粗配准后副图像回波信号矢量的自相关矩阵 求取粗配准后副图像回波信号矢量与主图像像素点s₁(i)的互相关矢量 上式中，E[]表示统计平均，H表示矩阵的共轭转置，*表示矩阵的共轭；

对粗配准后副图像回波信号矢量的自相关矩阵求逆，并与互相关矢量相乘，得到最佳滤波器参量ω_opt(i)，即其中-1表示求矩阵的逆；

步骤4、实现InSAR图像精配准

将步骤3得到的最佳滤波器参量ω_opt(i)与由步骤2构造的回波信号矢量相乘，得到副图像S₂中第i个像素点的精配准后的值s_opt2(i)

$s_{\mathrm{opt}2}\left(i\right)={\left({\mathrm{\ω}}_{\mathrm{opt}}\left(i\right)\right)}^H{\stackrel{\‾}{S}}_2\left(i\right)$

其中ω_opt(i)表示第i个的像素点最佳滤波器参量，是步骤2中副图像S₂中的第i个像素的3×3的矩阵构成的回波信号矢量，H表示矩阵的共轭转置；

步骤5、构造接收信号的协方差矩阵

利用步骤4精配准后副图像的相应像素点s_opt2(i)，构造联合最佳加权矢量s(i)，s(i)＝[s₁(i),s_opt2(i)]^T，这里T表示矩阵的转置；

利用最佳加权矢量s(i)，构造接收信号协方差矩阵C_s(i)为：

其中E[]表示统计平均，表示干涉相位导向矢量，表示主副图像中相对应的第i个像素点之间的干涉相位，x(i)表示加权矢量s(i)中的信号矢量，v(i)表示加权矢量s(i)中加性高斯噪声矢量，R(i)＝E[x(i)x(i)^H]表示信号矢量x(i)的自相关矩阵，σ²表示噪声矢量v(i)的方差，I₂表示二阶的单位矩阵，□表示Hadamard积，H表示矩阵的共轭转置；

步骤6、利用MUSIC算法分离信号和噪声子空间

将步骤5得到的接收信号的协方差矩阵C_s(i)，利用MUSIC算法分解得到信号和噪声子空间；具体步骤如下：对协方差矩阵C_s(i)进行特征值分解，求取它的最小特征值λ_min，与最小特征值λ_min对应的特征向量为噪声子空间记为u₂；

步骤7、构造谱函数估计干涉相位

将噪声子空间记为E_N＝u₂，利用干涉相位导向矢量与噪声子空间相乘，构造谱函数这里表示干涉相位估计值，T表示矩阵的转置，在[0,2π)内改变相位的值，P_MUSIC峰值位置对应的相位即为主副图像中第i个像素点的干涉相位值