高动态GPS 抗干扰算法仿真研究*

陈嘉鑫，黎柏春

（中国民航大学 航空工程学院，天津 300300）

0 引言

GPS（global positioning system）系统作为目前世界上应用最广泛的卫星导航系统，在国家经济建设和国防安全方面有着非常重要的地位。在高动态的环境下，GPS 接收机载体的速度和加速度会有较大的变化范围，再加上会有其他的干扰源会散发干扰信号，使得GPS 接收机对卫星信号的捕捉、跟踪和定位难度会进一步增大。因此，对于高动态环境下GPS性能研究仍然有非常大的研究空缺，需要不断完善。

目前，已有多种抗干扰技术，这些技术包含时域处理法、频域处理法、时频域处理法、空域处理法和空时域处理法等［1-3］。现阶段，美国在研究第3 代GPS，将抗干扰作为核心问题，对其进行了大量的研究和总结，许多研究成果都直接面向实际应用。在我国，抗干扰难题也成为国内许多高校和研究所的研究方向。例如西安电子科技大学提出的在微分约束下空域最小功率递推法［4-5］，中国民航大学的微分约束最小功率抗干扰算法，国防科大的GPS 空时抗干扰的仿真和研究成果等［6-10］。但是国内的抗干扰算法研究还不是太深入，许多成果还是基于国外提出方法的跟踪、改进和实现。

1 高动态环境下阵列接收信号模型

1.1 高动态环境的定义

早在1988 年，美国宇航局喷气推进实验室在一份关于高动态GPS 信号跟踪技术的报告中就明确给出了高动态的2 种定义。第1 种是在1 s 内有70 g的加速度斜升，也就是用70g/s 的加加速度持续1 s；第2 种比第1 种情况还要更剧烈，是0.5 s 内有50g的加速度斜升，也就是用100g/s 的加加速度持续0.5 s，而且50g/s 的加速度持续2 s，最后在0.5 s 产生50g/s 的加速度斜降。

1.2 高动态下阵列接收信号模型

设置M元等距离线阵，各阵元之间距离，其中λ为入射信号的载波波长，另设信号入射方向方位角为φ，是与阵列法线方向夹角所形成的方位角，φ∈[-π/2，π/2]。设最左侧的阵元为原点，可以推导出信号的入射导向矢量为

GPS 接收机接收到的信号中包含许多颗卫星信号，还有干扰信号和噪声，各信号及噪声之间相互独立。因为卫星离地球地面距离较远，可近似认为相对于稳定的接收机，接收到的卫星信号来向不变，而干扰相对于接收机距离较近，认为干扰信号来向会快速变化。设有W个卫星信号、Z个干扰信号，阵列接收到的信号模型表示为

式中：sw(t)表示第w个卫星信号(w= 1，2，...，W)；φw表示卫星信号的来向；jz(t)表示第z个干扰信号；φz表示干扰信号的来向；n(t)表示接收机热噪声矢量，一般是均值为0，方差为σ2n的高斯白噪声矢量。因为卫星距离地球远，信号到达接收机会损失很多能量，其电平会远低于噪声电平，所以可认为阵列接受信号的协方差矩阵主要由干扰信号与噪声信号决定，设协方差矩阵R为

式中：σ2z为第z个干扰信号的功率；σ2n为噪声功率；I为单位矩阵。因为自适应处理要求信号是稳定的，但式（3）中的协方差矩阵是不断变化的，所以只有少数可以用来估计自适应处理的协方差矩阵，式（3）转变为

式中：φz表示此刻的干扰来向；̂(t)表示此刻快拍数很少的采样数据，采样的协方差矩阵的第m行第n列元素则可以表示为

式中：δmn为Kroneckerδ函数，即当m=n时，δmn=1；当m≠n时，δmn=0。

2 Mailloux 算法原理及仿真

2.1 Mailloux 算法原理［11］

设M元等距离线阵，相邻阵元之间距离为d，接受了P个远距离的窄带信号，入射角设为φP，则阵列接受的快拍数据表示为

式中：x(t)=(x1(t)，x2(t)，...xM(t))T，x(t)为M× 1 维阵列数据向量；A为阵列流行矩阵，A=(a(φ1)，a(φ2)，...，a(φp))，矩阵中a(φp)表示第p个信号 源 的 导 向 矢 量表示为信号复包络向量，sP(t)是第p个信号源的复包络。n(t)=(n1(t)，n2(t)，...nM(t))T为M×1维阵列噪声向量。阵列的协方差矩阵定义为

式中：Rs=E[s(t)sH(t)]为信号的协方差矩阵；Ι为M维度单位矩阵；σ2n为噪声的功率。Mailloux 方法是假设在接收机接受的干扰信号源附近，存在K个相等强度的虚拟干扰，并按照一定的空间间隔排列，用来构建一个新的协方差矩阵Rˉ。用新的协方差矩阵进行自适应处理时相比于原协方差矩阵，会在接受的干扰方位附近产生较宽的零陷。新的协方差矩阵可以表示为

式中：R为没加虚拟干扰的接收信号的噪声协方差矩阵；符号⊙表示Hadamard 积；T为扩展矩阵［6］，表达式为

式 中：T(m，n) 表 示 矩 阵T的 第m行 第n列；Δmn=π(xm-xn)β/λ，xm和xn分别为等距线阵上第m，n个阵元的位置，β=B/(K- 1)，表示虚拟干扰源之间的距离，一般很小，B为虚拟的角度宽度拟干扰源。

根据式（9）可知，如果等距线阵用Mailloux 方法，不需要知道干扰信号的来向。当m=n时，扩展矩阵T的元素T(m，n) =K。新的协方差矩阵会使原协方差矩阵的噪声项功率增强K倍。

2.2 MVDR 算法仿真

仿真时，为了显示该方法形成的零陷较宽，考虑需要有对照方，所以选择MVDR（minimum variance distortionless response）算法作为常规方法，即对照方。设置阵元数为32 的等距线阵，快拍数为100，工作频率为3 × 109Hz，波长为光速除以工作频率，阵元间隔为波长的一半，干燥比为40 dB，信噪比为0 dB。因为高斯白噪声信号的功率是单位1，所以可以通过干燥比来反推得出干扰信号，通过信噪比来反推得出期望信号。设置2个互相独立的干扰信号入射至阵列接收模型，入射角度分别为-30°和30°，干扰频率 分 别 设 为5 × 109Hz 和6 × 109Hz，带 宽 为0.3 GHz，而期望信号的入射角在0°方向。将期望信号、干扰信号进行空域处理，将处理后的期望信号、干扰信号，和用randn 函数产生的噪声信号叠加，即可得阵列接收信号x(t)，通过x(t)根据式（3）得出协方差矩阵R。计算最优权值，进行空间角度扫描和计算空间导向矢量，最后画出方向图，得到MVDR 算法的方向图，如图1 所示。从图1 中可以明显看出MVDR 算法在2 处干扰方向上形成的零陷很窄。

图1 MVDR 算法的方向图仿真结果Fig.1 Simulation results of MVDR algorithm

2.3 Mailloux 算法仿真

Mailloux 方法则是需要在原协方差矩阵R上乘上扩展矩阵T，设2 个干扰信号有方向上的扰动，扰动的角度最大范围为6°，即B= 6∘，加入7 个虚拟干扰源。根据式（9）计算可得出扩展矩阵T，根据式（8）算出新协方差矩阵Rˉ。通过新的协方差矩阵来算最优权值，进行空间角度扫描和计算空间导向矢量，最后画出方向图，得到Mailloux 算法与MVDR 算法的方向图对比，如图2 所示。从图2 对比来看，Mailloux 方法能够明显地加宽干扰方向上形成的零陷。

图2 Mailloux 算法的方向图仿真结果Fig.2 Direction diagram simulation results of Mailloux algorithm

2.4 高动态下算法对比仿真

2.4.1 加入1 颗卫星信号仿真

以上是没有加入GPS 信号的仿真，本文研究GPS 抗干扰，所以以下是加入了GPS 信号的仿真。GPS 卫星信号可以用2 种方法产生，第1 种是直接调用已生成的GPS 数据文件，第2 种是用GPS 数据模拟器产生所需要的GPS 数据文件。本文采用的是第1 种方法，调用GPS 卫星信号产生程序，该程序调用已生成的GPS 卫星信号数据文件，先调用了1 颗卫星信号，卫星信号编号为1，读取0.1 s 的卫星信号数据，设置采样率为5 714 000 Hz，所以卫星信号数据为1 × 5 714 00 维矩阵。为保证程序能够运行，所以将快拍数也改为5 714 00。将产生的GPS 卫星信号矩阵替代原仿真程序中的期望信号。阵列接收信号模型重新设置为7 元阵列等距线阵，保留一个干扰信号，干扰来向设为-30°，干扰信号和噪声信号矩阵改变维度，其他如工作频率和波长都不变。

因为卫星轨道高度在20 000 km 以上，高动态GPS 接收机距离卫星很远，所以在短时间内，近似认为接收机接收到的GPS 卫星信号的来向不变，也就是φw不变，仿真中设为30∘。为了体现高动态，认为干扰信号距离接收机较近，所以认为干扰信号的波达方向会逐渐变化，但变化不大。仿真中为了直观方便，设置干扰信号的前N/2 个数据的波达方向为-30∘，后N/2 个数据的波达方向为-32∘。利用干扰信号的前N/2 个数据求解阵列加权向量，得到高动态下Mailloux 算法的方向图，如图3 所示。

图3 1 颗卫星的Mailloux 算法方向图仿真结果Fig.3 Direction diagram of one satellite based on Mailloux algorithm

从图3 中可看出，零陷明显加宽。

新建函数：

式中：WH为利用干扰信号的前N/2 个数据求解阵列加权向量；X为后N/2 个样本数据。

将新建函数作为输入参数带入捕获程序进行捕获。先计算每个C/A 码包含多少采样点，取相邻的2 个信号，信号长度为1 ms，去直流，计算采样时间间隔，计算1 ms 信号内的载波相位点，计算频率捕获的个数，产生一颗卫星的载波偏移初始向量，初始化32 颗卫星的载波频率，对已产生的当前通道的C/A 码做离散傅里叶变换，扫描频点，找出最大的值及其对应的频点，确认最大值对应的相位，获取与最大峰超过1 个码片的次最大峰，求与最大峰距离为1 个码片的区域边界，确定寻找次最大峰的搜索区域，保存各通道峰峰比值，判决，如果比值比预设阈值大，则表示能接收到该卫星信号。图4 为MVDR 算法1 号 卫星信号捕获图，图5 为Mailloux 算法1 号卫星信号捕获图。

图4 MVDR 算法1 号卫星信号捕获图Fig.4 Signal acquisition diagram of satellite No.1 base on MVDR algorithm

图5 Mailloux 算法1 号卫星信号捕获图Fig.5 Signal acquisition diagram of satellite No.1 based on Mailloux algorithm

2.4.2 加入4 颗卫星信号仿真

为验证算法可行度，加入多个卫星信号。调用4 颗卫星信号，分别是1，2，3，20 号卫星。图6 是加4颗卫星的Mailloux 算法的方向图，图7 是MVDR 算法4 颗卫星信号捕获图，图8 是Mailloux 算法4 颗卫星信号捕获图。当加入多颗卫星后，如果捕获图8中没有出现已加入的卫星序号，只需要将捕获门限值提高一点即可。从图7，8 的对比来看，Mailloux 算法在高动态条件下能够捕获卫星信号。

图6 4 颗卫星的Mailloux 算法的方向图Fig.6 Direction diagram of four satellites based on Mailloux algorithm

图7 MVDR 算法4 颗卫星信号捕获图Fig.7 Signal acquisition diagram of four satellites based on MVDR algorithm

图8 Mailloux 算法4 颗卫星信号捕获图Fig.8 Signal acquisition diagram of four satellites based on Mailloux algorithm

3 Zatman 算法原理及仿真

3.1 Zatman 算 法原理［12-14］

Zatman 方法是假设接收机接受的窄带干扰信号具有虚拟带宽Bw，而且在整条带宽内有平坦的功率谱，用来构建一个新的协方差矩阵Rˉ，用新的协方差矩阵进行自适应处理时相比于原协方差矩阵，会在接受的干扰方位附近产生较宽的零陷。新的协方差矩阵可以表示为

式中：T(m，n)表示矩阵T的第m行第n列；xm和xn分别为等距线阵上第m和n个阵元的位置；φz表示第z个干扰信号的来向；c为光速。

由式（11）可知，Zatman方法需要知道干扰信号的来向。当m=n时，扩展矩阵T的元素T(m，n) = 1，不会改变原协方差矩阵R中的噪声项功率。

3.2 Zatman 算法仿真

对照方依旧选择MVDR 算法。先不加卫星信号，使用期望信号。阵元数为32 的等距线阵，快拍数为100，使用2 个独立的干扰信号。工作频率为3 ×109Hz，波长等于光速除以工作频率，阵元间隔为波长的一半，干燥比为40 dB，信噪比为0 dB，设高斯白噪声的功率为1。通过干燥比反推得到干扰信号，通过信噪比反推得到期望信号。将各信号叠加得到x(t)，剩余的操作与第3 节中MVDR 算法仿真中的操作一致。MVDR 算法得到的方向图如图1 所示。

Zatman 算法是在MVDR 算法的基础矩阵上再乘上一个扩展矩阵T。假设干扰信号的频带扩展宽度为0.3 GHz，即Bw= 0.3 GHz。根据式（11）得到扩展矩阵，根据式（8）算出新协方差矩阵Rˉ。通过新的协方差矩阵来算最优权值，进行空间角度扫描和计算空间导向矢量，最后画出方向图，得到Zatman 算法与MVDR 算法的方向图对比，如图9 所示。

图9 Zatman 算法的方向图Fig.9 Direction diagram of Zatman algorithm

3.3 高动态下算法对比仿真

3.3.1 加入1 颗卫星信号仿真

以上是没有加卫星信号的仿真，以下是加入卫星信号的仿真。还是用第1 种方法直接调用已生成的GPS 数据文件。先调用1 号卫星信号，读取参数与第3 节里一致。

为保证程序能够运行，将快拍数也改为571 400。将产生的1 号卫星信号矩阵替代原仿真程序中的期望信号。设置7 元阵列等距线阵，保留一个干扰信号，干扰来向设为-30°。为了体现高动态，也将设置干扰信号的前N/2 个数据的波达方向为-30∘，后N/2个数据的波达方向为-32∘。用前N/2 个数据求解阵列加权向量，得到的加权向量处理后N/2 个数据，得到加宽的方向图，如图10 所示。从图10 中可看出，零陷加宽。

图10 1 号卫星的Zatman 算法的方向图Fig.10 Direction diagram of one satellite based on Zatman algorithm

根据式（10）得到新建函数，将新建函数送入捕获程序中捕获卫星信号。捕获程序流程与2.3中所写一致。MVDR 算法1 号卫星信号捕获图如图4 所示，图11 为Zatman 算法1 号卫星信号捕获图。从图11 中可知，Zatman 算法在高动态条件下可以捕获到1 号卫星。

图11 Zatman 算法1 号卫星信号捕获图Fig.11 Signal acquisition diagram of satellite No.1 based on Zatman algorithm

3.3.2 加入4 颗卫星信号仿真

为验证算法有效性，加入多个卫星信号进行仿真。调用4 颗卫星信号，分别是1，2，3，20 号卫星。图12 为加4 颗卫星的Zatman 算法的方向图，MVDR算法4 颗卫星信号捕获图如图7 所示，图13 为Zatman算法4 颗卫星信号捕获图。当加入多颗卫星后，捕获图也出现了Mailloux 算法类似的问题，同样只需要将捕获门限值提高一点即可。从图7，13 的对比来看，Zatman 算法在高动态下能够捕获卫星信号。

图12 4 颗卫星的Zatman 算法的方向图Fig.12 Direction diagram of four satellites based on Zatman algorithm

图13 Zatman 算法4 颗卫星信号捕获图Fig.13 Signal acquisition diagram of four satellites based on Zatman algorithm

4 Mailloux 算 法 和Zatman 算 法 的 优缺点对比分析

根据2 种算法的推导和仿真结果的对比可得，Mailloux 算法运算量适中，在线阵模型中不需要估计干扰来向，可以直接得到扩展矩阵，但在面阵模型中，该算法需提前给出预估得干扰信号来向。此外还有一个缺点，用T计算Rˉ时，协方差矩阵的噪声项功率会增强K倍。该算法在方向图干扰来向上形成的零陷比Zatman 算法较宽，比MVDR 算法宽很多。

Zatman 算法运算量不大，对于干扰噪声的协方差矩阵的功率没有增强，但无法保证角度的取值，不能进行角度的扩展，因此该算法不可应用于面阵，只可用于线阵。该算法在方向图干扰来向上形成的零陷比MVDR 算法宽，但相比Mailloux 算法要窄一些。

通过2 种算法得对比可知，Zatman 算法不能用在圆阵以及任意阵模型下，利用Mailloux 算法需要预估干扰的来向信息。但是，如果干扰方向随时间快速变化，难以得到干扰方向的来向信息，Mailloux算法便不能解决问题，与此同时也会增加很大的计算量。

5 结束语

本文首先研究了当阵列为等距线阵时接收信号模型。根据信号入射的导向矢量和信号中不同组成部分，得到阵列接收信号模型。并由信号模型求得协方差矩阵，来表示各个维度信号之间的相关性。在Matlab 上搭建了阵列接收信号模型。之后对Mailloux 算法和Zatman 算法原理进行了研究分析。仿真时在阵列接收信号模型的基础上，通过虚加干扰源和扩展频带来获得扩展矩阵，从而得到新的协方差矩阵。为更直观验证算法的有效性，选择MVDR 算法作为对照。最后为体现高动态的特点，仿真时改变干扰信号数据后半部分的来向。用前半部分数据计算新的协方差矩阵，得到的矩阵处理后半部分数据。

通过仿真验证可知，在高动态条件下，当阵元是等距线阵时，Mailloux 算法不需要知道干扰的来向信息，但Zatman 算法需要知道干扰的来向。相较于MVDR 算法，Mailloux 算法和Zatman 算法都能在方向图干扰来向上形成较宽的零陷，能够有效抑制掉干扰信号。根据需要，将3 种算法形成的新建函数送入捕获程序后，MVDR 算法不能捕获到卫星信号，而Mailloux 算法和Zatman 算法能捕获到卫星信号。

本文只验证了阵列接收信号模型是在等距线阵时算法的有效性。由于在飞机、导弹等一系列高速运动的载体上，圆阵列的应用比较广泛，后续将验证圆阵列和其他任意形模型下，Mailloux 算法和Zatman 算法在高动态条件下的抗干扰有效性。