基于FFT的微弱GPS信号频率精细估计

张洪伦 巴晓辉 陈 杰 周 航

1 引言

GPS接收机信号处理流程包括信号的捕获、跟踪和解算，为了尽快捕获到信号，频率搜索步长通常设置得比较大，捕获得到的频率误差在几十到几百赫兹，通常需要FLL（Frequency-Locked Loop）将频率牵引至一定范围后，再交给PLL（Phase-Locked Loop）进行跟踪。在弱信号情况下，环路等效噪声带宽设置比较小，使得环路锁定时间延长甚至无法正确牵引，影响接收机的TTFF（Time To First Fix）。同时，近年来高灵敏度GPS接收机设计[13]-受到广泛关注，文献[4,5]中提到了高灵敏度捕获方案，其中文献[4]能够捕获到20 dBHz的卫星信号，这样低的信号想要快速正确地转入跟踪状态，首先要取决于频率牵引能否成功。

已经有一些文献提出了针对频率精细估计的解决方案，文献[6,7]提到了基于相位测量的方法，文献[8]提出一种基于多项式拟合的方式，文献[9]提出用切比雪夫线性最小二乘曲线拟合法得到多普勒频移的精细化估计值，但是这些方法只适用于较强信号下的频率估计。文献[10-12]用线性调频 Z变换（CZT）对频率进行局部精细估计，提高了局部频率分辨率，但如果不克服比特跳变影响，也不适用于弱信号下频率估计。文献[6]在弱信号下用相干累加的同相分量作平方运算消除调制数据的影响，再作FFT运算，在理论上估计出的只是频率偏移的绝对值，并且忽略了正交分量的作用，性能无法达到最佳，而文献[13]中给出的基于FFT频率精细估计方法，充分利用相干累加能量值，对弱信号频率能够进行有效的牵引，文献[14]也提到这种方法，但是这两篇文章都没有对算法进行深入分析，并且只考虑了20 ms相干累加的情况，牵引范围只有[-12.5 Hz,12.5 Hz]。本文基于文献[13]，给出3种基于FFT的微弱GPS信号频率精细估计的建模方法，并且对这3种方法的性能进行分析对比。

文章由以下几部分构成，第1节为引言；第 2节为FFT频率估计原理；第3节为噪声性能分析；第4节为计算复杂度分析；第5节为实验结果；第6节对全文进行总结。

2 FFT频率估计原理

接收到的卫星信号可以表示为

P为接收的信号功率，D（ t）为导航数据，c（ t）为C/A码，fc为信号载波频率，n（ t）为高斯白噪声。

以fL表示捕获阶段得到的载波频率估计，cL（t）为本地产生的C/A码，且认为码相位与接收信号对齐，相干累加时间为 Tb，假设 Tb时间内没有比特反转， Ts为采样间隔， L = Tb/Ts，经过载波剥离和码剥离后，第k次相干累加结果为

其 中 Δ f = fc- fL, φ =π（fc- fL） （L - 1）Ts+φ,D（ k）表示第k次相干累加时的调制数据，NIk和 NQk是均值为零，方差为 σ2，并且相互独立的高斯随机变量。

从式（2）中得到，相干累加结果中还有调制的数据信息，并且每个比特周期为 20 ms，因此无法对式（2）的结果直接进行FFT运算获得频率偏移Δf的估计值。在现代化的GNSS信号中，为了获得更好的性能，一般都增加了导频通道，虽然导频通道没有调制数据信息，但是调制了二次码，二次码在频率牵引时还没有对齐，产生与数据信息相同的影响。

为了解决这个问题，需要对式（2）的结果进行处理，消除比特反转带来的影响，从而可以利用更长的数据进行FFT运算，以适应弱信号下频率估计。下面通过变量构造方法，消除比特反转影响，使得到的变量适用于FFT频率估计。

算法1 平方法构造变量

构造的新变量已经没有比特信息影响，频率和相位变为原来的两倍，噪声项为

算法2 差分法构造变量

构造的新变量中，还存在比特信息，后面将会对这个问题进行分析，噪声项为

算法 3 在比特边沿已知的情况下，每一个比特周期内构造一个变量，利用前10 ms相干累加和后10 ms相干累加，差分法构造变量

这种构造方法与算法2相同，但是比特信息影响被消除，代价是每20 ms的数据只利用了10 ms的能量，噪声项如算法2。

以算法1构造的变量 Xk为例，对FFT频率估计方法进行说明，分析过程忽略噪声影响。首先对Xk作N点FFT运算，并取模得

使（）f n幅值最大的n对应的谱线作为估计频率值，即

频率分辨率为 1 /2TbN，频率覆盖范围 - 1/4TbHz至 1 /4TbHz，由此得出，在FFT点数不变的情况下相干累加时间越长，频率分辨率越高，但是频率覆盖范围变小，在20 ms相干累加时频率覆盖范围仅为[-12.5 Hz, 12.5 Hz]。因此即使导频通道中不考虑二次码的影响，出于对频率覆盖范围的要求，也不可以任意增加相干累加时间长度，并且当载体相对卫星具有较大动态时，长时间相干累加会导致频谱扩展，计算出的峰值下降，使得发现概率减小。

不考虑比特信息的影响，算法2和算法3构造的变量，频率估计方法与分析方法相同。

3 噪声性能分析

通过变量构造，虽然消除了比特反转的影响，但是引入了更多的噪声项，不过构造的新变量，不受调制数据的影响，可以使用更长的数据长度，以提高检测弱信号的能力。下面对构造变量的噪声性能进行分析。由于 NIk, NIk-1,NQk和NQk-1是均值为零、方差为 σ2，并且相互独立的高斯随机变量，因此得到噪声项均值为 0。噪声方差的计算中用到一个结论，即若ξ1, ξ2, ξ3, ξ4是四元联合正态分布的随机变量，且均值都为0，则有

经过计算得到噪声方差为

在不考虑比特反转的情况下，并且相干累加时间相同，kX与kY的信号幅值相同，噪声对kY影响明显小于kX，理论上算法2具有更好的性能。

差分法构造变量的方法，由于前后两个相干结果的噪声项不相关，因此受到噪声影响要比平方法小，但是比特反转对算法2的影响相对算法1更严重。比特反转位置与相干累加时间段有两种关系，一种是相干累加没有跨越比特边界，另一种是相干累加时间跨越比特边界，以2 ms相干累加为例，分别如图1（a）和图1（b）所示，下面针对这两种情况分别进行讨论，分析这两种情况对算法1和算法2的影响。

图1 相干累加时间与比特反转关系示意图

（1）如图 1（a）所示，当所有相干累加时间段内不存在比特反转时，比特反转对算法1没有影响，而算法2构造的变量Yk在比特跳变处引入一个负号。考虑最坏情况，即每个比特边沿都有跳变，假设第n条线谱对应正确频率点，r为20 ms内相干累加次数，当2 ms相干累加时r等于10，则每r个相干累加结果中就会引入一个负号，则FFT运算峰值结果为

从式（16）得到，由于比特反转造成算法 2信号能量利用率并没有达到100%，对FFT峰值有贡献的能量仅为全部信号能量的（ 2）/r r- 。1 ms相干累加时，能量利用率为90%，而4 ms相干累加能量利用率为60%。

（2）如图1（b）所示，当有相干累加时间段跨越比特边界时，对算法1和算法2都会造成影响，在比特边界处构造的kX和kY会产生新的频率成分。对于Xk在20 ms内只会影响1个构造变量，对于 Yk则会影响连续的两个构造变量，受到影响的构造变量对FFT峰值不但没有贡献，还将产生不利影响。此时算法 1能量利用率为（r - 1 ）/r，算法 2还是（r - 2 ）/r。

由上面分析得出，算法1和算法2通过增加相干累加时间，可以提高性能，但是随着相干累加时间增加，由于比特边沿未知，能量利用率也随之下降，与增加相干累加时间产生相反的作用，并且增加相干累加时间，能够估计的频率范围减小。同时，幅值A随着频率偏移增大而减小，如果估计的频率刚好在两条谱线之间，也会造成频率泄漏，这些都对估计性能产生影响。

4 计算复杂度分析

各种频率估计方法，一般都是对相干累加结果进行处理，因此不考虑相干累加之前的各种运算。平方法和差分法构造变量都需要N次复数乘法运算，FFT 运算需要 N l og2N次复数加法和（N / 2）log2N次复数乘法运算，因此这两种构造变量方法用于FFT频率估计时运算量相同，平均对每一个相干累加结果进行处理的计算量为 lo g2N次复数加法和 1 + （log2N ）/2次复数乘法，表1列出N取1024, 2 ms相干累加与N取64, 20 ms相干累加的平均计算量，其中已将复数运算换算为实数运算。而传统的基于FLL的频率估计（牵引）方法，需要对每次相干累加结果进行鉴频和环路滤波，涉及到除法、反三角函数等运算，相位鉴别器输入为相干累加正交与同相分量比值，需要 1次除法，文献[15]给出一种快速除法算法，进行 3次迭代误差小于0.002，需要3次加法和6次乘法，文献[16]给出一种反正切函数实现方法，若逼近多项式取到五阶，则需要5次加法和9次乘法，计算频率差值需要2次加法和1次乘法，文献[17]给出的二阶环路滤波器需要3次加法和2次乘法，因此FLL的频率估计总运算量为13次加法和18次乘法。通过表1对比，可以得出本文给出的基于 FFT频率估计方法的运算量与传统方法相当。

表 1 FFT频率估计与FLL频率估计平均运算量对比

5 实验结果

为了验证前面方法的可行性，排除其它因素对FFT频率估计的影响，并且能够方便进行实验对比，产生的采样信号中频载波频率设定为1.25 MHz，采样率为5 MHz，采样点用浮点表示（没有量化），码相位已经对齐，各个实验所用的数据长度相同为1.28 s，当估计的频率值与实际频率偏差小于1.5 Hz时即认为频率估计正确。为了验证产生的中频率数据是正确的，将中频采样数据利用窄带-宽带方法[18]进行载噪比估计，所产生数据的信号强度与估计结果非常吻合。实验产生的中频数据如果存在调制信息时，考虑最坏情况，即每个比特都存在反转。

5.1 1 ms相干累加

设置频率偏移为 50 Hz，考虑没有数据调制和有数据调制两种情况，实验结果如图2所示。没有数据调制时，算法2在22 dBHz时频率正确估计概率在95%以上，比算法1高出1 dB，即使有调制数据，在最坏情况下，算法2依然好于算法1。1 ms相干累加，在有数据调制时，算法 1能量利用率100%，而算法2存在能量损失，但是其受噪声影响相对较小，即使10%的能量没有利用其性能依然比算法1高出1 dB。

5.2 2 ms相干累加

设置频率偏移为50 Hz，在调制数据后，2 ms相干累加，比特反转会产生两种影响，实验中考虑了所有情况。实验结果如图3所示，图中没有特别说明的认为一次相干累加时间内不存在比特反转，在相干累加时间内存在比特反转的，用偏移量来标

图2 1 ms相干累加时频率估计性能

记，例如图1（b）所示2 ms相干累加时，在第1 ms处出现比特反转，记1 ms偏移。实验结果显示，两种算法在没有数据调制情况下，算法2比算法1性能高出约2 dB，在有数据调制时，算法2性能有所下降，但是在21 dBHz时频率估计正确概率仍高于95%，比算法1高约1 dB。

5.3 4 ms相干累加

频率偏移设置为25 Hz，实验结果如图4所示。在不考虑数据调制时，算法2的性能在19 dBHz时能够达到95%以上的发现概率，但是由于比特反转，最坏情况下造成40%的能量损失，使得其性能开始低于算法1。算法1在21 dBHz时，最坏情况下的发现概率能够达到90%以上。

5.4 20 ms相干累加

在数据比特边沿已知的情况下，可以实现20 ms相干累加，实验结果如图5所示，频率偏移为10 Hz时，算法1和算法3在19 dBHz时发现概率在95%以上，频率偏移为5 Hz时，算法1在17 dBHz时仍然能够达到90%以上的发现概率。算法3构造变量的方法，在消除了比特边沿的同时，牺牲了一半的能量，因此性能相比算法1差一些。在5 Hz频率偏移的情况下，算法1发现概率比10 Hz频率偏移时高出1 dB，原因在于幅值A受到频率偏移的影响，频率偏移越大，幅值A越小，从噪声中检测出信号的概率越低。

通过前面的分析，结合实验结果，得到如下结论：

（1）在1 ms, 2 ms相干累加时，算法2比算法1性能高出1 dB, 1 ms相干累加时算法2在22 dBHz时发现概率能够超过 95%, 2 ms相干累加时在 21 dBHz时发现概率超过95%；

（2）4 ms相干累加时，算法1性能开始超过算法2，算法1在21 dBHz时发现概率在90%以上；

（3）20 ms相干累加时，算法 1与算法 3在 19 dBHz时，发现概率均超过95%。

图 3 2 ms相干累加时频率估计性能

图 4 4 ms相干累加时频率估计性能

图 5 20 ms相干累加时频率估计性能

由此得到，当要估计的频率偏差比较大时，可以使用算法2, 2 ms相干进行频率估计，当要估计的频率范围比较小时，并且信号非常弱，用算法 1，20 ms相干累加。实验中所用的数据长度是固定的，在实际应用中，可以通过增加数据长度的方式提高频率估计性能。实验中考虑比特反转时，假设比特反转每次都发生，实际中能量利用率要高于实验中的情况，因此得到的FFT峰值比实验中要高。

在没有比特边沿先验信息的情况下，又想得到比较高的性能时，需要20 ms相干累加，并且利用更长的数据，这时可以对20个可能的比特边界逐个进行检测，运用算法1或者算法3进行FFT运算，得到20个峰值，取最大峰值对应的估计频率作为最终结果。

本文给出的3种方法，能够在GPS信号强度较弱的情况下，根据不同情况和需求，经过一次FFT运算即可实现较大范围的频率估计。

6 结束语

本文针对传统 GPS接收机弱信号下频率牵引速度慢、成功率低等问题，给出 3种基于 FFT的GPS微弱信号频率精细估计的方法。算法经过了理论推导、噪声和性能分析，并且通过实验验证，结果表明本文给出的方法在没有比特边沿先验信息情况下，2 ms相干累加时差分法比平方法性能更好，前者能够实现21 dBHz信号强度下的频率估计，在有比特边沿先验信息情况下，差分法比平方法性能稍差，两种方法在信号强度低至19 dBHz时，频率估计正确率都能达到95%以上。

[1] 莫凡, 王新龙. 高灵敏度 GPS捕获技术研究综述[J]. 航空兵器, 2014（2）: 19-26.Mo Fan and Wang Xin-long. Review of acquisition techniques of high sensitivity GPS[J]. Aero Weaponry, 2014（2）: 19-26.

[2] Musumeci L, Dovis F, Silva P F, et al.. Design of a very high sensitivity acquisition system for a space GNSS receiver[C].IEEE/ION Position, Location and Navigation Symposium,California, 2014: 556-568.

[3] 李小捷. 弱信号环境下的GNSS信号捕获技术研究[D]. [博士论文], 西安电子科技大学, 2012.Li Xiao-jie. GNSS signal acquisition techniques in weak signal environment[D]. [Ph.D. dissertation], Xidian University,2012.

[4] 郑伦贵, 尤政, 张高飞, 等. 基于非相干积分的 GNSS弱信号捕获[J]. 清华大学学报（自然科学版）, 2014, 54（6）: 794-798.Zheng Lun-gui, You Zheng, Zhang Gao-fei, et al.. Acquisition of weak GNSS signals based on non-coherent integration[J].Journal of Tsinghua University（Science & Technology）, 2014,54（6）: 794-798.

[5] Wu Ling-juan, Lu Wei-jun, and Yu Dun-shan. Research of weak signal acquisition algorithms for high sensitivity GPS receivers[C]. Asia Pacific Conference on Postgraduate Research in Microelectronics & Electronics, Shanghai, 2009:173-176.

[6] 曾庆喜, 唐琳琳, 王庆, 等. GPS软件接收机宽载噪比载波频率精确捕获策略[J]. 中国惯性技术学报, 2013, 21（4）: 500-505.Zeng Qing-xi, Tang Lin-lin, Wang Qing, et al.. Acquisition of fine carrier frequency under wide carrier-to-noise ratio in GPS software receivers[J]. Journal of Chinese Inertial Technology, 2013, 21（4）: 500-505.

[7] 郑伦贵, 尤政, 张高飞, 等. GNSS接收机的一种快速捕获方法研究[J]. 仪器仪表学报, 2014, 35（4）: 807-813.Zheng Lun-gui, You Zheng, Zhang Gao-fei, et al.. Study on fast GNSS signal acquisition in GNSS receiver[J]. Chinese Journal of Scientific Instrument, 2014, 35（4）: 807-813.

[8] Zeng Dan-dan and Li Jin-hai. GPS signal fine acquisition algorithm[C]. International Conference on Information Science and Engineering （ICISE）, Hangzhou, 2010:3729-3732.

[9] 安毅, 蔡伯根, 宁滨, 等. 基于简化差分相干积累的北斗 B1频点信号精捕获算法[J]. 北京交通大学学报, 2014, 38（2）: 1-7.An Yi, Cai Bai-gen, Ning Bin, et al.. BeiDou B1 signals fine acquisition algorithm based on simplified differentially coherent integration[J]. Journal of Beijing Jiaotong University, 2014, 38（2）: 1-7.

[10] Aiello M, Cataliotti M and Nuccio S. A Chirp-Z transformbased synchronizer for power system measurements[J]. IEEE Transactions on Instrumentation and Measurement, 2005,54（3）: 1025-1032.

[11] 陈略, 唐歌实, 孟桥, 等. 一种卫星信号频率精确估计方法研究及应用[J]. 遥测遥控, 2010, 31（6）: 26-30.Chen Lue, Tang Ge-shi, Meng Qiao, et al.. Research and application of satellite signal frequency accurate estimation method[J]. Journal of Telemetry, Tracking and Command,2010, 31（6）: 26-30.

[12] 黄富彪, 何兵哲. 一种卫星信号载波频率精确估计算法[J]. 现代电子技术, 2012, 35（3）: 129-134.Huang Fu-biao and He Bing-zhe. A carrier frequency accurate estimation algorithm for satellite signal[J]. Modern Electronics Technique, 2012, 35（3）: 129-134.

[13] Ba Xiao-hui, Wang Yun, Chen Jie, et al.. An effective algorithm for pull-in and tracking processing of weak GNSS signal[C]. China Satellite Navigation Conference（CSNC）,Guangzhou, 2010: 601-605.

[14] 王勇松, 赵民建. 微弱 GPS信号精捕获算法[J]. 宇航学报,2010, 31（12）: 2671-2676.Wang Yong-song and Zhao Min-jian. Fine acquisition for weak GPS signal[J]. Journal of Astronautics, 2010, 31（12）:2671-2676.

[15] Baese U M. Digital Signal Processing with Field Programmable Gate Arrays[M]. Fourth Edition, New York:Springer, 2014: 102-103.

[16] 伍微, 刘小汇, 李峥嵘, 等. 实现定点DSP汇编层反正切函数的差分进化算法[J]. 系统工程与电子技术, 2005, 27（5）:926-928.Wu Wei, Liu Xiao-hui, Li Zheng-rong, et al.. Implementation of arc tangent function using assembly in fixed-point DSP based on differential evolution algorithm[J]. Systems Engineering and Electronics, 2005, 27（5）: 926-928.

[17] Kaplan D E and Christopher J H. Understanding GPS:Principles and Applications[M]. London: Artech House Inc,2006: 181.

[18] Parkinson B and Spilker J J. Global Positioning System:Theory and Applications[M]. Washington: AIAA, Inc, 1996:390-392.