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

具体实施方式

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

天波超视距雷达采用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}}{\Sigma }}{X}_0\left(m\right)e^{j2\pi \frac{m·0}{M}}---\left(8\right)$

$=\underset{m=(k-1)·M+1}{\overset{\mathrm{Mk}}{\Sigma }}{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}{\Sigma }}\underset{m=1}{\overset{M}{\Sigma }}{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}{\Sigma }}\underset{m=1}{\overset{M}{\Sigma }}{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)={\sigma }_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))·\lambda$

求代价函数关于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)·\lambda )---\left(17\right)$

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

$W^H\left(k\right)·X_c\left(k\right)=(R_{\mathrm{xd}}^H\left(k\right)-{\lambda }^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)$

进一步可得：

$\lambda ={(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}{\Sigma }}\underset{m=(k-1)m+1}{\overset{\mathrm{kM}}{\Sigma }}{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}{\Sigma }}\underset{m=(k-1)m+1}{\overset{\mathrm{kM}}{\Sigma }}{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(\omega \right)=\underset{m=1}{\overset{N_p}{\Sigma }}{y}_0\left(m\right)w\left(m\right)e^{j2\pi \frac{\mathrm{m\omega }}{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左右。