CN104182942A - Sar图像方位模糊抑制方法 - Google Patents

Abstract

本发明公开了一种SAR图像方位模糊抑制方法，包括步骤：获取参考图像、对原始图像结果进行方位向加窗、计算加窗后图像的模糊信号比、设定方位向窗长、更新方位向窗长。本发明的方法通过对图像中模糊区域的检测和模糊信号比的计算，自适应调整方位向窗函数长度，有效去除方位模糊，相对于背景中的现有技术，对场景没有特殊要求，同时能够实现模糊的自适应去除。

Description

SAR图像方位模糊抑制方法

技术领域

本发明属于雷达图像去模糊的领域，具体涉及SAR图像去除方位模糊抑制方法。

背景技术

方位模糊是合成孔径雷达成像领域所要解决的一个重要问题，其由于回波多普勒信号欠采样，使得高于脉冲重复频率的多普勒信号折叠到方位频谱的中心部分，从而产生了混淆的信号。由于方位模糊对SAR图像重点目标的检测和识别造成干扰，严重时甚至会形成假目标，极大地影响了图像的质量，因此抑制方位模糊具有重要的工程意义。

目前，针对方位模糊问题的解决方法主要可分为两类：一是从天线方向图加权的角度出发，尽量抑制天线旁瓣。但是，根据方位模糊比的计算公式，利用天线方向图加权与脉冲重复频率(PRF)的大小有很大关系，只有当PRF较大时加权对抑制方位模糊才是有效的。

另一类抑制方位模糊的方法是从信号处理角度出发的，如文献“Zhang Zhimin,OnSuppressing Azimuth Ambiguities of Synthetic Aperture Radar by Three Filters,IEEEInternational Conference on Radar,624-626,2001”中，首先构造3个滤波器，其中1个是主瓣信号的匹配滤波器，另外2个是用来获得第1模糊项的参考滤波器，然后从匹配滤波器后的信号中减去参考滤波器中的模糊信号，即可得到抑制了方位模糊的图像。然而，主瓣信号和模糊项信号与匹配滤波器函数进行卷积时会引入重叠项，也就是第1模糊项虽然被去除，但又会引入新的模糊项，从而使得该方法的模糊抑制效果受限。

文献“Moreira,Suppressing the Azimuth Ambiguities in Synthetic Aperture RadarImages,IEEE Transactions on Geoscience and Remote Sensing,vol.31,no.4,pp.885-895,1993”中，通过构造2个相位刚好相差180°的校正滤波器，在复数域相减构造出模糊像，然后从带有方位模糊的图像中减去模糊像实现方位模糊的抑制。然而该方法仅适用于点目标模糊的情况，对场景要求较高，当出现区域模糊的时候，模糊抑制效果不佳。

发明内容

本发明的目的是为了解决现有方法在方位模糊抑制存在的上述缺陷，提出了一种SAR图像方位模糊抑制方法。

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

术语1、SAR方位模糊几何模型：如图1所示，假设载机平台沿着直线飞行，波束中心对准波束照射范围的中心，其中，T、T^*分别为成像平面上的目标与目标的模糊部分。

术语2、方位模糊产生原因示意图：如图2所示，由于方位向采样频率(PRF)不足以及天线旁瓣的影响导致在方位向频谱产生混叠，导致模糊信号的产生，图中的阴影部分为混叠信号的部分。

本发明的技术方案为：一种SAR图像方位模糊抑制方法，具体流程如图4所示，包括步骤：

步骤S1：获取参考图像：

输入原始成像结果I_o，大小为M×N，其中，M、N分别为目标回波距离向和方位向采样点数，其解析表达式可表示为：

$I_0(\mathrm{\τ},t;x,y)=\mathrm{\σ}(x,y)\sin c\left\{B_{a}t\right\}\sin c\left\{B\right[\mathrm{\τ}-\frac{2R(t;x,y)}{c}\left]\right\}---\left(1\right)$

其中，σ(x,y)是反射系数，B_a为方位向多普勒带宽，B为发射信号带宽，λ为载波波长，目标平面点目标距离历史为V为平台速度，为天线俯仰角，R₀为方位时间为0时刻天线与点目标的初始斜距，θ₀为0时刻点目标对应的斜视角；

方位时间向量t＝[-PRI·N/2,-PRI·(N/2-1),···,PRI·(N/2-1)]，其中，PRI为脉冲重复间隔；距离时间向量τ＝[-1/f/_s·M/2,-1/f_s·(M/2-1),···,1/f_s·(M/2-1)]，其中，f_s为距离向采样率。

对I₀(τ,t；x,y)进行方位向逆傅里叶变换(IFFT)，然后乘以一个方位窗函数W_ac(f_t)，再进行方位向傅里叶变换，得到参考图像I_c可表示为：

$\begin{array}{c}{I}_c(\mathrm{\τ},t;x,y)=\mathrm{FFT}\left\{\mathrm{IFFT}\right\{I_0(\mathrm{\τ},t;x,y)\left\}{W}_{\mathrm{ac}}\left(f_t\right)\right\}\\ =\mathrm{\σ}(x,y)\sin c\left\{B_{a}t\right\}{w}_{\mathrm{ac}}\left(t\right)\sin c\left\{B\right[\mathrm{\τ}-\frac{2R(t;x,y)}{c}\left]\right\}\end{array}---\left(2\right)$

其中，w_ac(t)是W_ac(f_t)的时域表示，W_ac(f)为频域上的矩形窗函数，即 $W_{\mathrm{ac}}\left(f_t\right)=\left\{\begin{array}{cc}1,& f_0-\frac{f_{\mathrm{ac}}}{2}\≤f_t\≤f_0+\frac{f_{\mathrm{ac}}}{2}\\ 0,& \mathrm{otherwise}\end{array}\right.,$ f₀为方位向频率中心，频宽α_c为窗函数对应的窗长，一般取方位向采样点数N的五分之一；

对参考图像I_c进行边缘检测，得到边缘图像I_ce。这里，具体利用Sobel算子进行边缘检测。

然后，对边缘图像I_ce进行膨胀操作，得到膨胀后的图像I_ced；最后，对膨胀后结果中的孔隙进行填充，并对图像进行平滑处理，获得图像分割结果I_cs；这里，具体可以采用菱形结构元素对图像进行平滑处理。

步骤S2：对原始图像结果进行方位向加窗：

对原始成像结果I_o进行方位向逆傅里叶变换，然后乘以方位向窗函数W_a(f_t)，再进行方位向傅里叶变换，得到的方位向加窗图像I₁可表示为：

I₁(τ,t；x,y)＝FFT{IFFT{I₀(τ,t；x,y)}W_a(f_t)}  (3)

其中， $W_a\left(f_t\right)=\left\{\begin{array}{cc}1,& f_0-\frac{f}{2}\≤f_t\≤f_0+\frac{f}{2}\\ 0,& \mathrm{otherwise}\end{array}\right.,$ 频宽α为窗函数对应的窗长，若是第一次加窗，窗长α设为α₀＝N。

步骤S3、计算加窗后图像的模糊信号比(AASR)，具体如图5所示：

对方位加窗后的图像I₁进行边缘检测，得到边缘图像I_oe。这里，具体利用Sobel算子进行边缘检测。

然后，对边缘图像I_oe进行膨胀操作，得到膨胀后的图像I_oed。

最后，对膨胀后结果中的孔隙进行填充，并对图像进行平滑处理，获得图像分割结果I_os。这里，具体可以采用菱形结构元素对图像进行平滑处理。

将I₁与对应分割图I_os对应的图像矩阵相乘，得到I₁中目标区域的图像矩阵M_o，将M_o中每一个数值进行平方累加求和得到图像I₁目标区域的总能量

将I₁与参考图像分割图I_cs对应的图像矩阵相乘，得到I₁中不含方位模糊的图像矩阵M_a，对矩阵M_a中的每一个数值平方累加的和就是图像I₁中不含方位模糊的能量 $S_1={\left|\right|M_a\left|\right|}_2^2.$

根据下式计算图像I₁的方位模糊信号比A1：

$A_1=\frac{S_{\mathrm{ambiguity}}}{S-S_{\mathrm{ambiguity}}}=\frac{S-S_1}{S_1}=\frac{{\left|\right|M_0\left|\right|}_2^2-{\left|\right|M_a\left|\right|}_2^2}{{\left|\right|M_a\left|\right|}_2^2}---\left(4\right)$

式中，S_ambiguity是原始图像中方位模糊区域的能量。

判断是否是第一次获得加窗后图像的模糊信号比，如果是，则转入步骤4；如果此时不是第一次获得，则继续判断图像I₁的模糊信号比A₁是否大于AASR期望值A_e(A_e为满足成像质量要求的经验值，一般预先设定)。如果大于AASR期望值A_e，则转入步骤5；如果I₁的模糊信号比A₁小于或者等于AASR期望值A_e，则输出此时的图像I₁。

步骤S4、设定方位向窗长：

计算AASR差值e＝A₁-A_e，A_e为AASR期望值；

构造权值w，令w＝μ·e·α₀，μ为步长因子；

按下式设定窗长α：α＝α₀·w

然后转入步骤S2，将设定的α作为步骤S2中的方位向窗长，按(3)式对原始图像进行加窗处理。

步骤S5、更新方位向窗长：

计算AASR差值e＝A₁-A_e；

更新权值w，令w＝w+μ·e·α₀；

更新窗长值α＝α·w；

然后转入步骤S2，将更新后的α作为步骤S2中的方位向窗长，按(3)式对原始图像进行加窗处理。

本发明的有益效果：本发明的方法通过对图像中模糊区域的检测和模糊信号比(AASR)的计算，自适应调整方位向窗函数长度，有效去除方位模糊，相对于背景中的现有技术，本发明的方法对场景没有特殊要求，同时能够实现模糊的自适应去除。

附图说明

图1方位模糊几何模型。

其中，雷达波束中心指向场景中心，T*为目标T在成像平面上的方位模糊。

图2方位模糊产生原因示意图，由于方位向采样频率(PRF)不足以及天线旁瓣的影响导致在方位向频谱产生混叠，采样后导致模糊信号的产生，图中的阴影部分就是混叠信号的部分，也就是方位模糊。

图3方位模糊三维示意图，为方位模糊点目标仿真，可以清晰的看到实际信号和模糊信号在方位向的分布，一般在实际信号左右两边都会产生模糊，左、右模糊的能量接近并且低于实际信号的能量。仿真图中间的尖峰指的是实际信号，左右两侧的尖峰指的是模糊信号。

图4本发明的SAR图像方位模糊抑制方法流程示意图。

图5模糊信号比(AASR)计算流程示意图。

图6原始成像结果示意图。

图7参考图像结果示意图。

图8最终去模糊结果示意图。

具体实施方式

下面结合附图对本发明的实施例做进一步的说明。

本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2010b上验证正确。

本实施例提供的方法为：一种SAR图像方位模糊抑制方法，具体实施步骤如下：

步骤1、获取参考图像：

该方法需要进行初始化的参数如下：光速c＝300000Km/s；雷达发射脉冲的载频f₀＝300MHz；雷达扫描周期为T＝2s；雷达发射脉冲的带宽B＝150MHz；雷达扫描空间被划分成的方位向个数N＝256，雷达在距离向上的采样点数M＝512；雷达的脉冲重复频率为PRF＝100Hz。平台的运动速度为V＝150m/s，运动高度为H＝3000m。本发明利用温哥华海湾的一处场景进行仿真，通过面目标仿真，得到原始成像结果I_o。图6为原始图像仿真结果。

为了得到参考图像I_c，窗函数W_ac(f_t)的窗长α设为方位向采样点数的五分之一左右。由于方位向采样点数N＝256，则α＝56，方位加窗后获得参考图像I_c，并获得参考图像的分割图I_cs。图7为参考图像仿真结果。

步骤2、对原始图像结果进行方位向加窗：

对原始成像结果I_o进行方位向逆傅里叶变换并乘以方位向窗函数(窗函数W_a(f_t)的窗长α初始化值设为α₀＝256)，再进行方位向傅里叶变换，得到的方位向加窗后图像为I₁，转入步骤3计算模糊信号比。

步骤3、计算加窗后图像的模糊信号比(AASR,Azimuth Ambiguous Signal Ratio)

对加窗后图像I₁进行图像分割，获得分割图I_os。通过相关计算，得到成像结果I₁的模糊信号比A₁。

第一次获得成像结果的模糊信号比时，将转入步骤4设定方位向窗长。第一次之后得到成像结果的模糊信号比大于AASR期望值A_e＝-40dB，将转入步骤5中更新方位向窗长；成像结果的模糊信号比小于或者等于AASR期望值A_e，则输出此时的图像I₁。

步骤4、设定方位向窗长

计算AASR差值e＝A₁-A_e＝23.01dB；

计算此时窗长值α＝α₀·w，式中权值w＝μ·e·α₀，步长因子μ＝0.002，既保证了较快的收敛速度，又保证了可承受范围内的稳态误差。

将设定的α＝227转入步骤2，对原始图像进行方位向加窗。

步骤5、更新方位向窗长值：

计算AASR差值，并更新此时的权值w＝w+μ·e·α₀和窗长α＝α·w，将更新后的α作为步骤2中的方位向窗长对原始图像进行加窗处理。

当成像结果的模糊信号比A₁小于或者等于AASR期望值A_e＝-40dB，输出此时的方位向加窗后的图像结果，并得到对应方位向窗长α为124。

图8是最终成像结果示意图。从图中可以看出，本发明的方法能将场景中的方位模糊全部去除的同时，而且保证较好的分辨率。

由具体实施可以看出，本发明通过图像模糊信号比(AASR)的计算，自适应调整方位向窗函数长度，从而实现图像的方位模糊抑制。面目标仿真中，最终得到的图像方位模糊信号比(AASR)从-16.99dB降到-40dB左右，完成了方位模糊的去除。与现有的其他方位模糊抑制算法相比，本发明的方法从图像质量角度出发，自适应调整方位向窗函数长度，使得图像中的方位模糊得到快速有效地去除。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.一种SAR图像方位模糊抑制方法，具体包括步骤：

步骤S1：获取参考图像：

输入原始成像结果I_o，大小为M×N，其中，M、N分别为目标回波距离向和方位向采样点数，其解析表达式可表示为：

$I_0(\mathrm{\τ},t;x,y)=\mathrm{\σ}(x,y)\sin c\left\{B_{a}t\right\}\sin c\left\{B\right[\mathrm{\τ}-\frac{2R(t;x,y)}{c}\left]\right\}---\left(1\right)$

其中，σ(x,y)是反射系数，B_a为方位向多普勒带宽，B为发射信号带宽，λ为载波波长，目标平面点目标距离历史为V为平台速度，为天线俯仰角，R₀为方位时间为0时刻天线与点目标的初始斜距，θ₀为0时刻点目标对应的斜视角；

方位时间向量t＝[-PRI·N/2,-PRI·(N/2-1),···,PRI·(N/2-1)]，其中，PRI为脉冲重复间隔；距离时间向量τ＝[-1/f_s·M/2,-1/f_s·(M/2-1),···,1/f_s·(M/2-1)]，其中，f_s为距离向采样率；

对I₀(τ,t；x,y)进行方位向逆傅里叶变换，然后乘以一个方位窗函数W_ac(f_t)，再进行方位向傅里叶变换，得到参考图像I_c可表示为：

$\begin{array}{c}{I}_c(\mathrm{\τ},t;x,y)=\mathrm{FFT}\left\{\mathrm{IFFT}\right\{I_0(\mathrm{\τ},t;x,y)\left\}{W}_{\mathrm{ac}}\left(f_t\right)\right\}\\ =\mathrm{\σ}(x,y)\sin c\left\{B_{a}t\right\}{w}_{\mathrm{ac}}\left(t\right)\sin c\left\{B\right[\mathrm{\τ}-\frac{2R(t;x,y)}{c}\left]\right\}\end{array}---\left(2\right)$

其中，w_ac(t)是W_ac(f_t)的时域表示，W_ac(f)为频域上的矩形窗函数，即 $W_{\mathrm{ac}}\left(f_t\right)=\left\{\begin{array}{cc}1,& f_0-\frac{f_{\mathrm{ac}}}{2}\≤f_t\≤f_0+\frac{f_{\mathrm{ac}}}{2}\\ 0,& \mathrm{otherwise}\end{array}\right.,$ f₀为方位向频率中心，频宽α_c为窗函数对应的窗长；

对参考图像I_c进行边缘检测，得到边缘图像I_ce；

对边缘图像I_ce进行膨胀操作，得到膨胀后的图像I_ced；

对膨胀后结果中的孔隙进行填充，并对图像进行平滑处理，获得图像分割结果I_cs；

步骤S2：对原始图像结果进行方位向加窗：

对原始成像结果I_o进行方位向逆傅里叶变换，然后乘以方位向窗函数W_a(f_t)，再进行方位向傅里叶变换，得到的方位向加窗图像I₁可表示为：

I₁(τ,t；x,y)＝FFT{IFFT{I₀(τ,t；x,y)}W_a(f_t)}  (3)

其中， $W_a\left(f_t\right)=\left\{\begin{array}{cc}1,& f_0-\frac{f}{2}\≤f_t\≤f_0+\frac{f}{2}\\ 0,& \mathrm{otherwise}\end{array}\right.,$ 频宽α为窗函数对应的窗长；

步骤S3、计算加窗后图像的模糊信号比：

对方位加窗后的图像I₁进行边缘检测，得到边缘图像I_oe；

对边缘图像I_oe进行膨胀操作，得到膨胀后的图像I_oed；

对膨胀后结果中的孔隙进行填充，对图像进行平滑处理，获得图像分割结果I_os。

将I₁与对应分割图I_os对应的图像矩阵相乘，得到I₁中目标区域的图像矩阵M_o，将M_o中每一个数值进行平方累加求和得到图像I₁目标区域的总能量

将I₁与参考图像分割图I_cs对应的图像矩阵相乘，得到I₁中不含方位模糊的图像矩阵M_a，对矩阵M_a中的每一个数值平方累加的和就是图像I₁中不含方位模糊的能量 $S_1={\left|\right|M_a\left|\right|}_2^2;$

根据下式计算图像I₁的方位模糊信号比A1：

$A_1=\frac{S_{\mathrm{ambiguity}}}{S-S_{\mathrm{ambiguity}}}=\frac{S-S_1}{S_1}=\frac{{\left|\right|M_0\left|\right|}_2^2-{\left|\right|M_a\left|\right|}_2^2}{{\left|\right|M_a\left|\right|}_2^2}---\left(4\right)$

式中，S_ambiguity是原始图像中方位模糊区域的能量；

判断是否是第一次获得加窗后图像的模糊信号比，如果是，则转入步骤S4；如果此时不是第一次获得，则继续判断图像I₁的模糊信号比A₁是否大于AASR期望值A_e；如果大于AASR期望值A_e，则转入步骤S5；如果I₁的模糊信号比A₁小于或者等于AASR期望值A_e，则输出此时的图像I₁；

步骤S4、设定方位向窗长：

计算AASR差值e＝A₁-A_e，A_e为AASR期望值；

构造权值w，令w＝μ·e·α₀，μ为预先设定的步长因子；

按下式设定窗长α：α＝α₀·w，然后转入步骤S2，将设定的α作为步骤S2中的方位向窗长，按(3)式对原始图像进行加窗处理；

步骤S5、更新方位向窗长：

计算AASR差值e＝A₁-A_e；

更新权值w，令w＝w+μ·e·α₀；

更新窗长值α＝α·w；

然后转入步骤S2，将更新后的α作为步骤2中的方位向窗长对原始图像进行加窗处理。

2.根据权利要求1所述的SAR图像方位模糊抑制方法，其特征在于，步骤S1和步骤S3中所述的边缘检测具体利用Sobel算子进行边缘检测。

3.根据权利要求1所述的SAR图像方位模糊抑制方法，其特征在于，步骤S1和步骤S3中所述的对图像进行平滑处理具体利用菱形结构元素对图像进行平滑处理。

4.根据权利要求1所述的SAR图像方位模糊抑制方法，其特征在于，步骤S1所述的窗函数对应的窗长α_c取方位向采样点数N的五分之一。

5.根据权利要求1所述的SAR图像方位模糊抑制方法，其特征在于，第一次加窗时，步骤S2所述的窗函数对应的窗长α设为α₀＝N。