杜延磊, 高 帆, 刘 涛, 杨 健,*

(1. 中国科学院空天信息创新研究院, 北京 100101; 2. 清华大学电子工程系, 北京 100084;3. 海军工程大学电子工程学院, 湖北 武汉 430033)

0 引 言

海杂波是制约极化合成孔径雷达(synthetic aperture radar,SAR)海上目标检测和识别性能的主要因素之一,其物理机制复杂,且受雷达系统参数、观测几何参数及环境物理参数等诸多因素的影响,从而导致对其特征描述和抑制的难度大[1-2]。特别在中高海况的复杂海面条件下,对海杂波统计特性了解的欠缺成为在该背景下极化SAR海面目标检测、识别和跟踪困难的主要问题。因此,深入了解极化SAR海杂波的统计分布特性有利于改善杂波抑制效果,提升对海面弱小目标检测性能,这对基于极化SAR的海面目标探测预警、跟踪识别、海上救援等军事和民用都具有重要意义。

极化SAR海杂波是随机粗糙海面在各个极化通道对入射雷达电磁波的后向散射回波,因此,对极化SAR海杂波的统计建模和特性分析的前提是获得大量不同极化状态和参数条件下的海面回波数据。目前,获得海杂波的手段主要分为两类:一类是通过岸基或机载平台对海杂波进行实测获得;另一类则主要基于海面形态建模和电磁散射求解的理论方法对海杂波进行数值仿真[3]。前者的优点在于所获得的海杂波数据真实准确,基于此所构建的杂波模型对特定设备和特定应用场景的针对性更强,实用性更好。因此,为了提升各型海用雷达装备的海上目标探测能力,国内外开展了丰富多样的雷达海杂波测量试验,丁昊等对此进行了较为系统全面的介绍[1,4]。然而,海杂波的测量试验成本较高,且受限于试验所用的设备、平台性能以及所在的海域和时机,只能获得特定场景和参数条件下的海杂波数据[1]。而海杂波的数值仿真则具有较好的灵活性和普适性,可以一定程度上弥补实测试验的不足。这类方法首先利用海面形态建模方法生成随机粗糙海面,进而通过数值方法求解其散射方程(低频类方法)[5]或通过引入近似假设获得海面远场散射场的解析近似解(高频类方法)[6-7]。目前,在实际的海杂波研究中,上述两类海杂波获取手段可以实现相互验证、相互补充[1]。

海杂波的统计分布模型是极化SAR海上目标检测的重要先验知识[8-9]。对海杂波统计建模的研究经历了从高斯性到非高斯性、从单一极化到多极化、从纯统计模型到物理/统计结合模型的发展过程[2,10]。较早提出的几种适用于低分辨率单一极化海杂波的典型分布模型包括Rayleigh分布、Log-normal分布、Weibull分布、K分布等[2,11]。而随着SAR技术的发展,早期的杂波分布模型难以准确描述较高分辨率的极化SAR海杂波,尤其是其较长的拖尾表现出显著的非高斯性[2]。目前,复合高斯模型(又称乘积模型或纹理模型)是最为常用的一类极化SAR海杂波统计分布模型,其假设雷达杂波的散射向量可以分解为两部分的乘积,一是描述目标后向散射截面(radar cross section, RCS)空间变化的纹理特征,另一是服从复高斯分布的斑噪向量,并通过选择不同的纹理分布模型,形成各种具体的极化SAR海杂波模型[10]。

海杂波分布模型的参数估计也是其统计建模的关键问题之一。传统的参数估计方法主要包括最大似然估计方法和矩估计方法[11-12]。其中,利用最大似然估计参数时往往需要用到数值方法进行迭代求解,计算效率低[12]。而矩估计方法的特征函数基于Fourier变换获得,在分布模型较为复杂的时候难以甚至无法获得解析形式的特征函数,此外,基于矩的估计方法对数据敏感性也常常较高,参数估计的精度较低[13]。针对SAR海杂波的乘积模型,Nicolas通过利用Mellin变换替换Fourier变换来构建新的特征函数和累积量生成函数,从而在此理论框架下构建Mellin类统计量,并基于此进行海杂波模型参数估计,具有计算简单、精度高、通用性好的特点[14-15]。

目前,对海杂波特性的研究主要利用实测数据的建模分析,而基于具有一定普适性的数值仿真建模研究相对较少,有限的相关研究存在海面仿真尺寸较小、极化信息不足的问题,且对所获得杂波数据的统计建模和分析方法较为陈旧。鉴于此,本文利用成熟的Apel海浪谱模型,结合Monte Carlo方法仿真生成不同风速下的大尺寸二维粗糙海面,进而利用结合了物理光学(physical optics, PO)模型和弹跳射线法(shooting and bouncing ray,SBR)的电磁散射模型对多极化多角度的海杂波进行仿真,最后基于Mellin类统计量对仿真结果进行统计建模和特性分析。

1 极化SAR海杂波数值仿真模型

1.1 二维粗糙海面建模

粗糙海面的空间形态建模是开展极化SAR海杂波数值仿真和分析的基础。目前,海面形态建模最常用的方法是基于海浪谱模型的Monte Carlo方法,其假设随机粗糙海面由无数的谐波叠加而成,从而可以通过对海浪谱先滤波再进行逆Fourier变换获得海面的高度起伏,因此,该方法也被称为线性滤波方法[5,16]。基于此,二维粗糙海面高度f(x,y)可以表达为如下的Fourier级数展开:

(1)

式中:Lx和Ly分别为二维海面的长宽,bmn为高斯随机变量。由于实际的海面仿真中,海面的仿真尺寸和剖分密度是有限的,其决定了对海浪谱滤波的上下限,并具有如下关系[17-18]:

(2)

(3)

此外,根据式所示的海面高度Fourier级数展开特性,可以获得其系数bmn的周期性和共轭对称的特性,即

b(m+M)(n+N)=bmn

(4)

(5)

进而结合粗糙面协方差的定义可以推导获得bmn为[5]

(6)

图1所示为利用Apel谱和Monte Carlo方法生成的风速为6 m/s下的二维粗糙海面。

图1 利用Apel谱和Monte Carlo方法生成的6 m/s风速下二维粗糙海面Fig.1 Two-dimensional rough ocean surface at 6 m/s generated using Apel spectrum and Monte Carlo method

1.2 海面电磁散射计算

基于Monte Carlo方法生成的二维粗糙海面,采用理论的电磁散射方法对海面杂波信号进行仿真。现有常用的海面电磁散射计算方法大致可以分为低频方法和高频方法两大类[11]。低频方法主要是以矩量法[5]、时域有限差分法[22]、有限元法[23]等为代表的数值方法,其通过数值求解粗糙面问题的散射场,具有精度高、适用性广的优点,但由于其计算量巨大,实际中只能在较低的频段使用[24]。而高频方法在粗糙面散射问题的求解中通过引入一定的近似假设,从而获得散射场的解析形式表达,可以在特定应用范围内获得较为准确的结果,且求解速度快。

由于本文主要关注3～5级海况下海面X波段的极化SAR海杂波特性,在此条件下,二维粗糙海面的尺寸和起伏都远远大于入射波波长,高精度的低频数值方法须求解巨量的未知数,对计算资源的消耗巨大,难以满足三维问题下的海杂波仿真需求[24]。因此,本文采用高频方法PO方法和SBR方法以分别考虑粗糙海面的单次散射和多次散射,从而开展极化海杂波的仿真。此外,由于海水具有较大的介电常数,对海面模型采用理想电导体假设[25]。

考虑如图2所示的海面散射场景,入射电磁场Ei、Hi照射到海面一小面元s′,并形成散射场Es、Hs,并定义以海面模型质心为原点的三维坐标系O′X′Y′Z′。根据Maxwell方程组,散射场可由Stratton-Chu方程表达为[7]

图2 海面散射几何示意图Fig.2 Geometric configuration of ocean scattering

(7)

由于Stratton-Chu方程包含面积分项,无法获得精确的解析解,因此需要引入一定的近似假设,对方程进行简化求解。为求解Stratton-Chu方程,PO法引入3个近似假设[7]。① 高频近似:假设入射波波长远小于目标尺寸,从而近似认为入射波只具备粒子性,忽略其在粗糙面棱边及边缘绕射等电磁波波动性,且认为总的散射场仅由粗糙面上被入射波照射到的小面元散射场叠加而来,即,被遮挡的面元对总散射场无贡献。② 远场近似:假设散射场观测点远离目标粗糙面,且二者距离也远大于粗糙面尺寸,即|r|≫|r′|。③ 切平面近似:假设目标表面任意一点及其附近表面曲率半径比波长大得多,可以近似为平面。

基于上述三个近似假设,PO法可给出单站条件下(入射波方向i和散射波方向s相反,即i=-s)粗糙面元s′的后向散射场为

(8)

(9)

式中:面积分项可以利用斯托克斯公式进行计算,得到解析表达式为[26]

(10)

(11)

式中：l1,l2,l3分别表示三角面元s′的3条边矢量。

由式(9)可以发现,后向散射场的极化方向与入射场一致,因此,PO法无法给出交叉极化下的散射系数。鉴于此,本文进一步采用SBR方法来计算粗糙面的多次散射。

(12)

类似的,根据几何关系,可以获得第n次反射的入射波方向为

(13)

式中:I为单位矩阵。基于确定的反射路径和PO法,可以获得单站条件下,最终的SBR多次散射场为[27]

(14)

基于上述针对海面单次散射和多次散射计算获得单站散射场,最终接收天线端的总散射场表示为

(15)

(16)

式中：A表示仿真海面像元内雷达波照射面积；θi为雷达入射角。为了仿真分析极化SAR幅度图的海杂波统计特性[15,28]，本文后续主要采用如下定义的后向散射系数的幅度进行统计建模：

(17)

对于极化SAR系统而言,入射电场和散射电场也可由散射矩阵建立如下关系[8]:

(18)

式中:下标v和h表示水平极化和垂直极化;S为散射矩阵,也可以表达为向量形式为

s=vec(ST)=[Shh,Shv,Svh,Svv]T∈Cd

(19)

式中:vec(·)表示矩阵向量化操作;(·)T表示转置;d为向量维度。

2 基于Mellin类统计量的杂波建模

2.1 极化SAR海杂波的典型乘积模型

基于前文所介绍的海杂波数值仿真方法,可以计算获得成像单元内极化SAR散射向量s,其协方差矩阵定义为

Σ=E{ssH}

(20)

式中:E{·}表示期望;(·)H表示Hermitian(转置共轭)操作。对于视数为L的多视处理数据而言,多视极化协方差矩阵C为[15]

(21)

式中:Xf0表示矩阵正定;C表示复数集。

而随着雷达技术的发展,SAR图像的空间分辨率得到了显著的提升,传统的Gaussian统计模型越来越难以满足对高分辨率、非均匀海杂波的描述,尤其是大量数据表明海杂波分布的长拖尾表现出显著的非高斯特性[2]。因此,诸多学者提出了一类复合Gaussian海杂波统计分布模型,又称为乘积模型或纹理模型。这类模型假设雷达杂波的散射向量可以分解为两部分的乘积,一是描述目标后向散射截面空间变化的纹理特征,另一是服从复高斯分布的斑噪向量,并通过选择不同的纹理分布模型,形成各种具体的极化SAR海杂波模型[10]。

假设SAR海杂波各极化通道的纹理变化特性相同,则分辨单元内的散射向量k可以写作:

(22)

式中:τ为纹理参数,其均值为1;s与前文定义一致,表示斑噪向量。对于多视数据,其极化协方差矩阵表示为

(23)

对于给定的τ,散射向量k和极化协方差矩阵C的概率密度函数(probability density function, PDF)分别为[10]

(24)

(25)

式中:tr(·)表示矩阵的迹；Γd(L)为归一化因子,定义为

(26)

式中:Γ(·)表示gamma函数。

当纹理分布p(τ)已知时,海杂波分辨单元内散射向量k和极化协方差矩阵C的PDF为

(27)

(28)

目前,为了描述不同场景的极化SAR数据的杂波信号,学者们通过引入不同的纹理分布模型提出了多种极化SAR杂波的复合Gaussian统计模型。本文总结了其中较为常用的几种典型杂波分布统计模型,包括K分布、G0分布、Kummer-U分布、W分布、M分布等,表1给出了这几种分布模型纹理参量τ、散射向量k和极化协方差矩阵C的PDF。特别地,对于Wishart分布而言,可以认为其纹理模型为常数,即p(τ)=δ(τ-1)。

表1中,式(32)中的Kν表示ν阶第二类修正Bessel函数,式(34)中U(a,b,z)为第二类合流超几何函数,式(35)中W为Whittaker W函数,其为合流超几何函数U(a,b,z)的函数,式(36)中M为第一类合流超几何函数,亦即Kummer-M函数[10]。

2.2 Mellin类统计量的定义与计算

对于所获得的一组杂波数据,其统计建模一般来说包括两个部分,一是确定其分布模型,二是对所确定分布模型的参数进行估计。在极化SAR海杂波研究中,在经过早期利用简单的Rayleigh分布、Log-normal分布、Weibull分布等纯统计模型的发展之后,越来越多的研究认为中低分辨率条件下,极化SAR海杂波可以用K分布来描述,并基于此发展了大量K分布参数估计方法。其中,传统的参数估计方法主要包括最大似然估计方法和矩估计方法[11-12]。而利用最大似然估计参数时往往需要使用数值方法进行迭代求解,计算效率低。因此,基于矩的参数估计方法得到了广泛地使用。

矩估计方法的数学理论基础是:对于给定随机变量的分布函数,其特征函数可以表示为关于矩的级数展开形式,当该分布所有矩存在且有限,且其特征函数的幂级数展开在靠近原点处收敛,则该分布函数可由其特征函数唯一表示[15]。然而,矩估计方法的特征函数是基于Fourier变换获得,在分布模型较为复杂的时候难以甚至无法获得解析形式的特征函数,此外,基于矩的估计方法对数据敏感性也常常较高,参数估计的精度较低[28]。因此,本文采用基于Mellin类统计量的方法来估计杂波分布模型的参数。

不同于传统参数估计方法的特征函数基于Fourier变换,Mellin类统计量的定义与计算采用Mellin变换。对于复矩阵元函数f(Z),其Mellin变换定义为[15]

(29)

基于Mellin变换,矩阵变量Z的Mellin类特征函数(characteristic function, CF)φZ和Mellin类累积量生成函数(cumulant generation function, CGF)φZ定义为

φZ(χ)=E{Zχ-d}=M{pZ(Z)}(χ)=

φZ(χ)=ln[φZ(χ)]

(30)

进一步地,根据矩量和累积量的定义,Z的v阶Mellin类矩μv(Z)和v阶Mellin类累积量υv(Z)分别定义为

(31)

由(31)式可以发现μv(Z)=E{lnZ},因此Mellin类矩也被称为对数矩,类似的,Mellin类累积量也被称为对数累积量。当Z的所有对数矩和对数累积量存在时,有

(37)

对数矩量和对数累积量之间可以由Faadi Bruno公式的组合形式表示为[29]

(38)

从而可获得Z的前三阶Mellin类累积量为

(39)

2.3 基于Mellin类统计量的杂波模型参数估计

为了使用Mellin类统计量对海杂波模型参数进行估计,我们需要推到获得二者之间的严格关系。假设乘积模型的极化协方差矩阵写作:

(40)

(41)

(χ-d)ln|Σ|-(χ-d)dlnL

(42)

(43)

进一步地,对于纹理分量T,其Mellin类特征函数为

(44)

结合Mellin统计量定义和Faá di Bruno公式可得:

(45)

由于随机矩阵变量X=UV的Mellin类特征函数和统计量有如下性质[14]:

(46)

因此,乘积模型极化协方差矩阵C的Mellin类累积量为

(47)

代入式(43)最终有:

(48)

式(48)给出了极化SAR杂波的Mellin类累积量与乘积模型参数之间的关系,基于此,我们可以利用式(38)计算杂波数据的若干阶累积量,从而可以对杂波分布模型参数进行估算。对于式(48)中纹理模型的Mellin类统计量,表2给出了K分布、G0分布、Kummer-U分布、W分布、M分布等乘积模型的纹理分布对数累积量。

表2 典型纹理分布模型的对数累积量Table 2 Log-cumulants of typical texture distribution models

此外,在理论上,这5种典型乘积模型可以严格地划分二阶和三阶Mellin类累积量所形成的平面空间,如图3所示[9]。因此,利用样本数据的Mellin类统计量还可以判断杂波数据所服从的分布。

图3 典型杂波乘积模型在二阶/三阶Mellin类累积量空间的分布Fig.3 Typical product models in the space of 2nd/3rd Mellin-kind cumulants

3 仿真结果与分析

利用前文介绍的海面形态建模和雷达散射计算高频方法,我们对海面X波段全极化海杂波进行数值仿真,并基于Mellin类统计量方法开展海杂波统计建模和特性分析。

由于本文旨在对3～5级海况下海面X波段(9.6 GHz)高分辨率极化海杂波进行仿真建模和统计分析,因此,我们首先基于Apel海浪谱和Monte Carlo方法分别生成了6 m/s、9 m/s和12 m/s风速下的二维粗糙海面,如图4(a)所示。所生成海面的尺寸为128 m×128 m (4 096λ×4 096λ,λ为入射波波长),其中,为了仿真海面大尺度海浪,生成真实的中高海况海面,在海面仿真的过程中,本文使用了包含谱峰的海浪能量谱,这使得生成的海面包含了所有的大尺度风浪,同时也使得仿真海面尺寸范围较大,远远大于现有的低频数值仿真方法可以实现的大小。尤其在仿真12 m/s风速下海面时,为了包含海浪谱的谱峰,则生成的海面尺寸达到256 m×256 m。在后续的海面电磁散射计算中,为了与其他两个风速海面具有可比性,也为了减少运算量,我们对12 m/s风速海面进行裁剪,获得尺寸为128 m×128 m的海面模型。由图4(a)可以看出,随着风速增大,海面的粗糙起伏也逐渐变大,在本文的案例中,所生成的6 m/s、9 m/s和12 m/s风速下的海面均方根高度分别为6.13λ、12.80λ和19.88λ。此外,图4(a)中,风向由y轴正坐标指向负坐标。为了保证对海面刻画以及电磁散射仿真的精度,并考虑到海面小尺度粗糙度对雷达后向散射的影响,本文所仿真海面的剖分密度达到4个采样点/λ。

图4 不同风速下的二维粗糙海面和对应的雷达后向散射系数Fig.4 Rough ocean surfaces at various wind speeds and corresponding radar backscattering coefficients

基于所生成的不同风速二维粗糙海面,本文利用所建立的海面电磁散射高频仿真方法开展不同雷达入射角和方位角条件下X波段全极化海杂波数值仿真计算。在实际的计算中,根据SAR的成像空间分辨率概念,我们将海面剖分为64×64个大小为2 m×2 m的成像单元。图4(b)所示为入射角为20°、雷达方位角为逆风方向条件下,本文所计算的不同风速下二维海面X波段的VV极化雷达后向散射系数。由图可以发现,随着海面风速的增大,海面的雷达后向散射系数也增大,这与海面粗糙度增大,海面雷达散射增强的基本常识相符合。此外,考虑到雷达观测方向为逆风方向(沿y轴负坐标指向正坐标),由图4(a)和图4(b)可以发现,海面迎向入射波的面元后向散射系数较大,这与粗糙面雷达散射特性相一致。

3.1 仿真海杂波的后向散射系数与分布验证

在开展海杂波的建模和分析之前,我们需要对仿真海杂波的后向散射系数和统计分布准确性进行验证。本文选择X波段广泛使用的并由实测数据拟合获得的地球物理模式函数XMOD-2对仿真海杂波的同极化后向散射系数进行验证[30]。对于交叉极化,受限于目前较少的X波段交叉极化实测数据,我们采用Voronovich和Zavorotny在文献中给出的二阶小斜率近似模型(SSA-2)结果和测量数据对仿真海杂波的交叉极化后向散射系数进行验证[31-32]。图5所示为仿真海杂波后向散射系数与XMOD-2、SSA-2模型结果及实测数据的对比。

图5 利用其他模型和实测数据对仿真海杂波后向散射系数精度进行验证Fig.5 Accuracy validation of simulated sea clutter using other models and measurement

由图5可以看出,在同极化条件下,仿真海杂波的后向散射系数与XMOD-2模型有较好的吻合,最大误差不超过2 dB,这表明本文采用的海面雷达散射高频方法在同极化条件下对X波段海杂波具有较高的仿真精度。而对于交叉极化,尽管在高风速情况下仿真海杂波后向散射系数与SSA-2模型结果较为接近,但在低风速情况下,仿真结果则与参考数据有相对较大的误差,最大可达约10 dB。考虑到海面交叉极化后向散射能量远远低于同极化,且本文主要关注中高海况下海杂波的统计分布及定性的特性分析,这一误差大小对于交叉极化仿真是可以接受的。

为了验证数值仿真所得到海杂波的分布特性准确性,我们利用实测的海杂波分布与仿真海杂波分布进行对比。本文采用由岸基X波段固态功放监视/导航雷达于2021年1月6日15时6分利用凝视模式获得的海面杂波数据,该数据由海军航空大学的“雷达对海探测数据共享计划”发布,并包含了同步测量的气象水文等辅助资料。获得实测数据的雷达中心频率为9.4 GHz,极化方式为HH极化,距离分辨率为6 m,入射角为80°,雷达观测方位角为55.05°,其他的雷达系统参数在文献[3,33]具体给出。此外,根据与同步测量的海洋气象水文数据进行时空匹配,可知所测海杂波数据对应海面的风速为9.57 m/s,风向为61.73°,因此,考虑到雷达观测方位角,雷达观测方位相对于风向约为6.68°。

由于实测海杂波数据是未经辐射定标的散射强度结果,与仿真数据雷达后向散射系数无法进行直接对比,因此,本文对仿真和实测海杂波数据的进行了归一化处理,并基于归一化后的海杂波数据计算了其前三阶对数累积量,图6(a)所示为实测和仿真海杂波数据的二阶/三阶对数累积量对比。由图中散点分布可以看出实测和仿真海杂波具有非常相近的分布特征,且较为服从Wishart分布,这表明在此条件下,所获得的海杂波纹理特性并不显著,考虑到本案例较大的雷达入射角,这一特性是较为合理的。

图6 仿真海杂波分布与实测海杂波分布的对比Fig.6 Comparison of distributions of simulated sea clutter and measured sea clutter

基于图6(a)结果,我们采用K分布模型对归一化实测和仿真海杂波的雷达后向散射系数和散射场幅度的概率密度分布进行拟合,并利用所计算的二阶/三阶对数累积量对模型参数进行估计,最终仿真和实测海杂波分布的对比如图6(b)所示。可以看出,仿真海杂波与实测海杂波的分布具有良好的匹配,证明了利用前文的海面杂波高频数值仿真方法对海杂波分布仿真的准确性,以及利用Mellin类统计量开展海杂波统计分布和特性分析的可行性。

在本文的海杂波数值仿真中,尽管我们通过采用包含海浪谱谱峰考虑了所有大尺度风浪的贡献,但是小尺度波浪的贡献则由于剖分密度的限制进行了一定的截断。由式可知,4个采样点/λ的剖分密度表明海面波数大于402 rad/m的小尺度波浪在仿真海面中被舍去。因此,对于本文采用的电磁散射方法和海面剖分密度是否足以描述海浪细节对海杂波的贡献,尤其在5级高海况条件下,需要进一步进行验证。图7给出了5级海况下,分别利用4个采样点/λ和8个采样点/λ仿真海面获得的海杂波的二阶/三阶对数累积量及分布对比。特别地,这里海杂波的仿真像元大小为1 m×1 m。由图7(a)可以看出,采用不同剖分密度海面获得的极化SAR仿真海杂波具有相近的二阶/三阶对数累积量,且两组海杂波的后向散射系数平均差异在同极化下可以忽略,在交叉极化下仅约0.72 dB。在图7(b)中,由不同剖分密度海面获得的海杂波在同极化和交叉极化都分别具有几乎一致的分布。这表明海面较小尺度的波浪对海面雷达散射的贡献较小,这也与我们前期的研究结论相符合[17],本文采用4个采样点/λ的海面剖分密度可以满足海杂波仿真中对海浪细节的描述。

图7 不同剖分密度海面获得的仿真海杂波分布对比Fig.7 Comparison of distributions of simulated sea clutters from surfaces with different discretizing density

3.2 X波段高分辨率极化SAR海杂波建模与分析

在海杂波数值仿真的基础上,我们基于Mellin类统计量构建海杂波的分布模型,并分析了雷达入射角、海面风速风向对海杂波分布的影响。图8所示为不同风速下X波段同极化海杂波的分布,其中,仿真海杂波的雷达入射角为20°,风向为逆风方向。由图8(a)首先可以看出,VV极化和HH极化海杂波的二阶/三阶对数累积量非常接近,这主要是因为雷达入射角较小,VV和HH极化的雷达后向散射系数较为接近,表明VV和HH极化海杂波有着相近的分布,尤其是二者的纹理分布接近,因此在图8(b)只展示了VV极化结果。进一步地,我们可以发现高分辨率的同极化X波段海杂波并不服从传统的适用于中低分辨率海杂波的K分布模型,这说明在此条件下,海杂波的纹理特征不服从gamma分布。而在不同的风速下,X波段同极化海杂波表现出不同的纹理分布,并最终使得其服从不同的分布模型。具体来说,6 m/s和9 m/s风速下海杂波分别服从Kummer-U分布和M分布,但二者都比较接近G0分布,在12 m/s风速下,X波段同极化海杂波服从M分布。基于海杂波所服从的分布模型,我们利用所计算的对数累积量对相应的模型参数进行了估计,并在图8(b)中对比了不同风速下海杂波幅度的概率密度分布。由图8(b)可以看出,利用图8(a)所指示的分布模型可以很好地拟合海杂波数据。随着风速的增大,海杂波分布中后向散射系数幅度较大部分占比逐渐升高,海杂波分布由陡逐渐变缓,这与物理事实相符合。

图8 不同风速下X波段同极化海杂波的分布Fig.8 Distribution of X-band co-pol sea clutters at various wind speeds

图9展示了不同风速下X波段交叉极化海杂波的分布,相关的海杂波仿真参数与图8相同。类似于同极化,HV极化和VH极化海杂波的二阶、三阶对数累积量也非常接近,考虑到极化互易性原理,这是合理的,因此,在此处及后续的交叉极化海杂波分布模型中,我们仅展示HV极化。而与同极化不同的是,交叉极化海杂波的二阶、三阶对数累积量在不同风速下都集中在Wishart分布附近,表明交叉极化条件下高分辨率X波段海杂波没有显著的纹理特征,因此,在高分辨率条件下,极化SAR海杂波在四个极化通道不适宜使用相同的纹理模型来描述,这与传统针对中低分辨率海杂波在各极化通道采用相同纹理模型不一致。在图9(b)中我们利用K分布拟合了不同风速下交叉极化的海杂波分布,与同极化类似,随着风速的增大,海杂波分布由陡逐渐变缓。

图9 不同风速下X波段交叉极化海杂波的分布Fig.9 Distribution of X-band cross-pol sea clutters at various wind speeds

图10所示为不同入射角下X波段同极化海杂波的分布,其中,仿真海杂波的海面风速为9 m/s,风向为逆风方向。由图10(a)可以看出,雷达入射角对高分辨率X波段同极化海杂波的分布有着显著的影响,在入射角较小时,同极化海杂波有着显著的纹理特征,且VV和HH极化的二阶/三阶对数累积量非常接近。而在雷达入射角较大时,海杂波的纹理特征显著减弱,且二阶/三阶对数累积量在VV和HH极化间显示出差异。在图10(b)中我们分别利用M分布和K分布拟合雷达入射角为20°、40°、60°时的VV极化海杂波分布,可以看出所用模型可以很好地拟合海杂波分布,由于海面同极化后向散射能量随着入射角增大而减小,图中的分布模型也随着入射角增大而逐渐变陡。

图10 不同入射角下X波段同极化海杂波的分布Fig.10 Distribution of X-band co-pol sea clutters at various incidence angles

图11所示为不同风向下X波段同极化海杂波的分布,其中,仿真海杂波的雷达入射角为20°,海面风速为9 m/s。由图可以看出,该条件下海杂波的二阶/三阶对数累积量集中在G0分布附近,并在不同的风向下稍有区别。在图11(b)中,我们分别采用M分布、G0分布和Kummer-U分布对风向为逆风、侧风和45°时的海杂波分布进行拟合,可以发现海杂波分布随着风向由侧风向逆风变化逐渐变缓。

图11 不同风向下X波段同极化海杂波的分布Fig.11 Distribution of X-band co-pol sea clutters in various wind directions

图12展示了雷达视数对X波段同极化海杂波分布的影响,其中,仿真海杂波的雷达入射角为20°,海面风速为9 m/s,风向为侧风方向。由图12(a)可以看出,随着雷达视数的增大,海杂波的纹理特征逐渐减弱,考虑到雷达多视处理的基本原理及其空间滤波的效果,海杂波的纹理特征被平滑是合理的。进一步地,由图12(b)也可以看出,多视处理使得海杂波的极大值和极小值减少,海杂波的分布逐渐向中间移动。

图12 不同雷达视数条件下X波段同极化海杂波的分布Fig.12 Distribution of X-band co-pol sea clutters under conditions with various radar looks

4 结 论

中高海况下极化SAR海杂波的统计建模是海面舰船、导弹及相关飞行器检测预警的关键和难点问题。本文针对3～5级海况下X波段极化SAR海杂波开展数值仿真、统计建模与分析研究。首先利用成熟的Apel海浪谱模型和Monte Carlo方法建立了6 m/s、9 m/s和12 m/s风速下二维粗糙海面模型,进而采用高频方法开展不同风速、风向、雷达入射角、极化方式下的X波段海杂波仿真计算。为了获得准确的海面极化散射系数,我们采用PO方法和SBR方法以分别考虑粗糙海面的单次散射和多次散射。基于仿真海杂波数据,利用典型的乘积模型和Mellin类统计量开展高分辨率X波段极化SAR海杂波统计建模和特性分析。仿真结果表明:在2 m×2 m的高分辨率条件下,X波段同极化海杂波并不服从适用于中低分辨率海杂波的K分布,且在不同风速风向下呈现出不同的纹理特征;X波段交叉极化海杂波表现出较弱的纹理特征,服从Wishart分布。因此,X波段极化SAR海杂波在同极化和交叉极化条件下具有不同的纹理分布特征,这与传统针对中低分辨率海杂波在各极化通道采用相同纹理模型不一致。此外,随着入射角和雷达视数的增大,海杂波的纹理特征逐渐减弱。

尽管本文基于高频仿真方法和Mellin类统计量对X波段极化海杂波进行仿真建模,给出了海面风速、风向、雷达入射角等雷达系统参数和环境参数对极化海杂波影响的初步分析结果,但受限于模型在交叉极化的仿真精度,相关的分析仅在定性层面展开,未给出定量化的结果。在后续工作中,可进一步构建精细化的海面模型,并考虑海浪尖峰的绕射等现象,提高模型在交叉极化的仿真精度,并在本文工作的基础上增加仿真数据,建立极化SAR海杂波统计模型与雷达系统参数和海洋环境参数之间的经验模式函数,服务海面极化SAR海面目标检测应用。