CN103630910A - 一种gnss接收机设备的抗干扰方法 - Google Patents

Abstract

本发明提出一种GNSS接收机设备的抗干扰方法，本方法采用空间谱估计MUSIC算法，抗干扰效果性能良好。不同于传统的功率倒置算法，在干扰源附近有带状零陷分布；本方法能准确找到干扰的方向并形成有效抑制，当信干比降到20dB时，本方法依然能起到干扰抑制的效果。

Description

一种GNSS接收机设备的抗干扰方法

技术领域

本发明属于GNSS工程安全技术领域，特别涉及一种GNSS接收机设备的抗干扰方法；是一种波束形成的抗干扰波束形成算法，以防止GNSS接收端受到外界的连续干扰。

背景技术

卫星通信系统的信道非常脆弱，容易受到干扰和欺骗。GNSS抗干扰技术的目的是尽最大的努力抑制外来电磁能和电磁频谱对GNSS接收机的攻击，以有效的措施保障卫星通信的安全。总之，通信卫星承担着提供通信保障的任务，必须具有较强的抗干扰能力，因此深入广泛地研究抗干扰技术，提高它的抗干扰能力，具有非常重要的意义。

现有抗干扰技术有天线抗干扰技术、扩频技术、星上处理技术、自适应编码调制技术、扩展频段技术、无线光通信技术和限幅技术。

现有天线抗干扰技术中，天线自适应调零的原理是:在星上采用大型的、具有多波束的接收天线，组成一个赋形的天线照射到某一区域，当卫星检测到干扰时，自动将干扰方向的点波束关闭，从而达到抗干扰的目的。具体是将天线作为空间滤波手段，利用目标到卫星的已知几何关系，并根据已知的或可以推测的干扰源的位置，形成带指向的波束。

目前工程中使用功率倒置算法是在LCMV准则下的GPS接收机时空最小化方案，具体如下：采用线性约束最小方差（LCMV）准则，在保证期望信号不衰减的情况下，使阵列输出的方差（即输出功率）最小化。它的方程式为：

$\left\{\begin{array}{c}\underset{w}{\arg \min }E\left\{{\left|y\left(n\right)\right|}^2\right\}=w^H{R}_{\mathrm{xx}}w\\ s.t.w^{H}s=1\end{array}\right.$

式中，R_xx=E[x(n)x^H(n)]为接收数据协方差矩阵，w=[w₁,…,w_M]^T为加权矢量，s是一个给定的列向量为s=[1,0,…,0]^T。利用拉格朗日方程，可以很容易求得最优解为

$w_{\mathrm{opt}}={\left(s^H{R}_{\mathrm{xx}}^{-1}s\right)}^{-1}{R}_{\mathrm{xx}}^{-1}s$

该方法分辨率较低，并在干扰源方向的零陷周围形成大量的带状零陷，而且当信干噪比（SINR）下降时，性能劣化严重，甚至会造成干扰误判问题。

发明内容

针对背景技术存在的问题，本发明提供一种GNSS接收机设备的抗干扰方法，属于一种天线抗干扰技术。本发明采用MUSIC算法对传统的自适应天线零陷技术进行改进，应用到GPS接收机抗干扰中；由于MUSIC算法具备超分辨特性，利用该算法估计GNSS接收机的强干扰方向，及实现DOA的准确估计，然后在功率输出最小的方向上形成抗干扰的零陷。

为了规避在干扰源以外的方向上产生抑制多余的零陷，本发明采用如下技术方案：

一种GNSS接收机设备的抗干扰方法，将MUSIC算法应用到GNSS抗干扰领域，具体包括以下步骤，

步骤1、建立阵列干扰信号模型，接收阵列信号数据；

步骤2、由接收到的阵列信号数据得到数据协方差矩阵R_xx；

步骤3、对R_xx进行特征分解；

步骤4、对R_xx的特征值进行信号源数判断；

步骤5、确定信号子空间与噪声子空间；即，特征值小的形成噪声子空间，特征值大的是信号子空间；

步骤6、根据信号参数范围进行谱峰搜索；

信号信号参数范围：方位角0度-360度，仰角0度-90度；

步骤7、寻找信号分量的导向矢量与噪声子空间特征向量正交位置，即谱峰位置；

步骤8、仰角、方位角与干扰信号方向增益的关系式如下：

a(θ,φ)表示信号分量的导向矢量，V_n表示噪声子空间的特征向量

步骤9、在对应的位置生成零陷抑制，得到干扰抑制后的数据y(n)；

$y\left(n\right)=w_{\mathrm{opt}}^{H}x\left(n\right)\≈\frac{v}{{\mathrm{\σ}}^2}{s}^H{P}_{A_J}^{\⊥}x\left(n\right)$

其中，

$w_{\mathrm{opt}}=v{R}^{-1}s\≈\frac{v}{{\mathrm{\σ}}^2}[I-A_J{\left(A_J^H{A}_J\right)}^{-1}{A}_J^H]s=\frac{v}{{\mathrm{\σ}}^2}{P}_{A_J}^{\⊥}s,$

x(n)为输入的离散信号，A_J=[a₁(θ,φ),…,a_Q(θ,φ)]为a(θ,φ)的矢量表示，R_J=diag(V_L+1,...,V_Q)为V_n的矢量表示，即干扰的功率；σ²为高斯白噪声的功率，v=(s^HR^-1s)^-1为一常数，R=E[x(n)x^H(n)]为接收数据协方差矩阵；s=[1,0,…,0]^T为给定参考列向量；

表示与由A_J构造的干扰子空间正交的子空间投影矩阵；y(n)实际上为快拍数据投影后列向量

的第一个元素。

由于通过MUSIC算法获得干扰信号的来向，就可以利用正交投影的方法来抑制干扰。

所述步骤1的实现过程如下，

设定干扰信号模型满足以下四个条件；

1)阵元的间距要远大于阵元的尺寸，且不考虑阵元之间存在耦合情况；

2)把单位阵元当作点阵元来考虑；

3)系统噪声为均值为0，方差为δ²的高斯白噪声，并且每个阵元之间噪声是独立的；

4)每个阵元均为全方向天线，把接收信号当成是平面波；

依照上述设定，构造信号模型；假设空间中一组阵列天线的阵元数是M，且可以接收到的信号数量为N；设定参考阵元为第一个阵元，若第i个信号到天线阵列的入射角度为θ_i且信号定义为S_i(t)e^jωt，那么，第m个阵元所接收到的信号可以表示为：

$x_m\left(t\right)=a_{\mathrm{mi}}S(t+(m-1){\mathrm{\τ}}_i)e^{\mathrm{j\ω}(t+(m-1){\mathrm{\τ}}_i)},m=\mathrm{1,2},\·\·\·,M---\left(a1\right)$

式(a1)中， a_mi是第m个阵元对第i个信号源的幅度响应，每个阵元均为全向天线，a_mi为信号的振幅； τ_i是当入射信号以θ_i角度到达阵列相邻阵元的时延；ω是入射信号的角频率；假设辐射源S_i(t)为窄带信号：

S_i(t)≈S_i(t+(m-1)τ_i),m=1,2,…,M  (a2)

如果相邻阵元之间的相位差为

，可以得到阵列中每个阵元的输出矢量形式：

$\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \·\\ \·\\ \·\\ x_M\left(t\right)\end{array}\right]=\left[\begin{array}{c}{e}^{j(\mathrm{\ωt}+(1-1){\mathrm{\φ}}_i)}\\ e^{j(\mathrm{\ωt}+(2-1){\mathrm{\φ}}_i)}\\ \·\\ \·\\ \·\\ e^{j(\mathrm{\ωt}+(M-1){\mathrm{\φ}}_i)}\end{array}\right]*S_i\left(t\right)---\left(a3\right)$

可化简为：

$\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \·\\ \·\\ \·\\ x_M\left(t\right)\end{array}\right]=\left[\begin{array}{c}1\\ e^{j{\mathrm{\φ}}_i}\\ \·\\ \·\\ \·\\ e^{j(M-1){\mathrm{\φ}}_i}\end{array}\right]*S_i\left(t\right)e^{\mathrm{jwt}}---\left(a4\right)$

则式（a4）可以记为：

X(t)=a(φ_i)*S_i(t)e^jwt   (a5)

其中，X(t)=[x₁(t),x₂(t),…,x_M(t)]^Ｔ，

$a\left({\mathrm{\φ}}_i\right)={[1,e^{j{\mathrm{\φ}}_i},\·\·\·,e^{j(M-1){\mathrm{\φ}}_i}]}^T;$

a(φ_i)为导向向量，信源个数为N，接收到的阵列信号是：

X(t)=a(Φ₁)*S₁(t)e^jwt+a(Φ₂)*S₂(t)e^jwt+…+a(Φ_N)*S_N(t)e^jwt  (a6)

式（a6）表示为矩阵形式：

$X\left(t\right)=\left[\begin{array}{cccc}1& 1& \·\·\·& 1\\ e^{j{\mathrm{\Φ}}_1}& e^{j{\mathrm{\Φ}}_2}& \·\·\·& e^{j{\mathrm{\Φ}}_N}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ e^{j(M-1){\mathrm{\Φ}}_1}& e^{j(M-1){\mathrm{\Φ}}_2}& \·\·\·& e^{j(M-1){\mathrm{\Φ}}_N}\end{array}\right]*\left[\begin{array}{c}{S}_1\left(t\right)\\ S_2\left(t\right)\\ \·\\ \·\\ \·\\ S_N\left(t\right)\end{array}\right]e^{\mathrm{jwt}}---\left(a7\right).$

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

1、本发明采用的基于MUSIC算法的抗干扰防范在能够分辨GNSS接收机的强干扰方向，可以在GNSS抗干扰领域得到成功应用。

2、与传统的功率倒置算法相比，本发明采用的MUSIC算法的抑制能力更明显，

低信噪比的条件下优势明显，而且在干扰方向形成更深的零陷。

附图说明

图1为阵列天线接收信号几何关系；

图2为现有GPS接收机的原理框图；

图3为信干比为80dB时的平面均匀中心圆阵方向谱的立体图；

图4为信干比为80dB时的平面均匀中心圆阵方向谱的俯视图；

图5为信干比为20dB时的平面均匀中心圆阵方向谱的立体图；

图6为信干比为20dB时的平面均匀中心圆阵方向谱的俯视图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

本发明是在现有GPS接收机上实现的，现有接收机的原理如下：

GPS接收机是GPS导航卫星的用户设备，是实现GPS卫星导航定位的终端仪器。GPS接收机接收到的GPS卫星信号，经相关处理就能测定由卫星到接收机的信号传播时间延迟，并解算出接收机到卫星之间的距离（因含误差而被称为伪距），从而确定GPS接收机的位置和时间等信息。图2给出了一个通用GPS接收机的实现框图。它主要由天线、射频前端、中频数据采集和信号处理四部分组成。

硬件实现上面介绍的天线、射频前端，中频数据采集和信号处理四部分在一块基带板上进行处理。在基带板上，FPGA芯片是整个系统实现的核心器件，FPGA芯片包括数字混频模块、低通滤波模块、干扰抑制模块、GPS接收模块；FPGA芯片它除了实现数字混频、低通滤波、干扰抑制，以及GPS接收模块中的捕获、跟踪、解算外，还控制着A/D、D/A、串口等外围设备。

干扰抑制模块实现的功能主要包括数字下变频和干扰抑制。在射频前端的接收芯片GP2015的中频信号为4.309MHz，其内部本身的数字化处理为带通采样，采样率为5.714MHz，采样后的数字中频为1.405MHz。采用本发明的方法之前，先对模拟中频信号进行低通采样，采样率为22.857MHz，采样后，中频仍为4.309MHz，在干扰抑制之前首先需要进行数字下变频，并将单路的信号进行I/Q分离。

本发明的方法在FPGA内部实现，本发明的原理如下：

1、假设空间任意两个阵元，把其中一个作为参考阵元（位于坐标原点），另一个阵元位于坐标(x,y,z)，信号相对于坐标系的入射仰角为φ，方位角为θ，则两阵元间的波程差为，阵列天线接收信号几何关系如图1所示：

$\mathrm{\τ}=\frac{1}{c}(x\cos \mathrm{\θ}\cos \mathrm{\φ}+y\sin \mathrm{\θ}\cos \mathrm{\φ}+z\sin \mathrm{\φ})---\left(b1\right)$

τ是信号延时，c是波速。当入射信号为窄带信号时，忽略其包络的影响，并不考虑阵列中各个阵元的通道不一致和互耦等因素的影响，可以将M个阵元在n时刻的输出表示为如下的列矢量：

$\left[\begin{array}{c}{x}_1\left(n\right)\\ x_2\left(n\right)\\ \·\\ \·\\ \·\\ x_M\left(n\right)\end{array}\right]=\left[\begin{array}{cccc}{e}^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{11}}& e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{12}}& \·\·\·& e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{1N}}\\ e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{21}}& e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{22}}& \·\·\·& e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{2N}}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{M1}}& e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{M2}}& \·\·\·& e^{-j{\mathrm{\ω}}_0{\mathrm{\τ}}_{\mathrm{MN}}}\end{array}\right]\left[\begin{array}{c}{s}_1\left(n\right)\\ s_2\left(n\right)\\ \·\\ \·\\ \·\\ s_N\left(n\right)\end{array}\right]+\left[\begin{array}{c}{w}_1\left(n\right)\\ w_2\left(n\right)\\ \·\\ \·\\ \·\\ w_M\left(n\right)\end{array}\right]---\left(b2\right)$

其中，s_n(n)为第n个入射信号，w_m(n)为第m个阵元上的噪声，ω₀为入射信号的角频率，τ_mn为第n个信号到达第m个阵元时相对参考阵元的时延，上式可以简化为：

X(n)=AS(n)+W(n)   （b3）

式中，A称为导向矢量矩阵。

MUSIC算法是一种基于特征结构子空间的超分辨方法，是由Schmidt在1979年提出的，将其应用于GPS抗干扰是十分有效的。MUSIC算法可以做到角度分辨率的提高。令天线阵列接收协方差矩阵为：

$R_{\mathrm{xx}}=E\left\{X\left(n\right)X^H\left(n\right)\right\}={\mathrm{AR}}_{\mathrm{ss}}{A}^H+{\mathrm{\σ}}_n^{2}I---\left(b4\right)$

I是M×N单位阵，由于强干扰信号互不相关，所以AR_ssA^H是满秩矩阵，式中R_ss是干扰信号的自相关矩阵，定义为R_ss=E[s·s^H]，

为噪声方差。

令R_xx的特征值为{λ₀,…,λ_i,…,λ_M-1}，使得

|R_xx-λ_iI|=0   (b5)

利用式（b4），可以把它改写为

$|{\mathrm{AR}}_{\mathrm{ss}}{A}^H+{\mathrm{\σ}}_n^{2}I-{\mathrm{\λ}}_{i}I|=|{\mathrm{AR}}_{\mathrm{ss}}{A}^H-({\mathrm{\λ}}_i-{\mathrm{\σ}}_n^2)I|=0---\left(b6\right)$

则 AR_ssA^H的特征值V_i为 $v_i={\mathrm{\λ}}_i-{\mathrm{\σ}}_n^2---\left(b7\right)$

当入射干扰信号数D小于阵元数M时，M×M的矩阵AR_ssA^H是半正定的，且秩为D，这意味着AR_ssA^H的特征值ν_i中，有M-D个特征值为零，即：

$v_0,v_1,\·\·\·,v_{D-1}>v_i={\mathrm{\λ}}_i-{\mathrm{\σ}}_n^2=0---\left(b8\right)$

其中i=D,D+1,…,M-1

由此可知，R_xx中有N-D个特征值等于噪声方差

因此对于R_xx的特征值有：

${\mathrm{\λ}}_0,{\mathrm{\λ}}_1\·\·\·,{\mathrm{\λ}}_{D-1}>{\mathrm{\λ}}_D={\mathrm{\λ}}_{D+1}=\·\·\·={\mathrm{\λ}}_{M-1}={\mathrm{\σ}}_n^2---\left(b9\right)$

假设关于特征值 λ_i的特征向量为 q_i，满足：

（R_xx-λ_iI） q_i=0        （b10）

对于

个最小特征值

相关的特征向量，可推出：

$(R_{\mathrm{xx}}-{\mathrm{\λ}}_{i}I)q_i={\mathrm{AR}}_{\mathrm{ss}}{A}^H{q}_i+{\mathrm{\σ}}_n^{2}I{q}_i-{\mathrm{\σ}}_n^{2}I{q}_i=0---\left(b11\right)$

于是AR_ssA^H q_i=0     （b12）

因为A满秩，R_ss非奇异，所以

A^Hq_i=0      （b13）

其中

θ代表仰角，代表方位角。

这表明与

个最小特征值对应的特征向量张成的噪声子空间和

个导向矢量正交，即

构造一个包含噪声特征向量的矩阵

V_n={q_D,…,q_M-1}   （b15）

信号分量的导向矢量与噪声子空间特征向量正交，即

MUSIC谱中个最大谱峰为式（b16）的倒数，由下式实现波达方向的超分辨估计

当和V_n的正交性时输出功率最小，能够在干扰的来向形成零陷，此时式（b17）变为：

基于上述原理，本发明的步骤如下：

步骤1、建立阵列干扰信号模型，接收阵列信号数据；

本步骤的具体实现过程如下：

设定干扰信号模型满足以下四个条件；

1)阵元的间距要远大于阵元的尺寸，且不考虑阵元之间存在耦合情况；

2)把单位阵元当作点阵元来考虑；

3)系统噪声为均值为0，方差为δ²的高斯白噪声，并且每个阵元之间噪声是独立的；

4)每个阵元均为全方向天线，把接收信号当成是平面波；

依照上述设定，构造信号模型；假设空间中一组阵列天线的阵元数是M，且可以接收到的信号数量为N；设定参考阵元为第一个阵元，若第i个信号到天线阵列的入射角度为θ_i且信号定义为S_i(t)e^jωt，那么，第m个阵元所接收到的信号可以表示为：

$x_m\left(t\right)=a_{\mathrm{mi}}S(t+(m-1){\mathrm{\τ}}_i)e^{\mathrm{j\ω}(t+(m-1){\mathrm{\τ}}_i)},m=\mathrm{1,2},\·\·\·,M---\left(a1\right)$

式(a1)中， a_mi是第m个阵元对第i个信号源的幅度响应，每个阵元均为全向天线，a_mi为信号的振幅； τ_i是当入射信号以θ_i角度到达阵列相邻阵元的时延；ω是入射信号的角频率；假设辐射源S_i(t)为窄带信号：

S_i(t)≈S_i(t+(m-1)τ_i),m=1,2,…,M  (a2)

如果相邻阵元之间的相位差为

，可以得到阵列中每个阵元的输出矢量形式：

$\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \·\\ \·\\ \·\\ x_M\left(t\right)\end{array}\right]=\left[\begin{array}{c}{e}^{j(\mathrm{\ωt}+(1-1){\mathrm{\φ}}_i)}\\ e^{j(\mathrm{\ωt}+(2-1){\mathrm{\φ}}_i)}\\ \·\\ \·\\ \·\\ e^{j(\mathrm{\ωt}+(M-1){\mathrm{\φ}}_i)}\end{array}\right]*S_i\left(t\right)---\left(a3\right)$

可化简为：

$\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \·\\ \·\\ \·\\ x_M\left(t\right)\end{array}\right]=\left[\begin{array}{c}1\\ e^{j{\mathrm{\φ}}_i}\\ \·\\ \·\\ \·\\ e^{j(M-1){\mathrm{\φ}}_i}\end{array}\right]*S_i\left(t\right)e^{\mathrm{jwt}}---\left(a4\right)$

则式（a4）可以记为：

X(t)=a(φ_i)*S_i(t)e^jwt (a5)

其中，X(t)=[x₁(t),x₂(t),…,x_M(t)]^T，

$a\left({\mathrm{\φ}}_i\right)={1,e^{j{\mathrm{\φ}}_i},\·\·\·,e^{j(M-1){\mathrm{\φ}}_i}]}^T;$

a(φ_i)为导向向量，信源个数为N，接收到的阵列信号是：

X(t)=a(Φ₁)*S₁(t)e^jwt+a(Φ₂)*S₂(t)e^jwt+…+a(Φ_N)*S_N(t)e^jwt  (a6)

式（a6）表示为矩阵形式：

$X\left(t\right)=\left[\begin{array}{cccc}1& 1& \·\·\·& 1\\ e^{j{\mathrm{\Φ}}_1}& e^{j{\mathrm{\Φ}}_2}& \·\·\·& e^{j{\mathrm{\Φ}}_N}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ e^{j(M-1){\mathrm{\Φ}}_1}& e^{j(M-1){\mathrm{\Φ}}_2}& \·\·\·& e^{j(M-1){\mathrm{\Φ}}_N}\end{array}\right]*\left[\begin{array}{c}{S}_1\left(t\right)\\ S_2\left(t\right)\\ \·\\ \·\\ \·\\ S_N\left(t\right)\end{array}\right]e^{\mathrm{jwt}}---\left(a7\right).$

步骤2、由接收到的阵列信号数据得到数据协方差矩阵R_xx；

步骤3、对R_xx进行特征分解；

步骤4、对R_xx的特征值进行信号源数判断；

步骤5、确定信号子空间与噪声子空间；即，特征值小的形成噪声子空间，特征值大的是信号子空间；

步骤6、根据信号参数范围进行谱峰搜索；

信号信号参数范围：方位角0度-360度，仰角0度-90度；

步骤7、寻找信号分量的导向矢量与噪声子空间特征向量正交位置，即谱峰位置；

步骤8、仰角、方位角与干扰信号方向增益的关系式如下：

a(η,φ)表示信号分量的导向矢量，V_n表示噪声子空间的特征向量

步骤9、在对应的位置生成零陷抑制；

令A_J=[a_L+1(θ,φ),…,a_Q(θ,φ)]为a(θ,φ)的矢量表示，R_J=diag(V_L+1,...,V_Q)为V_n的矢量表示，即干扰的功率，σ²为高斯白噪声的功率，v=(s^HR^-1s)^-1为一常数，R=E[x(n)x^H(n)]为接收数据协方差矩阵；s=[1,0,…,0]^T为给定参考列向量，x(n)为输入的离散信号；

$w_{\mathrm{opt}}=v{R}^{-1}s\≈\frac{v}{{\mathrm{\σ}}^2}[I-A_J{\left(A_J^H{A}_J\right)}^{-1}{A}_J^H]s=\frac{v}{{\mathrm{\σ}}^2}{P}_{A_J}^{\⊥}s$

表示与由A_J构造的干扰子空间正交的子空间投影矩阵；x(n)为输入的离散信号；

干扰抑制后的数据y(n)实际上为快拍数据投影后列向量

的第一个元素。

由于通过MUSIC算法获得干扰信号的来向，就可以利用正交投影的方法来抑制干扰。

下面结合附图，详细说明本发明与现有的功率倒置算法技术相比，其具有的优点和有益效果：

1、为了验证本发明的正确性，进行了如下的仿真实验。期望卫星选取为C/A码，考虑天线阵的阵元数M为7，半径为半个波长。三个强窄带干扰为连续波干扰(CWI),入射仰角分别为30°、40°和70°，入射方位角分别为50°、220°、170°。此外假设白噪声为不相关高斯白噪声，平均功率为零。

2、实验1假设天线阵为平面均匀中心圆阵，取其中一个阵元放在圆心，其余6个均匀分布在半径为半个波长的圆上，输入信号的信噪比为-80dB，图3是平面均匀中心圆阵方向谱的立体图。x轴为仰角，y轴为方位角，z轴为方向谱增益。图的左边是采用MUSIC算法得到的方向谱增益，右图是采用功率倒置法得到的方向谱增益。从图中看出，它们都准确找到了对应的强窄带干扰的方向。同时MUSIC算法得到的抑制效果分别达到228.4dB、235.8dB、234.6dB。功率倒置算法得到的抑制效果分别达到139.1dB、138.1dB、142dB。这说明采用MUSIC算法，GNSS接收机的抗干扰性能得到显著提升。

3、图4是平面均匀中心圆阵方向谱的俯视图。图的左边采用MUSIC算法得到的零陷图，图中准确地找到三个强干扰信号的零陷，而图的右边是采用功率倒置算法得到的零陷图，图中除三个干扰的零陷点外，还出现了大量的随机带状分布的零陷。这说明MUSIC算法的滤波效果要优于功率倒置算法。

4、实验2同样假设天线阵为平面均匀中心圆阵，输入信号的信干比从原来的80dB降到20dB。图5的左边是采用MUSIC算法得到的方向谱增益，右图是采用功率倒置法得到的方向谱增益。从图5中看出，它们都准确找到了对应的强窄带干扰的方向。同时MUSIC算法得到的抑制效果分别达到103.1dB、121.9dB、87.33dB。而功率倒置算法得到的抑制分辨率低，无准确找到干扰的谱峰方向，误判率高。这说明采用MUSIC算法，GNSS接收机的抗干扰性能得到显著提升。

5、图6是平面均匀中心圆阵方向谱的俯视图。图的左边采用MUSIC算法得到的零陷图，图中依然能够准确地找到三个强干扰信号的零陷，与80dB条件下的零陷相比，零陷相对较宽，但深度变浅。而图的右边是采用功率倒置算法得到的零陷图，图中除三个干扰的零陷点外，还出现了大量的随机带状分布的零陷。并且大量的零陷连在一起，三个强干扰信号的零陷无法区别。

Claims (2)

1.一种GNSS接收机设备的抗干扰方法，其特征在于：将MUSIC算法应用到GNSS抗干扰领域，具体包括以下步骤，

步骤1、建立阵列干扰信号模型，接收阵列信号数据；

步骤2、由接收到的阵列信号数据得到数据协方差矩阵R_XX；

步骤3、对R_XX进行特征分解；

步骤4、对R_XX的特征值进行信号源数判断；

步骤5、确定信号子空间与噪声子空间；即，特征值小的形成噪声子空间，特征值大的是信号子空间；

步骤6、根据信号参数范围进行谱峰搜索；

信号参数范围：方位角0度--360度，仰角0度--90度；

步骤7、寻找信号分量的导向矢量与噪声子空间特征向量正交位置，即谱峰位置；

步骤8、仰角、方位角与干扰信号方向增益的关系式如下：

其中，a(θ,φ)表示信号分量的导向矢量，V_n表示噪声子空间的特征向量；

步骤9、在对应的位置生成零陷抑制，得到干扰抑制后的数据y(n)；

$y\left(n\right)=w_{\mathrm{opt}}^{H}x\left(n\right)\≈\frac{v}{{\mathrm{\σ}}^2}{s}^H{P}_{A_J}^{\⊥}x\left(n\right)$

其中，

$w_{\mathrm{opt}}=v{R}^{-1}s\≈\frac{v}{{\mathrm{\σ}}^2}[I-A_J{\left(A_J^H{A}_J\right)}^{-1}{A}_J^H]s=\frac{v}{{\mathrm{\σ}}^2}{P}_{A_J}^{\⊥}s,$

x(n)为输入的离散信号，A_J=[a₁(θ,φ),…,a_Q(θ,φ]为a(θ,φ)的矢量表示，R_J=diag(V_L+1,...,V_Q)为V_n的矢量表示，即干扰的功率；σ²为高斯白噪声的功率，v=(s^HR^-1s)^-1为一常数，R=E[x(n)x^H(n)]为接收数据协方差矩阵；s=[1,0,…,0]^T为给定参考列向量；

表示与由A_J构造的干扰子空间正交的子空间投影矩阵；y(n)实际上为快拍数据投影后列向量

的第一个元素。

2.根据权利要求1所述的一种GNSS接收机设备的抗干扰方法，其特征在于：

所述步骤1的实现过程如下，

设定干扰信号模型满足以下四个条件；

1)阵元的间距要远大于阵元的尺寸，且不考虑阵元之间存在耦合情况；

2)把单位阵元当作点阵元来考虑；

3)系统噪声为均值为0，方差为δ²的高斯白噪声，并且每个阵元之间噪声是独立的；

4)每个阵元均为全方向天线，把接收信号当成是平面波；

依照上述设定，构造信号模型；假设空间中一组阵列天线的阵元数是M，且可以接收到的信号数量为N；设定参考阵元为第一个阵元，若第i个信号到天线阵列的入射角度为θ_i且信号定义为S_i(t)e^jωt，那么，第m个阵元所接收到的信号可以表示为：

$x_m\left(t\right)=a_{\mathrm{mi}}S(t+(m-1){\mathrm{\τ}}_i)e^{\mathrm{j\ω}(t+(m-1){\mathrm{\τ}}_i)},m=\mathrm{1,2},\·\·\·,M---\left(a1\right)$

式(a1)中， a_mi是第m个阵元对第i个信号源的幅度响应，每个阵元均为全向天线，a_mi为信号的振幅； τ_i是当入射信号以θ_i角度到达阵列相邻阵元的时延；ω是入射信号的角频率；假设辐射源S_i(t)为窄带信号：

S_i(t)≈S_i(t+(m-1)τ_i),m=1,2,…,M   (a2)

如果相邻阵元之间的相位差为

，可以得到阵列中每个阵元的输出矢量形式：

$\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \·\\ \·\\ \·\\ x_M\left(t\right)\end{array}\right]=\left[\begin{array}{c}{e}^{j(\mathrm{\ωt}+(1-1){\mathrm{\φ}}_i)}\\ e^{j(\mathrm{\ωt}+(2-1){\mathrm{\φ}}_i)}\\ \·\\ \·\\ \·\\ e^{j(\mathrm{\ωt}+(M-1){\mathrm{\φ}}_i)}\end{array}\right]*S_i\left(t\right)---\left(a3\right)$

可化简为：

$\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \·\\ \·\\ \·\\ x_M\left(t\right)\end{array}\right]=\left[\begin{array}{c}1\\ e^{j{\mathrm{\φ}}_i}\\ \·\\ \·\\ \·\\ e^{j(M-1){\mathrm{\φ}}_i}\end{array}\right]*S_i\left(t\right)e^{\mathrm{jwt}}---\left(a4\right)$

则式（a4）可以记为：

X(t)=a(φ_i)*S_i(t)e^jwt   (a5)

其中，X(t)=[x₁(t),x₂(t),…,x_M(t)]^Ｔ，

a(φ_i)为导向向量，信源个数为N，接收到的阵列信号是：

X(t)=a(Φ₁)*S₁(t)e^jwt+a(Φ₂)*S₂(t)e^jwt+…+a(Φ_N)*S_N(t)e^jwt  (a6)

式（a6）表示为矩阵形式：

$X\left(t\right)=\left[\begin{array}{cccc}1& 1& \·\·\·& 1\\ e^{j{\mathrm{\Φ}}_1}& e^{j{\mathrm{\Φ}}_2}& \·\·\·& e^{j{\mathrm{\Φ}}_N}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ e^{j(M-1){\mathrm{\Φ}}_1}& e^{j(M-1){\mathrm{\Φ}}_2}& \·\·\·& e^{j(M-1){\mathrm{\Φ}}_N}\end{array}\right]*\left[\begin{array}{c}{S}_1\left(t\right)\\ S_2\left(t\right)\\ \·\\ \·\\ \·\\ S_N\left(t\right)\end{array}\right]e^{\mathrm{jwt}}---\left(a7\right).$

Images