SIRP法相干相关K分布雷达杂波的建模与仿真

SIRP 法相干相关K 分布雷达杂波的建模与仿真
gjj_hit@
所谓杂波仿真，实际上就是要生成一系列在幅度上服从特定的概率密度分布（pdf ）的相关随机序列，常见的杂波仿真方法有两种：零记忆非线性变换法（ZMNL ）和 球不变随机过程法（SIRP ）。

ZMNL 方法的基本思想是：首先产生相关的高斯随机过程，然后经过某种非线性变换得到所求的相关随机序列。

这种方法的缺点就是输入序列与输出序列间有复杂的非线性关系,因此必须寻找输入序列与输出序列的相关函数间的非线性对应关系。

SIRP 方法的基本思想是：产生一个相关的高斯随机过程，然后用具有所要求的单点概率密度函数的随机序列进行调制。

这种方法的缺点则是受所求的序列的阶数及自相关函数的限制,同时这种方法的计算量非常大,不易形成快速算法。

ISAR 是一种相干雷达，其海杂波必然是相干且时空相关的。

对于相干相关杂波，以往的方法都是将非相干的ZMNL 方法加以推广得到相干的ZMNL 模型。

这种方法得以应用的一个前提是已知非线性变换前后杂波相关系数的非线性关系，然而对于相干相关K 分布杂波却很难找到这样一种非线性变换，于是我们采取SIRP 方法来仿真ISAR 的海杂波。

K 分布适用于描述高分辨雷达的非均匀杂波，多用于对海杂波的模拟。

K 分布可以由一个均值是慢变化的瑞利分布来表示，其中这个慢变化的均值服从Γ分布。

K 分布的概率密度函数为：
()()()12;,K /,
(0,0)2x f x x x ν
νανανανα-⎛⎫
=∙∙>> ⎪Γ⎝⎭
(1)
其中，ν是形状参数，α是尺度函数，()Γ 是伽马函数，K ν是第二类修正贝赛尔函数。

杂波平均功率2σ，ν和α之间的关系可表示为：
2
2
2σαν
= (2)
对于大多数杂波来说，形状参数的取值范围是0ν<<∞，对于较小的ν的取值，如0.1ν→时，杂波有较长的托尾，ν→∞时的分布接近于瑞利分布。

图1给出了K 分布杂波序列的实现结构。

图1 相干相关K 分布杂波SIRP 方法
图中，1()w k 为一复高斯白噪声，线性滤波器1()H z 由()x k 的相关函数设计决定，2()w k 为一与1()w k 相互独立的实高斯噪声，线性滤波器2()H z 必须使得输出的高斯序列具有高度的相关性(相关函数接近于1)，ZMNL 变换使得输出的
()s k 的概率密度函数（pdf ）为杂波的特征pdf 。

对于K 分布来说，()s k 服从广
义χ分布，该分布的定义如下：
()()
()21
22exp ,
0X x f x x x ννννν-=-≥Γ
(3)
要用图1所示的模型产生K 分布杂波，需要产生符合广义K 分布的()s k 并设计线型滤波器1和线型滤波器2。

滤波器1的设计比较简单，它使输出()y k 具有所要产生杂波的功率谱，设计方法同ZMNL 法的滤波器设计。

由于我们对()s k 的相关函数不感兴趣，因此，可将滤波器2设计为一带宽很窄的低通滤波器，使得非线性变换随机序列的功率谱足够窄。

下面以一例MATLAB 仿真说明上述产生K 分布杂波的过程 [1]。

具体程序见“Matlab 程序”文件夹K_distribution.m 和nonline_eq_sirp.m 。

例：产生杂波的幅度概率密度函数的参数为 2.0ν=，0.5α=，功率谱密度为高斯谱，其3dB 带宽为40Hz 的K 分布杂波，滤波器1的设计采用傅立叶级数展开法，模拟的杂波的功率谱密度采用Burg 法估计得到。

一、 复高斯白噪声的产生
二、 滤波器1()H z 的设计——傅立叶级数展开法[2]
这种方法是通过将所希望的网络的频率特性展成傅立叶级数的方法求滤波器加权系数的，故称这种方法为傅立叶级数展开法。

众所周知，非递归滤波器可由以下差分方程来描述：
0,
0N
n i n i i y a x n N -==≤≤∑ (4)
式中：n i x -表示滤波器的第n i -个输入；n y 表示滤波器的第n 个输出；i a 为滤波器加权系数。

滤波器的传递函数可通过Z 变换求出：
()0
N
i i i H z a z -==∑
(5)
频率响应为：
0()N
jw
jwi i i H e a e -==∑
(6)
（6）式为数字滤波器的频率响应，令s w T =Ω，将（6）转化为模拟滤波器的频率响应：
()20
()s
s N
j T j fT i i i H f H e
a e πΩ-===∑
(7)
其中s T 为将模拟滤波器转化成数字滤波器时的抽样间隔（这里的s T 的单位也是频率单位，因为是在频域抽样），1/s s F T =为其抽样频率，s F 为模拟滤波器的频域周期。

又已知，杂波归一化的高斯谱密度为：
2
2
()exp 4f f S f σ⎛⎫
=- ⎪ ⎪⎝⎭
(8)
希望在输入白噪声时，有：
2
()()S f H f =
(9)
显然，所设计滤波器应有高斯响应：
2
2()exp()4f
f H f σ=-
(10)
将其展成傅立叶级数1：
()2s
N
j fnT n
n N
H f A e
π-=-=
∑ (11)
又由于()H f 为偶函数，所以：
()01
/2cos(2)N
n s n H f C C fnT π==+∑
(12)
其中：
2n n C A =
(13)
对式(7)取绝对值2，根据谱的偶函数特性知，式(12)中的n C 便等于式(7)中的i a ，即非递归滤波器频率响应的傅立叶级数展开式的系数，就是该滤波器的加权系数。

由于频率响应是给定的，于是使问题简单了。

为了求系数n C ，改变变量，将()()H f H t →的傅立叶变换写成：
2()()j ft F f H t e dt π∞
--∞
=⎰
(14)
将式(10)代入，得：
222
4()2f f
f F f σπσ-= (15)
当n 有限时，傅立叶级数的系数3
1 这里要知道
()H f 是以Fs 为周期的频谱函数。

2 对比式(11).
3
这里应用了周期信号的傅立叶变换与周期信号的傅立叶级数展开系数之间的关系[3]。

()2222
42f s
n T n s s f A T F nT T σπσ-== (16)
式中s T 为抽样间隔，至少应为杂波频谱宽度的倒数。

这样，在高斯谱已知的情况
下，非递归滤波器的加权系数i a 就由n C 完全确定了。

该滤波器的主要特点：
(1) 首先，它具备非递归滤波器的优点，结构简单，运算速度快。

(2) 便于硬件实现，特别适用于雷达模拟器。

(3) 要得到一个较好的响应，N 值应大于8。

(4) 这种方法对于输出序列的长度没有限制，取决于输入序列的长度。

这对雷达系统的性能测试具有重要意义，例如对虚警概率进行测试时，应给出足够长的序列，如虚警概率610F P -=时，其长度应大于810。

(5) 我们知道描述数字滤波器的差分方程时稳态情况下的差分方程，在输入
序列小于它的阶数时，输出序列仍处于暂态期，它们不满足给定的统计特性。

因此，在将其用于雷达模拟器时必须控制暂态输出。

下面设计一个非递归高斯滤波器，程序见“Matlab 程序”文件夹fseries.m 。

由例子知，20f Hz σ=，设采样频率512s F Hz =（杂波模拟时的脉冲重复频率）。

假定在傅立叶展开式中9N =，经计算，系数09C C 的数值给在表1中。

实际上，在表中也给出了1013C C 的数值，可以看出，1013C C 对频率响应的贡献已经是很小了。

根据此式所计算的功率谱密度曲线与理论值的差异如图2所示。

从图2可以看出，当N 取9时，说得到的功率谱曲线1与理论高斯谱模型曲线2重合在一起了；当N 取3时，功率谱密度曲线3要比理论高斯曲线2宽，并且在高端有小的起伏振荡。

计算表明，对于当N>8时，再增加谐波次数，效果并不明显，由表1看到，这时因为系数1013C C 的贡献太小的原因。

图2 滤波器输出的归一化的功率谱密度曲线
三、 滤波器2()H z 的设计以及广义χ分布变量()s k 的产生[4]
对于滤波器2的设计，由于通常非线性变换会使随机序列的功率谱展宽，所以，应使滤波器2的带宽足够窄，不致使非线性变换后输出的随机序列的功率谱展宽，以致()s k 序列的相关函数()s r m 不再接近于1，因而对2()H z 的要求主要是针对带宽的。

由于我们并不对控制()s k 的相关函数感兴趣，因此，可将滤波器2设计为一带宽很窄的低通滤波器，使得非线性变换随机序列的功率谱仍足够窄。

在本文例子的设计中，对于滤波器2，我们采用5阶butterworth 低通滤波器。

应用SIRP 的K 分布产生模型，必须产生广义χ分布随机变量()s k ，由于()s k 的平方即为伽马分布，所以可以产生伽马随机变量，再对它求平方根得到()s k 。

文献[5]4给出了伽马分布随机序列产生的方法，然而比较复杂。

对于图1所示的结构，文献[6]给出了非线性变换的表达式如下：
()
()()2(,)
1s s Q z v γνμ=-Γ
(17)
式中：(,)q p γ为不完全欧拉函数(Incomplete Eulerian Function)5。

文献[7]给出的
4 感谢大树文献服务网： 友情提供此篇文献的下载 5
这个函数的具体表达形式我还没有找到，谁找到的话请共享一下.
表达式如下：
()()222
,/()1g v E y s Q z απ⎡⎤=-⎣⎦
(18)
式中()()
1
01,a t b g a b e t dt b --=
Γ⎰为不完全伽马函数。

y 为滤波器1H 的输出，v 是K 分布形状参数，α为K 分布尺度参数。

由式(2)及()22E y σ=可知，式(18)可变化为：
()()2,2/1g v vs Q z π=-
(19)
()Q z 为标准正态随机变量的尾部面积，即有：
()2
2x
u Q z du ∞
⎛⎫=- ⎪⎝⎭⎰
(20)
将式(20)代入式(19)，并应用概率密度在全区间积分为1，得：
()211,2/
22g v vs erf π=
+ (21)
式中：()erf x 是误差函数，定义为：
()()2
exp x
erf x u du =
-⎰
(22)
因此，产生()s k 变量的问题转化为求式(21)的问题，这是一个非线性方程，可以
用二分法求解。

为了提到仿真时的速度，也可以提前对此非线性变换进行拟合，实际仿真时调用拟合得到的曲线关系即可。

后记
关于空域相关的海杂波仿真正在进一步研究之中。

回忆自己当初接触海杂波到现在写这个总结，中间经历了很多，有迷惑，也有失落，但在研学论坛的帮助下，在与我认识的同学(比如hahnc 和一位电子科技大学的同学（抱歉，不知道你研学上的id ）)的讨论中，我逐渐对海杂波有了初步的认识，只是最初步的认识。

为了感谢研学论坛，也为了总结一下以前的认识，特写了这篇总结(主要是摘抄，再组合， )，聊以共享。

希望各位研友再接再厉，在自己的研究道路上有新的成就。

参考文献
1.罗军辉，et al. Matlab 7.0在数字信号处理中的应用. 北京：机械工业出版社，
2005:180-184.
2.杨万海. 雷达系统建模与仿真. 西安：西安电子科技大学出版社，2007：
111-113.
3.郑君里，et al. 信号与系统（第二版）.北京：高等教育出版社，2000：145-149.
4.吕雁，史林. SIRP法相干相关K分布雷达杂波的建模与仿真. 现代雷达,
2002，24(2): 13-16.
5.Oliver, C. J.; Tough, R. J. A. On The Simulation of Correlated K-Distributed
Random Clutter. Opt. Acta, 1986,33(3):223-250.
6. E. Conte, M. Longo, and M. Lops. Modelling and simulation of non-Rayleigh
radar clutter. IEE Proceedings, Part F: Radar and Signal Processing, vol. 138, pp.
121-130, 1991.
7.杨俊岭, 吕韶昱, 万建伟.一种新的相干k分布模型及其在海杂波仿真中的应
用. 系统仿真学报，2007，19(2):250-253,260.。