CN101226236A - 基于旁瓣约束的天波超视距雷达自适应干扰抑制方法 - Google Patents

Abstract

本发明公开了一种基于旁瓣约束的天波超视距雷达自适应干扰抑制方法，主要解决分段自适应处理时海、地杂波谱展宽问题。具体过程是：对待检测单元和样本单元分别形成指向θ0的接收波束作为主天线；选取若干接收天线作为辅助天线；对样本单元主天线和辅助天线的数据分别进行3脉冲MTI处理；对待检测单元辅助天线各段数据分别进行多普勒滤波，保留零多普勒通道输出；利用样本单元数据求解给定约束条件下待检测单元辅助天线各段数据对应的自适应权值；对样本单元数据进行自适应滤波；对自适应滤波后的数据进行相干积累得到干扰抑制后的多普勒频谱。本发明具有干扰抑制性能高、运算量小的优点，可用于天波超视距雷达自适应干扰抑制。

Description

基于旁瓣约束的天波超视距雷达自适应干扰抑制方法

技术领域

本发明属于雷达技术领域，具体的说是一种基于旁瓣约束的天波超视距雷达自适应干扰抑制方法，用于在保证干扰抑制性能的条件下，消除时变自适应权值引起的海、地杂波谱展宽。

背景技术

天波超视距雷达SKW-OTHR以其大范围的监视能力、防低空突袭、抗隐身以及早期预警能力强等突出优点，受到了许多国家的重视。SKW-OTHR依靠电离层的折射作用进行远距离目标探测，其作用距离可达几千公里。虽然SKW-OTHR分辨率、精度不高，但作为早期远程预警具有重大价值。SKW-OTHR利用电离层进行下视探测，因此可以探测到远距离低空和超低空飞行的目标。此外，由于SKW-OTHR工作在高频波段，飞机、舰船等目标的尺寸与雷达的工作波长相近，处于谐振散射区，它们的雷达目标截面RCS取决于总尺寸而不是外形的细节，所以对微波雷达有效的控制外形使雷达截面减小技术，对超视距雷达几乎无效。同时SKW-OTHR作为一种远程警戒雷达能够提供大范围的早期预警且预警时间远远超过一般的雷达。

高频频段外部干扰和噪声比接收机内噪声高得多，通常可达20～40分贝，是限制SKW-OTHR检测性能的主要因素。流星余迹回波、大气噪声以及雷电干扰等经过电离层的传播可以传播很远，如果不对这些干扰进行有效抑制，SKW-OTHR的检测性能是难以保证的。此外，这些干扰普遍具有空变和时变的特性，给自适应处理带来相当多的困难。

对付空中定向干扰的有效方法是方向图置零，即采用加权的方法使接收阵列的合成方向图在干扰的来波方向形成零点，从而有效抑制干扰。但高频频段来源众多的各种干扰所普遍具有的空变和时变特性，要求形成方向图零点的权矢量也必须是时变的。由于SKW-OTHR工作的波长较长，为了获得较高的角度分辨率，其接收阵列孔径可达公里级，接收阵元可达上百个，如果对全部接收阵元做自适应处理，一方面用来估计阵列协方差矩阵的样本数严重不足；另一方面自适应处理的运算量也是十分巨大的。因此，通常采用子阵处理或广义旁瓣相消方式来降低运算量。为了使自适应处理的权矢量与干扰的空变、时变特性相适应，必须采用分段自适应处理的方法。在常规分段自适应处理中，通常对接收天线阵列合成方向图的主瓣方向增益有约束，虽然能够自适应的在干扰方向形成零点，但旁瓣电平随自适应权值的变化是时变的。此时，时变的旁瓣电平会对经由接收阵列合成方向图旁瓣进入的海、地杂波回波信号附加相位调制，导致海、地杂波回波信号时域相干性下降，在相干积累时会导致杂波谱的展宽，影响低速目标的检测。

澳大利亚研究人员G..A.Fabrizio在2004年IEEE的AES上发表论文，提出采用样本单元回波信号对自适应权值进行约束解决海地、杂波谱展宽。该方法隐含的假设条件是各段的自适应权值对与本段相邻的若干重复周期的干扰是有效的，即能够得到干扰完全抑制的仅含有海、地杂波的回波信号。这若干重复周期的海、地杂波回波信号将用于约束条件，由此得到的自适应权值将使相邻数据段海、地杂波自适应输出结果保持相干性，使相干积累后的海、地杂波谱不展宽。如果假设条件不完全成立，即干扰不能完全抑制，用于约束条件的数据不仅包含海、地杂波回波信号，还含有干扰剩余，将导致由此得到的自适应权值不能将干扰完全对消。由于干扰对消剩余信号进入了约束条件，所以干扰对消剩余在各段之间具有累积效应，导致总体干扰对消性能下降很多。此外，由于使用多个重复周期的信号用于约束条件，相对于常规的分段自适应处理，运算量增加了不少。

发明内容

本发明的目的是克服传统天波超视距雷达SKW-OTHR自适应干扰抑制方法的不足，提出一种基于旁瓣约束的天波超视距雷达自适应干扰抑制方法，以保证在干扰抑制性能的条件下，消除时变自适应权值引起的海、地杂波谱展宽。

为实现上述目的本发明基于来自接收阵列天线旁瓣的高频频段的外部干扰由于受电离层不稳定性的影响而表现出空变、时变的特点；以及外部干扰可看作点源的特点；而海、地杂波是由于发射信号引起的，在时域是相关的，在空域则是分布式的特点，提出如下自适应干扰抑制过程：

A、对待检测单元和样本单元的数据分别进行波束形成，得到指向待检测方位θ₀的接收波束，其输出作为主天线接收数据；

B、在接收阵列中选取N个阵元作为辅助天线，选定的N个辅助天线所对应的待检测单元和样本单元的数据作为辅助天线的接收数据；

C、对样本单元的主天线数据和各个辅助天线数据分别进行3脉冲MTI处理，获得仅含有干扰信息的样本数据；

D、对待检测单元和样本单元主天线以及各个辅助天线的数据分别进行相同的时域分段；

E、对待检测单元各个辅助天线各段数据分别进行多普勒滤波，并保留零多普勒通道输出；

F、采用常规自适应滤波算法求解待检测单元各个辅助天线第一段数据对应的自适应权值；

G、在给定约束条件下，利用样本单元数据分别求解待检测单元各个辅助天线第k段数据对应的自适应权值，k＝2，3，…，K；

H、使用待检测单元各个辅助天线各段数据对应的自适应权值分别对待测单元各段数据进行自适应滤波，获得干扰抑制后的输出结果；

I、对自适应滤波后的输出结果进行相干积累，获得干扰抑制后的多普勒频谱。

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

1、常规的分段自适应处理虽然能够抑制空变、时变干扰，但是会导致海、地杂波谱展宽，而采用本发明的方法不仅能有效抑制空变、时变干扰，而且能消除海、地杂波谱的展宽。

2、直接采用样本单元回波信号用于自适应权值约束的方法，如果存在干扰对消剩余，干扰对消剩余将具有累积效应，导致干扰对消性能下降很多；本发明的方法不仅能够有效抑制干扰，而且能够消除海、地杂波谱的展宽，干扰抑制性能优于采用样本单元回波信号用于自适应权值约束的方法。

3、本发明提出的基于旁瓣约束的天波超视距雷达自适应干扰抑制方法与常规分段自适应处理方法运算量基本相当，但比直接采用样本单元回波信号用于自适应权值约束的方法运算量少，复杂度低。

附图说明

图1是本发明的流程图；

图2是本发明中自适应旁瓣相消示意图；

图3是对待检测单元进行波束形成，当主波束指向为4°时，距离单元390的多普勒谱；

图4是对待检测单元进行波束形成，当主波束指向为4°时，经过常规分段自适应处理后，距离单元390的多普勒谱；

图5是对待检测单元进行波束形成，当主波束指向为4°时，采用样本单元回波信号用于约束条件，经过分段自适应处理后，距离单元390的多普勒谱；

图6是对待检测单元进行波束形成，当主波束指向为4°时，采用本发明的自适应处理方法进行自适应干扰抑制后，距离单元390的多普勒谱；

图7是对待检测单元进行波束形成，当主波束指向为9°时，距离单元390的多普勒谱；

图8是对待检测单元进行波束形成，当主波束指向为9°时，经过常规分段自适应处理后，距离单元390的多普勒谱；

图9是对待检测单元进行波束形成，当主波束指向为9°时，采用样本单元回波信号用于约束条件，经过分段自适应处理后，距离单元390的多普勒谱；

图10是对待检测单元进行波束形成，当主波束指向为9°时，采用本发明的自适应处理方法进行自适应干扰抑制后，距离单元390的多普勒谱。

具体实施方式

为了清楚的说明本发明的具体实施方法，首先简单介绍天波超视距雷达基本的接收数据格式及信号处理流程，以及本发明中所用到的一些术语的定义。

天波超视距雷达采用N_a个接收天线同时接收目标回波信号以及杂波、干扰等，在一个相干积累时间CIT内每个接收天线接收N_p个脉冲重复周期的信号；每个脉冲重复周期又包含若干距离单元，用N_r表示。在天波超视距雷达目标检测中，首先对同一个距离单元的来自N_a个接收天线的数据进行波束形成，以保证波束形成结果中来自特定方向的信号被正确接收，而抑制来自于其他方向的信号。然后对波束形成的输出结果进行相干积累，得到该距离单元在特定方向的多普勒谱。如果在多普勒域检测到目标，则目标的距离为该距离单元对应的距离，目标的方位为该波束的指向，目标的多普勒频移为该目标在多普勒谱中对应的值。

待检测单元：当前需要判别是否有目标存在的距离单元，有可能包含目标回波信号。

样本单元：不包含目标回波信号的距离单元。

待测单元接收数据：X(m)为N_a×1维矢量，m＝1，2，…，N_p；

样本单元接收数据：X_p(m，r)为N_a×1维矢量，m＝1，2，…，N_p，r＝1，2，…，N_r；

其中，m表示脉冲重复周期序号，r表示距离单元序号。

参照图1和图2，本发明的实施步骤如下：

1、对待检测单元和样本单元的数据分别进行波束形成，得到指向待检测方位θ₀的接收波束，其输出作为主天线接收数据。

对待检测单元和样本单元分别进行波束形成，得到指向θ₀的接收波束，波束形成权矢量为S(θ₀)，为降低接收波束的旁瓣电平，可对权矢量进行为加窗函数处理。

待检测单元波束形成的输出为：

d₀(m)＝S^H(θ₀)·X(m)，m＝1，2，…，N_p    (1)

样本单元波束形成的输出为：

d_p(m，r)＝S^H(θ₀)·X_p(m，r)，m＝1，2，…，N_p，r＝1，2，…，N_r    (2)

式中，H表示复共轭转置。

2、在接收阵列中选取N个阵元作为辅助天线，选定的N个辅助天线所对应的待检测单元和样本单元的数据作为辅助天线的接收数据，分别为：

X₀(m)，m＝1，2，…，N_p    (3)

X_0p(m，r)，m＝1，2，…，N_p，r＝1，2，…，N_r    (4)

式中，X₀(m)和X_0p(m，r)均为N×1维矢量，N为辅助天线个数。

3、对样本单元的主天线数据和辅助天线数据分别进行3脉冲MTI处理，获得仅含有干扰信息的样本数据。

通常认为天波超视距雷达存在“空阔区”，即存在不包含海、地杂波回波信号的距离段，但实际工作中，由于二次回波的出现，导致实际接收数据中几乎不存在“空阔区”。而在进行自适应处理时，样本单元应该仅含有干扰的信息，而不应该含有目标和杂波回波信号。

因此，为了得到仅含有干扰的样本数据，需要对样本单元数据进行MTI处理，滤除海、地杂波回波信号，本发明采用3脉冲MTI处理。如果杂波谱比较宽，也可以采用更多的脉冲进行MTI处理。经过3脉冲MTI处理后样本单元主天线和辅助天线的输出分别为：

d_s(m，r)＝d_p(m，r)-2d_p(m+1，r)+d_p(m+2，r)       (5)

X_s(m，r)＝X_0p(m，r)-2X_0p(m+1，r)+X_0p(m+2，r)    (6)

式中，m＝1，2，…，N_p-2，r＝1，2，…，N_r。

4、对待检测单元和样本单元主天线以及辅助天线的的数据分别进行相同的时域分段，分段依据高频频段外部干扰的时变特点，一般每段数据约为200～300毫秒的时间。

设每段数据包含M个脉冲的回波数据，共分为K段，即

N_p＝K×M    (7)

由于样本单元数据经过3脉冲MTI处理，所以时域总脉冲数N_p-2个，但是仍将样本单元数据分为K段，最后一段包含M-2个脉冲的数据，其余各段包含M个脉冲的数据。

5、对待检测单元各个辅助天线各段数据分别进行多普勒滤波，并保留零多普勒通道输出。保留零多普勒通道输出的过程可用式(8)表示，即

$X_c\left(k\right)=\underset{m=(k-1)\·M+1}{\overset{\mathrm{Mk}}{\mathrm{\Σ}}}{X}_0\left(m\right)e^{j2\mathrm{\π}\frac{m\·0}{M}}---\left(8\right)$

$=\underset{m=(k-1)\·M+1}{\overset{\mathrm{Mk}}{\mathrm{\Σ}}}{X}_0\left(m\right)$

式中，X_c(k)表示待检测单元各个辅助天线第k段数据的零多普勒通道输出；k＝1，2，…，K；X₀(m)是N×1维列矢量，表示待检测单元辅助天线的数据。

6、采用常规自适应滤波算法求解待检测单元各个辅助天线第一段数据对应的自适应权值，由图2给出的旁瓣相消示意图可知，自适应权值为下式的解：

$\underset{w\left(1\right)}{\arg \min }E\left[{|d_s(m,r)-W{\left(1\right)}^H\·X_s(m,r)\·W\left(1\right)|}^2\right]---\left(9\right)$

式中，m＝1，2，…，N_p-2，r＝1，2，…，N_r。

由式(9)可得

$W\left(1\right)=R_{\mathrm{xx}}^{-1}\left(1\right)\·R_{\mathrm{xd}}\left(1\right)---\left(10\right)$

式中，R_xx ^-1(1)表示对矩阵R_xx(1)的求逆运算。

在实际工作中，辅助天线的协方差矩阵R_xx(1)，辅助天线与主天线之间的互相关矢量R_xd(1)只能由样本单元的数据估计得到，即

${\hat{R}}_{\mathrm{xx}}\left(1\right)=\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{X}_s(m,r)\·X_s^H(m,r)---\left(11\right)$

${\hat{R}}_{\mathrm{xd}}\left(1\right)=\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{X}_s(m,r)\·d_s^H(m,r)---\left(12\right)$

将式(11)和式(12)代入式(10)可得第一段数据对应的自适应权值W(1)的估计值。

7、在给定约束条件下，利用样本单元数据分别求解待检测单元各个辅助天线第k段数据对应的自适应权值，k＝2，3，…，K。

为了保证经过分段自适应处理后的海、地杂波谱不展宽，对待检测单元辅助天线第k段数据对应的权值增加相应的约束条件，以保证经过空域自适应滤波后各段海、地杂波输出之间的相关性。根据图2所示的自适应旁瓣相消示意图可看出，上支路经由主瓣进入的海、地杂波已经做了约束，约束条件实际上是由S(θ₀)给出的。所以只需要对下支路的自适应权值增加约束，保证从接收天线阵列合成方向图旁瓣进入的海、地杂波回波信号在各段之间的相关性即可。

根据待检测单元当前段自适应滤波后海、地杂波输出与前一段自适应权值对本段海、地杂波滤波后的输出相等的条件，本发明中采用的具体约束条件如下：

    W^H(k-1)·X_c(k)＝W^H(k)·X_c(k)k＝2，3，…，K    (13)

此时，第k段的自适应权值为下式所给出的约束最小均方误差解：

$\left\{\begin{array}{c}\underset{w\left(k\right)}{\arg \min a}E\left[\right|{d_s(m,r)-W^H\left(k\right)\·X_s(m,r)\·W^H\left(k\right)|}^2]\\ s.t.W^H(k-1)\·X_c\left(k\right)=W^H\left(k\right)X_c\left(k\right)\end{array}\right.---\left(14\right)$

式中，m＝(k-1)·M+1，…，kM，r＝1，2，…，Nr，k＝2，3，…，K。

利用拉格阳乘数法，无约束代价函数可表示为：

$J\left(W\left(k\right)\right)={\mathrm{\σ}}_s^2-W^H\left(k\right)R_{\mathrm{xd}}\left(k\right)-R_{\mathrm{xd}}^H\left(k\right)\·W\left(k\right)+W^H\left(k\right)R_{\mathrm{xx}}\left(k\right)W\left(k\right)---\left(15\right)$

$+(W^H\left(k\right)\·X_c\left(k\right)-W^H(k-1)\·X_c\left(k\right))\·\mathrm{\λ}$

求代价函数关于W^H(k)的梯度得

令上式梯度等于0，得

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)(R_{\mathrm{xd}}\left(k\right)-X_c\left(k\right)\·\mathrm{\λ})---\left(17\right)$

将W(k)代入约束条件式(13)可得如下等式

$W^H\left(k\right)\·X_c\left(k\right)=(R_{\mathrm{xd}}^H\left(k\right)-{\mathrm{\λ}}^H\·X_c^H\left(k\right))R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)---\left(18\right)$

$=W^H(k-1)\·X_c\left(k\right)$

进一步可得：

$\mathrm{\λ}={(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}^{-1}[X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-X_c^H\left(k\right)\·W(k-1)]---\left(19\right)$

将式(19)代入式(17)，得到k＝2，3，…，K各段的自适应权值为

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right){(X_c^H\left(k\right)\·{R_{\mathrm{xx}}^{-1}\left(k\right)\·X}_c\left(k\right))}^{-1}---\left(20\right)$

$\·X_c^H\left(k\right)[R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-W(k-1)]$

在实际工作中，辅助天线的协方差矩阵R_xx(k)，辅助天线与主天线之间的互相关矢量R_xd(k)只能由样本单元的数据估计得到，即

${\hat{R}}_{\mathrm{xx}}\left(k\right)=\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{m=(k-1)m+1}{\overset{\mathrm{kM}}{\mathrm{\Σ}}}{X}_s(m,r)\·X_s^H(m,r)---\left(21\right)$

${\hat{R}}_{\mathrm{xd}}\left(k\right)=\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{m=(k-1)m+1}{\overset{\mathrm{kM}}{\mathrm{\Σ}}}{X}_s(m,r)\·d_s^H(m,r)---\left(22\right)$

式中，k＝2，3，…，K。

将

和

代入式(20)，即得到待检测单元辅助天线k＝2，3，…，K各段数据对应的自适应权值。

此外，式(20)所给出的第k段的自适应权值还有多种表达形式，如式(23)、式(24)和式(25)等。

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-\frac{R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}[R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-W(k-1)]---\left(23\right)$

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)[R_{\mathrm{xd}}\left(k\right)-\frac{X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}[R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-W(k-1)\left]\right]---\left(24\right)$

$W\left(k\right)=[I-\frac{R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}]R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)---\left(25\right)$

$+\frac{R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}W(k-1)$

需要说明的是，本发明所用的自适应权值计算公式即式(20)、(23)、(24)和(25)，与常规的自适应权值计算相比，计算量增加不大。因为在计算辅助天线协方差矩阵时，两者的运算量相当；而本发明中(X_c ^H·R_xx ^-1·X_c)^-1实际上是一个标量，相当于做除法，剩余部分只是简单的矢量乘法和减法，所以运算量增加不大。而直接采用样本单元回波信号进行约束的方法，也可以利用本发明所得到的计算公式，但此时(X_c ^H·R_xx ^-1·X_c)^-1项为一矩阵，需要做矩阵求逆，剩余部分还要做矩阵乘法和减法，所以运算量要增加不少。

8、使用待检测单元各个辅助天线各段数据对应的自适应权值分别对待测单元各段数据进行自适应滤波，获得干扰抑制后的输出结果。

自适应滤波的输出为：

y₀(m)＝d₀(m)-W^H(k)·X₀(m)    (26)

式中，m＝(k-1)·M+1，…，M，k＝1，2，3，…，K。

9、对自适应滤波后的输出结果进行相干积累，获得干扰抑制后的多普勒频谱。

对自适应滤波的输出y₀(m)进行相干积累，得到多普勒频谱即：

$Y\left(\mathrm{\ω}\right)=\underset{m=1}{\overset{N_p}{\mathrm{\Σ}}}{y}_0\left(m\right)w\left(m\right)e^{j2\mathrm{\π}\frac{\mathrm{m\ω}}{N_p}}---\left(27\right)$

式中，ω表示频率；w(m)表示窗函数，用于降低回波信号频谱的旁瓣电平。

本发明的优点可通过实测数据实验进一步说明。

1、实测数据录取及实验条件

(1)在本实验中，利用应答机发射的信号作为干扰源。干扰源大圆方位角为-29°，由于经过电离层的传播，干扰源对应的锥角约为-27.1°。本实验中接收阵列为等距线阵，所以其观测得到的方位角对应的是锥角，方位角相同，但仰角不同，对应的锥角是不同的。

(2)用全部接收阵元形成接收主波束，选取16个接收阵元作为辅助天线。距离单元390作为待检测单元，选取距离单元390附近100个距离单元的回波信号作为样本单元。在一个CIT中共包含512个脉冲的回波信号。

(3)使用同一组数据，同样的主波束指向，相同的辅助天线，分段自适应处理时，每段包含16个脉冲的回波数据，分别用常规分段自适应处理算法、直接采用样本单元回波信号约束的自适应处理算法以及本发明的自适应干扰抑制算法对数据进行处理。

2、实测数据处理结果

对待检测单元进行波束形成，当主波束指向为4°时，距离单元390的多普勒谱如图3；

对待检测单元进行波束形成，当主波束指向为4°时，经过常规分段自适应处理后，距离单元390的多普勒谱如图4；

对待检测单元进行波束形成，当主波束指向为4°时，采用样本单元回波信号用于约束条件，经过分段自适应处理后，距离单元390的多普勒谱如图5；

对待检测单元进行波束形成，当主波束指向为4°时，采用本发明的自适应处理方法进行自适应干扰抑制后，距离单元390的多普勒谱如图6；

对待检测单元进行波束形成，当主波束指向为9°时，距离单元390的多普勒谱如图7；

对待检测单元进行波束形成，当主波束指向为9°时，经过常规分段自适应处理后，距离单元390的多普勒谱如图8；

对待检测单元进行波束形成，当主波束指向为9°时，采用样本单元回波信号用于约束条件，经过分段自适应处理后，距离单元390的多普勒谱如图9；

对待检测单元进行波束形成，当主波束指向为9°时，采用本发明的自适应处理方法进行自适应干扰抑制后，距离单元390的多普勒谱如图10。

3、实测数据处理结果比较

从图3和图7中可以看出，海、地杂波谱集中在零频附近，而干扰信号集中在21赫兹附近。分别将图4与图3，图8与图7对比，可以看出，常规分段自适应处理算法虽然有效抑制了21赫兹附近的干扰，干扰对消比达30dB左右，但是在零频附近的海、地杂波谱却发生了展宽。

将图5分别与图4、图3对比，以及将图9分别与图8、图7对比，可以看出，直接采用样本单元回波信号用于自适应权值约束的方法，虽然能够消除海、地杂波谱的展宽，但是干扰对消性能下降很多，干扰对消比只有大约20dB左右，比常规分段自适应处理时干扰对消性能下降了10dB左右。

将图6分别与图5、图4、图3对比，以及将图10分别与图9、图8、图7对比，可以看出，本发明的基于旁瓣约束的天波超视距雷达干扰抑制方法，不仅能够消除海、地杂波谱的展宽，而且干扰对消性能与常规分段自适应处理基本相当，干扰对消比达30dB左右，比直接采用样本单元回波信号用于自适应权值约束的方法干扰对消性能提高了10dB左右。

Claims (8)

1.一种基于旁瓣约束的天波超视距雷达自适应干扰抑制方法，包括如下过程：

A、对待检测单元和样本单元的数据分别进行波束形成，得到指向待检测方位θ₀的接收波束，其输出作为主天线接收数据；

B、在接收阵列中选取N个阵元作为辅助天线，选定的N个辅助天线所对应的待检测单元和样本单元的数据作为辅助天线的接收数据；

C、对样本单元的主天线数据和各个辅助天线数据分别进行3脉冲MTI处理，获得仅含有干扰信息的样本数据；

D、对待检测单元和样本单元主天线以及各个辅助天线的数据分别进行相同的时域分段；

E、对待检测单元各个辅助天线各段数据分别进行多普勒滤波，并保留零多普勒通道输出；

F、采用常规自适应滤波算法求解待检测单元各个辅助天线第一段数据对应的自适应权值；

G、在给定约束条件下，利用样本单元数据分别求解待检测单元各个辅助天线第k段数据对应的自适应权值，k＝2，3，…，K；

H、使用待检测单元各个辅助天线各段数据对应的自适应权值分别对待测单元各段数据进行自适应滤波，获得干扰抑制后的输出结果；

I、对自适应滤波后的输出结果进行相干积累，获得干扰抑制后的多普勒频谱。

2.根据权利要求1所述的方法，其特征在于步骤C所述的对样本单元的数据进行3脉冲MTI处理，采用如下公式(5)和(6)进行：

d_s(m，r)＝d_p(m，r)-2d_p(m+1，r)+d_p(m+2，r)       (5)

X_s(m，r)＝X_0p(m，r)-2X_0p(m+1，r)+X_0p(m+2，r)    (6)

式中，d_p(m，r)是样本单元主天线的数据；d_s(m，r)是样本单元主天线的数据经过3脉冲MTI后的输出；X_0p(m，r)是样本单元辅助天线的数据；X_s(m，r)是样本单元辅助天线的数据经过3脉冲MTI后的输出；m＝1，2，…，N_p-2，r＝1，2，…，N_r，m表示脉冲重复周期序号，r表示距离单元序号，N_p和N_r分别表示样本单元总脉冲重复周期数和总距离单元数；d_p(m，r)和d_s(m，r)的维数均为1×1，X_p(m，r)和X_s(m，r)的维数均为N×1，N为辅助天线个数。

3.根据权利要求1所述的方法，其特征在于步骤D所述的时域分段，是按200～300毫秒的时间确定每段数据，设每段数据包含M个脉冲的回波数据，共分为K段，即

N_p＝K×M    (7)

样本单元数据的最后一段包含M-2个脉冲的数据，其余各段包含M个脉冲的数据。

4.根据权利要求1所述的方法，其特征在于步骤E所述的保留待检测单元辅助天线各段数据零多普勒通道输出，是按如下公式(8)进行：

$X_c\left(k\right)=\underset{m=(k-1)\·M+1}{\overset{\mathrm{Mk}}{\mathrm{\Σ}}}{X}_0\left(m\right)e^{j2\mathrm{\π}\frac{m\·0}{M}}---\left(8\right)$

$\underset{m=(k-1)\·M+1}{\overset{\mathrm{Mk}}{\mathrm{\Σ}}}{X}_0\left(m\right)$

式中，X_c(k)表示待检测单元各个辅助天线第k段数据的零多普勒通道输出；k＝1，2，…，K；X₀(m)是N×1维列矢量，表示待检测单元辅助天线的数据。

5.根据权利要求1所述的方法，其特征在于步骤G按如下过程进行：

G1、给定待检测单元各个辅助天线第k段数据对应的自适应权值的约束条件为：

W^H(k-1)·X_c(k)＝W^H(k)X_c(k)    (13)

式中，W(k)表示第k段数据对应的自适应权值；X_c(k)为检测单元辅助天线第k段数据的零多普勒通道输出；k＝2，3，…，K；H表示复共轭转置；

G2、利用下式所给出的约束最小均方误差，求解第k段数据对应的自适应权值

$\left\{\begin{array}{c}\underset{w\left(k\right)}{\arg \min }E\left[{|d_s(m,r)-W^H\left(k\right)\·X_s(m,r)\·W^H\left(k\right)|}^2\right]\\ s.t.W^H(k-1)\·X_c\left(k\right)=W^H\left(k\right)X_c\left(k\right)\end{array}\right.---\left(14\right)$

式中，m＝(k-1)·M+1，…，kM，r＝1，2，…，N_r，k＝2，3，…，K；X_c(k)是待检测单元第k段数据的零多普勒通道输出；

G3、利用拉格阳乘数法，求出自适应权值为：

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right){(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}^{-1}---\left(20\right)$

$\·X_c^H\left(k\right)\·[R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-W(k-1)]$

式中，k＝2，3，…，K；

R_xx(k)是第k段数据辅助天线的协方差矩阵，按照下式由样本单元数据估计得到

${\hat{R}}_{\mathrm{xx}}\left(k\right)=\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{m=(k-1)m+1}{\overset{\mathrm{kM}}{\mathrm{\Σ}}}{X}_s(m,r)\·X_s^H(m,r)---\left(21\right)$

R_xd(k)是第k段数据辅助天线与主天线之间的互相关矢量，按照下式由样本单元数据估计得到

${\hat{R}}_{\mathrm{xd}}\left(k\right)=\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{m=(k-1)m+1}{\overset{\mathrm{kM}}{\mathrm{\Σ}}}{X}_s(m,r)\·d_s^H(m,r)---\left(22\right)$

式中，d_s(m，r)和X_s(m，r)分别表示经过3脉冲MTI处理后样本单元主天线和辅助天线的输出。

6.根据权利要求5所述的方法，其特征在于步骤G3所述的利用拉格阳乘数法，求出自适应权值进一步用如下公式(23)表示：

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-\frac{R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}[R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-W(k-1)]---\left(23\right)$

7.根据权利要求5所述的方法，其特征在于步骤G4所述的利用拉格阳乘数法，求出自适应权值进一步用如下公式(24)表示：

$W\left(k\right)=R_{\mathrm{xx}}^{-1}\left(k\right)[R_{\mathrm{xd}}\left(k\right)-\frac{X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}[R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)-W(k-1)\left]\right]---\left(24\right)$

8.根据权利要求5所述的方法，其特征在于步骤G4所述的利用拉格阳乘数法，求出自适应权值进一步用如下公式(25)表示：

$W\left(k\right)=[I-\frac{R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}]R_{\mathrm{xx}}^{-1}\left(k\right)R_{\mathrm{xd}}\left(k\right)---\left(25\right)$

$+\frac{R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right)\·X_c^H\left(k\right)}{(X_c^H\left(k\right)\·R_{\mathrm{xx}}^{-1}\left(k\right)\·X_c\left(k\right))}W(k-1)$

Images