CN100545676C - 基于相关加权的干涉合成孔径雷达干涉相位估计方法 - Google Patents

Abstract

本发明公开了一种基于相关加权的干涉合成孔径雷达干涉相位估计方法。该方法的实现过程为：用传统的相关法以主图像1为参考对副图像2进行粗配准，得到粗配准图像3；利用主图像1和粗配准图像3，构造相关加权观测矢量si(i，w(i))；根据相关加权观测矢量估计协方差矩阵C_si(i，w(i))；由协方差矩阵得到代价函数J_BF；通过代价函数估计干涉相位，即用导向矢量和协方差矩阵进行波束扫描，并将波束形成的输出功率最大值所对应的干涉相位作为干涉相位估计结果。本发明具有运算量小和估计干涉相位准确之优点，可用于在合成孔径雷达SAR图像配准精度很差的条件下对其相应像素间的干涉相位进行准确地估计。

Description

基于相关加权的干涉合成孔径雷达干涉相位估计方法

技术领域

本发明涉及雷达探测技术领域，具体的说是一种基于相关加权联合单像素模型的干涉合成孔径雷达干涉相位估计方法，可以在合成孔径雷达SAR图像配准精度很差的条件下准确地估计相应像素间的干涉相位。

背景技术：

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

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

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

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

发明的内容

本发明的目的在于克服上述已有技术的不足，提供一种计算量小的基于相关加权的干涉合成孔径雷达干涉相位估计方法，以解决在SAR图像配准精度很差条件下可准确估计相应像素间的干涉相位问题。

实现本发明目的方法如下：

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

图像3；

B.利用主图像1和粗配准图像3，构造相关加权观测矢量si(i，w(i))，

si(i，w(i))＝[s₁(i)，s_3W(i，w(i))]^T＝[s₁(i)，w^T(i)s_3v(i)]^T

式中s_3W(i，w(i))＝w^T(i)s_3v(i)

w(i)＝[r₃₁(i-4，i)，r₃₁(i-3，i)，r₃₁(i-2，i)，r₃₁(i-1，i)，

           r₃₁(i，i)，r₃₁(i+1，i)，r₃₁(i+2，i)，r₃₁(i+3，i)，r₃₁(i+4，i)]^T

s_3v(i)＝[s₃(i-4)，s₃(i-3)，s₃(i-2)，s₃(i-1)，

              s₃(i)，s₃(i+1)，s₃(i+2)，s₃(i+3)，s₃(i+4)]^T

$r_{31}(m,i)=\frac{\left|\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{s}_3(m+k)s_1^{*}(i+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}}$

m＝i-4，i-3，i-2，i-1，i，i+1，i+2，i+3，i+4

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

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

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

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

C.根据相关加权观测矢量估计协方差矩阵为

C_si(i，w(i))＝E{si(i，w(i))si^H(i，w(i))}

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

D.由协方差矩阵得到代价函数：

J_BF＝a^H(φ_i)C_si(i，w(i))a(φ_i)

式中 $a\left({\mathrm{\φ}}_i\right)={[1,e^{{\mathrm{j\φ}}_i}]}^T$ 称为导向矢量，φ_i的取值范围是[-π，π]；

E.通过代价函数估计干涉相位，即用导向矢量和协方差矩阵进行波束扫描，并将波束形成的输出功率最大值所对应的干涉相位

作为干涉相位估计结果。

本发明与现有技术相比，具有以下优点：

(1)采用相关加权联合单像素模型，能够同时利用空域和复图像域的信息，比直接利用空域信息的处理结果更精确；

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

(3)可以在SAR图像配准精度很差的条件下准确地估计相应像素间的干涉相位

对本发明的目的、特征、优点可通过如下附图和实例详细描述。

附图说明

图1是本发明的方法流程图；

图2是本发明构造相关加权观测矢量的主图像1和粗配准图像3示意图；

图3(a)是精确配准时采用现有中值滤波方法对干涉相位图的处理结果；

图3(b)是配准误差为0.5像素时采用现有中值滤波方法对干涉相位图的处理结果；

图3(c)是配准误差为1个像素时采用现有中值滤波方法对干涉相位图的处理结果；

图4(a)是精确配准时采用现有均值滤波方法对干涉相位图的处理结果；

图4(b)是配准误差为0.5像素时采用现有均值滤波方法对干涉相位图的处理结果；

图4(c)是配准误差为1个像素时采用现有均值滤波方法对干涉相位图的处理结果；

图5(a)是精确配准时采用本发明方法得到的干涉相位图；

图5(b)是配准误差为0.5像素时采用本发明方法得到的干涉相位图；

图5(c)是配准误差为1个像素时采用本发明方法得到的干涉相位图；

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

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

具体实施方式

参照图1，实现本发明的过程如下

1.图像粗配准

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

2.构造相关加权观测矢量

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

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

$r_{31}(m,i)=\frac{\left|\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{s}_3(m+k)s_1^{*}(i+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(1\right)$

m＝i-4，i-3，i-2，i-1，i，i+1，i+2，i+3，i+4

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

i表示当前要计算的主图像1和粗配准图像3的像素标号，上标*表示共轭操作；

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

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

s_3W(i，w(i))＝w^T(i)s_3v(i)(2)

式中

s_3v(i)＝[s₃(i-4)，s₃(i-3)，s₃(i-2)，s₃(i-1)，

                                                        (3)

s₃(i)，s₃(i+1)，s₃(i+2)，s₃(i+3)，s₃(i+4)]^T

w(i)＝[r₃₁(i-4，i)，r₃₁(i-3，i)，r₃₁(i-2，i)，r₃₁(i-1，i)，

                                                        (4)

r₃₁(i，i)，r₃₁(i+1，i)，r₃₁(i+2，i)，r₃₁(i+3，i)，r₃₁(i+4，i)]^T

其中上标T表示转置操作；

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

si(i，w(i))＝[s₁(i)，s_3W(i，w(i))]^T＝[s₁(i)，w^T(i)s_3v(i)^T  (5)

3.估计协方差矩阵

用相关加权观测矢量si(i，w(i))估计相应的协方差矩阵C_si(i，w(i))，公式如下：

C_si(i，w(i))＝E{si(i，w(i))si^H(i，w(i))}(6)

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

4.估计干涉相位

a.由协方差矩阵得到代价函数：

J_BF＝a^H(φ_i)C_si(i，w(i))a(φ_i)(7)

式中 $a\left({\mathrm{\φ}}_i\right)={[1,e^{j{\mathrm{\φ}}_i}]}^T$ 称为导向矢量，φ_i的取值范围是[-π，π]；

b.通过代价函数估计干涉相位，即用导向矢量和协方差矩阵进行波束扫描，并将波束形成的输出功率最大值所对应的干涉相位

作为干涉相位估计结果，过程如下：

ba.以0.01弧度为步长，将φ_i在[-π，π]范围内的每一个取值代入代价函数J_BF中，分别计算其结果；

b_b.对所计算的结果进行排列，找出其中的最大值，即波束形成输出功率的最大值；

bc.将找到的最大值所对应的相位

作为干涉合成孔径雷达的干涉相位估计结果；

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

5.用所估计的干涉相位结果进行后续的相位展开，得到图像的数字高程图。

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

仿真数据描述：在两颗卫星对应的有效垂直基线为281.46米，卫星高度750公里，下视角45°的条件下，利用一幅图像分辨率为3米×3米的实测SAR图像作为地面场景背景来产生每一个SAR像素的回波功率，并利用两维汉明窗来模拟高程地形，生成一对InSAR图像。SAR图像中的信噪比SNR为23dB，其相干系数由垂直基线长度、局部地形坡度和SNR决定。

图3为在不同配准误差的情况下采用现有中值滤波方法对干涉相位图的处理结果：

其中图3(a)是精确配准时采用现有中值滤波方法对干涉相位图的处理结果，由图中可见，所得到的干涉相位图的条纹边缘很清晰，易知其估计的干涉相位结果很好。

其中图3(b)是配准误差为0.5像素时采用现有中值滤波方法对干涉相位图的处理结果，由图中可见，所得到的干涉相位图的条纹边缘很不清晰，易知其估计的干涉相位结果很差。

其中图3(c)是配准误差为1个像素时采用现有中值滤波方法对干涉相位图的处理结果，由图中可见，当配准误差达到1个像素时，已经无法得到干涉条纹，也就得不到所需要的干涉相位。

图4为在不同配准误差的情况下采用现有均值滤波方法对干涉相位图的处理结果：

其中图4(a)是精确配准时采用现有均值滤波方法对干涉相位图的处理结果，由图中可见，所得到的干涉相位图的条纹边缘很清晰，易知其估计的干涉相位结果很好。

其中图4(b)是配准误差为0.5像素时采用现有均值滤波方法对干涉相位图的处理结果，由图中可见，所得到的干涉相位图的条纹边缘很不清晰，易知其估计的干涉相位结果很差。

其中图4(c)是配准误差为1个像素时采用现有均值滤波方法对干涉相位图的处理结果，由图中可见，当配准误差达到1个像素时，已经无法得到干涉条纹，也就得不到所需要的干涉相位。

图3和图4表明，现有的中值滤波方法和均值滤波方法受图像配准误差的影响很大，当配准误差达到1个像素时，已经无法得到干涉条纹。

图5为在不同配准误差的情况下采用本发明方法得到的干涉相位图：

其中图5(a)是精确配准时采用本发明方法得到的干涉相位图，由图中可见，所得到的干涉相位图的条纹边缘很清晰，易知其估计的干涉相位结果很好。

其中图5(b)是配准误差为0.5像素时采用本发明方法得到的干涉相位图，由图中可见，所得到的干涉相位图的条纹边缘很清晰，易知其估计的干涉相位结果很好。

其中图5(c)是配准误差为1个像素时采用本发明方法得到的干涉相位图，由图中可见，所得到的干涉相位图的条纹边缘很清晰，易知其估计的干涉相位结果很好。

由图5可见，本发明方法对干涉相位的估计结果几乎不受图像配准误差的影响，即使配准误差达到1个像素，仍能获得几乎和精确配准一样的处理结果。

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

图6为实测数据ERS-1/2的处理结果：图6(a)是对两幅实测SAR图像经过粗配准后生成的干涉相位图，图6(b)是采用本发明方法对两幅实测SAR图像进行处理得到的干涉相位图。

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

Claims (1)

1.一种基于相关加权的干涉合成孔径雷达干涉相位估计方法，包括如下过程：

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

B.利用主图像(1)和粗配准图像(3)，按照如下过程构造相关加权观测矢量si(i，w(i))：

(B1).分别计算主图像(1)中像素i与粗配准图像(3)中像素m的相关系数：

$r_{31}(m,i)=\frac{\left|\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}{s}_3(m+k)s_1^{*}(i+k)\right|}{\sqrt{\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\left|s_3(m+k){|}^2\underset{k=-K}{\overset{K}{\mathrm{\Σ}}}\right|s_1(i+k){|}^2}}$

m＝i-4，i-3，i-2，i-1，i，i+1，i+2，i+3，i+4，

(B₂₎用所得到的相关系数r₃₁(m，i)分别对粗配准图像(3)中的像素m进行加权，得到元素：

s_3W(i，w(i))＝w^T(i)s_3v(i)；

(B3)用主图像(1)中的像素i和得到的元素s_3W(i，w(i))构造相关加权观测矢量为：

si(i，w(i))＝[s₁(i)，s_3W(i，w(i))]^T＝[s₁(i)，w^T(i)s_3v(i)]^T式中：

w(i)＝[r₃₁(i-4，i)，r₃₁(i-3，i)，r₃₁(i-2，i)，r₃₁(i-1，i)，

r₃₁(i，i)，r₃₁(i+1，i)，r₃₁(i+2，i)，r₃₁(i+3，i)，r₃₁(i+4，i)]^T

s_3v(i)＝[s₃(i-4)，s₃(i-3)，s₃(i-2)，s₃(i-1)，

s₃(i)，s₃(i+1)，s₃(i+2)，s₃(i+3)，s₃(i+4)]^T

s₁和s₃分别为主图像(1)和粗配准图像(3)所对应的复图像数据，

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

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

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

C.根据相关加权观测矢量估计协方差矩阵：

C_si(i，w(i))＝E{si(i，w(i))si^H(i，w(i))}

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

D.由协方差矩阵得到代价函数：

J_BF＝a^H(φ_i)C_si(i，w(i))a(φ_i)

式中， $a\left({\mathrm{\φ}}_i\right)=[1,e^{j{\mathrm{\φ}}_i}{]}^T$ 称为导向矢量，φ_i的取值范围是[-π，π]；

E.以0.01弧度为步长，将φ_i在[-π，π]范围内的每一个取值代入代价函数J_BF中，分别计算其结果，并将所计算的结果进行排列，找出其中最大值所对应的干涉相位

作为干涉相位估计结果。

Images