CN102176027B - 高动态干扰自跟踪稳健抑制方法 - Google Patents

Abstract

本发明公开了一种高动态干扰自跟踪稳健抑制方法，主要解决常规抗干扰方法对动态干扰抑制效果不佳的问题。其实现过程为：利用两重循环实现动态干扰抑制，内循环完成动态干扰的特征跟踪，用t-1时刻第i个干扰特征向量w_i(t-1)和t时刻接收数据X(：，t)获得t时刻特征波束输出y_i(t)，然后用t时刻特征波束输出y_i(t)和t-1时刻第i个动态干扰的特征值d_i(t-1)更新t时刻的动态干扰特征值d_i(t)和t时刻干扰特征向量w_i(t)，再用t时刻干扰特征向量w_i(t)和特征波束输出y_i(t)更新t时刻接收数据矢量X(：，t)；外循环用所有动态干扰特征向量w_i(t)组成加宽动态干扰特征矩阵F，然后估计出抗干扰权矢量h，再用抗干扰权矢量h对接收数据加权完成动态干扰抑制。本发明干扰抑制效果明显好，其输出信干噪比更接近理论最优值，可用于对空间动态干扰的抑制。

Description

高动态干扰自跟踪稳健抑制方法

技术领域

本发明属于信号处理技术领域，涉及阵列信号处理及自适应滤波技术，可以用于干扰在空间高动态运动时自跟踪抑制干扰同时保证信号稳定输出。

背景技术

卫星导航定位***可全天候、全天时提供位置、速度以及时间等信息，在民用和国防领域有着广阔的应用需求。但是，由于电磁环境日益恶劣并且到达接收机的导航信号非常微弱，使得导航接收机对干扰极其敏感，这种对干扰敏感的脆弱性严重影响导航的性能。实际中存在的高动态干扰抑制是自适应抗干扰面临的主要难题之一。这里高动态干扰是指干扰源与接收机存在相对运动，导致干扰源到达接收机的来波方向是快变的，例如干扰本身快速运动如安置于无人机上，而接收机不动，或者接收机高速运动如安置于飞机、导弹或卫星上而干扰不动，或者两者皆有运动。常规的抗干扰方法如采样协方差矩阵求逆SMI方法、最小均方误差LMS方法、递归最小二乘RLS方法等都适用于干扰的来波方向不随时间变化的静态干扰抑制情况。它们通过条件约束或迭代步骤去掉干扰对信号产生的影响，从而实现了在静态干扰背景下提取信号的目的。但是针对高动态干扰抑制，以上方法已不再适用，因为这些方法不能有效地跟踪动态干扰的特征，当干扰的来波方向快速变化时，使用这些方法进行干扰抑制无法实时在方向图的干扰位置形成零陷，即不具有零陷自跟踪性能，从而使阵列无法区分信号和动态干扰，导致高动态干扰抑制效果不佳。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种高动态干扰自跟踪稳健抑制方法，以快速跟踪干扰方向使方向图的凹口位置实时地对准干扰来波方向，并且保证主瓣对准信号来波方向，从而实现了高动态干扰抑制同时保证阵列稳定输出。

为实现上述目的，本发明的技术方案，包括如下步骤：

1.对天线阵列数字采样，得到一批包含信号和动态干扰的接收数据X，其中X为M×N维矩阵，M是阵列阵元个数，N是数据采样数；

2.采用内、外两重循环对接收数据X中动态干扰进行抑制：

1)设外循环变量t＝1，t表示采样时刻，初始化第i个动态干扰的特征值：d_i(t-1)＝0，i＝1，2，...，r，r是动态干扰数目，初始化第i个动态干扰的特征向量w_i(t-1)为除第i个元素为1以外其余元素为0的M×1维列矢量；

2)设内循环变量i＝1，i表示动态干扰序号；

3)用t-1时刻第i个干扰特征向量w_i(t-1)和t时刻接收数据X(：，t)获得t时刻特征波束输出：

其中H表示共轭转置；

4)用t时刻特征波束输出y_i(t)和t-1时刻第i个动态干扰的特征值d_i(t-1)更新t时刻的动态干扰特征值：d_i(t)＝(1-α+αβ)d_i(t-1)+α|y_i(t)|²，其中α为平滑因子，β为遗忘因子，用t时刻的动态干扰特征值d_i(t)更新t时刻动态干扰特征向量：

$w_i\left(t\right)=w_i(t-1)+\mathrm{\α}[X(:,t)-w_i(t-1)y_i\left(t\right)]\frac{y_i^{*}\left(t\right)}{d_i\left(t\right)},$ 其中*表示共轭；

5)用t时刻干扰特征向量w_i(t)更新t时刻接收数据矢量X(：，t)＝X(：，t)-w_i(t)y_i(t)，更新内循环变量i＝i+1，若更新后的内循环变量i小于等于动态干扰数目r，则返回至步骤3)继续执行内循环，否则执行步骤6)，完成t时刻对r个干扰的特征跟踪；

6)当外循环变量t能够整除分块处理的采样数T时，执行步骤7)至9)，否则跳至步骤10)；

7)利用r个动态干扰特征向量w_i(t)组成加宽动态干扰特征矩阵：F＝[w₁(t)，Bw₁(t)，...，w_r(t)，Bw_r(t)]，B是用各个阵元坐标形成的对角矩阵：

B＝diag(dx₁，dx₂，...，dx_M)，其中dx_m，m＝1，2，...，M，表示第m个阵元坐标；

8)根据加宽动态干扰特征矩阵估计出抗动态干扰权矢量：

$h={\mathrm{Pa}}_s/\sqrt{\left({a_s}^H{\mathrm{Pa}}_s\right)}$

其中P是由加宽干扰特征矩阵构造的正交补空间：P＝I_M-F(F^HF)^-1F^H，I_M是M×M维单位矩阵，a_s是M×1维信号导向矢量；

9)由抗动态干扰权矢量h得到抗动态干扰输出：

z＝h^HX(：，t+T+1：t+2T)；

10)更新外循环变量t＝t+1，若更新后的外循环变量t小于等于数据采样数N则返回至步骤2)，否则结束外循环，此时完成对接收数据X中动态干扰的抑制。

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

(a)本发明由于利用了平滑处理更新t时刻的动态干扰特征值d_i(t)和相应的特征向量w_i(t)，保证了方向图的主瓣对准信号的来波方向，凹口实时地对准动态干扰的来波方向，使之对干扰的特征跟踪更加准确。

(b)本发明由于用干扰特征向量w_i(t)和对角矩阵B组成动态干扰特征矩阵F，它加宽了干扰在空域位置上的零陷，使之对动态干扰抑制效果更加明显。

仿真结果表明：本发明与采样协方差矩阵求逆SMI方法相比，由于本发明的以上优点提高了阵列输出信干噪比，使其抑制动态干扰具有更加稳健的性能。

附图说明

图1是本发明实施例所用的天线阵列排布图；

图2是本发明方法的流程图；

图3是接收数据X中存在一个动态干扰时的方向图；

图4是用本发明与现有SMI方法分别对图3抑制后的信号频谱图；

图5是用本发明与现有SMI方法分别对图3抑制后的输出信干噪比图；

图6是接收数据X中存在两个动态干扰时的方向图；

图7是用本发明与现有SMI方法分别对图6抑制后的信号频谱图；

图8是用本发明与现有SMI方法分别对图6抑制后的输出信干噪比图。

具体实施方式

参照阵列模型图1和流程图2，本发明实现的步骤如下：

步骤1.对天线阵列数字采样，得到一批包含信号和动态干扰的接收数据X＝a_sS₀+a_iS_i+n，设信号的来波方向提前通过测向方法已知，其中X为M×N维矩阵，M是阵列阵元个数，N是数据采样数，S₀为接收数据中信号组成的矢量，a_s是信号导向矢量，S_i为接收数据中干扰组成的矢量，a_i为动态干扰导向矢量，n为接收数据中的噪声。

步骤2.采用内、外两重循环抑制动态干扰，内循环完成对接收数据X中r个干扰的特征跟踪，外循环利用这些干扰的特征进行动态干扰抑制，这里r是动态干扰数目：

(2.1)设外循环变量t＝1，t表示采样时刻，初始化第i个动态干扰的特征值：d_i(t-1)＝0，i＝1，2，...，r，初始化第i个动态干扰的特征向量w_i(t-1)为除第i个元素为1以外其余元素为0的M×1维列矢量；

(2.2)设内循环变量i＝1，i表示动态干扰序号；

(2.3)用t-1时刻第i个干扰特征向量w_i(t-1)和t时刻接收数据X(：，t)获得t时刻特征波束输出：其中H表示共轭转置；

(2.4)用t时刻特征波束输出y_i(t)和t-1时刻第i个动态干扰的特征值d_i(t-1)更新t时刻的动态干扰特征值：d_i(t)＝(1-α+αβ)d_i(t-1)+α|y_i(t)|²，其中α为平滑因子，β为遗忘因子，用t时刻的动态干扰特征值d_i(t)更新t时刻干扰特征向量：

其中*表示共轭；

(2.5)用t时刻干扰特征向量w_i(t)更新t时刻接收数据矢量X(：，t)＝X(：，t)-w_i(t)y_i(t)，更新内循环变量i＝i+1，若更新后的内循环变量i小于等于动态干扰数目r，则返回至步骤(2.3)继续内循环，否则执行步骤(2.6)，完成t时刻对r个干扰的特征跟踪；

(2.6)采用分快处理数据的方法，当外循环变量t能够整除分块采样数T时，执行步骤(2.7)至(2.9)，否则跳至步骤(2.10)；

(2.7)用r个动态干扰特征向量w_i(t)组成加宽动态干扰特征矩阵：F＝[w₁(t)，Bw₁(t)，...，w_r(t)，Bw_r(t)]，B是用各个阵元坐标形成的对角矩阵：B＝diag(dx₁，dx₂，...，dx_M)，其中dx_m，m＝1，2，...，M，表示第m个阵元坐标；

(2.8)根据加宽动态干扰特征矩阵估计出抗动态干扰权矢量：

$h={\mathrm{Pa}}_s/\sqrt{\left({a_s}^H{\mathrm{Pa}}_s\right)}$

其中P是由加宽干扰特征矩阵构造的正交补空间：P＝I_M-F(F^HF)^-1F^H，I_M是M×M维单位矩阵，a_s是信号导向矢量；

(2.9)由抗动态干扰权矢量h得到抗动态干扰输出Z：

Z＝h^HX(：，t+T+1：t+2T)

其中t是外循环变量，T是分块采样数；

(2.10)更新外循环变量t＝t+1，若更新后的外循环变量t小于等于数据采样数N则返回至步骤(2.2)，否则结束外循环，此时完成对接收数据X中动态干扰的抑制。

(1)实验环境

本发明采用n＝8的一个线性均匀传感器阵列，信号载频为F_c＝10⁹Hz，相邻传感器元素之间的距离是半波长λ/2，数据采样数N为1024，分块处理的采样数T为16，平滑因子α为0.5，遗忘因子β为0.92。

(2)实验内容与结果

实验1，对接收数据X中存在一个干噪比为40dB的动态干扰，其来波方向θ₁(t)＝20+0.01t°/s，接收数据X中信号的信噪比为0dB，其来波方向为0°的接收数据利用平滑处理更新t时刻的动态干扰特征值d_i(t)和相应的特征向量w_i(t)，得出该动态干扰的方向图并与采用SMI方法得出的相应方向图比较，其结果如图3所示；然后每隔15次数据采样次数估计1次抗干扰权矢量h，再用抗干扰权矢量h作用于阵列接收数据X得出对接收信号X干扰抑制后的信号频谱图，并与采用SMI方法抑制该干扰的信号频谱图对比，其结果如图4所示；最后分别用本发明方法和SMI方法统计抑制干扰后阵列的输出信干噪比，其结果如图5所示。

由图3可以看出，使用本发明可以使方向图的主瓣对准信号来波方向0°，并在干扰方向产生至少-60dB的零陷，而SMI方法使得方向图主瓣偏移，干扰的零陷最多达-40dB。这说明了本发明保证了方向图的主瓣对准信号的来波方向，并能在干扰方向上产生较深的零陷。

从图4可以看出，采用本发明方法进行动态干扰抑制之后，信号幅度为68dB，干扰幅度为31dB，二者之差为37dB，而采用SMI方法进行干扰抑制之后，信号幅度为60dB，干扰幅度为48dB，二者之差为12dB。这说明采用本发明方法抑制高动态干扰后信干比相对增大。从图4上还可以看出，用SMI方法进行干扰抑制之后干扰频谱周围出现未抗掉的频线分量，说明本发明方法抗干扰效果更好。

图5进一步说明了该方法的优越性。这里显示两种方法相应的输出信干噪比。可以看出本发明方法的输出信干噪比平均在8.5dB左右，接近理论最优值10×log10(8)≈9dB，而SMI方法的输出信干噪比平均在-1dB左右，这说明了本方法具有更加稳定的输出。

实验2，对接收数据X中存在两个动态干扰，其中一个是上述干扰θ₁(t)，另一个是来波方向为θ₂(t)＝5sin(πt/500)-30°/s、干噪比为30dB的动态干扰，接收数据X中信号的信噪比仍为0dB，其来波方向为0°的接收数据沿用上述处理方法先利用平滑处理更新t时刻的动态干扰特征值d_i(t)和相应的特征向量w_i(t)，得出这两个动态干扰的方向图并与采用SMI方法得出的相应方向图比较，其结果如图6所示；然后每隔15次数据采样次数估计1次抗干扰权矢量h，再用抗干扰权矢量h作用于阵列接收数据X得出对接收信号X干扰抑制后的信号频谱图，并与采用SMI方法抑制这两个干扰的信号频谱图对比，其结果如图7所示；最后分别用本发明方法和SMI方法统计抑制干扰后阵列的输出信干噪比，其结果如图8所示。

从图6可以看出，本发明方法能分别在各个干扰方向上形成-50dB和-60dB的零陷，同时使主瓣对准来波方向。

由图7可知，在抑制两个干扰时SMI方法抗干扰效果不及本发明方法，因为从频谱图上看出两个干扰都存在未被抑制掉的分量，而本方法则不存在该现象。

从图8中看出输出信干噪比平均在8dB左右，而SMI方法的输出信干噪比平均在-3dB。

综上所述，采用平滑处理和加宽零陷的方法抑制高动态运动干扰能够使方向图的凹口实时对准动态干扰的来波方向，同时保证主瓣对准信号的来波方向，相比SMI方法干扰抑制效果明显，输出信干噪比接近理论最优值。

Claims (1)

1.一种高动态干扰自跟踪稳健抑制方法，包括如下步骤：

（1）对天线阵列数字采样，得到一批包含信号和动态干扰的接收数据X，其中X为M×N维矩阵，M是阵列阵元个数，N是数据采样数；

（2）采用内、外两重循环对接收数据X中动态干扰进行抑制：

（2a）设外循环变量t＝1，t表示采样时刻，初始化第i个动态干扰的特征值：d_i(t-1)＝0，i＝1,2,...,r，r是动态干扰数目，初始化第i个动态干扰的特征向量w_i(t-1)为除第i个元素为1以外其余元素为0的M×1维列矢量；

（2b）设内循环变量i＝1，i表示动态干扰序号；

（2c）用t-1时刻第i个干扰特征向量w_i(t-1)和t时刻接收数据X(:,t)获得t时刻特征波束输出：其中H表示共轭转置；

（2d）用t时刻特征波束输出y_i(t)和t-1时刻第i个动态干扰的特征值d_i(t-1)更新t时刻的动态干扰特征值：d_i(t)＝(1-α+αβ)d_i(t-1)+α|y_i(t)|²，其中α为平滑因子，β为遗忘因子，用t时刻的动态干扰特征值d_i(t)更新t时刻干扰特征向量： $w_i\left(t\right)=w_i(t-1)+\mathrm{\α}[X(:,t)-w_i(t-1)y_i\left(t\right)]\frac{y_i^{*}\left(t\right)}{d_i\left(t\right)},$ 其中*表示共轭；

（2e）用t时刻干扰特征向量w_i(t)更新t时刻接收数据矢量X(:,t)＝X(:,t)-w_i(t)y_i(t)，更新内循环变量i＝i+1，若更新后的内循环变量i小于等于动态干扰数目r，则返回至步骤（2c）继续内循环，否则执行步骤（2f），完成t时刻对r个干扰的特征跟踪；

（2f）当外循环变量t能够整除分块处理的采样数T时，执行步骤（2g）-（2i），否则跳至步骤（2j）；

（2g）利用r个动态干扰特征向量w_i(t)组成加宽动态干扰特征矩阵：F=[w₁(t),Bw₁(t),...,w_r(t),Bw_r(t)],B是用各个阵元坐标形成的对角矩阵：B=diag(dx₁,dx₂,...,dx_M)，其中dx_m，m＝1,2,...,M，表示第m个阵元坐标；

（2h）根据加宽动态干扰特征矩阵估计出抗动态干扰权矢量：

$h={\mathrm{Pa}}_s/\√\sqrt{\left({a_s}^H{\mathrm{Pa}}_s\right)}$

其中P是由加宽干扰特征矩阵构造的正交补空间：P＝I_M-F(F^HF)^-1F^H，I_M是M×M维单位矩阵，a_s是M×1维信号导向矢量，(·)^-1表示矩阵求逆运算符；

（2i）由抗动态干扰权矢量h得到抗动态干扰输出：

z＝h^HX(:，t+T+1:t+2T)

（2j）更新外循环变量t＝t+1，若更新后的外循环变量t小于等于数据采样数N则返回至步骤（2b），否则结束外循环，此时完成对接收数据X中动态干扰的抑制。

Images