全球卫星导航系统(GPS)提供的导航服务已经广泛应用到全球各个领域，为了满足在城市峡谷、密林等信号衰减严重环境下的信号跟踪定位，对GPS接收机设计提出新的挑战^[1-2]。

由于载波跟踪环是独立GPS接收机最脆弱的环节^[3]，一般情况下先于码跟踪环失锁，因此载波跟踪环在某种程度上决定了未受辅助的GPS接收机的性能，成为设计的关键。由此，文献[4-8]针对如何提高载波跟踪灵敏度进行了研究，其中文献[4-7]将卡尔曼滤波及其扩展方法应用于载波跟踪，文献[8]给出一种基于FFT的高灵敏度GPS跟踪算法，这些方法都在一定程度上改善了跟踪灵敏度。

在单点定位时，接收端在同一时刻为所有可视卫星所共有，不同卫星信号之间具有潜在相关性，同时信号跟踪与导航解算之间也具有相关性，跟踪信息可用于用户状态估计，用户状态与卫星状态决定伪码相位与载波频率偏移^[9-10]，但文献[4-8]都忽略了这些相关性。文献[10]引入矢量跟踪的概念，而传统方法称为标量跟踪，提出矢量延迟锁定环(vector delay lock loop，VDLL)，通过一个滤波器同时对所有通道码跟踪环进行处理，代替分立的延迟锁定环(delay lock loop，DLL)，充分利用不同通道间的相关性。文献[11]使用经验跟踪门限的方法对VDLL和矢量延迟/频率锁定环(vector delay/frequency lock loop，VDFLL)的热噪声性能进行分析，通过理论推导证明矢量跟踪性能优于标量跟踪。文献[12-15]针对VDFLL及其在弱信号环境下的性能进行分析研究，但针对VFLL的研究很少。由于载波环比码环更容易失锁，并且考虑到VFLL实现较VDFLL简单，本文针对VFLL及其在微弱GPS环境下的信号跟踪进行研究，相比现有文献详细的给出了VFLL理论推导及实现过程，并以最小二乘估计方法证明VFLL在载波跟踪性能上优于FLL。

1 VFLL基本原理

基于位置解算完成的前提条件下，给出方向余弦矩阵，同时根据获得的星历数据计算当前时刻卫星的速度，就可以实现VFLL，进行速度估计和频率跟踪。下面给出具体实现过程。

在ECEF(earth-centered, earth-fined)坐标系下，用户速度v_u=$\left[ {{{\dot x}_{\rm{u}}}\;{{\dot y}_{u}}\;{{\dot z}_{\rm{u}}}} \right]$，估计值记为${{\mathit{\boldsymbol{\hat v}}}_{\rm{u}}}$=$\left[ {{{\hat{\dot{x}}}}_{u}}\ {{{\hat{\dot{y}}}}_{u}}\ {{{\hat{\dot{z}}}}_{\rm{u}}} \right]$。第i颗卫星速度v_{s, i}=$\left[ {{{\dot{x}}}_{\rm{s},i}}\ {{{\dot{y}}}_{s,i}}\ {{{\dot{z}}}_{s,i}} \right]$。第i颗卫星指向用户的径向单位向量a_{su, i}=[a_{x, i} a_{y, i} a_{z, i}]，i=1, 2, …, N，用户相对第i颗卫星的径向速度(即卫星至用户的速度矢量在径向投影)为v_{r, i}=(v_u-v_{s, i})a_{su, i}，径向速度估计值${{\mathit{\boldsymbol{\hat v}}}_{{\rm{r,}}i}} = $(${{{\mathit{\boldsymbol{\hat{v}}}}}_{\rm{u}}}$-${{{\mathit{\boldsymbol{\hat{v}}}}}_{\rm{s},i}}$)a_{su, i}。钟漂为${{{\dot{t}}}_{\rm{u}}}$，估计值为${{{\hat{\dot{t}}}}_{\rm{u}}}$。

标称载波频率为f_c，接收机接收到的第i颗卫星信号的载波频率为f_{r, i}，本地估计载波频率为f_{l, i}，先不考虑噪声，则鉴频器输出为

$ {\rm{F}}{{\rm{D}}_i} = {f_{{\rm{r}},i}} - {f_{{\rm{l}},i}} $

(1)

考虑用户与卫星相对运动，载波频率为

$ {{f'}_{{\rm{r}},i}} = \left( {1 - \frac{{{v_{{\rm{r}},i}}}}{{\rm{c}}}} \right){f_{\rm{c}}} $

(2)

但是用户时钟频率与系统时钟频率存在偏移，时钟漂移${{{\dot{t}}}_{\rm{u}}}$、f_{r, i}以及f′_{r, i}的关系^[3]为

$ {{f'}_{{\rm{r}},i}} = \left( {1 + {{\hat t}_{\rm{u}}}} \right){f_{{\rm{r}},i}} $

(3)

式中：${{{\hat{t}}}_{\rm{u}}}$的单位为s/s，为计算方便，将光速c与${{{\dot{t}}}_{\rm{u}}}$相乘，记c${{{\hat{t}}}_{\rm{u}}}$为${{{\dot{t}}}_{\rm{u}}}$，其单位变为m/s，将${{{\hat{t}}}_{\rm{u}}}={{{\dot{t}}}_{\rm{u}}}/c$代入式(3)变为

$ {{f'}_{{\rm{r}},i}} = {f_{{\rm{r}},i}} + \frac{{{f_{{\rm{r}},i}}}}{c} \cdot {{\dot t}_{\rm{u}}} \approx {f_{{\rm{r}},i}} + \frac{{{f_{\rm{c}}}}}{c} \cdot {{\dot t}_{\rm{u}}} $

(4)

式中：f_c与f_{r, i}相差很小，可近似认为相等^[3]。

由式(2)、(4)整理得

$ {f_{{\rm{r}},i}} = {{f'}_{{\rm{r}},i}} - \frac{{{f_{\rm{c}}}}}{c} \cdot {{\dot t}_{\rm{u}}} = \left( {1 - \frac{{{v_{{\rm{r}},i}} + {{\dot t}_{\rm{u}}}}}{c}} \right){f_{\rm{c}}} = {f_{\rm{c}}} + {f_{\rm{d}}} $

(5)

根据式(5)，若计算出径向速度${{{\mathit{\boldsymbol{\hat{v}}}}}_{\rm{r,}i}}$与钟漂${{{\hat{\dot{t}}}}_{\rm{u}}}$，本地载波频率f_{l, i}估计值为

$ {f_{{\rm{l}},i}} = \left( {1 - \frac{{{{\mathit{\boldsymbol{\hat v}}}_{{\rm{r}},i}} + {{\hat {\dot t}}_{\rm{u}}}}}{c}} \right){f_{\rm{c}}} = {f_{\rm{c}}} + {{\hat f}_{\rm{d}}} $

(6)

用f_{l, i}更新本地载波发生器，从而由式(1)、(5)、(6)可推出当前时刻第i通道频率鉴别器输出为

$ {\rm{F}}{{\rm{D}}_i} = \left[ {\left( {{{\mathit{\boldsymbol{\hat v}}}_{{\rm{r}},i}} - {\mathit{\boldsymbol{v}}_{{\rm{r}},i}}} \right) + \left( {{{\mathit{\boldsymbol{\hat {\dot t}}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\dot t}}}_{\rm{u}}}} \right)} \right]\frac{{{f_{\rm{c}}}}}{c} $

(7)

式中：${{{\mathit{\boldsymbol{\hat{v}}}}}_{\rm{r,}i}}$、${{{\hat{\dot{t}}}}_{\rm{u}}}$和f_{l, i}为上一时刻估计值，进一步整理得到

$ \begin{array}{*{20}{c}} {{\rm{F}}{{\rm{D}}_i} \cdot \frac{c}{{{f_{\rm{c}}}}} = \left[ {\left( {{{\mathit{\boldsymbol{\hat v}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\hat v}}}_{{\rm{s,}}i}}} \right){\mathit{\boldsymbol{a}}_{{\rm{su}},i}} - \left( {{\mathit{\boldsymbol{v}}_{\rm{u}}} - {\mathit{\boldsymbol{v}}_{{\rm{s,}}i}}} \right){\mathit{\boldsymbol{a}}_{{\rm{su}},i}}} \right] + }\\ {\left( {{{\mathit{\boldsymbol{\hat {\dot t}}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\dot t}}}_{\rm{u}}}} \right) = \left( {{{\mathit{\boldsymbol{\hat {\dot x}}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\dot x}}}_{\rm{u}}}} \right){a_{x,i}} + \left( {{{\mathit{\boldsymbol{\hat {\dot y}}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\dot y}}}_{\rm{u}}}} \right){\mathit{\boldsymbol{a}}_{y,i}} + }\\ {\left( {{{\mathit{\boldsymbol{\hat {\dot z}}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\dot z}}}_{\rm{u}}}} \right){a_{z,i}} + \left( {{{\mathit{\boldsymbol{\hat {\dot t}}}}_{\rm{u}}} - {{\mathit{\boldsymbol{\dot t}}}_{\rm{u}}}} \right)} \end{array} $

(8)

式(8)右边为用户径向分量速度和钟漂的估计误差之和，称为伪距率误差，用δ${{{\tilde{\rho }}}_{i}}$表示。令$\delta {{{\dot{x}}}_{\rm{u}}}={{{\hat{\dot{x}}}}_{\rm{u}}}-{{{\dot{x}}}_{\rm{u}}},\delta {{{\dot{y}}}_{\rm{u}}}={{{\hat{\dot{y}}}}_{\rm{u}}}-{{{\dot{y}}}_{\rm{u}}},\delta {{{\dot{z}}}_{\rm{u}}}={{{\hat{\dot{z}}}}_{\rm{u}}}-{{{\dot{z}}}_{\rm{u}}},\delta {{{\dot{t}}}_{\rm{u}}}={{{\hat{\dot{t}}}}_{\rm{u}}}-{{{\dot{t}}}_{\rm{u}}}$，并且考虑噪声影响，则式(8)简写为

$ \delta {{\tilde {\dot \rho} }_i} = \delta {{\dot x}_{\rm{u}}}{a_{x,i}} + \delta {{\dot y}_{\rm{u}}}{a_{y,i}} + \delta {{\dot z}_{\rm{u}}}{a_{z,i}} + \delta {{\dot t}_{\rm{u}}} + {v_{{{\dot \rho }_i}}} $

(9)

式中: ${{\mathit{\boldsymbol{v}}}_{{{{\dot{\rho }}}_{i}}}}$为第i颗卫星所在跟踪通道的噪声项，满足${{\mathit{\boldsymbol{v}}}_{{{{\dot{\rho }}}_{i}}}}$~N(0, σ²_vi)，E(${{\mathit{\boldsymbol{v}}}_{{{{\dot{\rho }}}_{i}}}}{{\mathit{\boldsymbol{v}}}_{{{{\dot{\rho }}}_{j}}}}$)=σ²_viδ(i-j)。

1.1 扩展卡尔曼滤波模型

取状态变量$\mathit{\boldsymbol{\dot{X}=}}\left[ {{{\dot{x}}}_{\rm{u}}}\ {{{\dot{y}}}_{\rm{u}}}\ {{{\dot{z}}}_{\rm{u}}}\ {{{\dot{t}}}_{\rm{u}}} \right]$，状态误差为$\mathit{\boldsymbol{\dot{X}}}=\left[ \delta {{{\dot{x}}}_{\rm{u}}}\ \delta {{{\dot{y}}}_{\rm{u}}}\ \delta {{{\dot{z}}}_{\rm{u}}}\ \delta {{{\dot{t}}}_{\rm{u}}} \right]$，则观测方程为

$ \left[ {\begin{array}{*{20}{c}} {\delta {{\tilde {\dot \rho} }_1}}\\ {\delta {{\tilde {\dot \rho} }_2}}\\ \vdots \\ {\delta {{\tilde {\dot \rho} }_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{a_{x,1}}}&{{a_{y,1}}}&{{a_{z,1}}}&1\\ {{a_{x,2}}}&{{a_{y,2}}}&{{a_{z,2}}}&1\\ \vdots&\vdots&\vdots&\vdots \\ {{a_{x,N}}}&{{a_{y,N}}}&{{a_{z,N}}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\delta {{\dot x}_{\rm{u}}}}\\ {\delta {{\dot y}_{\rm{u}}}}\\ {\delta {{\dot z}_{\rm{u}}}}\\ {\delta {{\dot t}_{\rm{u}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{v}}_{{{\dot \rho }_1}}}}\\ {{\mathit{\boldsymbol{v}}_{{{\dot \rho }_2}}}}\\ \vdots \\ {{\mathit{\boldsymbol{v}}_{{{\dot \rho }_N}}}} \end{array}} \right] $

(10)

其中观测矩阵为

$ {\mathit{\boldsymbol{H}}_{N \times 4}} = \left[ {\begin{array}{*{20}{c}} {{a_{x,1}}}&{{a_{y,1}}}&{{a_{z,1}}}&1\\ {{a_{x,2}}}&{{a_{y,2}}}&{{a_{z,2}}}&1\\ \vdots&\vdots&\vdots&\vdots \\ {{a_{x,N}}}&{{a_{y,N}}}&{{a_{z,N}}}&1 \end{array}} \right] $

(11)

令Δ${\mathit{\boldsymbol{\tilde {\dot \rho} }}}$为伪距率误差观测向量，真值为Δ${\mathit{\boldsymbol{\dot \rho }}}$，V为噪声项，则观测方程简写为

$ \mathit{\Delta }\mathit{\boldsymbol{\tilde {\dot \rho} }} = {\mathit{\boldsymbol{H}}_{N \times 4}} \cdot \mathit{\Delta }\mathit{\boldsymbol{\dot X}} + \mathit{\boldsymbol{V}} $

(12)

观测噪声协方差矩阵为

$ {\mathit{\boldsymbol{R}}_v} = {\rm{E}}\left( {\mathit{\boldsymbol{V}} \cdot {\mathit{\boldsymbol{V}}^{\rm{T}}}} \right) = \left[ {\begin{array}{*{20}{c}} {\sigma _{{v_1}}^2}&0&0& \cdots &0\\ 0&{\sigma _{{v_2}}^2}&0& \cdots &0\\ \vdots&\vdots&\vdots&\vdots&\vdots \\ 0&0&0& \cdots &{\sigma _{{v_N}}^2} \end{array}} \right] $

(13)

状态方程表示为

$ \left[ {\begin{array}{*{20}{c}} {\delta {{\dot x}_{{\rm{u,}}k + 1}}}\\ {\delta {{\dot y}_{{\rm{u,}}k + 1}}}\\ {\delta {{\dot z}_{{\rm{u,}}k + 1}}}\\ {\delta {{\dot t}_{{\rm{u,}}k + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\delta {{\dot x}_{{\rm{u,}}k}}}\\ {\delta {{\dot y}_{{\rm{u,}}k}}}\\ {\delta {{\dot z}_{{\rm{u,}}k}}}\\ {\delta {{\dot t}_{{\rm{u,}}k}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{w_x}}\\ {{w_y}}\\ {{w_z}}\\ {{w_{{{\dot t}_{\rm{u}}}}}} \end{array}} \right] $

(14)

其中状态转移矩阵为

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right] $

(15)

过程噪声协方差矩阵为

$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {\sigma _{\delta {{\dot x}_{\rm{u}}}}^2}&0&0&0\\ 0&{\sigma _{\delta {{\dot y}_{\rm{u}}}}^2}&0&0\\ 0&0&{\sigma _{\delta {{\dot z}_{\rm{u}}}}^2}&0\\ 0&0&0&{\sigma _{\delta {{\dot t}_{\rm{u}}}}^2} \end{array}} \right] $

(16)

1.2 卡尔曼滤波具体过程

当接收机稳定跟踪并且得到PVT解算之后，就可以初始化和启动VFLL。

1) VFLL状态的每一次更新，需要根据位置解算计算出方向余弦矩阵H_N×4；

2) 对系统状态进行时间更新$\mathit{\boldsymbol{\dot X}}_{k + 1}^ - = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k}\mathit{\boldsymbol{\dot X}}_k^ + $；

3) 对误差协方差矩阵进行时间更新P_k+1^-=Φ_kP_k⁺Φ_k^T+Q_k；

4) 由于用上一时刻状态误差ΔX_k⁺的估计值对状态变量进行修正，状态误差ΔX_k⁺时间更新在步骤2中已经完成，所以步骤8中需要对状态误差置0，其时间更新结果ΔX_k+1^-=Φ_kΔ$\mathit{\boldsymbol{\dot{X}}}_{k}^{+}$=0；

5) 计算观测量残差，第i通道鉴频器输出的伪距率误差观测量为δ${{{\tilde{\dot{\rho }}}}_{i,k+1}}$，其预测值为δ${{{\hat{\dot{\rho }}}}_{i,k+1}}$，则观测量残差为

$ {z_{i,k + 1}} = \delta {{\tilde {\dot \rho} }_{i,k + 1}} - \delta {{\hat {\dot \rho} }_{i,k + 1}} $

(17)

由式(1)、(5)、(6)、(8)，观测量可以表示为δ${{{\tilde{\dot{\rho }}}}_{i,k+1}}$=(f_d-${{{\tilde{f}}}_{\rm{d}}}$)·c/f_c，从而预测量表示为δ${{{\hat{\dot{\rho }}}}_{i,k+1}}=\left( {{{\hat{f}}}_{\rm{d}}}-{{{\tilde{f}}}_{\rm{d}}} \right)$·c/f_c，其中f_d为频率偏移真值，${{{\tilde{f}}}_{\rm{d}}}$为上一时刻估计的频率偏移，可由本地载波发生器直接得到，${{{\tilde{f}}}_{\rm{d}}}$为当前时刻频率偏移预估计值，利用步骤2)结果通过式(6)计算得到，所有观测量残差用向量z_k+1表示；

6) 卡尔曼增益矩阵为

$ {\mathit{\boldsymbol{K}}_{k + 1}} = \mathit{\boldsymbol{P}}_{k + 1}^ - \mathit{\boldsymbol{H}}_{k + 1}^{\rm{T}}{\left[ {{\mathit{\boldsymbol{H}}_{k + 1}}\mathit{\boldsymbol{P}}_{k + 1}^ - \mathit{\boldsymbol{H}}_{k + 1}^{\rm{T}} + {\mathit{\boldsymbol{R}}_{k + 1}}} \right]^{ - 1}} $

(18)

7) 状态误差观测量更新结果$\Delta \mathit{\boldsymbol{\dot{X}}}_{k+1}^{+}=\Delta \mathit{\boldsymbol{\dot{X}}}_{k+1}^{-}+{{\mathit{\boldsymbol{K}}}_{k+1}}{{\mathit{\boldsymbol{Z}}}_{k+1}}={{\mathit{\boldsymbol{K}}}_{k+1}}{{\mathit{\boldsymbol{Z}}}_{k+1}}$；

8) 更新系统状态$\mathit{\boldsymbol{\dot{X}}}_{k+1}^{+}=\Delta \mathit{\boldsymbol{\dot{X}}}_{k+1}^{-}+\Delta \mathit{\boldsymbol{\dot{X}}}_{k+1}^{+}$，并对$\Delta \mathit{\boldsymbol{\dot{X}}}_{k+1}^{+}$置0；

9) 更新误差协方差矩阵P_k+1⁺=(I-K_k+1H_k+1)P_k+1^-；

10)利用式(2)或式(5)对各个通道载波频率进行估计。

2 VFLL与FLL伪距率误差方差比较

弱信号跟踪中，通常使用FLL代替相位锁定环(phase lock loop, PLL)进行频率跟踪，假设有N个通道同时跟踪，鉴频器输出N个伪距率误差观测量。位置的最大似然估计需要将所有的观测信号同时在一个矢量环路中进行处理^[10]，同理，对所有的伪距率误差观测量同时进行处理才有可能获得最优的速度估计。当单颗卫星信号强度不足以进行频率跟踪时，来自所有卫星的信号可能能够进行速度估计，通过速度估计值反过来对频率进行跟踪。引言部分已经指出接收机中一些潜在相关性，并且当同时跟踪多于4颗卫星信号时，矢量跟踪需要估计的变量较标量跟踪少^[11]，估计方程为超定方程，都可能使跟踪性能提高，但是超定方程并不一定意味着性能更好，观测量与接收机状态的关联方式至关重要。

文献[11]通过加权最小二乘的方法给出VDLL跟踪性能优于DLL的证明，但是在证明过程中假设所有通道的噪声独立同分布为高斯噪声，并且未给出W_N×N对角线元素小于1的证明。本文基于文献[11]，利用加权最小二乘方法，在更一般情况下证明VFLL跟踪性能优于FLL。

假设提供N个伪距率误差观测量，数据模型为

$ \left[ {\begin{array}{*{20}{c}} {\delta {{\tilde {\dot \rho} }_1}}\\ {\delta {{\tilde {\dot \rho} }_2}}\\ \vdots \\ {\delta {{\tilde {\dot \rho} }_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0&0& \cdots &0\\ 0&1&0& \cdots &0\\ \vdots&\ddots&\ddots&\ddots&\vdots \\ 0&0&0& \cdots &1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\delta {{\dot \rho }_1}}\\ {\delta {{\dot \rho }_2}}\\ \vdots \\ {\delta {{\dot \rho }_N}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{v_{{{\dot \rho }_1}}}}\\ {{v_{{{\dot \rho }_2}}}}\\ \vdots \\ {{v_{{{\dot \rho }_N}}}} \end{array}} \right] $

(19)

简记为

$ \mathit{\Delta }\mathit{\boldsymbol{\tilde {\dot \rho} }} = {\mathit{\boldsymbol{I}}_{N \times N}} \cdot \mathit{\Delta }\mathit{\boldsymbol{\dot \rho }} + \mathit{\boldsymbol{V}} $

(20)

其加权最小二乘估计为

$ \mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho }}} = {\left( {\mathit{\boldsymbol{I}}_{N \times N}^{\rm{T}}\mathit{\boldsymbol{R}}_v^{ - 1}{\mathit{\boldsymbol{I}}_{N \times N}}} \right)^{ - 1}} \cdot \mathit{\boldsymbol{I}}_{N \times N}^{\rm{T}} \cdot \mathit{\boldsymbol{R}}_v^{ - 1} \cdot \mathit{\Delta }\mathit{\boldsymbol{\tilde {\dot \rho} }} = \mathit{\Delta }\mathit{\boldsymbol{\tilde {\dot \rho} }} $

(21)

伪距率误差估计的期望E{Δ${\mathit{\boldsymbol{\hat {\dot \rho} }}}$}=Δ${\mathit{\boldsymbol{\dot \rho }}}$，为无偏估计，其协方差矩阵为

$ \begin{array}{*{20}{c}} {{\mathop{\rm var}} \left\{ {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho }}}} \right\} = {\mathop{\rm cov}} \left( {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho } }},\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho } }}} \right) = }\\ {{\rm{E}}\left\{ {\left( {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho } }} - {\rm{E}}\left\{ {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho } }}} \right\}} \right) \cdot {{\left( {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho } }} - {\rm{E}}\left\{ {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot \rho } }}} \right\}} \right)}^{\rm{T}}}} \right\} = {\mathit{\boldsymbol{R}}_v}} \end{array} $

(22)

由式(12)可得到状态误差最小二乘估计为

$ \mathit{\Delta }\hat {\dot X} = {\left( {\mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}}\mathit{\boldsymbol{R}}_v^{ - 1}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)^{ - 1}} \cdot \mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}} \cdot \mathit{\boldsymbol{R}}_v^{ - 1} \cdot \mathit{\Delta }\mathit{\boldsymbol{\tilde {\dot \rho} }} $

(23)

状态误差估计的期望E{Δ${\mathit{\boldsymbol{\hat{\dot{X}}}}}$}=Δ ${\mathit{\boldsymbol{\dot{X}}}}$，为无偏估计，其协方差矩阵为

$ {\mathop{\rm var}} \left( {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot X}}}} \right) = {\mathop{\rm cov}} \left( {\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot X}}},\mathit{\Delta }\mathit{\boldsymbol{\hat {\dot X}}}} \right) = {\left( {\mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}}\mathit{\boldsymbol{R}}_v^{ - 1}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)^{ - 1}} $

(24)

为了能够和FLL跟踪进行比较，用估计的状态误差对伪距率误差进行估计为

$ \mathit{\Delta }\mathit{\boldsymbol{\tilde {\dot \rho} = }}{\mathit{\boldsymbol{H}}_{N \times 4}} \cdot \mathit{\Delta }\mathit{\boldsymbol{\hat {\dot X}}} $

(25)

其期望E{Δ${\hat{\dot{\rho }}}$}= H_N×4·Δ ${\mathit{\boldsymbol{\dot{X}}}}$=Δ${\mathit{\boldsymbol{\dot \rho }}}$为无偏估计，协方差矩阵为

$ \begin{array}{*{20}{c}} {{\mathop{\rm var}} \left\{ {\mathit{\Delta }\hat {\dot \rho} } \right\} = {\mathop{\rm cov}} \left( {\mathit{\Delta }\hat {\dot \rho} ,\mathit{\Delta }\hat {\dot \rho} } \right) = }\\ {{\mathit{\boldsymbol{H}}_{N \times 4}}{{\left( {\mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}}\mathit{\boldsymbol{R}}_v^{ - 1}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)}^{ - 1}}\mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}} = {\mathit{\boldsymbol{W}}_{N \times N}}} \end{array} $

(26)

以上所有推导假设(H_N×4^TR_v^-1H_N×4)非奇异，对比式(22)、(26)，当N=4时，VFLL与FLL伪距率误差估计的方差相同，下面考虑当N>4的情况。

A、B为任意两个n×n实矩阵，当A-B为正定矩阵时记为A>B，当A-B为非负定时记为A≥B。由Schwarz不等式^[16]得到：

$ {\mathit{\boldsymbol{Q}}^{\rm{T}}}\mathit{\boldsymbol{Q}} \ge {\left( {{\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{Q}}} \right)^{\rm{T}}}{\left( {{\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{P}}} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{Q}}} \right) $

(27)

式中：(P^TP)为非奇异，当且仅当存在某个矩阵S使Q=PS时，取等号。

R_v为对角矩阵且非奇异，将R_v作如下分解：

$ {\mathit{\boldsymbol{R}}_v} = {\left( {{{\mathit{\boldsymbol{R'}}}_v}} \right)^{\rm{T}}}{{\mathit{\boldsymbol{R'}}}_v} $

(28)

根据式(27) (28)推导如下：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}_v} = {{\left( {{{\mathit{\boldsymbol{R'}}}_v}} \right)}^{\rm{T}}}{{\mathit{\boldsymbol{R'}}}_v} \ge {{\left( {{{\left( {{{\left( {{{\mathit{\boldsymbol{R'}}}_v}} \right)}^{ - 1}}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)}^{\rm{T}}}{{\mathit{\boldsymbol{R'}}}_v}} \right)}^{\rm{T}}} \cdot }\\ {{{\left( {{{\left( {{{\left( {{{\mathit{\boldsymbol{R'}}}_v}} \right)}^{ - 1}}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)}^{\rm{T}}}\left( {{{\left( {{{\mathit{\boldsymbol{R'}}}_v}} \right)}^{ - 1}}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)} \right)}^{ - 1}} \cdot }\\ {\left( {{{\left( {{{\left( {{{\mathit{\boldsymbol{R'}}}_v}} \right)}^{ - 1}}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)}^{\rm{T}}}{{\mathit{\boldsymbol{R'}}}_v}} \right) = {\mathit{\boldsymbol{H}}_{N \times 4}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}}\mathit{\boldsymbol{R}}_v^{ - 1}{\mathit{\boldsymbol{H}}_{N \times 4}}} \right)}^{ - 1}}\mathit{\boldsymbol{H}}_{N \times 4}^{\rm{T}} = {\mathit{\boldsymbol{W}}_{N \times N}}} \end{array} $

(29)

当N>4时，R′_v秩为N，((R′_v)^-1H_N×4)秩为4，不可能取到等号，此时R_v-W_N×N为正定矩阵，其对角线上元素都大于0，从而证明W_N×N对角元素小于R_v相应的对角线元素。

通过对比VFLL和FLL对伪距率误差估计的协方差分析，可得当可见卫星超过4颗时，矢量跟踪方法估计的伪距率误差方差变小，性能优于标量跟踪方法。同时，由W_N×N的表达式可以得到矢量跟踪的性能与可见卫星数目及卫星的分布情况有关^[11]。

3 基于扩展卡尔曼滤波器的VFLL实验结果与分析

实验之前，需要考虑下面几个方面：

1) 考虑低动态情况下的应用，卡尔曼滤波更新速度不需要很长，取0.5 s更新周期；

2) 为了得到高的跟踪灵敏度，针对GPS信号，使用20 ms数据进行鉴频；

3) 每个跟踪通道频率鉴别时刻不同，但在低动态情况下，可以忽略这种差异，使用距离卡尔曼滤波更新时刻最近的鉴别时刻对应的鉴频结果；

4) 为了进一步提高载波跟踪灵敏度，可以使用更新时刻之前的几个鉴频结果进行平均操作，再作为观测量输入VFLL。

分别采用仿真程序和思博伦信号发生器产生GPS数据，选择静止场景、9颗卫星，信号强度逐渐递减，如表 1所示，并且为了对比，仿真程序产生一个载噪比(C/N₀)恒定为40 dB/Hz的信号，其他参数设置与前面提到的用仿真程序产生的数据一致。

码跟踪环，作20 ms相干累加，4次非相干，等效噪声带宽0.25 Hz，位置解算使用最小二乘方法，将位置信息输入VFLL模块。

采用VFLL对仿真程序产生的两个GPS数据分别进行处理，图 1为某颗卫星载波跟踪结果，并且给出理论C/N₀，随着信号强度降低，估计频率在理论值附近波动，当信号强度低至13 dB/Hz时，VFLL无法对信号进行跟踪。图 2为Spirent STR4500信号发生器产生GPS数据中某颗卫星信号跟踪结果，并给出相应C/N₀估计结果，从图中可以得到在540~660 s时C/N₀估计值在14 dB/Hz附近波动，输出频率偏移没有出现大的波动，即信号低至14 dB/Hz时VFLL依然可以进行载波跟踪。文献[2]在没有比特先验信息情况下能够跟踪15.5 dB/Hz卫星信号，而本文给出的VFLL能够跟踪至14 dB/Hz。

4 结论

1) 通过对比VFLL和FLL对伪距率误差估计的协方差分析，可得当可见卫星超过4颗时，矢量跟踪方法估计的伪距率误差方差变小，性能优于标量跟踪方法。

2) 采用仿真程序和Spirent信号发生器产生GPS数据，选择静止场景、9颗卫星，信号低至14 dB/Hz时基于扩展卡尔曼滤波器的VFLL依然可以进行载波跟踪，与文献[2]相比跟踪灵敏度提高了1.5 dB左右。