基于变迹滤波的SAR图像旁瓣抑制技术

徐光耀, 樊霞良

(解放军91336部队, 河北 秦皇岛　066326)

在SAR/ISAR成像过程中,雷达传感器所获得的相位历史数据经过运动补偿和距离徙动校正之后,目标的散射中心分布和SAR图像就成为了一对傅里叶变换对(FT)。直接采用基于傅里叶变换的算法进行成像,所得到的图像具有较高的旁瓣。常用的SAR图像旁瓣抑制方法是通过数据域加窗来实现的,加窗处理的优点是几乎没有增加成像算法的复杂度和运行时间,然而加窗是对有限序列的截断,会产生吉布斯效应(数学上,即sinc函数插值),这会导致图像主瓣的展宽,从而降低图像的分辨率[1]。

空间变迹滤波技术[2-5](Spatially Variant Apodization,SVA),可以看作多变迹滤波的极限形式,其通过在图像域对每个像素点寻求最优的窗函数权值,使得输出的图像具有最窄的主瓣(sinc函数主瓣)和最低的旁瓣。需要特别注意的是SVA算法要求图像的点扩展函数(Point Spread Function,PSF)必须满足sinc核函数假设,这在实际中基本都是满足的。

1　变迹滤波算法

1.1　双变迹(DA)和复数双变迹(CDA)算法

经过加窗处理后,副瓣会得到较好的抑制,但会导致主瓣的展宽(约2倍)。常见的基于余弦底座的窗函数有汉宁窗和汉明窗,汉宁窗的第一副瓣约为-31.5dB,汉明窗的第一副瓣约为-43.8dB。

通过一种非线性双变迹(DA)处理所得到的图像,既具有较好的主瓣分辨率,又具有较低的副瓣。具体的处理过程如下:

Step1:分别计算两种不同窗函数得到脉冲响应函数,一种是不加窗,一种是加窗,比如汉宁窗;

Step2:在每个像素点处,选择模值最小者作为输出的像素值。

实际上,将上述的窗类型扩展到三个,便可以得到三变迹。双变迹和三变迹的脉冲响应函数如图1所示。

图1　不同变迹处理的脉冲响应

从图1可以看出,随着变迹次数(窗函数个数)的增多图像的副瓣,特别是第二副瓣以后的副瓣有明显的降低,两次变迹第二副瓣降到-31dB,三次变迹第二副瓣已经降到-45dB以下。

由于DA算法实际上只应用了图像的模值信息,而蕴含在复图像中的相位信息却被忽略了。因此为了充分利用复图像的实部(I路信号或者同相信号)以及虚部(Q路信号或者正交信号)(其对应的窗函数的脉冲响应的实部和虚部如图2所示),可以采用一种CDA算法[10]。

图2　不同窗函数的实部和虚部

通过以上的分析,可以推测随着变迹次数的进一步增多(比如每个像素点对应不同窗函数),那么就可以在保证主瓣宽度的基础上使得副瓣最低——这就是空间变迹算法(SVA)的核心思想。

1.2　1-D SVA算法

对于一个由很多复杂的散射体构成的目标而言,DA和CDA算法都可以获得比较好的视觉效果。在不考虑运算量的情况下,三次、四次或者更高阶的变迹处理,可以进一步提高副瓣的抑制效果。SVA算法可以保证主瓣具有较好的分辨率,同时在不改变杂波分布强度的基础上几乎完全消除了副瓣。

具有余弦底座的窗函数系可表示为

A(n)=1+2wcos(2πn/N)0≤w≤0.5

(1)

式中,w=0时,表示矩形窗;当w=0.43时表示汉明窗;当w=0.5时表示汉宁窗。对式(1)进行快速傅里叶变换便可得到奈奎斯特采样的脉冲响应:

a(m)=δm,0-w(δm,-1+δm,1)

(2)

式中,δm,n表示Kronecker delta函数:

(3)

从式(3)可以看出,时域加窗相当于在图像域的三点循环卷积。假设未加窗时的一维距离像可表示为

g(m)=I(m)+iQ(m)

(4)

那么,利用式(2)的三点卷积函数得到加了余弦底座窗后的图像变为

g′(m)=g(m)-w(m)[g(m-1)+g(m+1)],

0≤w(m)≤0.5

(5)

按照DA和CDA算法的同样的处理模式,SVA算法可以归结为如下的约束优化问题:

(6)

对于上述问题的求解可以从两个不同的方面来考虑,一是I、Q联合SVA算法,另一种是I、Q分离SVA算法,共有四种不同的组合求解,这里只给出I、Q分离SVA算法的原理。

I、Q分离SVA算法是分别对未加窗图像的I路(实部)和Q路(虚部)进行单独处理,以I路信号为例,令式(5)为0，可得

(7)

将式(7)直接代入式(5)可得:

g′(m)=

(8)

从以上推导可知,I、Q分离SVA实际上是分别使I2和Q2的输出最小。如果令y=0.5[g(m-1)+g(m+1)],那么I、Q分离SVA算法可以重新表示为

(9)

以上的分析是建立在整数倍奈奎斯特采样的基础上,如果是两倍、三倍或者更大的整数k倍时,SVA算法只需要通过将式(9)中的像素点m、m-1、m+1,变为m、m-k、m+k即可。

1.3　2-D SVA算法

对于类似点目标而言,其在二维图像域表现为二维的sinc函数,那么2-D SVA算法就可以通过两种途径来实现,一是径向和横向分别作1-D SVA,二是采用直接的2-D SVA算法,下面介绍下直接2-D SVA算法的原理。

直接2-D SVA算法选取当前像素周围的8个像素进行加权,具体的窗函数可表示为

(10)

根据式(10)中的wm和wn是否相等以及I、Q的联合和分离可以有四种不同的组合方式,在此只对比较常用的wm≠wm、I、Q分离直接2-D SVA算法进行推导。

令

Qm=g(m-1,n)+g(m+1,n)

(11)

Qn=g(m,n-1)+g(m,n+1)

(12)

P=g(m-1,n-1)+g(m+1,n+1)+g(m-1,n+1)+g(m+1,n-1)

(13)

则有

g′(m,n)=g(m,n)-wmwnP-wmQm-wnQn

(14)

由式(14)可以看出,当分别固定wm或者wn时,g′(m,n)与wm或wn均呈线性关系,因此极值处位于wm和wn的两端处。那么二维I、Q分离SVA算法可描述为:

Step 1:分别计算(wm,wn)=(0,0.5),(0.5,0),(0.5,0.5)处的g′(m,n)值,如果三者中任一个与g(m,n)的符号相反,则令g′(m,n)=0;

Step 2:计算(wm,wn)=(0,0),(0,0.5),(0.5,0),(0.5,0.5)四个g′(m,n)值,取模值最小者作为g′(m,n)。

1.4　非整型奈奎斯特采样SVA算法

上述对于SVA算法的讨论都是基于整数倍的奈奎斯特采样,当不满足整数倍的奈奎斯特采样时,频域加权函数可表示为:

(15)

式中,fs表示采样频率,a、α是SVA算法的两个参数。对于式(15)的求解可归结为给定α(m)时的约束最小化问题,那么非整型奈奎斯特采样的SVA算法的结果可表示为

gs(m)=ag(m)+α(m)[g(m-M)+g(m-M)]

(16)

2　仿真及结果分析

2.1　一维变迹滤波处理结果

图3(a)～(d)给出了一维变迹滤波的结果。图3(a)是FFT的结果,可以看出图像的旁瓣比较高;图3(b)是加汉宁窗(Hanning)FFT的结果,虽然旁瓣得到了抑制,但是分辨率有所下降,特别是前两根谱线已经难以区分;图3(c)是CDA的结果,可以看出主瓣分辨率没有降低,旁瓣得到了较好的抑制;图3(d)是SVA处理结果,具有和图3(c)类似的结论,相对于加窗处理,这两种方法都在主瓣不展宽的基础上有效抑制了旁瓣,从而提高了谱线分辨率(大约3Hz)。

图3　一维变迹滤波处理结果

2.2　二维变迹滤波处理结果

为了进一步分析非整型奈奎斯特采样情况下变迹滤波对高分辨率图像的处理效果,这里采用某型船只的二维实测成像数据进行处理,奈奎斯特采样取为2.5倍。处理结果如图4所示。

图4　舰船目标的二维变迹滤波结果

图4(a)是FFT结果,具有很强的旁瓣;图4(b)加窗(Hanning窗)FFT的结果,旁瓣得到了抑制,主瓣展宽比较明显;图4(c)是CDA的结果,主瓣基本没有展宽,旁瓣得到了较好的抑制;图4(d)是I、Q分离SVA的处理结果,可以看出主瓣没有任何展宽,而且旁瓣得到了较好的抑制,而且SVA的结果优于CDA的结果。图5对不同算法的距离向切片和方位向切片进行了对比。

图5　舰船目标的两维切片

图5分别给出了图像中心处的方位向切片和距离向切片,由图5可知SVA算法所得的结果具有最低的旁瓣和最窄的主瓣,这与上述对于二维图像的分析是一致的。

3　结束语

本文针对常规加窗方法对SAR图像旁瓣进行抑制会导致主瓣展宽,图像分辨率和动态范围均会降低的问题,在采用变迹滤波方法进行旁瓣抑制的基础上,对非整数倍奈奎斯特采样下的SVA算法进行了改进,一维仿真结果以及二维实测数据处理结果表明,本文的方法能够保证主瓣不展宽的情况下有效抑制旁瓣,提高图像的动态范围。

参考文献:

[1]许小剑, 黄培康. 雷达系统及其信息处理[M]. 北京: 电子工业出版社, 2010: 10-46.

[2]Ni C, Wangle Y F, Xu X H, et al. A SAR sidelobe suppression algorithm based on modified spatially variant apodization[J]. Sci China Tech Sci, 2010, 53:2542-2551.

[3]LIU Y Z, XU X J, Azimuth Sidelobe Suppression Technique for Near Field MIMO Radar Imaging[C]∥Proc of SPIE on Remote Sensing XXI: Image and Signal Processing. Lorenzo Bruzzone: SPIE, 2015, 96431E:1-8.

[4]R. Goodman, S. Tummala, and W. Carrara, Issues in Ultra-wideband, Widebeam SAR Image Formation[C]∥Proc of IEEE International Radar Conference, Alexandria, VA, May 1995:479-485.

[5]Vu V T, Sjögren T K, and Pettersson M I, On Apodization Techniques for Ultra-Wideband SAR Imaging[C]. Proc of EURAD 2009,Rome, Italy, Sep. 2009:529-532.

[6]Stoica P, Moses R. Introduction to Spectral Analysis[M]. Upper Saddle River, NJ: Prentice Hall[M], 1997:23-65.

[7]XU X J, Narayanan R M. Enhanced Resolution in SAR/ISAR imaging Using Iterative Sidelobe Apodization[J]. IEEE Transactions on Image Processing, 2005, 14(4):537-546.

[8]XU X J, Narayanan R M. SAR Image Enhancement Using Noninteger Nyquist SVA Technique[J]. IEEE Letters, 2002: 298-301.

[9]Stankwitz H C, Kosek M R. Super-resolution for SAR/ISAR RCS measurement using spatially variant apodization (super-SVA)[C]∥Proc of 17th Meeting Symp. Antenna Measurement Techniques Association, Williamsburg, VA, 1995:251-256.

[10] Stankwitz H C, Dallaire R J, Fienup J R. Nonlinear Apodization for Sidelobe Control in SAR Imagery[J]. IEEE Transactions on Aerospace and Electronic System, 1995, 31(1):267-278.

[11] Smith B H. Generalization of Spatially Variant Apodization to Noninteger Nyquist Sampling Rates[J]. IEEE Transactions on Image Processing, 2000, 19(6):1088-1093.

[12] Sparse Aperture Fill for SAR Using Super-SVA[C].in Proc. 1996 IEEE Radar Conf., Ann Arbor, MI, May 1996:70-75.