本文主要介绍 GPS 是如何扩频与解扩的。
GPS 信号发送
简单来讲,GPS 数据的发送与接收分为以下几个步骤:
![image-20220524171208992](https://i-blog.csdnimg.cn/blog_migrate/7bb8a15eb12a81dbbd642d6441e5fa12.png)
如图所示,在发送时首先进行扩频,然后调制(通常是 BPSK)。接收时正好相反,先解调再解扩。
调制
在调制中主要使用了三种载波(L1、L2、L5)。
L1、L2、L5 的中心频率均为原子钟核心频率 f 0 = 10.23 MHz f_0=10.23\text{ MHz} f0=10.23 MHz 的整数倍。具体的,有
- f 1 = 154 f 0 = 1575.42 MHz f_1=154f_0=1575.42\text{ MHz} f1=154f0=1575.42 MHz
- f 2 = 120 f 0 = 1227.60 MHz f_2=120f_0=1227.60\text{ MHz} f2=120f0=1227.60 MHz
- f 5 = 115 f 0 = 1176.45 MHz f_5=115f_0=1176.45\text{ MHz} f5=115f0=1176.45 MHz
这样做的好处在于更加容易生成载波信号。
扩频
GPS 主要使用了三种扩频码:P 码、C/A 码、M 码。其中,M 码为军用,缺少相关资料。
P 码(Precise code)频率为 10.23 Mbps,运行在 L1 和 L2 上。P 码总长度较长,共有 2 × 1 0 14 2\times 10^{14} 2×1014 bit,发送需要 37 周。每个卫星(1-32)和地面站(33-37)都被分配了其中的一部分,该部分周期为 7 天,对应的是星历更新的周期。
P 码采用的是相位调制,如图所示:
![image-20220524172709793](https://i-blog.csdnimg.cn/blog_migrate/9fd7e7c0aef6da26b65a132eac7708d4.png)
加密过后的 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 码为例。
![image-20220524173604302](https://i-blog.csdnimg.cn/blog_migrate/9247715e4d497893a7f104452598c102.png)
如图是一个典型的 C/A 码生成器结构。其中,G1 和 G2 为 LSFR。输出时,通过相位选择器选择 G2 中的两个位置做异或,所得结果再与 G1 最后一位异或,即可得到扩频码。相位选择器选择的位置共有 37 种,每个卫星或地面站均不相同。
G1 和 G2 的具体原理如下所示:
![image-20220524173918234](https://i-blog.csdnimg.cn/blog_migrate/219b5bee7c415043b2d442b706ca2148.png)
![image-20220524173924815](https://i-blog.csdnimg.cn/blog_migrate/c80b1c04aeced5e906905cb794066210.png)
他们都利用的是有限域中的运算。其中,
- G1 为 1 + x 3 + x 10 1+x^3+x^{10} 1+x3+x10
- G2 为 1 + x 2 + x 3 + x 6 + x 8 + x 9 + x 10 1+x^2+x^3+x^6+x^8+x^9+x^{10} 1+x2+x3+x6+x8+x9+x10
根据文档,G2 寄存器设置有初始值和时延。具体的值如下所示:
![image-20220524174506382](https://i-blog.csdnimg.cn/blog_migrate/bcb29a652df8197327dab4ae11bc32e5.png)
不过,初始值和时延不是必须的。在没有时延的情况下,如果将初始值全部置 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 码:
![](https://i-blog.csdnimg.cn/blog_migrate/8cfdb86d291ddf4255e11805dfbd7195.png)
![](https://i-blog.csdnimg.cn/blog_migrate/aa6d9935d13ae8f39cc02fc13467aa98.png)
当然,C/A 码具体的值不是我们所关心的。我们需要关心的是它的相关性。
![](https://i-blog.csdnimg.cn/blog_migrate/caa06cefb2832fdda8a5b22d4944a2cf.png)
![](https://i-blog.csdnimg.cn/blog_migrate/277a4480fbc2d6613d8dda896e1f33d8.png)
左图为第 1 个 C/A code 自相关,右图为第 1 和第 2 个 C/A code 互相关。它们的明显区别在于:自相关有着非常突出的峰值,而互相关则没有。因此,我们可以通过识别峰值来找到正确的 C/A 码。
我们将第 1 个 C/A code 做相移(即,将尾部的一部分移动到头部),再与原来的 C/A 码做互相关,得到下图:
![](https://i-blog.csdnimg.cn/blog_migrate/caa06cefb2832fdda8a5b22d4944a2cf.png)
![](https://i-blog.csdnimg.cn/blog_migrate/e92c53804f1cd6a3dfd51a3e4540fae3.png)
通过对比可以发现,尽管相移后仍然有明显的峰值,但峰值处的相关值降低了。我们做出了峰值与相移大小的关系图:
![](https://i-blog.csdnimg.cn/blog_migrate/e19c4478e1f3dbb183d37ef1e61e48aa.png)
由图可以看出,它是一个完美的等腰三角形。在没有相移的情况下,相关值最大。
GPS 信号接收
接收机信号处理分为 4 个步骤,如下所示:
![image-20220524182221824](https://i-blog.csdnimg.cn/blog_migrate/e13588373d41db927fa52a74a68d5dda.png)
其中,跟踪分为载波环(用于解调)和码环(用于解扩)。这里,我们只讨论与主题有关的码环。
根据上文得到的结论,在解扩时,我们尝试使用不同的 C/A 码解扩,如果看到明显的峰值,则表明选对了 C/A 码。选择正确的 C/A 码后,需要使用峰值与相移大小的关系来使得相位对齐。具体做法如下图所示:
![image-20220524182913117](https://i-blog.csdnimg.cn/blog_migrate/c282d376319f5998771b5ff663a64ce2.png)
假设当前随意选取的相位为 P,其与最大值对应的相位差为 δ c p \delta_{cp} δcp。我们将其向前移动 d 个相位到点 E,再向后移动 d 个相位到点 L。根据等腰三角形的性质可以得出,如果 E 和 L 的值相等,则 P 点在最大值,如下图所示。
![image-20220524183313551](https://i-blog.csdnimg.cn/blog_migrate/8225db6e60681c2258f49715902f9124.png)
而如果 E 和 L 不等,则需要进行调整,向左或向右偏移。然而,由于噪声等因素的存在,实际情况并不是完美的三角形(如下图所示),故无法直接计算出需要偏移的量,只能一点点尝试。
![image-20220524183734467](https://i-blog.csdnimg.cn/blog_migrate/1096fdbe5587342887f7dc017ee0ef4a.png)
具体的流程如下所示:
![image-20220524183805587](https://i-blog.csdnimg.cn/blog_migrate/a6bdd3c7d8cf05e24b1e7890b8ce6f5d.png)
右下角的环路滤波器接收计算出的偏移量,并生成实际的偏移量。具体的,计算方法为
δ
c
p
^
(
n
)
=
(
1
−
α
)
δ
c
p
^
(
n
−
1
)
+
α
δ
c
p
(
n
)
\hat{\delta_{cp}}(n)=(1-\alpha)\hat{\delta_{cp}}(n-1)+\alpha\delta_{cp}(n)
δcp^(n)=(1−α)δcp^(n−1)+αδcp(n)
其中,
δ
c
p
^
(
n
)
\hat{\delta_{cp}}(n)
δ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
I_E,I_L,Q_E,Q_L
IE,IL,QE,QL。然后根据
E
=
I
E
2
+
Q
E
2
L
=
I
L
2
+
Q
L
2
E=\sqrt{ {I_E}^2+{Q_E}^2}\\ L=\sqrt{ {I_L}^2+{Q_L}^2}
E=IE2+QE2L=IL2+QL2
得到幅值。
当然,前面的步骤是相干的。我们也可以使用非相干的办法,这在低信噪比的情况下表现会更好。
再计算出幅值后,进入鉴别器,计算期望的相位调整值 δ c p \delta_{cp} δcp。具体的计算方法有三种:
-
幅值法
δ c p = 1 2 E − L E + L \delta_{cp}=\frac{1}{2}\frac{E-L}{E+L} δcp=21E+LE−L -
功率法
δ c p = 1 2 E 2 − L 2 E 2 + L 2 \delta_{cp}=\frac{1}{2}\frac{E^2-L^2}{E^2+L^2} δcp=21E2+L2E2−L2 -
点积功率法
δ c p = 1 2 E − L P \delta_{cp}=\frac{1}{2}\frac{E-L}{P} δcp=21PE−L
我们分别作出了这三种方法的图像:
![](https://i-blog.csdnimg.cn/blog_migrate/4e168f1094b4e52b10a5ee6b5fb04382.png)
可以看到,三条线共有三个交点。如果当前在最左侧的交点左侧,则说明 E、P、L 三点均在等腰三角形的左侧,无法继续计算;最右侧交点的右侧同理。
同时,我们可以看到幅值法是线性的——这是非常好的特性。因此,幅值法也是用得最多的方法。