基于频谱的SAR图像低散射区域散射值的估计方法

摘要

一种基于频谱的SAR图像低散射区域散射值的估计方法，包括步骤如下：根据SAR接收的原始数据生成并保存单视复图像；将单视复图像分成n个子图像块，对每个子图像块进行多普勒中心频率估计并将多普勒中心频率移至零频位置；计算每个子图像块的散射率初值和多普勒中心频率移至零频的多普勒频谱；根据每个子图像块的散射率初值和多普勒频谱，运用最大似然估计和牛顿迭代法估计出每个子图像块的相对散射值；计算每个子图像块的绝对散射值。本发明从频域推导出目标散射值和频谱的关系，并推得相关方程式，经过最大似然估计方法和牛顿迭代法估计方程中的散射值，且牛顿迭代法的初值考虑了系统噪声和方位模糊因素，因此估计结果更加准确。

说明书

技术领域

本发明属于合成孔径雷达(Synthetic Aperture Radar，SAR)图像处理领域，更具体地涉及一种基于频谱的SAR图像低散射区域散射值的估计方法。

背景技术

在合成孔径雷达(SAR)图像中，低散射区域的归一化雷达截面散射值在图像中呈现为暗区域。海洋中SAR图像上暗区域出现的频率较高，例如海面溢油、有机薄膜、低风速区、锋面、上升流、内波的暗条带和涌浪；陆地SAR图像中也有很多暗区域情况，例如山的阴面、光滑的机场跑道等都是经典的低散射区域。

一方面，SAR图像中低散射区域的信号强度接近甚至低于SAR系统的噪底，从而影响低散射区域的图像强度。以海洋表面为例，在中低风速和低入射角下，对于L、C和X波段海洋表面的平均散射值范围从-15dB到-25dB，而海洋表面低散射区域的散射值要比海洋表面的平均散射值更低，海洋表面低散射区域的散射值通常要低于-30dB，但大多数星载SAR系统的等效噪声系数(NESZ)的范围是-20dB到-30dB。另外，海洋SAR图像中低散射区域的后向散射信号强度通常要低于SAR系统的噪底。

另一方面，强散射区域的方位模糊会影响低散射区域的图像强度，这种影响主要是因为天线方位向旁瓣照射的区域反射回来的信号的多普勒频率超过了雷达脉冲重复频率(PRF)。目标的方位模糊信号位置和真实位置有一定偏移，这个偏移取决于PRF、平台速度和SAR系统的多普勒中心频率。经典星载SAR的方位模糊值大约是-15dB到-20dB。如果强散射区域在低散射区域的模糊信号强度比低散射区域的信号强度高出15dB到20dB，方位模糊信号就会严重影响低散射区域的图像强度。方位模糊信号在陆地和水面的交界处出现较为频繁，因为陆地目标的散射值比水面要强得多。

以上分析显示，要估计低散射区域的散射值的准确值，就要考虑系统噪声和方位模糊的影响。而SAR图像的标准辐射定标算法只考虑了系统噪声，却忽视了方位模糊的影响。现有提高判读图像质量的方法是简单的去除系统噪声或者是从图像域优化图像质量。当进行校准时，低散射区域的散射值会不可避免地被系统噪声和方位模糊所影响，这就会导致对后期图像的处理和判读带来困惑。

发明内容

(一)要解决的技术问题

基于以上问题，本发明的目的在于提出一种基于频谱的SAR图像低散射区域散射值的估计方法，用于解决以上技术问题中的至少之一。

(二)技术方案

为了实现上述目的，本发明提出了一种基于频谱的SAR图像低散射区域散射值的估计方法，其包括步骤如下：

步骤1、根据SAR接收的原始数据生成并保存单视复图像；

步骤2、将单视复图像分成n个子图像块，对每个子图像块进行多普勒中心频率估计，并将多普勒中心频率移至零频位置；其中，n为大于等于1的正整数；

步骤3、计算每个子图像块的散射率初值和多普勒中心频率移至零频的多普勒频谱；

步骤4、根据每个子图像块的散射率初值和多普勒频谱，运用最大似然估计和牛顿迭代法估计出每个子图像块的相对散射值；

步骤5、计算每个子图像块的绝对散射值。

进一步地，上述每个子图像块的散射率初值表示为：

其中，为图像强度的相对值，x₀、y₀分别是方位向和距离向的坐标，N₀是SAR系统的固有噪声，是像素位置[x₀+D_x-L/2，y₀+D_y]和[x₀+D_x+L/2，y₀+D_y]之间的平均相对散射值，是像素位置[x₀-D_x-L/2，y₀-D_y]和[x₀-D_x+L/2，y₀-D_y]之间的平均相对散射值，L是计算多普勒频谱的数据长度；D_x和D_y分别是方位模糊信号位置和真实目标位置之间在方位向和距离向的偏移量，σ_min为绝对散射值的最小值，A为方位模糊因子。

进一步地，上述牛顿迭代法的迭代公式为：

其中，M是每个子图像块的多普勒谱的点数，F_r是SAR系统的雷达脉冲重复频率，P_c(f_i)、P_L(f_i)、P_R(f_i)分别为多普勒频谱主瓣、左旁瓣和右旁瓣的第i个点数的能量谱，p_n-X(f_i)、p_n(f_i)、p_n+X(f_i)分别是第n-X个、第n个和第n+X个子图像块的多普勒频谱的第i个点数的能量谱，α表示天线俯仰角，分别是第n-2X个、第n-X个、第n个、第n+X个和第n+2X个子图像块的相对散射值，i为大于等于1的自然数。

进一步地，上述最大似然估计的公式表示为：

其中，上式中的g(p_n(f₁)>n(f₂) …>n(f_M))、g(p_n-X(f₁)>n-X(f₂) …>n-X(f_M))、g(p_n+X(f₁)>n+X(f₂) …>n+X(f_M))分别是第n个、第n-X个和第n+X个子图像块的概率密度函数，p_n(f_i)是第n个子图像块的第i个频率点的能量谱。

进一步地，根据以下公式计算每个图像块的绝对散射值：

其中，R、α和G分别表示目标斜距、天线俯仰角和特定像素的系统增益，g(α)为俯仰角为α时的两路天线增益，K是定标常数，R_ref、α_ref和G_ref分别表示参考位置处的斜距、俯仰角和系统增益。

进一步地，上述牛顿迭代法的初值为每个子图像块的散射率初值。

进一步地，上述步骤1中生成单视复图像时，方位向的压缩采用无权匹配滤波器。

进一步地，上述单视复图像的方位向能量谱与SAR原始数据具有相同的模型特征：

p_s(f，x₀，y₀)＝p_r(f，x₀，y₀)|H(f)|²＝p_r(f，x₀，y₀)；

其中，是方位向无权匹配滤波器，p_s(f，x₀，y₀)和p_r(f，x₀，y₀)分别是单视复图像和SAR原始数据的能量谱，f为多普勒中心频率，λ为SAR系统波长，V是SAR平台的运动速度。

(三)有益效果

本发明提出的一种基于频谱的SAR图像低散射区域散射值的估计方法，具有以下有益效果：

1、本发明从频域推导出目标散射值和频谱的关系，并推得相关方程式，经过最大似然估计方法和牛顿迭代法估计方程中的散射值，且牛顿迭代法的初值考虑了系统噪声和方位模糊因素，因此估计结果更加准确；

2、本发明从原始数据进行处理生成单视复图像，且在方位向的压缩采用无权匹配滤波器，只会改变多普勒频谱的相位，而不会改变多普勒普的幅值，因此从本质上提取数据的有用信息，避免了多普勒普强度的衰减；

3、对于低散射区域，本发明比传统去除系统噪声的方法效果更佳，经过本发明处理后会有更加明显的纹理特征。

附图说明

图1是本发明一实施例提出的SAR原始数据和其他各种能量成分的图解；

图2是本发明一实施例提出的基于频谱的SAR图像低散射区域散射值的估计方法的方法流程示意图；

图3是本发明一实施例提出的基于频谱的SAR图像低散射区域散射值的估计方法的详细流程示意图；

图4是本发明一实施例提出的选取的ERS-2卫星数据通过无权值匹配滤波成像后单视复图像；

图5是本发明一实施例提出的多普勒中心频率移动前后的对比图；

图6(a)是本发明一实施例提出的估计散射值方法的处理结果；

图6(b)是传统估计散射值方法的处理结果；

图7是本发明一实施例提出的估计散射值方法的估计结果图6(a)和图6(b)中白线位置的剖面图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

一种SAR图像的标准辐射定标算法如下所示：

其中，I、R、α和G分别是图像强度、目标斜距、天线俯仰角和特定图像像素的系统增益；g(α)是在俯仰角为α时两路天线的增益，N₀是系统噪声，K是定标常数，R_ref、α_ref和G_ref分别是参考位置处的斜距、俯仰角和系统增益。

然而，在商业星载SAR数据产品中几乎不会提供一个精确的系统噪声N₀。尽管提供了比较精确的系统噪声N₀，低散射区域的散射值估计依然可能会获得小于等于零的毫无意义的散射值，因为低散射区域的图像强度本身就是一个很可能低于系统噪声N₀的随机变量。况且，在大多数实用的散射值定标应用中，系统噪声也是被忽略的，因此，等式(1)又被简化为：

图1给出了实际SAR数据中包含的各种能量的成分图解，从图中可以看出，实际SAR数据中包括SAR原始数据、目标左侧和右侧强目标造成的方位模糊、系统噪声。因此需要去除系统噪声和方位模糊信号，才能更精确的得到SAR原始数据信息。

基于以上原因，本发明公开了一种基于频谱的SAR图像低散射区域散射值的估计方法，其包括步骤如下：

步骤1、根据SAR接收的原始数据生成并保存单视复图像；

步骤2、将单视复图像分成n个子图像块，对每个子图像块进行多普勒中心频率估计，并将多普勒中心频率移至零频位置；

步骤3、计算每个子图像块的散射率初值和多普勒中心频率移至零频的多普勒频谱；

步骤4、根据每个子图像块的散射率初值和多普勒频谱，运用最大似然估计和牛顿迭代法估计出每个子图像块的相对散射值；

步骤5、计算每个子图像块的绝对散射值。

其中步骤1中生成单视复图像时，方位向的压缩成像采用无权匹配滤波器，因为在大多数应用中，SAR原始信号的方位分辨率太低，为了提高方位分辨率，方位向匹配滤波器必须要应用在SAR原始信号上，并且将其转换成单视复图像。选择使用方位向无权匹配滤波器，则方位向匹配滤波器只会改变多普勒频谱的相位，而不会改变多普勒频谱的幅值。所以，单视复图像的方位向能量谱与SAR原始数据有相同的模型特征。单视复图像和SAR原始数据之间的关系式可以表示为：

p_s(f，x₀，y₀)＝p_r(f，x₀，y₀)|H(f)|²＝p_r(f，x₀，y₀)；(3)

其中，是方位向无权匹配滤波器，p_s(f，x₀，y₀)和p_r(f，x₀，y₀)分别是单视复图像和SAR原始信号的能量谱，f为多普勒频率，λ为SAR系统波长，V是SAR平台的运动速度。

对每个子图像块进行多普勒中心频率估计，并将多普勒中心频率移至零频位置，其详细过程为对每个子图像块的时域图像进行傅里叶变换得到每个子图像块的多普勒频谱，然后将多普勒谱的中心频率移至零频位置，再进行方位向的逆傅里叶变换，即可得到多普勒谱中心频率移动的单视复图像的每个子图像块。

计算每个子图像块的散射率初值和多普勒中心频率移至零频的多普勒频谱：

此处先对SAR原始信号的多普勒频谱进行分析，根据雷达回波录取原理，其形式可以表示为：

方位向的多普勒频谱是由天线方向图的主瓣和其他位置落在该区域的旁瓣与散射率的加权得到的；其中，(x₀，y₀)是进行傅里叶变换区域的中心位置，x₀和y₀分别是方位向和距离向的坐标。E[·]是数学期望，p_r(f)是SAR原始数据的方位向能量谱。P_a(f)是散射值为0dB的理想点目标的能量谱，它的模型取决于双向天线方向图。另外，f₀是多普勒中心频率，F_r是SAR系统的脉冲重复频率，是像素位置[x₀-nD_x-L/2>0-nD_y]和[x₀-nD_x+L/2>0-nD_y](L是计算多普勒频谱的数据长度)之间的平均散射值，D_x和D_y分别是方位模糊信号位置和真实目标位置之间在方位向和距离向的偏移量，它们可以被写为：

在公式(4)中，n＝0对应于天线主瓣反射回来的信号，n≠0对应于方位模糊信号的影响。通常情况下，方位模糊信号只考虑n＝-1和1的情况，也就是方位天线的第一旁瓣。因此，等式(4)可以简化为：

单视复图像的一小块的真实能量谱实际上是一个随机过程。因为单视复图像的信号是一个复高斯过程，所以能量谱的每个样本的概率密度函数服从著名的指数分布

由等式(6)可知，特定区域的后向散射信号由三个频谱构成：p_s(f，x₀，y₀)，p_s(f，x₀-D_x，y₀-D_y)，p_s(f，x₀+D_x，y₀+D_y)。因此所有频率点的联合概率密度函数是：

每个子图像块的散射率初值为根据SAR的工作原理分析推导出的值，可以表示为：

其中，为图像强度的相对值，是像素位置[x₀+D_x-L/2，y₀+D_y]和[x₀+D_x+L/2，y₀+D_y]之间的平均相对散射值，是像素位置[x₀-D_x-L/2，y₀-D_y]和[x₀-D_x+L/2，y₀-D_y]之间的平均相对散射值，σ_min为绝对散射值的最小值，通过对图像的散射值分析得到；A为方位模糊因子，其可表示为

每个子图像块多普勒中心频率移至零频的多普勒频谱即：将多普勒谱中心频率移动的单视复图像的每个子图像块的时域图像，再次经过傅里叶变换得到多普勒频谱。

根据每个子图像块的散射率初值和多普勒频谱，估计出每个子图像块的相对散射值：根据多普勒频谱估计散射值是一个经典的贝叶斯估计问题，整个单复视图像的估计方程可以表示为

为清楚表示，此处将g函数分成两行表示，实则应为p(f_{i-X_1})，...，p(f_{i_M})....

对于单个子图像块，其估计方程为

然后运用最大似然估计和牛顿迭代法通过求解下面的方程，就可以得到等式(10)的估计结果，最大似然估计方法可表示为

其中，上式中的g(p_n(f₁)>n(f₂) …>n(f_M))、g(p_n-X(f₁)>n-X(f₂) …>n-X(f_M))、g(p_n+X(f₁)>n+X(f₂) …>n+X(f_M))分别是第n个、第n-X个和第n+X个子图像块的概率密度函数，p_n(f_i)是第n个子图像块的第i个频率点的能量谱。

牛顿迭代公式为

其中，M是每个子图像块的多普勒谱的点数，F_r是SAR系统的雷达脉冲重复频率，P_c(f_i)、P_L(f_i)、P_R(f_i)分别为多普勒频谱主瓣、左旁瓣和右旁瓣的第i个点数的能量谱，p_n-X(f_i)、p_n(f_i)、p_n+X(f_i)分别是第n-X个、第n个和第n+X个子图像块的多普勒频谱的第i个点数的能量谱，α表示天线俯仰角，分别是第n-2X个、第n-X个、第n个、第n+X个和第n+2X个子图像块的相对散射值，i为大于等于1的自然数。

公式(12)是将公式(8)带入到公式(11)中，再运用约束求解方法，使得估计结果不出现负值，所以有公式(12)的最后一项公式。的迭代初值如下即为每个子图像块的散射率初值见公式(9)。其迭代初值为每个子图像块的散射率初值，直至迭代收敛得到准确散射率的相对值

计算每个图像块的绝对散射值：由于迭代得到的是一个相对的后向散射值，而不是绝对散射值。因此如果已知辐射定标的K常数，用迭代得到的代替公式(1)中的I-N₀，即可将估计的相对后向散射值转换为绝对散射值，表达式如下：

以下通过具体实施例对本发明提出的基于频谱的SAR图像低散射区域散射值的估计方法进行详细描述。

实施例

如图2所示，本实施例提出一种基于频谱的SAR图像低散射区域散射值的估计方法，其包括步骤如下：

步骤1、根据SAR接收的原始数据生成并保存单视复图像；

步骤2、将单视复图像分成n个子图像块，对每个子图像块进行多普勒中心频率估计，并将多普勒中心频率移至零频位置；

步骤3、计算每个子图像块的散射率初值和多普勒中心频率移至零频的多普勒频谱；

步骤4、根据每个子图像块的散射率初值和多普勒频谱，运用最大似然估计和牛顿迭代法估计出每个子图像块的相对散射值；

步骤5、计算每个子图像块的绝对散射值。

该方法的详细流程图如图3所示：

首先，根据SAR接收的原始数据，在方位向采用无权匹配滤波器生成单视复图像并保存此图像，本实施例选取ERS-2卫星的数据进行绝对散射值的求解，通过无权值匹配滤波成像后的单视复图像如图4所示；

然后将单视复图像分成n个子图像块，图4中白色框部分为其中的一个子图像块，采用傅里叶变换将n个子图像块的时域图像转化成多普勒频谱，将得到的多普勒频谱的中心频率移至零频位置后，采用方位向逆傅里叶变换得到中心频率移动的每个子图像块的时域图像，中心频率未移动和移动后的多普勒频谱的对比图如图5所示；

根据SAR的工作原理，通过公式(9)计算出每个子图像块的散射率初值，同时采用傅里叶变换将中心频率移动的每个子图像块的时域图像转换得到每个子图像块的多普勒频谱；

根据每个子图像块的散射率初值和多普勒频谱，采用公式(10)-公式(12)，通过最大似然估计和牛顿迭代法得到每个子图像块的相对散射值，迭代过程中，每迭代一次得到的值应该与上一次迭代得到的值相比较，直到循环收敛为止。

根据迭代得到的每个子图像块的相对散射率值，采用公式(13)得到每个子图像块的绝对散射值。

图6(a)和(b)分别是图4中白色矩形框部分通过本实施例处理方法和传统减去噪声的方法得到的两幅图像；

图7中的两条曲线分别是图6两幅图中白线部分的截面图。白线经过的三条线是海洋内波。其中虚线是图6(b)通过传统减去噪声的方法处理后的结果，实线是图6(a)通过本实施例提出的方法处理后得到的结果。从图7中的对比可以看出不论是传统处理方法还是本实施例的处理方法，对内波波峰(相对强散射区域)处理结果是相同的，而对低散射区域的内波波谷地区，传统方法只能得到大约10db的能量强度，而本实施例处理后可以将波谷位置的能量强度恢复至-25到-30db左右，这说明本实施例恢复暗区域纹理特征的能力更接近内波真实值。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。