CN102565798A - 基于相关加权联合子空间投影的InSAR干涉相位估计方法 - Google Patents

Abstract

一种基于相关加权联合子空间投影的InSAR干涉相位估计方法。其包括用传统的相关法以主图像作为参考对副图像进行粗配准，得到粗配准图像；利用主图像和粗配准图像构造相关加权观测矢量；计算相关加权联合协方差矩阵和加权联合相关函数矩阵；计算噪声子空间和信号子空间；用噪声子空间投影来构造代价函数，并通过代价函数估计干涉相位。本发明提供的基于相关加权联合子空间投影的InSAR干涉相位估计方法可用于在合成孔径雷达SAR图像配准精度很差的条件下对其相应像素间的干涉相位进行准确地估计。

Description

基于相关加权联合子空间投影的InSAR干涉相位估计方法

技术领域

本发明属于雷达探测技术领域，特别是涉及一种基于相关加权联合子空间投影的InSAR干涉相位估计方法，其可以在合成孔径雷达SAR图像配准精度很差的条件下准确地估计出相应像素间的干涉相位。

背景技术

干涉合成孔径雷达InSAR是一种获取地面数字高程图的重要遥感技术。InSAR技术中的三大关键处理步骤是SAR图像配准、干涉相位估计即干涉相位噪声滤波和相位展开。若图像配准的精度较差，将使后面的相位展开难以获得令人满意的结果。因为传统的相位展开方法，例如枝切法、区域增长法和最小二乘LS算法的性能受图像配准精度的影响很大，都要求图像配准的精度达到1/10～1/100个分辨单元，否则将会严重影响相位展开的性能。然而，对于相干性较低和不同几何形变的两幅图像，精确配准会存在困难。因此，研究对图像配准误差具有强稳健的干涉相位估计方法具有重要的实用价值。

目前，对图像配准误差稳健的干涉相位估计方法主要有以下几种：

1.李真芳等人在2006年的IEEE Trans.On GRS上发表的文章《ImageAuto-Coregistration and InSAR Interferogram Estimation Using Joint Subspace Projection》中，采用了相邻像素的相干信息和空间投影技术，虽然该方法可以在存在配准误差时获得满意的干涉相位估计结果，但是在估计干涉相位时首先要确定噪声子空间维数，若噪声子空间维数估计不准，必定影响干涉相位的估计结果。

2.李海等人在2007年电子学报第三期上发表的文章《InSAR自适应图像配准的干涉相位估计方法》中，提出了基于加权联合单像素模型的干涉相位估计方法，虽然在估计干涉相位时不需要确定噪声子空间维数，但是在估计干涉相位时需要首先确定配准误差方向，并且确定最优权值时需要搜索，因此计算量较大。

发明内容

为了解决上述问题，本发明的目的在于提供一种基于相关加权联合子空间投影的InSAR干涉相位估计方法，以解决在SAR图像配准精度很差条件下无法准确估计出相应像素间的干涉相位问题。

为了达到上述目的，本发明提供的基于相关加权联合子空间投影的InSAR干涉相位估计方法包括按顺序进行的下列步骤：

1)用传统的相关法以主图像作为参考对副图像进行粗配准，得到粗配准图像；

2)利用主图像和粗配准图像构造相关加权观测矢量；

3)计算相关加权联合协方差矩阵和加权联合相关函数矩阵；

4)计算噪声子空间和信号子空间；

5)用噪声子空间投影来构造代价函数，并通过代价函数估计干涉相位。

在步骤2)中，所述的构造相关加权观测矢量的方法是：

a.分别计算主图像中像素a与粗配准图像中像素m的相关系数：

$r_{31}(m,a)=\frac{\left|\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{s}_3(m+k)s_1^{*}(a+k)\right|}{\sqrt{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{\left|s_3(m+k)\right|}^2\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{\left|s_1(i+k)\right|}^2}}\left(\begin{array}{c}m=a-5,a-4,a-3,a-1,a,a+1,a+3,a+4,a+5\\ a=i-1,i,i+3,i+4\end{array}\right)$

b.用所得到的相关系数r₃₁(m，a)分别对粗配准图像中的像素m进行加权，得到元素s_3W(a，w)，公式如下：

$s_{3W}(a,w)=w_a^T{V}_a$

其中：

V_a＝[s₂(a-5)，s₂(a-4)，s₂(a-3)，s₂(a-1)，s₂(a)，

                    s₂(a+1)，s₂(a+3)，s₂(a+4)，s₂(a+5)]^T

                                        (a＝i-1，i，i+3，i+4)

w_a＝[r₃₁(a-5，a)，r₃₁(a-4，a)，r₃₁(a-3，a)，r₃₁(a-1，a)，

        r₃₁(a，a)，r₃₁(a+1，a)，r₃₁(a+3，a)，r₃₁(a+4，a)，r₃₁(a+5，a)]^T

                                                (a＝i-1，i，i+3，i+4)

c.用主图像中的像素a(a＝i-1，i，i+3，i+4)和得到的元素s_3W(a，w)构造相关加权观测矢量：

$\mathrm{js}(i,w)={[s_1(i-1),s_{2W}(i-1),s_1\left(i\right),s_{2W}\left(i\right),s_1(i+3),s_{2W}(i+3),s_1(i+4),s_{2W}(i+4)]}^T$

$={[s_1(i-1),w_{i-1}^T{V}_{i-1},s_1\left(i\right),w_i^T{V}_i,s_1(i+3),w_{i+3}^T{V}_{i+3},s_1(i+4),w_{i+4}^T{V}_{i+4}]}^T$

其中：s₁和s₃分别为主图像和粗配准图像所对应的复图像数据，

i表示当前要计算的主图像和粗配准图像的像素标号，

2K+1为从相邻像素中获得的独立同分布样本数，

上标T表示转置操作，上标*表示共轭操作。

在步骤3)中，所述的相关加权联合协方差矩阵用下式得出：

C_js(i，w)＝E{js(i，w)js^H(i，w)}

所述的加权联合相关函数矩阵用下式得出：

$R_{\mathrm{js}}(i,w)=|C_{\mathrm{js}}(i,w)-{\mathrm{\σ}}_n^{2}I|$

式中，js(i，w)为相关加权观测矢量，E表示统计平均，上标H表示共轭操作。

在步骤4)中，所述的噪声子空间是由相关加权联合协方差矩阵C_js(i，w)(维数为8×8)的小特征值所对应的特征矢量构成；信号子空间是由对加权联合相关函数矩阵R_js(i，w)进行特征值分解而获得的大特征值矢量

(k＝1，2，3，4)和导向矢量的Hadamard积所构成。

在步骤5)中，所述的干涉相位估计方法为：

用噪声子空间投影来构造代价函数：

式中

为导向矢量，φ_i的取值范围是[-π，π]，上面的代价函数达到最小值时所对应的干涉相位

即为干涉相位估计结果。

本发明提供的基于相关加权联合子空间投影的InSAR干涉相位估计方法与现有技术相比具有以下优点：

(1)采用相关加权联合观测矢量，使其所对应的协方差矩阵的噪声子空间维数不受配准误差的影响，即噪声子空间维数在任意配准误差(不超过一个像素)的情况下都和精确配准时的噪声子空间维数一样，因此不需要在估计干涉相位时首先确定噪声子空间维数，避免了由于噪声子空间维数估计不准影响干涉相位估计结果的问题；

(2)直接用相关系数对观测矢量进行加权，不需要搜索，大大降低运算量；

(3)同时采用联合子空间投影技术，因而可以在SAR图像配准精度较差的情况下得到准确的干涉相位估计结果。

附图说明

图1为本发明提供的基于相关加权联合子空间投影的InSAR干涉相位估计方法流程图。

图2为配准误差方向未知时本发明构造相关加权观测矢量时主图像和粗配准图像示意图。

图3(a)为精确配准时现有均值滤波方法的处理结果。

图3(b)为精确配准时现有中值滤波方法的处理结果。

图3(c)为精确配准时现有Lee自适应滤波方法的处理结果。

图3(d)为精确配准时现有自适应contoured滤波方法的处理结果。

图3(e)为精确配准时本发明方法的处理结果。

图4(a)为配准误差为0.5像素时现有均值滤波方法的处理结果。

图4(b)为配准误差为0.5像素时现有中值滤波方法的处理结果。

图4(c)为配准误差为0.5像素时现有Lee自适应滤波方法的处理结果。

图4(d)为配准误差为0.5像素时现有自适应contoured滤波方法的处理结果。

图4(e)为配准误差为0.5像素时本发明方法的处理结果。

图5(a)为配准误差为1像素时现有均值滤波方法的处理结果。

图5(b)为配准误差为1像素时现有中值滤波方法的处理结果。

图5(c)为配准误差为1像素时现有Lee自适应滤波方法的处理结果。

图5(d)为配准误差为1像素时现有自适应contoured滤波方法的处理结果。

图5(e)为配准误差为1像素时本发明方法的处理结果。

图6(a)为对两幅实测SAR图像经过粗配准后生成的干涉相位图。

图6(b)为采用本发明方法对两幅实测SAR图像进行处理得到的干涉相位图。

具体实施方式

下面结合附图和具体实施例对本发明提供的基于相关加权联合子空间投影的InSAR干涉相位估计方法进行详细说明。

如图1所示，本发明提供的基于相关加权联合子空间投影的InSAR干涉相位估计方法包括按顺序进行的下列步骤：

1)用传统的相关法以主图像作为参考对副图像进行粗配准，得到粗配准图像；

在实测得到的主、副两幅SAR图像中，以主图像作为参考，利用传统的相关法对副图像进行图像粗配准处理，得到粗配准图像，在该图像粗配准处理中，不要求像传统的干涉相位估计方法那样要求图像配准的精度必须达到亚像素级，即配准精度要达到1/10到1/100像素，只要求图像配准精度达到像素级，即配准精度允许达到一个分辨单元就够了，因此大大减轻了图像配准的难度。

2)利用主图像和粗配准图像构造相关加权观测矢量；

利用主图像和粗配准图像构造相关加权观测矢量js(i，w)，主图像和粗配准图像如图2所示，具体构造过程如下：

a.分别计算主图像中像素a与粗配准图像中像素m的相关系数：

$r_{31}(m,a)=\frac{\left|\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{s}_3(m+k)s_1^{*}(a+k)\right|}{\sqrt{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{\left|s_3(m+k)\right|}^2\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{\left|s_1(i+k)\right|}^2}}\left(\begin{array}{c}m=a-5,a-4,a-3,a-1,a,a+1,a+3,a+4,a+5\\ a=i-1,i,i+3,i+4\end{array}\right)---\left(3\right)$

b.用所得到的相关系数r₃₁(m，a)分别对粗配准图像中的像素m进行加权，得到元素s_3W(a，w)，公式如下：

$s_{3W}(a,w)=w_a^T{V}_a---\left(2\right)$

其中：

V_a＝[s₂(a-5)，s₂(a-4)，s₂(a-3)，s₂(a-1)，s₂(a)，

        s₂(a+1)，s₂(a+3)，s₂(a+4)，s₂(a+5)]^T

                        (a＝i-1，i，i+3，i+4)  (3)

w_a＝[r₃₁(a-5，a)，r₃₁(a-4，a)，r₃₁(a-3，a)，r₃₁(a-1，a)，

        r₃₁(a，a)，r₃₁(a+1，a)，r₃₁(a+3，a)，r₃₁(a+4，a)，r₃₁(a+5，a)]^T

                                    (a＝i-1，i，i+3，i+4)  (4)

c.用主图像中的像素a(a＝i-1，i，i+3，i+4)和得到的元素s_3W(a，w)构造相关加权观测矢量：

$\mathrm{js}(i,w)={[s_1(i-1),s_{2W}(i-1),s_1\left(i\right),s_{2W}\left(i\right),s_1(i+3),s_{2W}(i+3),s_1(i+4),s_{2W}(i+4)]}^T$ (5)

$={[s_1(i-1),w_{i-1}^T{V}_{i-1},s_1\left(i\right),w_i^T{V}_i,s_1(i+3),w_{i+3}^T{V}_{i+3},s_1(i+4),w_{i+4}^T{V}_{i+4}]}^T$

其中：s₁和s₃分别为主图像和粗配准图像所对应的复图像数据，

i表示当前要计算的主图像和粗配准图像的像素标号，

2K+1为从相邻像素中获得的独立同分布样本数，

上标T表示转置操作，上标*表示共轭操作；

3)计算相关加权联合协方差矩阵和加权联合相关函数矩阵；

用相关加权观测矢量js(i，w)估计相应的相关加权联合协方差矩阵：

C_js(i，w)＝E{js(i，w)js^H(i，w)}                                (6)

加权联合相关函数矩阵用下式进行估计：

$R_{\mathrm{js}}(i,w)=|C_{\mathrm{js}}(i,w)-{\mathrm{\σ}}_n^{2}I|---\left(7\right)$

式中，E表示统计平均，上标H表示共轭操作；

4)计算噪声子空间和信号子空间；

相关加权联合协方差矩阵C_js(i，w)(维数为8×8)的特征分解有如下特点：特征值λ₁＞λ₂＞λ₃＞λ₄＞＞λ₅＞λ₆＞λ₇＞λ₈，(l＝1，…，4)为小特征值λ_a(a＝5，6，…，8)所对应的特征矢量. $N_c=[{\mathrm{\β}}_{\mathrm{njs}}^{\left(1\right)},{\mathrm{\β}}_{\mathrm{njs}}^{\left(2\right)},{\mathrm{\β}}_{\mathrm{njs}}^{\left(3\right)},{\mathrm{\β}}_{\mathrm{njs}}^{\left(4\right)}]$ 构成C_js(i，w)的噪声子空间；

对R_js(i，w)进行特征值分解来获得大特征值矢量

(k＝1，2，3，4)，

(k＝1，2，3，4)位于C_js(i，w)的信号子空间；

5)用噪声子空间投影来构造代价函数，并通过代价函数估计干涉相位；

用噪声子空间投影来构造代价函数：

式中：

为导向矢量，φ_i的取值范围是[-π，π]，上面的代价函数达到最小值时所对应的干涉相位

即为干涉相位估计结果；

对主图像和粗配准图像中的每一像素分别进行上述五步操作，就可以得到整个地形的干涉相位图；

本发明的效果可以通过以下仿真结果进一步说明。

仿真数据描述：仿真假定3颗卫星按Cartwheel构形进行编队飞行，我们选择其中某一轨道时刻的两颗卫星位置进行仿真，对应的有效垂直基线为562.93米，卫星高度750公里，入射角50°，利用一幅实测SAR图像作为地面场景背景来产生每一个SAR像素的回波功率，并利用两维窗来模拟高程地形，生成一对InSAR图像。SAR图像中的信噪比(SNR)为16dB。相干系数由垂直基线长度、局部地形坡度和SNR决定。2颗卫星接收的SAR图像数据按统计模型产生。

图3为在精确配准情况下采用不同方法对干涉相位图的处理结果：

其中图3(a)为精确配准时现有均值滤波方法的处理结果。

图3(b)为精确配准时现有中值滤波方法的处理结果。

图3(c)为精确配准时现有Lee自适应滤波方法的处理结果。

图3(d)为精确配准时现有自适应contoured滤波方法的处理结果。

图3(e)为精确配准时本发明方法的处理结果。

由图中可见，精确配准时上述不同滤波方法所得到的干涉相位图的条纹边缘很清晰。

图4为在配准误差为0.5像素情况下采用不同方法对干涉相位图的处理结果：

其中图4(a)为配准误差为0.5像素时现有均值滤波方法的处理结果。

图4(b)为配准误差为0.5像素时现有中值滤波方法的处理结果。

图4(c)为配准误差为0.5像素时现有Lee自适应滤波方法的处理结果。

图4(d)为配准误差为0.5像素时现有自适应contoured滤波方法的处理结果。

图4(e)为配准误差为0.5像素时本发明方法的处理结果。

由图中可见，配准误差为0.5像素时现有方法所得到的干涉相位图的条纹边缘很不清晰，易知其估计的干涉相位结果很差，但本发明方法所得到的干涉相位图的条纹边缘很清晰。

图5为在配准误差为1像素情况下采用不同方法对干涉相位图的处理结果：

其中图5(a)为配准误差为1像素时现有均值滤波方法的处理结果。

图5(b)为配准误差为1像素时现有中值滤波方法的处理结果。

图5(c)为配准误差为1像素时现有Lee自适应滤波方法的处理结果。

图5(d)为配准误差为1像素时现有自适应contoured滤波方法的处理结果。

图5(e)为配准误差为1像素时本发明方法的处理结果。

由图中可见，配准误差为1像素时现有方法已经无法得到干涉条纹，也就得不到所需要的干涉相位，但本发明方法所得到的干涉相位图的条纹边缘很清晰。

由图3-图5的仿真结果可以看出：现有的均值滤波、中值滤波、Lee自适应滤波和自适应contoured滤波方法受图像配准误差的影响很大，配准误差越大，滤波效果越差，当配准误差达到1个像素时，已经无法得到干涉条纹。相反本发明方法几乎不受配准误差的影响，在SAR图像配准精度较差的情况下仍能获得和精确配准一样的处理结果。

下面利用实测数据来验证本发明方法的处理结果。所用的实测数据为欧洲空间局ESA的地球资源卫星ERS-1和ERS-2录取的河北张义地区的数据。

图6(a)为对两幅实测SAR图像经过粗配准后生成的干涉相位图。

图6(b)为采用本发明方法对两幅实测SAR图像进行处理得到的干涉相位图。

由实测数据的处理结果可以看出，SAR图像在没有经过精确配准处理时，能够获得满意的干涉相位估计结果。

Claims (5)

1.一种基于相关加权联合子空间投影的InSAR干涉相位估计方法，其特征在于：所述的InSAR干涉相位估计方法包括按顺序进行的下列步骤：

1)用传统的相关法以主图像作为参考对副图像进行粗配准，得到粗配准图像；

2)利用主图像和粗配准图像构造相关加权观测矢量；

3)计算相关加权联合协方差矩阵和加权联合相关函数矩阵；

4)计算噪声子空间和信号子空间；

5)用噪声子空间投影来构造代价函数，并通过代价函数估计干涉相位。

2.根据权利要求1所述的基于相关加权联合子空间投影的InSAR干涉相位估计方法，其特征在于：在步骤2)中，所述的构造相关加权观测矢量的方法是：

a.分别计算主图像中像素a与粗配准图像中像素m的相关系数：

$r_{31}(m,a)=\frac{\left|\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{s}_3(m+k)s_1^{*}(a+k)\right|}{\sqrt{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{\left|s_3(m+k)\right|}^2\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{\left|s_1(i+k)\right|}^2}}\left(\begin{array}{c}m=a-5,a-4,a-3,a-1,a,a+1,a+3,a+4,a+5\\ a=i-1,i,i+3,i+4\end{array}\right)$

b.用所得到的相关系数r₃₁(m，a)分别对粗配准图像中的像素m进行加权，得到元素s_3W(a，w)，公式如下：

$s_{3W}(a,w)=w_a^T{V}_a$

其中：

V_a＝[s₂(a-5)，s₂(a-4)，s₂(a-3)，s₂(a-1)，s₂(a)，

        s₂(a+1)，s₂(a+3)，s₂(a+4)，s₂(a+5)]^T

                          (a＝i-1，i，i+3，i+4)

w_a＝[r₃₁(a-5，a)，r₃₁(a-4，a)，r₃₁(a-3，a)，r₃₁(a-1，a)，

       r₃₁(a，a)，r₃₁(a+1，a)，r₃₁(a+3，a)，r₃₁(a+4，a)，r₃₁(a+5，a)]^T

                          (a＝i-1，i，i+3，i+4)

c.用主图像中的像素a(a＝i-1，i，i+3，i+4)和得到的元素s_3W(a，w)构造相关加权观测矢量：

$\mathrm{js}(i,w)={[s_1(i-1),s_{2W}(i-1),s_1\left(i\right),s_{2W}\left(i\right),s_1(i+3),s_{2W}(i+3),s_1(i+4),s_{2W}(i+4)]}^T$

$={[s_1(i-1),w_{i-1}^T{V}_{i-1},s_1\left(i\right),w_i^T{V}_i,s_1(i+3),w_{i+3}^T{V}_{i+3},s_1(i+4),w_{i+4}^T{V}_{i+4}]}^T$

其中：s₁和s₃分别为主图像和粗配准图像所对应的复图像数据，

i表示当前要计算的主图像和粗配准图像的像素标号，

2K+1为从相邻像素中获得的独立同分布样本数，

上标T表示转置操作，上标*表示共轭操作。

3.根据权利要求1所述的基于相关加权联合子空间投影的InSAR干涉相位估计方法，其特征在于：在步骤3)中，所述的相关加权联合协方差矩阵用下式得出：

C_js(i，w)＝E{js(i，w)js^H(i，w)}

所述的加权联合相关函数矩阵用下式得出：

$R_{\mathrm{js}}(i,w)=|C_{\mathrm{js}}(i,w)-{\mathrm{\σ}}_n^{2}I|$

式中，js(i，w)为相关加权观测矢量，E表示统计平均，上标H表示共轭操作。

4.根据权利要求1所述的基于相关加权联合子空间投影的InSAR干涉相位估计方法，其特征在于：在步骤4)中，所述的噪声子空间是由相关加权联合协方差矩阵C_js(i，w)(维数为8×8)的小特征值所对应的特征矢量构成；信号子空间是由对加权联合相关函数矩阵R_js(i，w)进行特征值分解而获得的大特征值矢量

(k＝1，2，3，4)和导向矢量的Hadamard积所构成。

5.根据权利要求1所述的基于相关加权联合子空间投影的InSAR干涉相位估计方法，其特征在于：在步骤5)中，所述的干涉相位估计方法为：

用噪声子空间投影来构造代价函数：

式中

为导向矢量，φ_i的取值范围是[-π，π]，上面的代价函数达到最小值时所对应的干涉相位

即为干涉相位估计结果。

Images