GPS扩频
本文主要介绍 GPS 是如何扩频与解扩的。
GPS 信号发送
简单来讲,GPS 数据的发送与接收分为以下几个步骤:
如图所示,在发送时首先进行扩频,然后调制(通常是 BPSK)。接收时正好相反,先解调再解扩。
调制
在调制中主要使用了三种载波(L1、L2、L5)。
L1、L2、L5 的中心频率均为原子钟核心频率 \(f_0=10.23\text{ MHz}\) 的整数倍。具体的,有
- \(f_1=154f_0=1575.42\text{ MHz}\);
- \(f_2=120f_0=1227.60\text{ MHz}\);
- \(f_5=115f_0=1176.45\text{ MHz}\).
这样做的好处在于更加容易生成载波信号。
扩频
GPS 主要使用了三种扩频码:P 码、C/A 码、M 码。其中,M 码为军用,缺少相关资料。
P 码(Precise code)频率为 10.23 Mbps,运行在 L1 和 L2 上。P 码总长度较长,共有 \(2\times 10^{14}\) bit,发送需要 37 周。每个卫星(1-32)和地面站(33-37)都被分配了其中的一部分,该部分周期为 7 天,对应的是星历更新的周期。
P 码采用的是相位调制,如图所示:
加密过后的 P 码被称为 P(Y) 码。其区别是,P(Y) 码在第四个子帧中有一个标志告诉接收器扩频码已被加密,一次来防止欺骗。
C/A 码(Civilian Acquisition code)顾名思义为民用码。他的速率为 P 码的十分之一,即 1.023 Mbps。每颗卫星的 C/A 码有 1023 位,每毫秒即可广播一次。它运行在 L1 上,旨在提供标准定位服务(SPS)。作为对比,P(Y) 码则提供的是精准定位服务(PPS)。
在设计之初,美国低估了 C/A 码的定位精准度,于是通过选择可用性 (SA) 故意降低到 95% 的时间为水平 ±100 米、垂直 ±175 米精度。这一措施直到 21 世纪初才取消。
P 码和 C/A 码的生成都利用的是线性移位反馈寄存器(LSFR)。我们以 C/A 码为例。
如图是一个典型的 C/A 码生成器结构。其中,G1 和 G2 为 LSFR。输出时,通过相位选择器选择 G2 中的两个位置做异或,所得结果再与 G1 最后一位异或,即可得到扩频码。相位选择器选择的位置共有 37 种,每个卫星或地面站均不相同。
G1 和 G2 的具体原理如下所示:
他们都利用的是有限域中的运算。其中,
- G1 为 \(1+x^3+x^{10}\)
- G2 为 \(1+x^2+x^3+x^6+x^8+x^9+x^{10}\)
根据文档,G2 寄存器设置有初始值和时延。具体的值如下所示:
不过,初始值和时延不是必须的。在没有时延的情况下,如果将初始值全部置 1,最终的结果相同。
MATLAB 仿真代码如下:
function g = C_A_Code(num)
taps = [2 ,6 ; 3 ,7 ; 4 ,8 ; 5 ,9 ; 1 ,9 ;
2 ,10; 1 ,8 ; 2 ,9 ; 3 ,10; 2 ,3 ;
3 ,4 ; 5 ,6 ; 6 ,7 ; 7 ,8 ; 8 ,9 ;
9 ,10; 1 ,4 ; 2 ,5 ; 3 ,6 ; 4 ,7 ;
5 ,8 ; 6 ,9 ; 1 ,3 ; 4 ,6 ; 5 ,7 ;
6 ,8 ; 7 ,9 ; 8 ,10; 1 ,6 ; 2 ,7 ;
3 ,8 ; 4 ,9 ; 5 ,10; 4 ,10; 1 ,7 ;
2 ,8 ; 4 ,10
];
g2s = [1, 1, 0, 0, 1, 0, 0, 0, 0, 0; % 1440
1, 1, 1, 0, 0, 1, 0, 0, 0, 0; % 1620
1, 1, 1, 1, 0, 0, 1, 0, 0, 0; % 1710
1, 1, 1, 1, 1, 0, 0, 1, 0, 0; % 1744
1, 0, 0, 1, 0, 1, 1, 0, 1, 1; % 1133
1, 1, 0, 0, 1, 0, 1, 1, 0, 1; % 1455
1, 0, 0, 1, 0, 1, 1, 0, 0, 1; % 1131
1, 1, 0, 0, 1, 0, 1, 1, 0, 0; % 1454
1, 1, 1, 0, 0, 1, 0, 1, 1, 0; % 1626
1, 1, 0, 1, 0, 0, 0, 1, 0, 0; % 1504
1, 1, 1, 0, 1, 0, 0, 0, 1, 0; % 1642
1, 1, 1, 1, 1, 0, 1, 0, 0, 0; % 1750
1, 1, 1, 1, 1, 1, 0, 1, 0, 0; % 1764
1, 1, 1, 1, 1, 1, 1, 0, 1, 0; % 1772
1, 1, 1, 1, 1, 1, 1, 1, 0, 1; % 1775
1, 1, 1, 1, 1, 1, 1, 1, 1, 0; % 1776
1, 0, 0, 1, 1, 0, 1, 1, 1, 0; % 1156
1, 1, 0, 0, 1, 1, 0, 1, 1, 1; % 1467
1, 1, 1, 0, 0, 1, 1, 0, 1, 1; % 1633
1, 1, 1, 1, 0, 0, 1, 1, 0, 1; % 1715
1, 1, 1, 1, 1, 0, 0, 1, 1, 0; % 1746
1, 1, 1, 1, 1, 1, 0, 0, 1, 1; % 1763
1, 0, 0, 0, 1, 1, 0, 0, 1, 1; % 1063
1, 1, 1, 1, 0, 0, 0, 1, 1, 0; % 1706
1, 1, 1, 1, 1, 0, 0, 0, 1, 1; % 1743
1, 1, 1, 1, 1, 1, 0, 0, 0, 1; % 1761
1, 1, 1, 1, 1, 1, 1, 0, 0, 0; % 1770
1, 1, 1, 1, 1, 1, 1, 1, 0, 0; % 1774
1, 0, 0, 1, 0, 1, 0, 1, 1, 1; % 1127
1, 1, 0, 0, 1, 0, 1, 0, 1, 1; % 1453
1, 1, 1, 0, 0, 1, 0, 1, 0, 1; % 1625
1, 1, 1, 1, 0, 0, 1, 0, 1, 0; % 1712
1, 1, 1, 1, 1, 0, 0, 1, 0, 1; % 1745
1, 1, 1, 1, 0, 0, 1, 0, 1, 1; % 1713
1, 0, 0, 1, 0, 1, 1, 1, 0, 0; % 1134
1, 1, 0, 0, 1, 0, 1, 1, 1, 0; % 1456
1, 1, 1, 1, 0, 0, 1, 0, 1, 1; % 1713
];
delay = [5 , 6 , 7 , 8 , 17 , ...
18 , 139, 140, 141, 251, ...
252, 254, 255, 256, 257, ...
258, 469, 470, 471, 472, ...
473, 474, 509, 512, 513, ...
514, 515, 516, 859, 860, ...
861, 862, 863, 950, 947, ...
948, 950
];
n = 10; % size of LFSR
L = 2^n - 1; % number of chips in a code
g1_sel = [0 0 1 0 0 0 0 0 0 1]; % G1 generator
g1 = ones(1, n); % G1 vector
g2_sel = [0 1 1 0 0 1 0 1 1 1]; % G2 generator
g2 = g2s(num, :); % G2 vector
tap = taps(num, :);
for i = 1 : (delay(num) - 1)
g2 = [mod(sum(g2 .* g2_sel), 2), g2(1:n - 1)];
end
for i = 1 : L
g(i) = mod(g1(n) + mod(sum(g2(tap)), 2), 2);
g1 = [mod(sum(g1 .* g1_sel), 2), g1(1 : n - 1)];
g2 = [mod(sum(g2 .* g2_sel), 2), g2(1 : n - 1)];
end
利用程序,我们可以生成第一个和第二个 C/A 码:
当然,C/A 码具体的值不是我们所关心的。我们需要关心的是它的相关性。
左图为第 1 个 C/A code 自相关,右图为第 1 和第 2 个 C/A code 互相关。它们的明显区别在于:自相关有着非常突出的峰值,而互相关则没有。因此,我们可以通过识别峰值来找到正确的 C/A 码。
我们将第 1 个 C/A code 做相移(即,将尾部的一部分移动到头部),再与原来的 C/A 码做互相关,得到下图:
通过对比可以发现,尽管相移后仍然有明显的峰值,但峰值处的相关值降低了。我们做出了峰值与相移大小的关系图:
由图可以看出,它是一个完美的等腰三角形。在没有相移的情况下,相关值最大。
GPS 信号接收
接收机信号处理分为 4 个步骤,如下所示:
其中,跟踪分为载波环(用于解调)和码环(用于解扩)。这里,我们只讨论与主题有关的码环。
根据上文得到的结论,在解扩时,我们尝试使用不同的 C/A 码解扩,如果看到明显的峰值,则表明选对了 C/A 码。选择正确的 C/A 码后,需要使用峰值与相移大小的关系来使得相位对齐。具体做法如下图所示:
假设当前随意选取的相位为 P,其与最大值对应的相位差为 \(\delta_{cp}\)。我们将其向前移动 d 个相位到点 E,再向后移动 d 个相位到点 L。根据等腰三角形的性质可以得出,如果 E 和 L 的值相等,则 P 点在最大值,如下图所示。
而如果 E 和 L 不等,则需要进行调整,向左或向右偏移。然而,由于噪声等因素的存在,实际情况并不是完美的三角形(如下图所示),故无法直接计算出需要偏移的量,只能一点点尝试。
具体的流程如下所示:
右下角的环路滤波器接收计算出的偏移量,并生成实际的偏移量。具体的,计算方法为
\[\hat{\delta_{cp}}(n)=(1-\alpha)\hat{\delta_{cp}}(n-1)+\alpha\delta_{cp}(n)\]其中,\(\hat{\delta_{cp}}(n)\) 为实际采用的偏移量,\(\alpha\) 为滤波系数,通常取 0.05。
经过滤波器后,C/A 码晶体振荡器根据偏移量生成 C/A 码,并通过移位寄存器生成上文提到的 E 和 L。E 和 L 分别与原信号的 I 路和 Q 路相乘,并计算 1 ms 内的相关值,分别得到 \(I_E,I_L,Q_E,Q_L\)。然后根据
\[E=\sqrt{ {I_E}^2+{Q_E}^2}\\ L=\sqrt{ {I_L}^2+{Q_L}^2}\]得到幅值。
当然,前面的步骤是相干的。我们也可以使用非相干的办法,这在低信噪比的情况下表现会更好。
再计算出幅值后,进入鉴别器,计算期望的相位调整值 \(\delta_{cp}\)。具体的计算方法有三种:
-
幅值法
\[\delta_{cp}=\frac{1}{2}\frac{E-L}{E+L}\] -
功率法
\[\delta_{cp}=\frac{1}{2}\frac{E^2-L^2}{E^2+L^2}\] -
点积功率法
\[\delta_{cp}=\frac{1}{2}\frac{E-L}{P}\]
我们分别作出了这三种方法的图像:
可以看到,三条线共有三个交点。如果当前在最左侧的交点左侧,则说明 E、P、L 三点均在等腰三角形的左侧,无法继续计算;最右侧交点的右侧同理。
同时,我们可以看到幅值法是线性的——这是非常好的特性。因此,幅值法也是用得最多的方法。
Comments