CN102472815B - 对从在相同区域上采集的sar图像获得的干涉图进行滤波的方法 - Google Patents

Abstract

对从SAR图像得到的干涉图进行滤波的方法，SAR图像由合成孔径雷达在相同区域上采集，该方法包括以下步骤：a)利用采集观测系统，通过在相同区域上的SAR传感器，采集N个系列雷达图像(Al..AN)，以便允许对公共栅格上的数据进行重新采样；b)在公共栅格上进行重新采样之后，从公共栅格选择像素；c)计算所选择的像素的相干矩阵，即对每个可能的可用图像对，估算复数相干值；d)对于源矢量θ，使这里未知的元素最大化，泛函如下，其中，

为提取复数实部的运算符，γ_nm为相干矩阵元素(n，m)的模数，k为正实数，φ_nm为相干矩阵元素(n，m)的相位，θ_n和θ_m为未知矢量θ的元素n和m元素。假设只在泛函T中出现相位差，则位置矢量的值将较少被估算为加常数，例如，可以通过设定θ₁＝0来将其固定。由此得到的相位值可以构成相位值经过滤波的矢量。

Description

对从在相同区域上采集的SAR图像获得的干涉图进行滤波的方法

本发明涉及对从在相同区域上采集的SAR图像获得的干涉图(interferogram)进行滤波的方法。

众所周知，合成孔径雷达(synthetic aperture radar，SAR)或SAR系统产生二维图像。二维图像中的一维称为距离(range)，为雷达与被照射目标之间的直线距离。二维图像中的另一维称为方位角(azimuth)，并且与“距离”垂直。

SAR型雷达一般在400Mhz与10Ghz之间的频率下运行，并且通常安装在飞机中或者安装在250到800公里高度的轨道上运行的人造卫星平台上。雷达天线以垂直于(飞机或人造卫星)平台运动的方向，即以垂直于地面的方向，瞄准地面，相对于天底方向(nadir direction)具有在20到80度之间的天底偏角(off-nadir angle)。

利用所述系统，通过沿着方位角方向，利用文献中众所周知的算法，对远大于实际尺寸的天线进行合成(其传感器由此得名)，可以获得空间分辨率为几米的地球表面图像。

SAR最重要的特性在于，其传感器为相干型的，因此，图像为复数矩阵，其中，值的大小取决于来自被照射物体的反射功率(即，到它们的雷达截面)，而相位决定于目标的属性以及到雷达的距离。至于雷达图像I的每个像素，用一定距离坐标r和方位角坐标a表示，由此形成一个复数：

$I(r,a)=x+i\·y=A\·e^{\mathrm{i\ψ}}$

式中，x和y表示该数的实部和虚部，A为其幅值，为ψ相位值，i为虚数单位或-1的平方根。

假设不论晴天以及多云天气都能够获得图像，那么，SAR成像适合于各种应用；其中，首要的是与目标识别和分类以及变化检测相关的应用和干涉测量应用。后者通常用于得到数字高程模型和/或根据多时相(multi-temporal)SAR数据集，对地面的表面变形进行分析。

给定在相同的目标区域上获得的两个一般SAR图像，标示为I_n和I_m，干涉图Φ_nm被定义为将一个图像与另一个图像的共轭复数进行复数相乘的结果：

${\mathrm{\Φ}}_{\mathrm{nm}}=I_n\·I_m^{*}=A_n\·e^{i{\mathrm{\ψ}}_n}\·A_m{e}^{-{\mathrm{i\ψ}}_m}=A_n{A}_m{e}^{i({\mathrm{\ψ}}_n-{\mathrm{\ψ}}_m)}$

式中，I^*表示通过改变相位值的符号而被变换成共轭复数的图像。它跟随干涉图的相位，也被称为干涉相位(interferometric phase)，由两个图像的相位差决定。

(用距离坐标r和方位角坐标a表示的)SAR图像的每个像素的相位可以被看作与被照射物体的属性相关的，称为“反射率相位(reflectivityphase)”的贡献ξ，以及与电磁波的光程(optical path)相关，由此与传播装置的特性和传感器目标距离相关的贡献d之和：

ψ(r，a)＝ζ(r，a)+d(r，a)

如果目标的电磁特性保持为不随时间变化，那么，在与反射率相位相关的项(ζ)中将没有变化，因此，考虑采样数量而揭示的任何相位变化将与光程的可能变化相关。

特别感兴趣的是，导致对运动的历史序列进行重构的应用，即对地面上反射率保持不变并且能够分离传输装置产生的相位贡献的目标随时间的任何运动进行监控的能力。

可操作地，假定在不同时间、在相同区域上采集的一系列N个SAR图像，但是在公共栅格上重复采样，人们希望对图像的N个相位值的历史序列的每个像素(即，对作用在感兴趣的区域上的每个采样值)进行计算，对其应用合适的算法，以便估算雷达目标沿着由传感器视线指定的方向的运动。

估算的质量很大程度上取决于一定分辨率的元素的反射相位随时间保持不变的事实。如果这个假设被证实，则通过计算各个采样相对于例如第一图像的相位差，可以得到仅由于光程变化引起的贡献。因此，应该注意，仅是两个采样之间的相位差而不是单个图像的相位提供信息，假定不知道反射率相位值并且假定从像素到像素的变化：因此，使任何表面变形图能够可视化的是干涉测量相位。

在所有可用采样中，真实雷达场景的反射率值通常不是常数，除了数量有限的称为永久散射体的目标以外，可以对这些目标应用特殊算法(欧洲专利EP-1183551，于2007年11月27日提交的第MO2007A000363号意大利专利申请)。

但是，有许多其他像素，可以仅在某些干涉图上提取与光程相关的信息(即感兴趣的信号)，或者，更一般地说，信噪比依据所考虑的图像对而显著变化。这意味着，反射率相位项(ζ)恒定不变的假说可以被部分满足。

导致反射率变化的主要机制有两个：(a)时间去相关(temporaldecorrelation)，即目标的电磁特性随时间变化；(b)几何或空间去相关，由采样几何学方面的变化导致。第一种机制取决于干涉图所谓的“时间基线”，或者，在不同时间获得了一组图像的情况下，用于生成干涉图的两个图像之间的时间差。第二种机制却取决于所谓的“几何基线”，即在进行两次采样期间，传感器跟踪的弹道轨迹之间的距离。

为了测量由一般图像n和m生成的干涉图的信噪比(即质量)，通常的做法是利用文献已知的参数，称为干涉相干性(interferometric coherence)ρ_nm，或者，更简单地说，相干性。一般情况下，像素与像素之间的相干性如下式：

${\mathrm{\ρ}}_{\mathrm{nm}}=\frac{E(I_n\·I_m^{*})}{\sqrt{E(I_n\·I_n^{*})\·E(I_m\·I_m^{*})}}$

式中，E(.)表示称为“期望”统计运算符。可操作地，期望运算符由在以当前像素为中心的适当窗口F上计算的空间平均值代替。在选择估算窗口时，应该尽可能选择统计学上同类的样本群，以便得到可靠的统计估算，并且，估算中使用的像素数量可以因此从像素到像素而变化(2009年4月3日提交的专利申请n.MI2009A000535)。

因此，按照下式计算估算相干性(c_nm)：

$c_{\mathrm{nm}}={\mathrm{\γ}}_{\mathrm{nm}}\·e^{i{\mathrm{\φ}}_{\mathrm{nm}}}=\frac{\underset{p\∈F}{\mathrm{\Σ}}{x}_n\left(p\right)x_m^{*}\left(p\right)}{\sqrt{\underset{p\∈F}{\mathrm{\Σ}}{\left|x_n\left(p\right)\right|}^2\underset{p\∈F}{\mathrm{\Σ}}{\left|x_m\left(p\right)\right|}^2}}$

已经用γ_nm、其相位φ_nm和估算窗口F的第p个元素x(p)表示出估算的相干性模数c_nm。由此计算的相干性是一个复数，它根据像素而变化，由于归一化，其模数的范围在0和1之间(分别为最小和最大相关性，即，信噪比为零或无穷大)，并且，其相位为在估算窗口中使用的像素相位的平均值。

假设在公共栅格上重复采样的N个SAR图像，相干性c_nm可以被看作称为相干矩阵的N×N矩阵的元素，对于采样场景的每个像素，这个矩阵可以描述全部可用SAR图像的干涉测量质量。即，假设在相同区域获得了一组N个SAR图像，并且可以在公共栅格上为这些图像进行重复数据采样，对于每个像素，可以使它与一个N×N元素的矩阵相关，该矩阵中的一般元素为可用图像组中的图像n与m之间的复数相干性的估算值。

对于所谓的永久散射体，有一个模数值将趋向于常数并接近于一的相干矩阵，表示这种类型的目标的信噪比很高，对于所考虑的每个干涉测量对，保持高信噪比。如前所述，永久散射体仅是真实场景中的少部分像素。受去相关性现象影响的绝大多数像素的特征在于，可以证明，相干矩阵的幅值显著变化。换句话说，在图像之间并由此在干涉图之间，相同像素的信噪比明显变化。

在永久散射体的情况下，提取关于一般像素光程的历史序列不是特别困难的操作：对它们而言，可以生成N-1个相干图，所有相干图包含相同场景，例如第一个场景，并且，具有足够高的信噪比，从而保证应用后续算法估算目标运动。

让我们通过例子考虑，一组在公共栅格上、在不同时间(t₁..t₅)重复采样的5个SAR图像，从第一次采样的时刻t₁开始，哪一个趋向于用来重构一般像素的光程的历史序列(图1)。如果选择的像素为永久散射体，则可以简单地得到四个相干图(t₂，t₁)，(t₃，t₁)，(t₄，t₁)和(t₅，t₁)，这使得能够对关于参考时间t₁的(关于当前像素的)目标光程进行估算(图1中的图A)。一般情况下，希望创建N个元素的历史序列，(对应于时刻t₁的)第一个值将被设为零。

在选择的像素不是永久散射体的情况下，情况非常困难，并且，因此不能讲前面段落中考虑的一个或多个干涉图的相干性显示到可接受水平。第一种解决方案可以是通过反复试验，找到允许重构光程值的完整历史序列的质量良好的干涉测量对，并且将得到的结果合成到各个干涉图中(例如，如图1中图B和图C所示，这是与图A所示不同的N-1个干涉图的两种配置，但是能够得到历史序列的N个值)。但是，有理由假设，如果在对与所讨论的像素有关的相干矩阵进行分析的基础上，这样的操作将更有效，这种操作通过构成给出了数据集的所有可能干涉测量对的概要图。

在这一点应该注意，一般相干矩阵的元素不仅使得能够通过计算矩阵的参模估算干涉图的信噪比，而且能够利用相位值，为图像的每个可能对提供经过滤波的干涉测量相位值。由于可以根据前面段落中陈述的相干性的定义进行推断，因此，由在适当估算窗口F上计算的干涉测量相位值的空间平均值，给出相干矩阵的一半元素的相位φ_nm：至少在同类统计群的情况下，对于非零信噪比的干涉图，这个操作能够显著降低噪声水平，同类统计群的特征在于有相同的光程值。一方面，求平均值的过程使得噪声水平下降，然而另一方面，这意味着三角形关系将不再满足： ${\mathrm{\φ}}_{\mathrm{nm}}\≠\∠\left(e^{i{\mathrm{\φ}}_{\mathrm{nl}}}{e}^{i{\mathrm{\φ}}_{\mathrm{lm}}}\right)$

即一般情况下，相位是不连续的(例如，像永久反射器的情况那样，φ₂₁与φ₃₂相加等于φ₃₁将不再成立)。重构N个相位值的历史序列，由于反射率相位导致的贡献被最佳补偿，因而使信噪比最大，这就要求开发一种合适的算法。

由此可以将问题概括如下：给出关于一般图像像素的相干矩阵，希望得到N个相位值矢量θ＝{θ_l...θ_N}，以适当方式考虑是有可以数据，即，与是有可能干涉图(数量为N(N-1)/2)相关的相位值以及它们的质量。本发明提出了得到这个矢量的方法。

出于实现这个目标的目的，按照本发明的方法提供了以下步骤：

a)通过SAR传感器，利用采集观测系统(acquisition geometry)采集一系列N个图像(A1......AN)，使得能够在公共栅格上对数据进行重复采样；

b)在公共栅格上进行了重复采样之后，从公共栅格选择一个像素；

c)计算选择的像素的相干矩阵，即估算每对可用图像的复数相干值；

d)相对于源矢量θ，使这里的一个未知元素最大化，泛函：

这里，

为提取复数实部的运算符，γ_nm为相干矩阵元素(n，m)的模数，k为正实数，φ_nm为相干矩阵元素(n，m)的相位，θ_n和θ_m为未知矢量θ的元素n和m元素。

假设在泛函中只出现相位差，位置矢量的值将较少被估算为加常数，例如，可以通过设定θ₁＝0来将其固定。由此得到的相位值θ_n可以构成相位值经过滤波的矢量。

为了提高相干模数而对指数k进行的选择取决于如何考虑相位的权重，并取决于估算的相干值的可能极化。可操作地，将它设置为等于1或2已经得到了良好结果。重要的是指出，使泛函最大化(不考虑k值)的极好启动点是自动矢量的相位值矢量，这个自动矢量与关于相干矩阵的主自动值相关。

应该注意，尽管在很强的非线性的基础上，但所提出的最优化不要求对相干矩阵求逆；假如相干矩阵是病态的，那么这一点在操作上相当重要。还应该注意，所提出的泛函实际上是加权求和，权重与相干矩阵的模数相关：因此，将希望将重点放在具有高信噪比特征的相位项：由此得到的矢量θ将具有这样的元素，它们必须相对于特征在于高相干值的相干矩阵的元素的相位具有较大角度，即，模数项有较大值。

重要的是注意，场景的每个像素的矢量θ一旦已知，将可以用为图像的各个像素计算的第n和第m个矢量θ之间的相位差代替第n个图像与第m图像之间的一般干涉图的相位，由此创建事实上考虑了可用图像的全部数据集的、经过滤波的干涉图。

相同的方法还可用于与相位值有关的参数化估算：如果希望这些值先天已知(例如，作为各个干涉图的时间基线和几何基线的函数的多项式规律)，可以再次利用所提出的泛函并且不再使相位值t最优化来估算所述参数，但是直接使用未知参数。例如当希望估算运动的平均速度，并且一旦各个干涉图的时间基线和几何基线已知而提升相干矩阵的雷达目标时，就是这种情况。

出于说明理由本发明能够得到什么结果的目的，由RADARSAT卫星在199年5月4日和2008年1月5日在飞跃夏威夷岛期间，在相同名字的群岛中，拍摄的85张SAR图像被进行处理。在已经将图像置于公共栅格上之后(图2，它通过图形表示有N个与感兴趣的区域相关的图像可用，并且，在各次采集中同类的像素对应于系统的地面分辨单元)，将与所谓的永久散射体相关的算法(专利EP-1183551)应用于这些图像，以提取用于场景像素运动的时间系列。相同的程序在通过将本发明的处理应用于原始SAR图像而获得的数据上(在泛函中设定k＝2以最优化)，即，用根据各个矢量θ获得的相位值代替原始相位值。如2009年4月3日提交的专利申请MI2009A000535中描述的，对图像的每个像素计算这些值，利用泛函T的指数k为1，估算相干矩阵的估算下游。举例来说，图3中示出了与感兴趣的区域的像素相关的相干矩阵的幅值(应该注意与数量N＝85个可用图像对应的矩阵的维数并且其值为0到1)。图4示出了在启动数据为未经滤波的干涉图的情况下(历史序列A)与在启动数据为通过源矢量重构干涉图的情况下(下面的历史序列B)，与像素有关的运动的历史序列之间的比较。证明噪声减小。测量的时间轴(按天计量时间)示出为沿着图中的x轴，在y轴上示出了地面上的目标范围在-30到+30mm之间的估算运动。

作为第二个例子，仅为了便于图的可视化，被处理的场景的区域被剪裁(图5中高亮部分)，以便示出原始干涉图与通过优化相位矢量重构的干涉图之间的两个直接比较(图6和图7，其中示出了干涉图的相位值)。在这些特征在于高空间基线的干涉图中，某些区域的特征在于低信噪比。在图中，左侧示出了原始干涉图，在右侧可以欣赏到按照本发明的处理并且由此用矢量θ的元素的相位差代替原始干涉测量相位而重构的相同的干涉图。可操作地，一旦得到了场景的每个像素的矢量θ，则用从与当前像素有关的矢量θ提取的相位差θ_n-θ_m代替在第n与第m图像之间得到的干涉图的每个像素的相位。该效果值得注意：噪声阻碍了干涉条纹的识别能力，这里提出的技术使得这种情况显著减少，由此使得能够清楚地分辨感兴趣的信号。

Claims (7)

1.一种对从SAR图像得到的干涉图进行滤波的方法，所述SAR图像由合成孔径雷达在相同区域上采集，所述方法包括以下步骤： 

a)利用采集观测系统，通过在相同区域上的SAR传感器，采集N个系列雷达图像(A1..AN)，以便允许对公共栅格上的数据进行重新采样； 

b)在公共栅格上进行重新采样之后，从公共栅格选择像素； 

c)计算所选择的像素的相干矩阵，即对每个可能的可用图像对，估算复数相干值； 

d)对于源矢量θ，使这里未知的元素最大化，泛函： 

这里，

为提取复数实部的运算符，γ_nm为相干矩阵元素c_nm的模数，k为正实数，φ_nm为相干矩阵元素c_nm的相位，θ_n和θ_m为未知矢量θ的元素n和元素m，除了加常数以外，对未知因数的值进行估算，并且，由此得到的相位值θ_n构成了经过滤波的相位值的矢量。 

2.如权利要求1所述的方法，其特征在于，通过设定θ₁＝0来固定所述加常数。 

3.如权利要求1所述的方法，其特征在于，利用以下公式估算相干矩阵的每个元素： 

已经用围绕所选择的像素的合适估算窗口F进行表示，x(p)为估算窗口F的第p个元素，n和m为属于在公共栅格上重新采样的N个SAR图像的第n个和第m个图像。 

4.如权利要求1所述的方法，其特征在于，用多项式关系将源矢 量θ的N个元素联系起来，可归纳成以下形式： 

θ_n＝g(t_n，B_n) 

式中，t_n为第一次采集的第n个图像的采集时间，B_n为第一次采集的第n个图像的法向基线，关于多项式g的系数，泛函被最大化。 

5.如权利要求1所述的方法，其特征在于，用线性关系将源矢量θ的N个元素联系起来，可归纳成以下形式： 

θ_n＝C_v·v·t_n+C_h·h·B_n

式中，C_v和C_h为取决于采集观测系统的并且在使用的传感器上的已知参数，t_n为第一次采集的第n个图像的采集时间，B_n为第一次采集的第n个图像的法向基线，在最优化处理中，t_n和B_n也已知，仅估计v和h的值，v和h的值分别与平均位移速度以及占据所选择的像素的目标的高度有关。 

6.如权利要求1所述的方法，其特征在于，利用SAR采集的N个图像是在不同时间采集的。 

7.如权利要求1所述的方法，其特征在于，利用SAR采集的N个图像是以不同视角采集的。 

Images