CN103592662B - 一种gps信号接收机的载波跟踪方法及环路 - Google Patents

Abstract

本发明提出了一种GPS信号接收机的载波跟踪方法及环路，所述方法包含：步骤101，经过天线接收GPS卫星信号，然后对接收的GPS卫星信号进行下变频和相关处理，得到相关结果；步骤102，将相关结果输入鉴别器计算接收信号与本地信号的频率差或相位差，其中，当载波环为锁频环时，所述鉴别器为鉴频器；当载波环为锁相环时，所述鉴别器是鉴相器；步骤103，将鉴别结果采用卡尔曼滤波器进行滤波；步骤104，将滤波更新后的输出结果再输入环路滤波器进行滤波，进而输出载波多普勒频率，依据载波多普勒频率和载波频率产生本地载波。本发明利用改进的卡尔曼滤波对跟踪环中的载波环进行改进和优化，因此在简化计算量的同时，大大提高了对高灵敏度GPS信号的跟踪精度。

Description

一种GPS信号接收机的载波跟踪方法及环路

技术领域

本发明涉及卫星导航接收机领域，特别是涉及一种GPS信号接收机的载波跟踪方法及环路。

背景技术

GPS（Global Positioning System）具有全能性、全球性、全天候、连续性和实时性的导航、定位和定时功能，能为用户提供精密的三维坐标、速度和时间。GPS技术在民用、军用领域都应用广泛，产生了极好的经济、社会效益，成为信息时代不可缺少的重要技术手段。卫星导航接收机技术是GPS技术的重要组成部分，直接影响GPS的定位、授时等作用。

传统的导航接收机一般应用于地表、户外等导航信号强度较大的环境，在这些环境中导航信号没有受到遮蔽、多径等影响。在室内、隧道、密林等类似的工作环境中，导航信号受到20dB-30dB的衰减，这大大低于普通卫星导航接收机的跟踪门限。因此很有必要研究一种对高灵敏度GPS信号的跟踪方法。

不少文献中采用基于卡尔曼滤波算法的跟踪环路以实现对高灵敏度信号的跟踪。基于卡尔曼滤波算法的跟踪环路结构主要有以下两种形式：1、卡尔曼滤波代替原有的鉴别器；2、将卡尔曼滤波加在鉴别器和环路滤波器之间。所用到的卡尔曼滤波主要有线性卡尔曼滤波、扩展卡尔曼滤波（Extended Kalman Filter,EKF）和平淡卡尔曼滤波（Unscented Kalman Filter,UKF）三种，其中EKF和UKF是非线性的。以运动模型对所用卡尔曼滤波分类：常速度模型对应CV模型，常加速度模型对应CA模型。

虽然现有技术存在将卡尔曼滤波加在鉴别器和环路滤波器之间的情况，但是大多数情况是采用非线性卡尔曼滤波，采用这种非线性卡尔曼滤波虽然效果相对较好但是计算量很大；当采用线性卡尔曼滤波时，由于现有技术选取的观测量大于1个，计算量复杂；且没有采取抑制滤波器发散等算法，使得采用现有技术的卡尔曼滤波精度低且计算量很大。

发明内容

本发明的目的在于，为克服现有技术中的滤波算法复杂精度低的问题，从而提供一种GPS信号接收机的载波跟踪方法及环路。

为了实现上述目的，本发明提供一种GPS信号接收机的载波跟踪环路，包含：鉴别器和环路滤波器，其特征在于，所述鉴别器和环路滤波器之间还包含一个采用线性卡尔曼滤波算法的卡尔曼滤波器，且所述卡尔曼滤波器与环路滤波器之间还设置一条反馈通路，用于将环路滤波器的输出信号反馈至卡尔曼滤波器；

其中，所述线性卡尔曼滤波算法所选取的量测量和状态量之间的关系是线性的；

所述量测量为即观测量，指的是鉴频器的输出，具体为是接收信号与本地复制信号的频率差；

所述状态量为：基于CV模型的FLL中的状态量有2个，分别为接收信号与本地复制信号的频率差、接收信号的多普勒频率变化率；基于CA模型的FLL中的状态量有3个状态量，包括：接收信号与本地复制信号的频率差、接收信号的多普勒频率变化率及多普勒频率加速度；

所述量测量和状态量之间的关系式为：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}$

其中，为在k-1时刻对k时刻观测值的预测，为在k-1时刻对k时刻状态值的预测，H为观测矩阵。

可选的，上述卡尔曼滤波器包含：

预测模块，用于构造如下预测方程，基于构造的预测方程得到当前状态预测值和预测误差协方差矩阵；

状态预测方程：

$\begin{array}{c}{\hat{X}}_{(k/k-1)}={\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k/k-1)}=F\×{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\\ =F{\×X}_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\end{array}---\left(1\right)$

上述方程最左边为“k-1”时刻的状态预测值；

预测误差协方差矩阵：

$\begin{array}{c}{\hat{P}}_{(k/k-1)}={\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k/k-1)}=F\×{\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k-1/k-1)}\×F^T+D\×D^T\\ =F\×P_{(k-1/k-1)}\×F^T+Q\end{array}---\left(2\right)$

上述方程最左边的是“k-1”时刻的预测误差协方差矩阵；

量测预测方程：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}---\left(3\right)$

新息协方差矩阵：

$S_{\left(k\right)}=H\×{\hat{P}}_{(k/k-1)}\×H^T+R_{\left(k\right)}---\left(4\right)$

上述方程中最左边的S_(k)是k时刻的新息协方差矩阵，上述方程中的R_(k)是观测噪声v_(k)的协方差矩阵；

滤波模块，用于基于预测方程得到如下滤波方程，再基于滤波方程对观测量进行卡尔曼滤波；

增益矩阵：

$K_{\left(k\right)}={\hat{P}}_{(k/k-1)}\×H^T\×{S^{-1}}_{\left(k\right)}---\left(5\right)$

上述方程最左边的K_(k)是k时刻的增益矩阵；

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})---\left(6\right)$

上述方程最左边的X_(k/k)为k时刻的状态滤波值，Z_(k)为k时刻的量测值；

误差协方差矩阵：

$P_{(k/k)}={\hat{P}}_{(k/k-1)}-K_{\left(k\right)}\×S_{\left(k\right)}\×{K_{\left(k\right)}}^T;---\left(7\right)$

上述方程最左边的P_(k/k)为k时刻的误差协方差矩阵。

可选的，上述卡尔曼滤波器还包含：噪声统计模块，野值判别模块，发散判别和抑制模块；

噪声统计模块，用于预测观测噪声；

野值判别模块，用于基于预测的观测噪声检测野值，并对检测得到的野值进行剔除操作；

发散判别和抑制模块，用于当野值判别模块判断得出某观测量不是野值时，再进行滤波器是否发散的判断；具体判断策略为：将第k次最终的卡尔曼滤波结果与第k次观测量进行比较，如果两者差值的绝对值小于设定的门限，则滤波器工作正常，否则判定滤波器发散，当滤波器发散时，将第k次观测值作为第k次滤波值输出，并初始化卡尔曼滤波器。

可选的，上述噪声统计模块进一步包含：

参考量选取单元，用于将观测量的前若干个点做为参考量；

均值获取单元，用于对选定的参考量求取平均值；

再将参考点分别与所求平均值做差，得到一组由若干个差值组成的数据集合，将这个集合的方差作为该观测量的观测噪声方差；

其中，所述观测量为鉴别器的输出值。

可选的，上述野值判别模块进一步包含：

第一处理单元，用于获得实际观测值与速度差预测值的差值的绝对值；

第二处理单元，用于获得新息正交性预测得到的观测值的若干倍与速度差预测值的差值的绝对值，其中所述若干倍为2-6倍；

比较判别单元，用于将第一处理单元和第二处理单元的输出结果进行比较；如果第一处理单元的处理结果大于第二处理单元的处理结果，则判定该观测值为野值；和

校正处理单元，用如下的活化函数f_(k)对判定为野值的观测量进行校正，避免滤波器的发散。

此外，本发明还提供了一种GPS信号接收机的载波跟踪方法，所述方法包含：

步骤101，经过天线接收GPS卫星信号，然后对接收的GPS卫星信号进行下变频和相关处理，得到相关结果；

步骤102，将相关结果输入鉴别器计算接收信号与本地信号的频率差或相位差，其中，当载波环为锁频环时，所述鉴别器为鉴频器；当载波环为锁相环时，所述鉴别器是鉴相器；

步骤103，将鉴别结果采用卡尔曼滤波器进行滤波，其中，所述的卡尔曼滤波采用如下状态滤波方程进行滤波：

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})$

其中，为采用状态预测方程得到的函数值，K_(k)为卡尔曼滤波器的增益矩阵，为采用量测预测方程的函数值，Z_(k)为量测值且所有量测值组成一个一维矩阵，k为卡尔曼滤波次数；

步骤104，将滤波更新后的输出结果再输入环路滤波器进行滤波，进而输出载波多普勒频率，依据载波多普勒频率和载波频率产生本地载波。

上述步骤101进一步包含：

步骤101-1，经过天线接收GPS卫星信号，然后对接收的GPS卫星信号进行射频下变频处理，得到数字中频信号；

步骤101-2，将数字中频信号分成两路进行混频，即，将其中一路和本地载波发生器生成的同相载波相乘，将另一路和本地载波发生器生成的正交载波相乘；

步骤101-3，将混频信号分别与超前码、即时码和滞后码相乘累加完成相关运算，得到相关结果。

可选的，当载波环为锁频环时，基于匀速模型的卡尔曼滤波过程采用如下步骤和公式得到状态滤波方程：

步骤201）选取状态量和观测量，所述状态量为接收信号与本地信号的频率差和接收信号的多普勒变化率；观测量为鉴频器的鉴频输出结果；

步骤202）建立如下的状态方程和量测方程如下：

状态方程：

$\begin{array}{c}{X}_{\left(k\right)}={\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{\left(k\right)}=\left(\begin{array}{cc}1& \mathrm{\ΔT}\\ 0& 1\end{array}\right)\×X_{(k-1)}+\left[\begin{array}{c}-1\\ 0\end{array}\right]\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}+\left(\begin{array}{c}{\mathrm{\ΔT}}^2/2\\ \mathrm{\ΔT}\end{array}\right)\×w_{(k-1)}\\ =F\×X_{(k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}+D\×w_{(k-1)}\end{array}$

其中，x_f表示接收信号与本地信号的频率差，x_a为多普勒频率变化率；w_(k-1)为过程噪声，表示多普勒频率加速度；ΔT为积分清除时间；F是从t_k-1到t_k的状态转移矩阵；G为输入矩阵；D为过程噪声的转移矩阵；为环路滤波器的输出值先经过数据处理模块再反馈为卡尔曼滤波模块的输入值，为前后两次的差值；

量测方程：

y_(k)＝[1 ΔT]×X_(k)+v_(k)

＝H×X_(k)+v_(k)

其中，H为观测矩阵；y_(k)为观测量，是鉴频器的输出；v_(k)为k时刻的观测值与实际信号值的误差；

步骤203）基于状态方程和量测方程得到如下预测方程：

状态预测方程：

$\begin{array}{c}{\hat{X}}_{(k/k-1)}{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k/k-1)}=F\×{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\\ =F\×X_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\end{array}$

预测协方差矩阵：

$\begin{array}{c}{\hat{P}}_{(k/k-1)}={\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k/k-1)}=F\×{\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k-1/k-1)}\×F^T+D\×D^T\\ =F\×P_{(k-1/k-1)}\×F^T+Q\end{array}$

量测预测方程：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}$

新息协方差矩阵：

$S_{\left(k\right)}=H\×{\hat{P}}_{(k/k-1)}\×H^T+R_{\left(k\right)};$

步骤204）基于预测方程得到增益矩阵，进而得到滤波方程和误差协方差方程：

增益矩阵：

$K_{\left(k\right)}={\hat{P}}_{(k/k-1)}\×H^T\×{S^{-1}}_{\left(k\right)}$

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})$

误差协方差矩阵：

$P_{(k/k)}={\hat{P}}_{(k/k-1)}-K_{\left(k\right)}\×S_{\left(k\right)}\×{K_{\left(k\right)}}^T$

其中，为在“k-1”时刻对k时刻状态值的预测，是“k-1”时刻的预测协方差矩阵，S_(k)是k时刻的新息协方差矩阵，R_(k)是观测噪声v_(k)的协方差矩阵，K_(k)是k时刻的增益矩阵，Z_(k)为k时刻的量测值。本发明有益效果如下：

本发明利用改进的卡尔曼滤波对跟踪环中的载波环进行改进和优化，因此在简化计算量的同时，大大提高了对高灵敏度GPS信号的跟踪精度。总之本发明提出基于线性卡尔曼滤波的载波环路能够滤除载波环路中鉴别结果中的噪声，提高信噪比。

附图说明

图1是传统载波跟踪环路逻辑框图；

图2是本发明实施例中卡尔曼滤波的FLL逻辑框图；

图3是本发明实施例中卡尔曼滤波的PLL逻辑框图；

图4是本发明实施例中卡尔曼滤波工作原理简述框图；

图5是本发明实施例中采用改进算法的卡尔曼滤波的实现流程图。

具体实施方式

下面结合附图对本发明的技术方案进行详细阐述。

本发明解决的技术问题是改进接收机跟踪环路以实现对高灵敏度GPS信号的跟踪。

现有技术的跟踪环路由载波跟踪环路和码环(DLL)组成。载波跟踪环路复制载波与接收信号进行混频以实现载波剥离及将信号下变频到基带。载波环又分为锁相环(PLL)和锁频环(FLL)两种，在动态应力较大的情况下适合用锁频环FLL，而在对精度要求高的时候则应该选择用PLL。为了在不同的环境中能达到更好地跟踪效果，现有的跟踪技术中往往是将FLL和PLL以某种方式结合起来使用。码环复制伪码与接收信号相乘以实现伪码剥离和信号解扩。

为了实现高灵敏度GPS信号跟踪的技术问题，本发明的技术方案是采用线性卡尔曼滤波对载波环路进行改进。现有大部分文献资料是将卡尔曼滤波应用在PLL中，而本发明不仅将卡尔曼滤波应用在PLL中且还应用在FLL中。此外，本发明采用的是线性卡尔曼滤波而不是非线性的EKF和UKF，大大减少了计算量；不仅如此，在线性卡尔曼滤波中采用的是一维观测量，这又进一步减少了计算量。为了防止在某些特殊情况下卡尔曼滤波造成环路失锁，在卡尔曼滤波中加入了三点算法，分别为噪声统计算法、野值的判别与修正算法和抑制发散算法。

本发明所述方法包括以下几个步骤：

步骤S101，经过天线接收GPS卫星信号，然后对接收的GPS卫星信号进行射频下变频处理，得到数字中频信号。

步骤S102，将数字中频信号分成两路，一路和本地载波发生器生成的同相载波相乘，另一路和本地载波发生器生成的正交载波相乘。

步骤S103，混频信号分别与超前码(E)、即时码(P)和滞后码(L)相乘累加完成相关运算。

步骤S104，相关结果分别输出给码环鉴别器和载波环鉴别器，FLL对应的是鉴频器，PLL对应的是鉴相器。

步骤S105，码环鉴别器的鉴别结果输出给码环环路滤波器；载波环鉴别器的鉴别结果输出给卡尔曼滤波模块，经过卡尔曼滤波更新后的鉴别结果再输出给载波环环路滤波器。另外，载波环中还存在一条从环路滤波器到卡尔曼滤波模块的反馈回路，这是一条必不可少、至关重要的反馈回路。

步骤S106，码环环路滤波器的输出为码多普勒频率，码多普勒频率加上码频率一起输出给码数控振荡器产生本地码；载波环环路滤波器的输出为载波多普勒频率，载波多普勒频率加上载波频率一起输出给载波环数控振荡器产生本地载波。

重复步骤S101-S106，直到鉴别器输出值无限接近零。

进一步，步骤S105中载波环FLL的卡尔曼滤波过程包括以下步骤：

步骤1，为卡尔曼滤波选取合适的状态量和观测量。基于CV模型卡尔曼滤波的FLL中有两个状态量，分别为接收信号与本地信号的频率差、接收信号的多普勒变化率；观测量为鉴频器的鉴频结果。基于CA模型卡尔曼滤波的FLL有三个状态量，比CV模型多了一个多普勒频率加速度状态量。以下步骤均基于CV模型，CA模型与CV模型的逻辑步骤基本保持一致，不再赘述。

步骤2，建立状态方程和量测方程如下：

状态方程：

$\begin{array}{c}{X}_{\left(k\right)}={\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{\left(k\right)}=\left(\begin{array}{cc}1& \mathrm{\ΔT}\\ 0& 1\end{array}\right)\×X_{(k-1)}+\left[\begin{array}{c}-1\\ 0\end{array}\right]\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}+\left(\begin{array}{c}{\mathrm{\ΔT}}^2/2\\ \mathrm{\ΔT}\end{array}\right)\×w_{(k-1)}\\ =F\×X_{(k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}+D\×w_{(k-1)}\end{array}$

其中，x_f表示接收信号与本地信号的频率差，x_a为多普勒频率变化率；w_(k-1)为过程噪声，表示多普勒频率加速度；ΔT为积分清除时间。F是从t_k-1到t_k的状态转移矩阵；G为输入矩阵；D为过程噪声的转移矩阵；为环路滤波器的输出值先经过数据处理模块再反馈为卡尔曼滤波模块的输入值，为前后两次的差值。

量测方程：

y_(k)＝[1 ΔT]×X_(k)+v_(k)

＝H×X_(k)+v_(k)

其中，H为观测矩阵；y_(k)为观测量，是鉴频器的输出；v_(k-1)为观测噪声，即叠加在鉴频器输出值上的噪声。

步骤3，卡尔曼滤波可以分为两个部分：预测过程和滤波过程。预测过程负责及时向前推算当前状态变量（即，状态量）和误差协方差估计值，以便为下一个时间状态构造先验估计；滤波过程负责反馈预测值，将先验估计和新的观测量结合以构造改进的后验估计。因此步骤3根据状态方程和量测方程推写出预测方程并加以实现。预测方程如下：

状态预测方程：

$\begin{array}{c}{\hat{X}}_{(k/k-1)}{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k/k-1)}=F\×{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\\ =F\×X_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\end{array}$

预测协方差矩阵：

$\begin{array}{c}{\hat{P}}_{(k/k-1)}={\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k/k-1)}=F\×{\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k-1/k-1)}\×F^T+D\×D^T\\ =F\×P_{(k-1/k-1)}\×F^T+Q\end{array}$

量测预测方程：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}$

新息协方差矩阵：

$S_{\left(k\right)}=H\×{\hat{P}}_{(k/k-1)}\×H^T+R_{\left(k\right)}$

步骤4，推写出滤波方程并加以实现，滤波方程如下：

增益矩阵：

$K_{\left(k\right)}={\hat{P}}_{(k/k-1)}\×H^T\×{S^{-1}}_{\left(k\right)}$

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})$

误差协方差矩阵：

$P_{(k/k)}={\hat{P}}_{(k/k-1)}-K_{\left(k\right)}\×S_{\left(k\right)}\×{K_{\left(k\right)}}^T$

进一步，步骤S105中载波环PLL的卡尔曼滤波过程和FLL的基本相似，CV模型卡尔曼滤波在PLL和FLL中有以下4点区别：

1、状态量的不同，PLL中卡尔曼滤波的中状态量为接收信号与本地复制信号的相位差、接收信号的频率；

2、输入矩阵不一样，PLL中卡尔曼滤波的输入矩阵G的第一行第一列为时间ΔT，而FLL中卡尔曼滤波的输入矩阵G的第一行第一列为-1；

3、PLL中卡尔曼滤波观测量为鉴相器输出，即接收信号与本地复制信号的实际相位差；

4、PLL中卡尔曼滤波的输入值直接等于环路滤波器的输出频率。

PLL中的卡尔曼滤波CA模型和FLL中也存在以上4点区别，其它的逻辑步骤和FLL的卡尔曼滤波CA模型一致。

实施例

图2所示为本发明设计的加入了卡尔曼滤波的FLL，通过对比图1的传统载波跟踪环路结构可知，本发明改进的FLL在鉴频器和环路滤波器之间加入了卡尔曼滤波，且存在一条从环路滤波器经过数据处理模块到卡尔曼滤波模块的反馈回路，这是一条必不可少、至关重要的反馈回路。图3所示为本发明设计的加入了卡尔曼滤波的PLL，在鉴相器和环路滤波器之间加入了卡尔曼滤波，同样存在一条反馈回路，与FLL不同的是中间没有数据处理环节。

如图4所示，卡尔曼滤波分为预测和滤波两个过程。在估算出状态变量和误差协方差矩阵初值的前提下，向前推算出状态变量和误差协方差矩阵，利用推算出的预测量计算卡尔曼滤波增益，并且利用观测量更新估计值和误差协方差。

整个卡尔曼滤波的具体实现过程如图5所示，将实现过程划分为了六个模块，分别为顶层控制模块、噪声统计模块、预测模块、野值判别模块、滤波模块和发散判别、抑制模块。

模块一顶层控制模块

顶层控制模块的功能主要有两个：连接各子模块之间的信号，控制各子模块的启止时间。

模块二噪声统计模块

在卡尔曼滤波的实现过程中，如何预测观测噪声是非常关键的一步，在很大程度上决定了所设计的卡尔曼滤波性能。

将观测量(鉴别值)的前15个点做为参考量来对该观测量进行噪声估计，这15个点的观测值则直接作为卡尔曼滤波输出值。具体方法是先对参考点求平均值，再将参考点分别与所求平均值做差，得到一组由15个差值组成的数据，将这组数据的方差作为该观测量的观测噪声方差。

模块三预测模块

预测过程负责及时向前推算当前状态变量和误差协方差估计值，以便为下一个时间状态构造先验估计。逐一实现以下预测方程。

$\begin{array}{c}{\hat{X}}_{(k/k-1)}{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k/k-1)}=F\×{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\\ =F\×X_{(k-1/k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}\end{array}$

$\begin{array}{c}{\hat{P}}_{(k/k-1)}={\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k/k-1)}=F\×{\left(\begin{array}{cc}{p}_{\mathrm{ff}}& p_{\mathrm{fa}}\\ p_{\mathrm{af}}& p_{\mathrm{aa}}\end{array}\right)}_{(k-1/k-1)}\×F^T+D\×D^T\\ =F\×P_{(k-1/k-1)}\×F^T+Q\end{array}$

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}$

$S_{\left(k\right)}=H\×{\hat{P}}_{(k/k-1)}\×H^T+R_{\left(k\right)};$

模块四野值判别模块

在观测值中，如果某个观测值与预测值相差太大，则称此观测值为野值。野值的出现将导致滤波的不准确，甚至导致滤波器的发散，因此如何有效地剔除野值是卡尔曼滤波的关键。

在本算法中，对野值的处理主要是用活化函数对观测值加权，具体的算法如下：

${Z\_\mathrm{Quadrature}}_{\left(k\right)}=\sqrt{S_{\left(k\right)}+H\×{\stackrel{\widehat{\→}}{X}}_{(k/k-1)}\×{{\stackrel{\widehat{\→}}{X}}_{(k/k-1)}}^T\×H^T}---\left(10\right)$

其中，Z_Quadrature_(k)为由新息正交性即信号与噪声的正交性预测得到的观测值。

f_(k)＝Z_Quadrature_(k)/Z_(k) (11)

将实际观测值Z_(k)与速度差预测值（的第一个值）的差值的绝对值和由新息正交性预测得到的观测值Z_Quadrature_(k)的若干倍与速度差预测值的差值的绝对值进行比较；如果前者大于后者，则判定该观测量为野值，用上述的活化函数f_(k)对实际观测量进行校正，避免滤波器的发散。

模块五滤波模块

滤波过程负责反馈预测值，将先验估计和新的观测量结合以构造改进的后验估计。逐一实现以下滤波方程。

$K_{\left(k\right)}={\hat{P}}_{(k/k-1)}\×H^T\×{S^{-1}}_{\left(k\right)}$

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})$

$P_{(k/k)}={\hat{P}}_{(k/k-1)}-K_{\left(k\right)}\×S_{\left(k\right)}\×{K_{\left(k\right)}}^T$

模块六发散判别、抑制模块

卡尔曼滤波的发散会带来非常严重的后果，不仅卡尔曼滤波失效，且整个跟踪环路也会随着卡尔曼滤波器的发散而失锁。因此如何抑制卡尔曼滤波器的发散是成败的关键。

如果某观测量被判为不是野值则要进行滤波器是否发散的判断，否则不进行判断。将第k次最终的卡尔曼滤波结果与第k次观测量比较大小，如果两者差值的绝对值小于所设定的门限，则认为滤波器工作正常，没有发散，否则认为发散。若判为发散，将第k次观测值作为第k次滤波值输出，并初始化卡尔曼滤波并且重新启动。

由上述实施例可以看出，本发明通过改进接收机跟踪环路，大大提高了对高灵敏度GPS信号的跟踪精度。

尽管为示例目的，已经公开了本发明的优选实施例，本领域的技术人员将意识到各种改进、增加和取代也是可能的，因此，本发明的范围应当不限于上述实施例。

Claims (7)

1.一种GPS信号接收机的载波跟踪环路，包含：鉴频器和环路滤波器，其特征在于，所述鉴频器和环路滤波器之间还包含一个采用线性卡尔曼滤波算法的卡尔曼滤波器，且所述卡尔曼滤波器与环路滤波器之间还设置一条反馈通路，用于将环路滤波器的输出信号反馈至卡尔曼滤波器；

其中，所述线性卡尔曼滤波算法所选取的量测量和状态量之间的关系是线性的；

所述量测量即为观测量，指的是鉴频器的输出，具体为接收信号与本地复制信号的频率差；

所述状态量为：基于CV模型的FLL中的状态量有2个，分别为接收信号与本地复制信号的频率差、接收信号的多普勒频率变化率；基于CA模型的FLL中的状态量有3个，包括：接收信号与本地复制信号的频率差、接收信号的多普勒频率变化率及多普勒频率加速度；

所述量测量和状态量之间的关系式为：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}$

其中，为在k-1时刻对k时刻观测值的预测，为在k-1时刻对k时刻状态值的预测，H为观测矩阵。

2.根据权利要求1所述的GPS信号接收机的载波跟踪环路，其特征在于，所述卡尔曼滤波器包含：

预测模块，用于构造如下预测方程，基于构造的预测方程得到当前状态预测值和预测误差协方差矩阵；

状态预测方程：

$\begin{array}{c}{\hat{X}}_{(k/k-1)}={\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k/k-1)}=F\×{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k-1/k-1)}+G\×{\mathrm{\ωre}}_{(k-1)}\\ =F\×X_{(k-1/k-1)}+G\×{\mathrm{\ωre}}_{(k-1)}\end{array}---\left(1\right)$

上述方程最左边为“k-1”时刻对k时刻状态值的预测；

预测误差协方差矩阵：

$\begin{array}{c}{\hat{P}}_{(k/k-1)}={\left(\begin{array}{cc}{p}_{ff}& p_{fa}\\ p_{af}& p_{aa}\end{array}\right)}_{(k/k-1)}=F\×{\left(\begin{array}{cc}{p}_{ff}& p_{fa}\\ p_{af}& p_{aa}\end{array}\right)}_{(k-1/k-1)}\×F^T+D\×D^T\\ =F\×P_{(k-1/k-1)}+F^T+Q\end{array}---\left(2\right)$

上述方程最左边的是“k-1”时刻的预测误差协方差矩阵；

量测预测方程：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}---\left(3\right)$

新息协方差矩阵：

$S_{\left(k\right)}=H\×{\hat{P}}_{(k/k-1)}\×H^T+R_{\left(k\right)}---\left(4\right)$

上述方程中最左边的S_(k)是k时刻的新息协方差矩阵，上述方程中的R_(k)是观测噪声v_(k)的协方差矩阵；

滤波模块，用于基于预测方程得到如下状态滤波方程，再基于滤波方程对观测量进行卡尔曼滤波；

增益矩阵：

$K_{\left(k\right)}={\hat{P}}_{(k/k-1)}\×H^T\×{S^{-1}}_{\left(k\right)}---\left(5\right)$

上述方程最左边的K_(k)是k时刻的增益矩阵；

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})---\left(6\right)$

上述方程最左边的X_(k/k)为k时刻的状态滤波值，Z_(k)为k时刻的量测值_；

误差协方差矩阵：

$P_{(k/k)}={\hat{P}}_{(k/k-1)}-K_{\left(k\right)}\×S_{\left(k\right)}\×{K_{\left(k\right)}}^T;---\left(7\right)$

上述方程最左边的P_(k/k)为k时刻的误差协方差矩阵；

其中，x_f表示接收信号与本地复制信号的频率差，x_a为多普勒频率变化率，F是从t_k-1到t_k的状态转移矩阵；G为输入矩阵，为环路滤波器的输出值先经过数据处理模块再反馈为卡尔曼滤波模块的输入值，为前后两次的差值，D为过程噪声的转移矩阵。

3.根据权利要求1所述的GPS信号接收机的载波跟踪环路，其特征在于，所述卡尔曼滤波器还包含：噪声统计模块，野值判别模块，发散判别和抑制模块；

噪声统计模块，用于预测观测噪声；

野值判别模块，用于基于预测的观测噪声检测野值，并对检测得到的野值进行剔除操作；

发散判别和抑制模块，用于当野值判别模块判断得出某观测量不是野值时，再进行滤波器是否发散的判断；具体判断策略为：将第k次最终的卡尔曼滤波结果与第k次观测量进行比较，如果两者差值的绝对值小于设定的门限，则滤波器工作正常，否则判定滤波器发散，当滤波器发散时，将第k次观测量作为第k次滤波值输出，并初始化卡尔曼滤波器。

4.根据权利要求3所述的GPS信号接收机的载波跟踪环路，其特征在于，所述噪声统计模块进一步包含：

参考量选取单元，用于将观测量的前若干个点做为参考量；

均值获取单元，用于对选定的参考量求取平均值；

再将参考量分别与所求平均值做差，得到一组由若干个差值组成的数据集合，将这个集合的方差作为该观测量的观测噪声方差；

其中，所述观测量为鉴频器的输出值。

5.根据权利要求3所述的GPS信号接收机的载波跟踪环路，其特征在于，所述野值判别模块进一步包含：

第一处理单元，用于获得实际观测量与速度差预测值的差值的绝对值；

第二处理单元，用于获得新息正交性预测得到的观测量的若干倍与速度差预测值的差值的绝对值，其中所述若干倍为2-6倍；

比较判别单元，用于将第一处理单元和第二处理单元的输出结果进行比较；如果第一处理单元的处理结果大于第二处理单元的处理结果，则判定该观测量为野值；和

校正处理单元，用活化函数f_(k)对判定为野值的观测量进行校正，避免滤波器的发散。

6.一种GPS信号接收机的载波跟踪方法，所述方法包含：

步骤101，经过天线接收GPS卫星信号，然后对接收的GPS卫星信号进行下变频和相关处理，得到相关结果；

步骤102，将相关结果输入鉴别器计算接收信号与本地信号的频率差或相位差，其中，当载波环为锁频环时，所述鉴别器为鉴频器；当载波环为锁相环时，所述鉴别器是鉴相器；

步骤103，将鉴别结果采用卡尔曼滤波器进行滤波，其中，所述的卡尔曼滤波采用如下状态滤波方程进行滤波：

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})$

其中，为采用状态预测方程得到的函数值，K_(k)为卡尔曼滤波器的增益矩阵，为采用量测预测方程得到的函数值，Z_(k)为量测值且所有量测值组成一个一维矩阵，k为卡尔曼滤波次数；

步骤104，将滤波更新后的输出结果再输入环路滤波器进行滤波，进而输出载波多普勒频率，依据载波多普勒频率和载波频率产生本地载波；

其中，

当载波环为锁频环时，基于匀速模型的卡尔曼滤波过程采用如下步骤和公式得到状态滤波方程：

步骤201)选取状态量和观测量，所述状态量为接收信号与本地信号的频率差和接收信号的多普勒变化率；观测量为鉴频器的鉴频输出结果；

步骤202)建立如下的状态方程和量测方程：

状态方程：

$\begin{array}{c}{X}_{\left(k\right)}={\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{\left(k\right)}=\left(\begin{array}{cc}1& \mathrm{\Δ}T\\ 0& 1\end{array}\right)\×X_{(k-1)}+\left[\begin{array}{c}-1\\ 0\end{array}\right]\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}+\left(\begin{array}{c}\mathrm{\Δ}{T}^2\mathrm{/2}\\ \mathrm{\Δ}T\end{array}\right)\×w_{(k-1)}\\ =F\×X_{(k-1)}+G\×{\mathrm{\ω}}_{{\mathrm{re}}_{(k-1)}}+D\×w_{(k-1)}\end{array}$

其中，x_f表示接收信号与本地信号的频率差，x_a为多普勒频率变化率；w_(k-1)为过程噪声，表示多普勒频率加速度；ΔT为积分清除时间；F是从t_k-1到t_k的状态转移矩阵；G为输入矩阵；D为过程噪声的转移矩阵；为环路滤波器的输出值先经过数据处理模块再反馈为卡尔曼滤波模块的输入值，为前后两次的差值；

量测方程：

y_(k)＝[1 ΔT]×X_(k)+v_(k)

＝H×X_(k)+v_(k)

其中，H为观测矩阵；y_(k)为观测量，是鉴频器的输出；v_(k)为k时刻的观测量与实际信号值的误差；

步骤203)基于状态方程和量测方程得到如下预测方程：

状态预测方程：

$\begin{array}{c}{\hat{X}}_{(k/k-1)}={\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k/k-1)}=F\×{\left[\begin{array}{c}{x}_f\\ x_a\end{array}\right]}_{(k-1/k-1)}+G\×{\mathrm{\ωre}}_{(k-1)}\\ =F\×X_{(k-1/k-1)}+G\×{\mathrm{\ωre}}_{(k-1)}\end{array}$

预测误差协方差矩阵：

$\begin{array}{c}{\hat{P}}_{(k/k-1)}={\left(\begin{array}{cc}{p}_{ff}& p_{fa}\\ p_{af}& p_{aa}\end{array}\right)}_{(k/k-1)}=F\×{\left(\begin{array}{cc}{p}_{ff}& p_{fa}\\ p_{af}& p_{aa}\end{array}\right)}_{(k-1/k-1)}\×F^T+D\×D^T\\ =F\×P_{(k-1/k-1)}+F^T+Q\end{array}$

量测预测方程：

${\hat{Z}}_{(k/k-1)}=H\×{\hat{X}}_{(k/k-1)}$

新息协方差矩阵：

$S_{\left(k\right)}=H\×{\hat{P}}_{(k/k-1)}\×H^T+R_{\left(k\right)};$

步骤204)基于预测方程得到增益矩阵，进而得到状态滤波方程和误差协方差矩阵：

增益矩阵：

$K_{\left(k\right)}={\hat{P}}_{(k/k-1)}\×H^T\×{S^{-1}}_{\left(k\right)}$

状态滤波方程：

$X_{(k/k)}={\hat{X}}_{(k/k-1)}+K_{\left(k\right)}\×(Z_{\left(k\right)}-{\hat{Z}}_{(k/k-1)})$

误差协方差矩阵：

$P_{(k/k)}={\hat{P}}_{(k/k-1)}-K_{\left(k\right)}\×S_{\left(k\right)}\×{K_{\left(k\right)}}^T$

其中，为在“k-1”时刻对k时刻状态值的预测，是“k-1”时刻的预测误差协方差矩阵，S_(k)是k时刻的新息协方差矩阵，R_(k)是观测噪声v_(k)的协方差矩阵，K_(k)是k时刻的增益矩阵，Z_(k)为k时刻的量测值。

7.根据权利要求6所述的GPS信号接收机的载波跟踪方法，其特征在于，所述步骤101进一步包含：

步骤101-1，经过天线接收GPS卫星信号，然后对接收的GPS卫星信号进行射频下变频处理，得到数字中频信号；

步骤101-2，将数字中频信号分成两路进行混频，即，将其中一路和本地载波发生器生成的同相载波相乘，将另一路和本地载波发生器生成的正交载波相乘；

步骤101-3，将混频信号分别与超前码、即时码和滞后码相乘累加完成相关运算，得到相关结果。