CN105929371A - 一种基于协方差矩阵估计的机载雷达杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于协方差矩阵估计的机载雷达杂波抑制方法，思路为：计算快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl和N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl，并计算降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl和N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵R_kl；计算第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵和N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵并计算空时自适应滤波处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波，进而计算阵元—多普勒域机载雷达回波对应的距离‑多普勒谱。

Description

一种基于协方差矩阵估计的机载雷达杂波抑制方法

技术领域

本发明属于雷达杂波抑制技术领域，特别涉及基于协方差矩阵估计的机载雷达杂波抑制方法，并且是基于非均匀杂波条件下级联二维空时处理方法，适用于斜侧阵机载雷达杂波抑制或非正侧阵机载雷达杂波抑制。

背景技术

机载预警雷达已经成为现代战争的预警情报、战场监视和作战指挥系统的重要组成部分，它可以显著增加对低空突防目标的探测距离，大大提高武器系统的效能，在现在战争中扮演着越来越重要的角色，受到各军事大国的高度重视。机载预警雷达的波束指向一般是水平或稍下方向，使得机载预警雷达不可避免地会接收到地杂波。由于技术水平限制，使得阵列机载雷达天线的旁瓣电平较高，且阵列机载雷达回波中地杂波较强，可达90dB，远强于目标回波；此外，由于飞机平台的运动，地杂波的多普勒范围大大展宽，使得可供检测的多普勒范围大大降低，有很大一部分目标都“湮没”在地杂波中。因此，为了可靠检测远距离目标，首先要解决机载雷达的杂波抑制问题。

空时自适应处理(STAP)技术的诞生追溯到20世纪70年代，在杂波协方差矩阵和目标信号分别确知的条件下，Brennan和Reed等在期刊Aerospace and ElectronicSystems(AES)上，提出了全空时自适应处理(STAP)的概念和理论，其思想是将阵列信号处理的基本原理推广到由脉冲和阵元采样的两维场中；全空时自适应处理(STAP)能够取得比较理想的机载雷达杂波抑制效果，但是，机载雷达杂波协方差矩阵精确已知的条件在工程实现中难以满足，通常情况下是由参考单元估计得到；根据Reed，Mallett和Brennan三人提出的理论(即RMB准则)可知，为了使空时自适应处理(STAP)器的性能损失不超过3dB，空时自适应处理(STAP)器训练样本个数至少应为系统自由度的两倍，而且，机载雷达杂波协方差矩阵须从独立同分布的训练样本估计而来，所以全空时自适应处理(STAP)的运算量和设备复杂度令人难以接受。

为了降低运算量，德国的R.Klemm博士于1987年在期刊Signal Processing上提出了辅助通道法，将空时自适应处理(STAP)器维数从NM降至N+M-1，此处N表示机载雷达天线阵元个数，M表示一个相干处理间隔内机载雷达接收的脉冲个数；在没有误差的理想情况下，该方法能够取得比较理想的杂波抑制效果，性能接近最优空时自适应处理(STAP)处理，并且降低了估计杂波协方差矩阵所需要的独立同分布参考单元数目，但在有误差情况下，该方法的杂波抑制性能并不理想。

国内的西安电子科技大学在降维空时自适应处理(STAP)研究方面也展开了大量工作，并提出了局域化处理法，即因子化方法(1DT)和扩展因子化方法(mDT)；该局域化处理法首先利用低旁瓣的多普勒滤波器对机载雷达回波进行局域化处理，然后在空域自适应处理或空时域联合自适应处理将雷达杂波抑制掉。

1994年，H.Wang等人提出了局域化联合处理(JDL)方法，该方法利用两维傅里叶变换将接收到的雷达杂波由阵元-脉冲域变换到角度-多普勒域，然后分别在相邻的空域维和时域维选取3个波束进行联合自适应处理，该方法的维数由1DT(或3DT)的N(或3N)进一步降到9，显著减少计算量。

发明内容

针对上述已有技术在非均匀杂波环境下协方差矩阵估计不够准确的问题，本发明的目的在于充分利用机载雷达阵元和机载雷达脉冲回波，并提出一种基于协方差矩阵估计的机载雷达杂波抑制方法，该种基于协方差矩阵估计的机载雷达杂波抑制方法基于迭代方式的协方差矩阵估计形式，从而进行机载雷达杂波抑制，旨在非均匀杂波环境下提高机载雷达杂波的抑制性能。

为达到上述技术目的，本发明采用如下技术方案予以实现。

一种基于协方差矩阵估计的机载雷达杂波抑制方法，包括以下步骤：

步骤1，分别设定机载雷达的天线阵面在方位向均匀排列N个阵元，机载雷达回波中包含的距离门个数为L,一个相干处理间隔内机载雷达接收的脉冲个数为M，假设x_nlm表示第n个阵元、第l个距离门、第m个脉冲的机载雷达回波，得到第n个阵元、第l个距离门接收的机载雷达脉冲记为x_nl，并对x_nl进行快速傅里叶变换，得到快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl，进而依次计算得到第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达回波z_nkl和N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl；

其中，k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门总个数；

步骤2，初始化：设k表示当前多普勒通道索引，l表示当前的距离门索引，k和l初始值均为1；设定p为迭代次数，且p的初始值为1；

步骤3，获取N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练样本集Q_kl，Q_kl＝[z_k(l-MN) … z_k(l-1) z_k(l+1) … z_k(l+MN) S]，并计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集Q′_kl，Q′_kl＝[q₁ q₂ … q_v … q_V]，Q_kl中包含的元素分别与Q′_kl中包含的元素一一对应，q_v表示Q′_kl中第v个优化训练样本；q_V对应Q_kl中z_l的相位差导向矢量S，且下标V＝2MN+1；

其中，z_l表示与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量，z_k(l-1)表示N个阵元的第k个多普勒通道在第l-1个距离门接收的机载雷达回波，z_k(l+1)表示N个阵元的第k个多普勒通道在第l+1个距离门接收的机载雷达回波，z_k(l-MN)表示N个阵元的第k个多普勒通道在第l-MN个距离门接收的机载雷达回波,z_k(l+MN)表示N个阵元的第k个多普勒通道在第l+MN个距离门接收的机载雷达回波；

步骤4，对N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl进行降维处理，计算得到降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl，进而计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵R_kl；其中，如果当前多普勒通道索引k＝1时，N个阵元的第k-1个多普勒通道在第l个距离门接收的机载雷达回波z_(k-1)l为N个阵元的第K个多普勒通道在第l个距离门接收的机载雷达回波z_Kl；如果当前多普勒通道索引k＝K时，N个阵元的第K+1个多普勒通道在第l个距离门接收的机载雷达回波z_(K+1)l为N个阵元的第1个多普勒通道在第l个距离门接收的机载雷达回波z_1l；

步骤5，计算第p次迭代后第v个优化训练样本对N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数进而计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵其中，为V×1维列向量，V表示中包含的元素个数；

步骤6，对第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵中前V-1个元素依次进行平方归一化，得到第p次迭代后中第v′个元素对应的加权归一化系数然后计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵其中，v′∈{1,2,…V-1}，V表示中包含的元素个数；

步骤7，判断第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵是否达到稳态，假定如果则令p加1，返回步骤5；

若迭代终止，此时得到的第p次迭代后的加权系数矩阵达到稳态，进而得到的第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵为N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵

步骤8，利用N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵对降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z′k_l进行空时自适应处理，得到空时自适应滤波处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波y_kl；

步骤9，令l加1，重复步骤3—步骤8，直到l＝L时，第k个多普勒通道的L个距离门各自对应的机载雷达回波杂波抑制完成后，再令k加1，重复步骤3步骤8，直到k＝K时，此时K个多普勒通道各自对应的L个距离门机载雷达回波杂波全部抑制完成，并得到阵元—多普勒域机载雷达回波对应的距离-多普勒谱。

本发明与现有技术相比具有以下优点：

第一，本发明中提出的方法相比常规脉冲多普勒方法，充分使用导向矢量信息，使得处理过程更具有稳健性；

第二，使用常规扩展因子化法3DT估计机载雷达杂波协方差矩阵时，通常采用广义内积(GIP)方法作为非均匀检测器剔除训练样本的样本奇点，该方法在均匀杂波条件下表现较好，但在非均匀杂波环境的条件下表现极差；本发明方法在非均匀杂波环境下依然能够估计出机载雷达杂波协方差矩阵，从而更好地抑制机载雷达杂波；

第三，使用常规扩展因子化法3DT估计机载雷达杂波协方差矩阵时，各训练样本分别使用相同的加权系数，导致在非均匀条件下，不能更好地补偿机载雷达杂波起伏；而本发明方法修正了这一缺陷。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细说明。

图1是本发明方法的实现流程图；

图2a是使用脉冲多普勒方法进行杂波抑制后得到的距离多普勒图，其中横坐标为多普勒通道，纵坐标为距离门；

图2b是使用本发明方法进行杂波抑制后得到的距离多普勒图，其中横坐标为多普勒通道，纵坐标为距离门；

图3a是分别使用脉冲多普勒处理方法和本发明方法进行杂波抑制后的杂波剩余中1-300号距离门对应的杂波剩余图，其中横坐标表示多普勒通道，纵坐标表示杂波剩余，单位为dB；

图3b是分别使用脉冲多普勒处理方法和本发明方法进行杂波抑制后的杂波剩余中300-530号距离门对应的杂波剩余图。

具体实施方式

参考图1，为本发明方法的实现流程图；本发明的一种基于协方差矩阵估计的机载雷达杂波抑制方法，包括以下步骤：

步骤1，分别设定机载雷达的天线阵面在方位向均匀排列N个阵元，机载雷达回波中包含的距离门个数为L,一个相干处理间隔内机载雷达接收的脉冲个数为M，假设x_nlm表示第n个阵元、第l个距离门、第m个脉冲的机载雷达回波，得到第n个阵元、第l个距离门接收的机载雷达脉冲记为x_nl，并对x_nl进行快速傅里叶变换(FFT)，得到快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl，进而依次计算得到第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达回波z_nkl和N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl。

其中，k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门个数。

具体地，分别设定机载雷达的天线阵面在方位向均匀排列N个阵元，机载雷达回波中包含的距离门个数为L,一个相干处理间隔内机载雷达接收的脉冲个数为M，同时N也表示空时处理的空域自由度，M也表示空时处理的时域自由度；假设x_nlm表示第n个阵元、第l个距离门、第m个脉冲的机载雷达回波，m∈{1,2,…,M}，那么第n个阵元、第l个距离门接收的机载雷达脉冲记为x_nl，其表达式为：

x_nl＝[x_nl1 x_nl2 … x_nlm … x_nlM]^T

进而得到N个阵元在第l个距离门各自接收的机载雷达脉冲，并得到与所述N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快排矢量x_l，其表达式为：

$x_l={\left[\begin{array}{cccc}{x}_{1l}^T& x_{2l}^T& \mathrm{...}& x_{Nl}^T\end{array}\right]}^T$

因此，对第n个阵元在第l个距离门接收的机载雷达脉冲x_nl进行快速傅里叶变换(FFT)，得到快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl，其表达式为：

z_nl＝T^Hx_nl＝[z_nl1 z_nl2 … z_nlk … z_nlK]^T

进而得到快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲，并得到与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量z_l，其表达式为：

$z_l={\left[\begin{array}{cccccc}{z}_{1l}^T& z_{2l}^T& \mathrm{...}& z_{nl}^T& \mathrm{...}& z_{Nl}^T\end{array}\right]}^T$

其中，上标T表示转置，x_nlm表示第n个阵元、第l个距离门、第m个脉冲的机载雷达回波，z_nl表示第n个阵元在第l个距离门接收的机载雷达脉冲，z_nlk表示第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达脉冲，T表示加权快速傅里叶变换(FFT)矩阵，x_nl表示第n个阵元、第l个距离门接收的机载雷达脉冲，上标H表示共轭转置；如果第n个阵元在第l个距离门接收的机载雷达脉冲x_nl做FFT处理时进行了fftshift处理，那么加权快速傅里叶变换(FFT)矩阵T的第k列矢量T_k为：

T_k＝[a₀ a₁e^{jπ(k-K/2-1)/(K/2)} … a_k′e^{jπ(k-K/2-1)/(K/2)(k′)} … a_K-1e^{jπ(k-K/2-1)/(K/2)(K-1)}]^T

其中，k′∈{0,1,2,…,K-1}，K表示快速傅里叶变换(FFT)的点数，相当于对N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快排矢量x_l做K点快速傅里叶变换(FFT)，K为2的幂次，且K值大于等于一个相干处理间隔内机载雷达接收的脉冲个数M；同时，K也表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；a_k'表示加权快速傅里叶变换(FFT)矩阵中第k'个加权系数；然后依次计算得到第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达回波z_nkl和N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl。

所述第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达回波z_nkl和所述N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl，其表达式分别为：

$z_{nkl}=T_k^H{x}_{nl}=\underset{m=0}{\overset{M-1}{\Σ}}{a}_m{x}_{nl(m+1)}{e}^{-j\mathrm{\π}(k-K/2-1)/(K/2)m}$

z_kl＝[z_1kl z_2kl … z_nkl … z_Nkl]^T

其中，m∈{0,1,2,…,M-1}，M表示一个相干处理间隔内机载雷达接收的脉冲个数，a_m表示加权快速傅里叶变换(FFT)矩阵中第m个加权系数；x_nl(m+1)表示第n个阵元、第l个距离门、第m+1个脉冲的机载雷达回波，x_nl表示第n个阵元、第l个距离门接收的机载雷达脉冲，上标H表示共轭转置，T_k表示加权快速傅里叶变换(FFT)矩阵T的第k列矢量；K表示快速傅里叶变换(FFT)的点数，相当于对N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快排矢量x_l做K点快速傅里叶变换(FFT)，K一般选择2的幂次，且K值大于等于一个相干处理间隔内机载雷达接收的脉冲个数M；同时，K也表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门个数，x_nl(m+1)表示第n个阵元、第l个距离门、第m+1个脉冲的机载雷达回波。

步骤2，初始化：包括两层循环，第一层为多普勒通道，第二层为距离门，设k表示当前多普勒通道索引，l表示当前的距离门索引，k和l初始值均为1；设定p为迭代次数，且p的初始值为1。

步骤3，获取N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练样本集Q_kl，Q_kl＝[z_k(l-MN) … z_k(l-1) z_k(l+1) … z_k(l+MN)S]，并计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集Q′_kl，Q′_kl＝[q₁ q₂ … q_v … q_V]，Q_kl中包含的元素分别与Q′_kl中包含的元素一一对应，q_v表示Q′_kl中第v个优化训练样本；q_V对应Q_kl中z_l的相位差导向矢量S，且下标V＝2MN+1。

其中，z_l表示与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量，z_k(l-1)表示N个阵元的第k个多普勒通道在第l-1个距离门接收的机载雷达回波，z_k(l+1)表示N个阵元的第k个多普勒通道在第l+1个距离门接收的机载雷达回波，z_k(l-MN)表示N个阵元的第k个多普勒通道在第l-MN个距离门接收的机载雷达回波,z_k(l+MN)表示N个阵元的第k个多普勒通道在第l+MN个距离门接收的机载雷达回波。

具体地，为了保证空时自适应处理的性能，在此对称选取第l个距离门前后的2MN个距离门样本(不包括自身)，以及z_l的相位差导向矢量S，且z_l表示与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量，表示矩阵的Kronecker乘积，S_t表示z'_kl空时处理过程中对应的时域导向矢量，S_s表示z'_kl空时处理过程中对应的空域导向矢量，z'_kl表示降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波，且j表示虚数运算单位，表示x_nlm相对x_nl1的相位偏移，f_d表示归一化的时域多普勒频率， 表示x_nlm相对x_1lm的相位偏移，f_s表示空域频率，N表示机载雷达天线阵面方位向均匀包含的阵元个数，N也表示z'_kl空时处理过程中对应的空域自由度；然后获取N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练样本集Q_kl，Q_kl＝[z_k(l-MN) … z_k(l-1) z_k(l+1) … z_k(l+MN) S]；其中，z_k(l-1)表示N个阵元的第k个多普勒通道在第l-1个距离门接收的机载雷达回波，z_k(l+1)表示N个阵元的第k个多普勒通道在第l+1个距离门接收的机载雷达回波，z_k(l-MN)表示N个阵元的第k个多普勒通道在第l-MN个距离门接收的机载雷达回波,z_k(l+MN)表示N个阵元的第k个多普勒通道在第l+MN个距离门接收的机载雷达回波，M表示一个相干处理间隔内机载雷达接收的脉冲个数，N表示机载雷达天线阵面方位向包含的阵元个数。

为了下标方便管理，构造N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集Q′_kl，Q′_kl＝[q₁ q₂ … q_v … q_V]，Q_kl中包含的元素分别与Q′_kl中包含的元素一一对应，q_v表示Q′_kl中第v个优化训练样本；q_V对应Q_kl中z_l的相位差导向矢量S，且下标V＝2MN+1。

步骤4，对N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl进行降维处理，计算得到降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl，进而计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵R_kl；其中，如果当前多普勒通道索引k＝1时，N个阵元的第k-1个多普勒通道在第l个距离门接收的机载雷达回波z_(k-1)l为N个阵元的第K个多普勒通道在第l个距离门接收的机载雷达回波z_Kl；如果当前多普勒通道索引k＝K时，N个阵元的第K+1个多普勒通道在第l个距离门接收的机载雷达回波z_(K+1)l为N个阵元的第1个多普勒通道在第l个距离门接收的机载雷达回波z_1l。

具体地，本发明对N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl进行3DT降维处理，降维处理后的N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl自由度降为3N，该方法利用当前多普勒通道k以及与其相邻的两个多普勒通道k-1和k+1进行降维处理，且如果当前多普勒通道索引k＝1时，N个阵元的第k-1个多普勒通道在第l个距离门接收的机载雷达回波z_(k-1)l为N个阵元的第K个多普勒通道在第l个距离门接收的机载雷达回波z_Kl；如果当前多普勒通道索引k＝K时，N个阵元的第K+1个多普勒通道在第l个距离门接收的机载雷达回波z_(K+1)l为N个阵元的第1个多普勒通道在第l个距离门接收的机载雷达回波z_1l，进而计算得到降维处理后的第k个多普勒通道、第l个距离门的机载雷达回波的向量形式表达式为：

经过整理，得到降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl，其表达式为：

${z^{\′}}_{kl}={(I_N\⊗B_k)}^H{x}_l$

其中，I_N表示N×N的单位阵，B_k＝[T_k-1 T_k T_k+1]，x_l表示N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快排矢量，T_k表示加权快速傅里叶变换(FFT)矩阵T的第k列矢量，B_k表示加权快速傅里叶变换(FFT)矩阵T的组合矩阵，称为降维矩阵，表示矩阵的Kronecker乘积，上标H表示共轭转置。

依据上面给出的数据矢量表示形式，获得降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl，并根据RMB准则，估计阵元—多普勒域机载雷达回波的协方差矩阵所需样本数必须大于2MN-3，M表示z'_kl空时处理过程中对应的时域自由度，降维处理后变为3，亦即M＝3；N表示z'_kl空时处理过程中对应的空域自由度，N也表示机载雷达天线阵面方位向包含的阵元个数；z'_kl表示降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波。

然后计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵R_kl，然后将R_kl的初值设置为MN维单位矩阵I_kl，N表示z'_kl空时处理过程中对应的空域自由度，N也表示机载雷达天线阵面方位向包含的阵元个数；M表示一个相干处理间隔内机载雷达接收的脉冲个数。

步骤5，计算第p次迭代后第v个优化训练样本对N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数进而计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵其中，为V×1维列向量，V表示中包含的元素个数。

具体地，所述第p次迭代后第v个优化训练样本对N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数和所述第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵其表达式分别为：

${\mathrm{\β}}_{vkl}^p=\frac{q_v^H{\left(R_{kl}^{p-1}\right)}^{-1}{z^{\′}}_{kl}}{q_v^H{\left(R_{kl}^{p-1}\right)}^{-1}{q}_v}$

${\mathrm{\β}}_{kl}^p={\[\left|{\mathrm{\β}}_{1kl}^p\right|,\left|{\mathrm{\β}}_{2kl}^p\right|,...,\left|{\mathrm{\β}}_{vkl}^p\right|,...,\left|{\mathrm{\β}}_{(V-1)kl}^p\right|,\left|{\mathrm{\β}}_{Vkl}^p\right|\]}^T$

其中，为V×1维列向量，V表示中包含的元素个数；表示第p-1次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵，表示第p次迭代后z_l的相位差导向矢量S对的加权系数，z_l表示与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量。

步骤6，对第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵中前V-1个元素依次进行平方归一化，得到第p次迭代后中第v′个元素对应的加权归一化系数然后计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵其中，v′∈{1,2,…V-1}，V表示中包含的元素个数。

具体地，所述第p次迭代后中第v′个元素对应的加权归一化系数其表达式为：

$|{\widehat{\mathrm{\β}}}_{v^{\′}kl}^p{|}^2=\frac{|{\mathrm{\β}}_{v^{\′}kl}^p{|}^2}{\underset{v^{\′}=1}{\overset{V-1}{\Σ}}|{\mathrm{\β}}_{v^{\′}kl}^p{|}^2},v^{\′}\∈\{1,2,...V-1\}$

其中，v′∈{1,2,…V-1}，V表示中包含的元素个数，表示第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵，表示第p次迭代后第v′个训练样本对第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数，q_v'表示Q′_kl中第v'个优化训练样本，Q′_kl表示N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集。

利用第p次迭代后中第v′个元素对应的加权系数归一化系数计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵 q_v'表示Q′_kl中第v'个优化训练样本，Q′_kl表示N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集。

步骤7，判断第p次迭代后的加权系数矩阵是否达到稳态，假定如果则令p加1，返回步骤5；

若迭代终止，此时得到的第p次迭代后的加权系数矩阵达到稳态，进而得到的第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵为N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵

步骤8，利用N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵对降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z′_kl进行空时自适应(STAP)处理，得到空时自适应(STAP)滤波处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波y_kl。

具体地，利用N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵通过N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波空时自适应权值计算得到空时自适应(STAP)滤波处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波y_kl。

$y_{kl}=w_{kl}^H{z}_{kl}^{\′}$

其中，w_kl表示N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波空时自适应权值，z′_kl表示降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波，k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门个数。

步骤9，令l加1，重复步骤3—步骤8，直到l＝L时，第k个多普勒通道的L个距离门各自对应的机载雷达回波杂波抑制完成后，再令k加1，重复步骤3—步骤8，直到k＝K时，此时K个多普勒通道各自对应的L个距离门机载雷达回波杂波全部抑制完成，并得到阵元—多普勒域机载雷达回波对应的距离-多普勒谱。

通过以下仿真试验对本发明效果作进一步验证说明。

(一)实验参数及实验条件

机载雷达的天线阵面结构采用斜侧阵面阵，各个阵元均匀排列在椭圆天线阵面上，天线阵面的方位向和俯仰向分别均匀排列56个阵元和10个阵元，阵元间距为d＝λ/2，λ为载波波长，仿真实验的杂波数据脉冲个数为500，距离门个数为530。

(二)实验内容及结果分析

分别采用常规脉冲多普勒滤波方法和本发明方法雷达回波进行杂波抑制，得到的两幅距离多普勒图分别如图2a、图2b所示；图2a是使用脉冲多普勒方法进行杂波抑制后得到的距离多普勒图，图2b是使用本发明方法进行杂波抑制后得到的距离多普勒图，其中，图2a和图2b中的横坐标分别为多普勒通道，纵坐标分别为距离门，并且白色区域分别表示杂波抑制后的杂波剩余功率分布，黑色区域分别表示杂波抑制后的噪声剩余功率分布。

从图2a中可以看出，采用常规脉冲多普勒滤波方法进行杂波抑制后，杂波剩余功率较强，主瓣所占的多普勒宽度较大；从图2b中可以看出，本发明方法能够有效降低杂波剩余，对旁瓣杂波的抑制性能较好，而且主瓣所占的多普勒范围也变窄，有利于后续对目标进行检测。

分别常规脉冲多普勒滤波方法和本发明方法对回波数据进行杂波抑制，绘制两种方法的杂波剩余比较图，如图3a和图3b所示；图3a是分别使用脉冲多普勒处理方法和本发明方法进行杂波抑制后的杂波剩余中1-300号距离门对应的杂波剩余图，其中横坐标表示多普勒通道，纵坐标表示杂波剩余，单位为dB；图3b是分别使用脉冲多普勒处理方法和本发明方法进行杂波抑制后的杂波剩余中300-530号距离门对应的杂波剩余图。

从图3a和图3b可以看出，与常规脉冲多普勒滤波方法相比，本发明方法使主瓣杂波变窄且使旁瓣杂波功率降低，而且，与常规脉冲多普勒滤波方法相比，本发明方法的旁瓣杂波的抑制能力更佳。

综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.一种基于协方差矩阵估计的机载雷达杂波抑制方法，其特征在于，包括以下步骤：

步骤1，分别设定机载雷达的天线阵面在方位向均匀排列N个阵元，机载雷达回波中包含的距离门个数为L,一个相干处理间隔内机载雷达接收的脉冲个数为M，假设x_nlm表示第n个阵元、第l个距离门、第m个脉冲的机载雷达回波，得到第n个阵元、第l个距离门接收的机载雷达脉冲记为x_nl，并对x_nl进行快速傅里叶变换，得到快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl，进而依次计算得到第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达回波z_nkl和N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl；

其中，k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门总个数；

步骤2，初始化：设k表示当前多普勒通道索引，l表示当前的距离门索引，k和l初始值均为1；设定p为迭代次数，且p的初始值为1；

步骤3，获取N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练样本集Q_kl，Q_kl＝[z_k(l-MN) … z_k(l-1) z_k(l+1) … z_k(l+MN) S]，并计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集Q′_kl，Q′_kl＝[q₁ q₂ … q_v … q_V]，Q_kl中包含的元素分别与Q′_kl中包含的元素一一对应，q_v表示Q′_kl中第v个优化训练样本；q_V对应Q_kl中z_l的相位差导向矢量S，且下标V＝2MN+1；

其中，z_l表示与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量，z_k(l-1)表示N个阵元的第k个多普勒通道在第l-1个距离门接收的机载雷达回波，z_k(l+1)表示N个阵元的第k个多普勒通道在第l+1个距离门接收的机载雷达回波，z_k(l-MN)表示N个阵元的第k个多普勒通道在第l-MN个距离门接收的机载雷达回波,z_k(l+MN)表示N个阵元的第k个多普勒通道在第l+MN个距离门接收的机载雷达回波；

步骤4，对N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl进行降维处理，计算得到降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl，进而计算得到N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵R_kl；其中，如果当前多普勒通道索引k＝1时，N个阵元的第k-1个多普勒通道在第l个距离门接收的机载雷达回波z_(k-1)l为N个阵元的第K个多普勒通道在第l个距离门接收的机载雷达回波z_Kl；如果当前多普勒通道索引k＝K时，N个阵元的第K+1个多普勒通道在第l个距离门接收的机载雷达回波z_(K+1)l为N个阵元的第1个多普勒通道在第l个距离门接收的机载雷达回波z_1l；

步骤5，计算第p次迭代后第v个优化训练样本对N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数进而计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵其中，为V×1维列向量，V表示中包含的元素个数；

步骤6，对第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵中前V-1个元素依次进行平方归一化，得到第p次迭代后中第v′个元素对应的加权归一化系数然后计算得到第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵其中，v′∈{1,2,…V-1}，V表示中包含的元素个数；

步骤7，判断第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵是否达到稳态，假定如果则令p加1，返回步骤5；

若迭代终止，此时得到的第p次迭代后的加权系数矩阵达到稳态，进而得到的第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵为N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵

步骤8，利用N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波的最优协方差矩阵对降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z′_kl进行空时自适应处理，得到空时自适应滤波处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波y_kl；

步骤9，令l加1，重复步骤3—步骤8，直到l＝L时，第k个多普勒通道的L个距离门各自对应的机载雷达回波杂波抑制完成后，再令k加1，重复步骤3步骤8，直到k＝K时，此时K个多普勒通道各自对应的L个距离门机载雷达回波杂波全部抑制完成，并得到阵元—多普勒域机载雷达回波对应的距离-多普勒谱。

2.如权利要求1所述的一种基于协方差矩阵估计的机载雷达杂波抑制方法，其特征在于，在步骤1中，所述第n个阵元、第l个距离门接收的机载雷达脉冲记为x_nl，其表达式为：

x_nl＝[x_nl1 x_nl2 …x_nlm… x_nlM]^T

所述快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl，以及所述第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达回波z_nkl和N个阵元的第k个多普勒通道在第l个距离门接收的机载雷达回波z_kl，其表达式分别为：

z_nl＝T^Hx_nl＝[z_nl1 z_nl2 …z_nlk… z_nlK]^T

$z_{nkl}=T_k^H{x}_{nl}=\underset{m=0}{\overset{M-1}{\Σ}}{a}_m{x}_{nl(m+1)}{e}^{-j\mathrm{\π}(k-K/2-1)/(K/2)m}$

z_kl＝[z_1kl z_2kl … z_nkl … z_Nkl]^T

其中，m∈{0,1,2,…,M-1}，M表示一个相干处理间隔内机载雷达接收的脉冲个数，a_m表示加权快速傅里叶变换矩阵中第m个加权系数；x_nl(m+1)表示第n个阵元、第l个距离门、第m+1个脉冲的机载雷达回波，x_nl表示第n个阵元、第l个距离门接收的机载雷达脉冲，上标T表示转置，x_nlm表示第n个阵元、第l个距离门、第m个脉冲的机载雷达回波，z_nl表示第n个阵元在第l个距离门接收的机载雷达脉冲，z_nlk表示第n个阵元、第k个多普勒通道在第l个距离门接收的机载雷达脉冲，T表示加权快速傅里叶变换矩阵，T_k表示加权快速傅里叶变换矩阵T的第k列矢量，上标H表示共轭转置，k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门个数，M表示一个相干处理间隔内机载雷达接收的脉冲个数。

3.如权利要求1所述的一种基于协方差矩阵估计的机载雷达杂波抑制方法，其特征在于，在步骤4中，所述降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波z'_kl和所述N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵R_kl，其表达式分别为：

${Z^{\′}}_{kl}={(I_N\⊗B_k)}^H{x}_l$

$R_{kl}=\frac{1}{2MN}\underset{v=1}{\overset{2MN}{\Σ}}{q}_v{q}_v^H$

其中，I_N表示N×N的单位阵，N表示z'_kl空时处理过程中对应的空域自由度，N也表示机载雷达天线阵面方位向包含的阵元个数；M表示一个相干处理间隔内机载雷达接收的脉冲个数，B_k＝[T_k-1 T_k T_k+1]，x_l表示N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快排矢量，T_k表示加权快速傅里叶变换矩阵T的第k列矢量，B_k表示加权快速傅里叶变换矩阵T的组合矩阵，表示矩阵的Kronecker乘积，上标H表示共轭转置，q_v表示Q′_kl中第v个优化训练样本，Q′_kl表示N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集；k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数。

4.如权利要求1所述的一种基于协方差矩阵估计的机载雷达杂波抑制方法，其特征在于，在步骤5中，所述第p次迭代后第v个优化训练样本对N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数和所述第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵其表达式为：

${\mathrm{\β}}_{vkl}^p=\frac{q_v^H{\left(R_{kl}^{p-1}\right)}^{-1}{z^{\′}}_{kl}}{q_v^H{\left(R_{kl}^{p-1}\right)}^{-1}{q}_v}$

${\mathrm{\β}}_{kl}^p={\[\left|{\mathrm{\β}}_{1kl}^p\right|,\left|{\mathrm{\β}}_{2kl}^p\right|,...,\left|{\mathrm{\β}}_{vkl}^p\right|,...,\left|{\mathrm{\β}}_{(V-1)kl}^p\right|,\left|{\mathrm{\β}}_{Vkl}^p\right|\]}^T$

其中，为V×1维列向量，表示第p-1次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵，表示第p次迭代后z_l的相位差导向矢量S对的加权系数，z_l表示与快速傅里叶变换后的N个阵元在第l个距离门各自接收的机载雷达脉冲相对应的空时数据快拍矢量。

5.如权利要求1所述的一种基于协方差矩阵估计的机载雷达杂波抑制方法，其特征在于，在步骤6中，所述第p次迭代后中第v′个元素对应的加权归一化系数和所述第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差归一化矩阵其表达式为：

$\begin{array}{cc}|{\widehat{\mathrm{\β}}}_{v^{\′}kl}^p{|}^2=\frac{|{\mathrm{\β}}_{v^{\′}kl}^p{|}^2}{\underset{v^{\′}=1}{\overset{V-1}{\Σ}}|{\mathrm{\β}}_{v^{\′}kl}^p{|}^2}& v^{\′}\∈\{1,2,...V-1\}\end{array}$

${\hat{R}}_{kl}^p=\underset{v^{\′}=1}{\overset{V-1}{\Σ}}|{\widehat{\mathrm{\β}}}_{v^{\′}kl}^p{|}^2{q}_{v^{\′}}{q}_{v^{\′}}^H$

其中，v′∈{1,2,…V-1}，V表示中包含的元素个数，表示第p次迭代后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波协方差加权系数矩阵，表示第p次迭代后第v′个训练样本对第k个多普勒通道、第l个距离门的机载雷达回波杂波协方差矩阵的加权系数，q_v'表示Q′_kl中第v'个优化训练样本，Q′_kl表示N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波扩展训练优化样本集。

6.如权利要求1所述的一种基于协方差矩阵估计的机载雷达杂波抑制方法，其特征在于，在步骤8中，所述空时自适应滤波处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波y_kl，其表达式为：其中，上标H表示共轭转置，w_kl表示N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波杂波空时自适应权值，z′_kl表示降维处理后N个阵元的第k个多普勒通道、第l个距离门的机载雷达回波；

k∈{1,2,…,K}，K表示快速傅里叶变换后的第n个阵元在第l个距离门接收的机载雷达脉冲z_nl包含的多普勒通道总个数；l∈{1,2,…L}，L表示机载雷达回波中包含的距离门个数。