基于解重扩算法的GNSS矢量接收机抗干扰方法与流程

本发明属于卫星导航抗干扰技术领域，特别是涉及一种基于解重扩算法的gnss(全球导航卫星系统)矢量接收机抗干扰方法。

背景技术：

基于矢量接收机的抗干扰研究目前主要分为信号处理与矢量跟踪环路简单级联和深度结合两大类，2012年plan小组提出了根据相关器输出幅度的变化来检测欺骗式干扰，然后利用矢量跟踪环路能够桥接被遮挡真实信号的特点，保证接收机不能跟踪上干扰信号。2012年dlr将特征波束形成技术与矢量跟踪环路紧耦合来提高矢量接收机的稳定性；2016年又利用多波束矢量接收机进行欺骗式干扰的检测，但是该方法复杂度较大；2015年美国的伊利诺伊大学的研究人员利用多个矢量接收机的辅助位置信息进行欺骗式干扰的检测，但是该方法需要多个接收机，并且需要保证接收机之间的通信协议。2014年哈尔滨工程大学的李强等人提出了一种利用矢量跟踪环路预测的信号作为cross-score算法中的参考信号形成权矢量来抑制干扰，2016年徐明山等人将盲波束形成技术与矢量跟踪环路结合来提高矢量跟踪环路的鲁棒性，但以上两种方法在复杂的干扰环境下矢量跟踪环路的性能会下降。为了保障gnss的安全应用，迫切需要研究能够提高接收机在复杂电磁环境下性能并且对阵列流形误差和干扰功率不敏感的抗干扰技术。这不仅在理论研究方面有重要的意义，而对推进gnss在军事应用等领域的应用具有重要的意义。

技术实现要素：

本发明所要解决的技术问题是，提供一种不需要已知阵列流形，对卫星信号提高增益，提高矢量接收机的定位性能，抗压制性干扰性能以及在复杂电磁环境下的稳健性的基于解重扩算法的gnss矢量接收机抗干扰方法。

为了达到上述目的，本发明提供的基于解重扩算法的gnss矢量接收机抗干扰方法包括按顺序进行的下列步骤：

1)将阵列天线接收到的数字中频信号向干扰正交补空间投影来抑制压制性干扰，并将抑制干扰后的数据送入接收机，由标量跟踪环路完成定位解算后给出接收机状态矢量以及跟踪环路参数初值；

2)对步骤1)得到的接收机状态矢量和跟踪环路参数进行预测，利用预测得到的跟踪环路参数构造参考信号，将接收信号和参考信号进行相干累加求解得到权矢量；

4)利用步骤2)获得的权矢量对接收信号加权，得到抑制干扰后的相关结果；

4)将上述步骤3)获得的抑制干扰后的相关结果送入鉴别器中而对接收机状态矢量和跟踪环路参数进行更新校正，以对接收机状态矢量和跟踪环路参数进行误差校正和预测值的补偿。

在步骤1)中，所述的将阵列天线接收到的数字中频信号向干扰正交补空间投影来抑制压制性干扰，并将抑制干扰后的数据送入接收机，由标量跟踪环路完成定位解算后给出接收机状态矢量以及跟踪环路参数初值的方法是：首先利用接收到的数字中频信号xs,j(t)构造协方差矩阵然后利用协方差矩阵构造正交补空间得到抑制压制性干扰的数据将上述数据x(t)送入普通标量接收机进行捕获得到接收信号的粗略载波频率、粗略码相位并对其进行逐步跟踪得到较为精确的跟踪环路参数初值，即载波频率、载波相位、码频率、码相位；将只包含数据码的跟踪处理后的信号进行解码，利用卫星信号到达接收机的时间来计算伪距ρk＝cτ，其中c为光速，τ为时间，最后基于线性迭代技术计算得到用户的位置之后得到接收机的状态矢量初值xk，即三维位置xkykzk，三维速度vx,kvy,kvz,k以及包括钟差tb,k、钟漂td,k在内的时间。

在步骤2)中，所述的对步骤1)得到的接收机状态矢量和跟踪环路参数进行预测，利用预测得到的跟踪环路参数构造参考信号，将接收信号和参考信号进行相干累加求解得到权矢量的方法是：选取时间间隔为1ms,一般认为接收机在这段极短的时间内是匀速运动的，此时可以根据k时刻接收机的状态矢量xk对k+1时刻的状态矢量进行预测得到即三维位置、三维速度、钟差、钟漂；然后根据预测到的状态矢量预测得到k+1时刻接收机和卫星之间的伪距ρk+1，将伪距对时间求导得到预测的伪距率利用伪距和伪距率得到k+1时刻跟踪环路参数的预测值，即码频率载波频率码相位利用跟踪环路的预测值得到参考信号以采样率fs取nms的数据样本将抑制压制性干扰后的接收信号与参考信号相干累加得到即时码的相关输出对即时码的相关输出进行协方差矩阵求解并取最大特征值对应的特征向量得到最优权矢量w。

在步骤3)中，所述的利用步骤2)获得的权矢量对接收信号加权，得到抑制干扰后的相关结果的方法是：将权矢量w对抑制压制性干扰后的接收信号进行加权得到加权信号y(t)＝w^hx(t)，对加权信号进行逐步跟踪处理得到抑制干扰后的相关结果，即早码相关值e(n)、即时码相关值p(n)、晚码相关值l(n)。

在步骤4)中，所述的将上述步骤3)获得的抑制干扰后的相关结果送入鉴别器中而对接收机状态矢量和跟踪环路参数进行更新校正，以对接收机状态矢量和跟踪环路参数进行误差校正和预测值的补偿的方法是：将上一步骤得到的早码相关值e(n)、即时码相关值p(n)、晚码相关值l(n)送入鉴别器中得到码相位误差△τk+1、载波频率误差△fk+1，将码相位误差和载波频率误差分别乘以相应的波长得到伪距误差△ρk+1、伪距率误差此时根据误差状态量建立扩展kalman滤波器状态方程和观测方程；得到矢量跟踪环路的状态方程和测量方程之后，根据k+1时刻的状态矢量误差将k+1时刻的状态矢量更新为xk+1＝xk+1+δxk+1，根据接收机更新校正后的三维位置、三维速度、时间对k+1时刻的码相位、码频率和载波频率进行更新校正；经过上面的校正和补偿后得到k+1时刻接收机状态矢量和跟踪环路参数。

本发明提供的基于解重扩算法的gnss矢量接收机抗干扰方法能产生高增益的同时可以提高接收机的定位性能，并且具有抗压制性干扰能力，无需知道卫星信号和干扰来向信息，对阵列流形误差不敏感，稳健性较好。

附图说明

图1是接收机工作过程中受到来向较远的大功率压制性干扰的情况下，采用本发明方法得到的相应卫星信号各个时间的波束图。

图2是接收机工作过程中受到来向较远的大功率压制性干扰的情况下，采用本发明方法得到的各颗卫星信号各个时间的载噪比。

图3是接收机工作过程中受到来向较远的大功率压制性干扰的情况下，将本发明方法与普通矢量跟踪环路的载噪比进行对比的结果。

图4是接收机工作过程中受到来向较远的大功率压制性干扰的情况下，将本发明方法与普通矢量跟踪环路的鉴相器输出进行对比的结果。

图5是接收机工作过程中受到来向较远的大功率压制性干扰的情况下，将本发明方法与普通矢量跟踪环路的定位误差进行对比的结果。

图6a和图6b分别是两个卫星在接收机工作过程中受到来向较近的大功率压制性干扰的情况下，本发明方法得到的波束图。

图7是在接收机工作过程中受到来向较近的大功率压制性干扰的情况下，本发明方法与普通矢量跟踪环路的定位误差进行对比的结果。

具体实施方式

下面参照附图和具体实施例对本发明提供的基于解重扩算法的gnss矢量接收机抗干扰方法给予详细说明。为了方便起见，下面以gps(globalpositioningsystem)系统和m个阵元的阵列天线为例加以说明，实际上，该方法适用于任意阵列天线，且不需要知道阵列流形，并且对阵列流形误差稳健。

本发明提供的基于解重扩算法的gnss矢量接收机抗干扰方法包括按顺序进行的下列步骤：

1)将阵列天线接收到的数字中频信号向干扰正交补空间投影来抑制压制性干扰，并将抑制干扰后的数据送入接收机，由标量跟踪环路完成定位解算后给出接收机状态矢量以及跟踪环路参数初值；

考虑一个由m个阵元组成的均匀线阵，阵元间距为d，入射信号的波长为λ，取第一个阵元为参考点，如果有g个有用卫星信号入射到均匀线阵上，到达角度分别为θ^g(g＝1,2...,g)，在远场窄带假设条件下，阵列天线接收的数字中频信号可表示为：

其中xs,j(t)＝[x1(t),x2(t),...,xm(t)]是m×1维阵列天线数据矢量；e(t)代表m×1维阵列天线接收噪声矢量，通常假设服从高斯分布；是第g颗卫星信号的复包络；为压制性干扰，a(θ^q)为第q个压制性干扰的导向矢量；a(θ^g)是第g颗卫星信号的导向矢量，且有：

gnss中频信号的数学模型可表示为：

其中为第g颗卫星信号的幅度，为数据位，为扩频码，τ^g为码相位，f0为载波频率，为载波初始相位。将公式(3)代入公式(1)中得到：

对上述数字中频信号xs,j(t)向干扰正交补空间投影来消除压制性干扰，即：

将上述抑制干扰后的数据送入接收机，由标量跟踪环路完成定位解后得到接收机状态矢量和跟踪环路参数初值，接收机状态矢量可表示为x＝[xyzctbvxvyvzctd]^t，其中p＝[xyz]^t表示接收机在ecef坐标系中的三维位置坐标，单位是m，v＝[vxvyvz]^t表示接收机在ecef坐标系中的三维速度，单位为m/s，ctb表示接收机的钟差，前面乘以c(光速)使得其单位为m,ctd表示接收机的钟漂，前面乘以c使得其单位为m/s。跟踪环路参数包括码相位、码频率、载波频率和载波相位。

2)对步骤1)得到的接收机状态矢量和跟踪环路参数进行预测，利用预测得到的跟踪环路参数构造参考信号，将接收信号和参考信号进行相干累加求解得到权矢量；

首先对接收机状态矢量进行预测。当选取的两个时刻之间足够短，比如1ms时，一般情况下可以认为接收机在这段极短的时间内匀速运动。根据k时刻接收机状态矢量xk可以预测k+1时刻的接收机状态矢量其中k+1时刻的位置分量和速度分量的预测值为：

同理k+1时刻的钟差和钟漂的预测值为：

其中tk,k+1表示k时刻到k+1时刻的时间间隔，这样就得到了k+1时刻的接收机状态矢量的预测值。

然后对跟踪环路参数进行预测。经过卫星时钟和误差校正后第g颗卫星的伪距ρ^g和伪距率的计算公式为：

其中表示第g颗卫星与接收机的几何距离(x^g是第g颗卫星x方向的坐标，其他类似)，表示第g颗卫星与接收机的相对速度在视线方向的变化率，其中l^g表示第g颗卫星与接收机连线方向的单位向量，计算公式如下：

k+1时刻接收机状态矢量的预测值可由公式(6)、(7)算出，卫星在k+1时刻的位置可由星历信息得到,那么k+1时刻伪距和伪距率就可由公式(8)预测得到，因此k+1时刻码频率和载波频率的预测值为：

其中fcode表示c/a码的频率(1.023mhz)，fif表示下变频后的中频频率。可以根据k+1时刻的卫星位置和接收机的预测位置由公式(9)得到。同样码相位是与距离相互联系起来的，k+1时刻的码相位可以写成：

以第g颗卫星信号为例，利用公式(10-11)可得到码频率载波频率和码相位(由于载波相位是同步的，因此本发明中不需预测)，再结合矢量跟踪环路工作前的标量跟踪环路给出的相应卫星信号的导航电文重构一个参考信号：

取长度为nms的数据样本利用公式(12)得到的参考信号通过以下方式得到权矢量：

首先将接收信号和参考信号进行相干累加得到即时码的相关输出：

其次对公式(13)进行协方差矩阵求解：

由于每个跟踪环路只处理一颗卫星信号，因此对协方差矩阵r进行特征值分解并取最大特征值对应的导向矢量u求得权矢量，计算公式如下：

w^g＝u/|u|(15)

其中为第g颗卫星信号c/a码自相关函数，为第g颗卫星信号与其余几颗卫星信号c/a码的互相关函数，为码相位误差和载波频率误差，tcoh为相干累积时间，本发明中为1ms，为载波相位误差；为第g颗卫星信号与相应即时码的相关结果，由于不同卫星信号的c/a码的互相关系数不为零，对于第g颗卫星来说其余g-1颗卫星信号在相关的过程中相当于干扰信号为跟踪环路中的测量噪声矢量，p^g(n)为第g颗卫星跟踪环路的相关输出。因此权矢量w^g正比于导向矢量且导向矢量和真实卫星信号导向矢量相关性接近1。

3)利用步骤2)获得的权矢量对接收信号加权，得到抑制干扰后的相关结果；

利用步骤(2)中预测到的权矢量对接收信号进行加权得到：

对加权后的信号进行跟踪，结果如下：

其中e(n)为早码相关值、p(n)为即时码相关值、l(n)为晚码相关值。

对于第g颗卫星信号来说得到的权矢量w^g和其余几颗卫星信号导向矢量a(θq)相关性较弱，利用该权矢量可以减小干扰信号在相关过程中的贡献，该权矢量对干扰在相关输出的贡献作用为：

4)将上述步骤3)获得的抑制干扰后的相关结果送入鉴别器中而对接收机状态矢量和跟踪环路参数进行更新校正，以对接收机状态矢量和跟踪环路参数进行误差校正和预测值的补偿；

将步骤(3)中获得的抑制干扰后的相关结果送入鉴别器，输出码相位误差和载波频率误差即：

其中表示每个跟踪环路中鉴别器的测量噪声，将码相位误差和载波频率误差分别乘以相应的波长可转化成伪距误差和伪距率误差将伪距和伪距率计算公式(8)进行一阶泰勒展开就可以得到公式：

这样就将伪距误差和伪距率误差和接收机的三维位置误差和速度误差联系起来，并将上式作为扩展kalman滤波器的测量方程。选取误差状态量δx作为扩展kalman滤波器的状态量，选取测量值z(鉴相器和鉴频器的输出转化成的伪距和伪距率误差)作为扩展kalman滤波器的观测量，此时扩展kalman滤波器的状态方程和测量方程为：

其中fk,k+1表示一步状态转移矩阵，如式(23)所示，wk表示扩展kalman滤波器的处理噪声矩阵，vk表示扩展kalman滤波器的测量方程的噪项，hk+1是由公式(21)得出的观测量与状态量之间的投影矩阵。

得到矢量跟踪环路的状态方程和测量方程之后，根据k+1时刻的状态矢量误差可以将k+1时刻的状态矢量更新为xk+1＝xk+1+δxk+1，k+1时刻的码相位、码频率和载波频率可以更新为：

式中的是由校正后的接收机位置得到的，不再是前面的预测值。经过上面的校正和补偿后得到了k+1时刻接收机状态矢量和跟踪环路参数，并将扩展kalman滤波器的状态量δx置0作为下一次更新的初值，从而形成一个闭合的回路，完成一次完整的基于波束形成的gnss矢量接收机的矢量跟踪过程。

本发明提供的基于解重扩算法的gnss矢量接收机抗干扰方法可以通过以下仿真结果进一步说明。

图1是采用本发明方法得到的权矢量波束图。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入干扰来向为-75度的压制性干扰，干噪比为30db。图1中各个时间得到的波束图指向真实卫星信号来向。从图1可以看出，本发明方法得到的权矢量较稳定。

图2是采用本发明方法对所有卫星信号跟踪处理后得到的载噪比。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入干扰来向为-75度的压制性干扰，干噪比为30db。图2中的深红色、蓝色、绿色、浅红色、黑色的线条分别代表prn1，prn3，prn14，prn20，prn22五个卫星信号的载噪比随时间的变化。从图2可以看出，本发明方法可以对每颗卫星信号形成较高的增益，而且载噪比水平较平稳。

图3是将本发明方法的载噪比输出与普通矢量跟踪环路方法的载噪比进行对比的结果。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入干扰来向为-75度的压制性干扰，干噪比为30db。从图3可以看出，本发明方法可以使卫星信号载噪比提高9db左右。

图4是将本发明方法的鉴相器输出与普通矢量跟踪环路方法的鉴相器输出进行对比的结果。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入干扰来向为-75度的压制性干扰，干噪比为30db。从图4可以看出，本发明方法提高了跟踪环路性能，减小了鉴相器输出误差。

图5是将本发明的定位误差与普通矢量跟踪环路方法的定位误差进行对比的结果。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入干扰来向为-75度的压制性干扰，干噪比为30db。从图5可以看出，本发明方法提高了跟踪环路性能，减小了定位误差，表1为定位误差统计结果，可见本发明方法的定位误差比普通矢量跟踪环路方法整体减小0.5m左右。

图6是在矢量跟踪环路工作的过程中受到与真实卫星信号来向较近的大功率压制性干扰的情况下，本发明方法25s后在prn3(图6a)，prn14(图6b)跟踪环路形成的权矢量。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入两个来向分别为9°,33°的压制性干扰，干噪比为30db。通过图6发现25s之后prn3(图6a)，prn14(图6b)跟踪环路利用本发明方法提供的高增益可以在被削弱卫星信号的来向上形成一定的增益，以保证卫星信号被持续跟踪。

图7是在矢量跟踪环路工作的过程中受到与真实卫星信号来向较近的大功率压制性干扰的情况下，本发明方法的定位误差与普通矢量跟踪环路方法的定位误差进行对比的结果。在本实施例中，阵列天线为七阵元均匀线阵，五个gps卫星信号prn1，prn3，prn14，prn20，prn22分别从50°,35°,10°,-10°,-35°方向入射到阵列天线上，信噪比为-20db，信号长度为70s，25s开始加入两个来向分别为9°,33°的压制性干扰，干噪比为30db。从图7可以看出，本发明方法具有较好的抗压制性干扰的性能，普通矢量跟踪环路方法受到压制性干扰之后会失锁，进而导致定位误差增大，然而本发明方法在复杂电磁环境下仍然具有较好的稳健性。

表1定位误差