1. INTRODUCTION
GNSS has been widely applied to the aeronautics, astronautics, land and ocean fields (Li et al., 2010). Autonomous navigation, enabled by GNSS, has been considered to be an enabling technology for the purpose of increasing reliability and accuracy and reducing operating costs. Due to the good geometry coverage, GNSS has been successfully applied to low-orbit satellites. However, the receiver meets a challenging environment in geostationary or high Earth orbit satellites (Lorga et al., Reference Lorga, Silva, Dovis, Cintio, Kowaltschek, Jimenez and Jansson2011), named high-orbit spacecraft, which play a crucial role in meteorological, television transmission and wide area augmentation systems (Axelrad et al., Reference Axelrad, Bradley, Tombasco, Mohiuddin and Donna2010).
For effective GNSS-based navigation for high-orbit spacecraft, the specifically designed on-board receiver has to overcome several major challenges. Due to occultation of the Earth, the poor geometric coverage of the GNSS signal is always present, causing about 65% of the main lobe signal to be blocked (Moreau et al., Reference Moreau, Axelrad, James, Garrison and Long2000). Moreover, the large transmission distance brings a remarkable attenuation in the signal level, even lower than −185 dBW at the apogee (Capuano et al., Reference Capuano, Botteron and Farine2013). At the same time, the signal continuity of the high sensitivity receiver turns to exacerbation with the motion and rotation of the vehicle (William et al., Reference William, Bo and Michael2006).
Generally, the receiver adopts a scalar tracking loop (STL). It exploits a multi-channel independent parallel structure to estimate the code phase and Doppler. Although the STL performs well under a high C/N0 environment, the tracking errors will increase significantly and even lose lock in case of low C/N0 (Lashley et al., Reference Lashley, Bevly and Hung2009). By contrast, the VTL, which is based on the common centre phase position of the receiver antenna, relies on a centralised filter to process the information from all channels (Kanwal et al., Reference Kanwal, Hurskainen and Nurmi2010). It makes full use of the redundant information and achieves the mutual assistance of each channel (Hsu et al., Reference Hsu, Jan, Groves and Kubo2015). However, the errors in the contaminated channel will diffuse along with the propagation of useful information. This results in the decrease of tracking accuracy of normal channels, even having a serious effect on the reliability of the navigation system.
In order to improve the adaptive ability of VTL for time-varying noise, the most common method is to allocate an impact factor for each channel based on the measurement information so as to adjust the gain matrix and noise covariance matrix of the filter (Kim et al., Reference Kim, Jee and Im2011). Furthermore, the impact factors can be determined more accurately by the noise models of the vector delay lock loop (VDLL) and the vector frequency lock loop (VFLL) (Dardin et al., Reference Dardin, Calmettes, Priot and Toumeret2013), which can be referred to together as VDFLL. However, it tends to diverge quickly if the outputs of VDLL and VFLL become abnormal. Therefore, the Bayesian filtering algorithm (Schaffrin and Kwon, Reference Schaffrin and Kwon2002; Garcia and Hausotte, Reference Garcia and Hausotte2013) and the neural network algorithm (Jwo et al., Reference Jwo, Wen and Lee2015) have been introduced to enhance the reliability of VTL. In the traditional Bayesian algorithm, the noise parameters can be estimated by a large quantity of sample data instead of the latest measurements. However, the accurate integral results of the probability density function in the Bayesian algorithm cannot be obtained easily. The neural network algorithm will face the problem of computational complexity and slow convergence, and can easily fall into a local minimum.
To enhance the VTL reliability in the case of signal outage, plenty of solutions have been provided. According to the C/N0, the VTL can detect the states of each channel and isolate the outage channels automatically (Zou et al., Reference Zou, Lian, Wu, Xu and Zhang2017). However, due to the weak strength of high-orbit GNSS signal, the estimation accuracy of C/N0 is limited. Another solution is to employ artificial intelligent modules, such as radial basis function neural networks (Kurban and Beždok, Reference Kurban and Beždok2009) and fuzzy neural networks (Li et al., Reference Li, Wang, Li, Gao and Tan2010), to enhance the traditional methods. Meanwhile, the support vector machine (SVM) approach has also been explored to avoid local minimisation and over-fitting problems in a neural network (Frangos et al., Reference Frangos, Kealy, Gikas and Hasnur2010). It achieves the interpolation for outage data by setting up a sliding window function (Bhatt et al., Reference Bhatt, Aggarwal, Devabhaktuni and Bhattacharya2014). However, the accuracy will decrease with the increase of signal unavailable time.
These facts represent the main motivation of the present research, that is, an adaptive vector tracking scheme applied to the high-orbit GNSS receiver is proposed. The second section develops the crucial problems of traditional VTL faced in the high-orbit environment. The following section provides a mathematical description of the VTL filter and a deep analysis of the measurement noise. In the fourth section, the implementation steps of adaptive algorithm are described in detail as well as the proposed VTL scheme. Finally, the simulation results and conclusions are presented.
2. PROBLEM DESCRIPTION
The main target of a tracking loop is to estimate the code phase and Doppler of a received signal accurately. The specific structure of traditional VTL is shown in Figure 1. The red lines represent the paths of error propagation and diffusion among channels.
Specific structure of traditional VTL.
As shown in Figure 1, the local carrier and code are correlated respectively with the received signal. Based on the DLL and FLL, the measurements of centralised Kalman filter (KF) can be provided. They are utilised to estimate the states (Δ P, Δ V, Δ t). With the assistance of ephemeris, the receiver's PVT and the tracking parameters of numeric control oscillators can be further calculated.
However, the weak GNSS signal can further increase the probability of signal unavailability. If it happens in a certain channel, the outputs of other discriminators, i.e., the code phase error observations, would also tend to diverge because of the information diffusion among channels. Besides, the measurement noise of KF mostly depends on the received signal and discriminator function. Since the signal strength is always changing, the measurement noise holds a time-varying characteristic, causing decrease in the estimated accuracy of states (Chen et al., Reference Chen, Wang and Xu2014). In fact, both of them limit the generalisation ability of the VTL in the high-orbit receiver.
3. MODELLING AND ANALYSIS FOR VTL FILTER
The centralised filter is an indispensable structure to realise the closed-loop feedback of VTL. In this section, the centralised unscented Kalman filter (UKF) (Gustafsson and Hendeby, Reference Gustafsson and Hendeby2012) is employed to deal with the nonlinearity of discriminator in the weak signal environment. Further, the factors that affect the measurement noise have been analysed in theory.
3.1. Centralised UKF model
UKF models include a linear system equation and a nonlinear measurement equation. The former can be modelled as
$$\begin{align}\mathop{\textit{\textbf{X}}}^{\bullet} ={\textit{\textbf{FX}}}+{\textit{\textbf{GW}}} \end{align}$$
where 
${\textbf{X}}=[ \delta \rho _{1} , \; \; \mathop{\delta \rho _{1} }\limits^{\bullet} , \; \; \; \cdots \, , \; \; \delta \rho _{n} \; , \; \; \mathop {\delta \rho _{n} }\limits^{\bullet}] $ is the state vector, 
${\textit{\textbf{F}}}$ is the system matrix, 
${\textit{\textbf{G}}}$ is the process noise driven matrix, and 
${\textit{{\textbf{W}}}}$ is the process noise. They can be given by
$$\begin{align} {{\textbf{F}}}&=\left[\begin{array}{@{}cccc@{}} {{{\textbf{F}}}_1 } & {0} & \cdots & {0} \\ {0} & {{\textit{\textbf{F}}}_2 } & \cdots & {0} \\ \cdots & \cdots & \cdots & \cdots \\ {0} & {0} & \cdots & {{{\textbf{F}}}_n } \end{array} \right], \quad {{\textbf G}}=\left[\begin{array}{@{}ccccc@{}} 0 & 0 & 0 & 1 & 0 \\ {e_x^{(1)} } & {e_y^{(1)} } & {e_z^{(1)} } & 0& 1 \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & 0 & 1 & 0 \\ {e_x^{(n)} } & {e_y^{(n)} } & {e_z^{(n)} } & 0& 1 \\ \end{array} \right]\label{eq2} \end{align}$$
$$\begin{align} {\textit{\textbf{W}}}&=\left[\begin{array}{@{}ccccc@{}} {w_x } & {w_y } & {w_z } & {w_b } & {w_d } \end{array}\right]\label{eq3} \end{align}$$
where 
$\delta \rho _{i} $ and 
$\mathop {\delta \rho _{i} }\limits^{\bullet} $ represent the residuals contained in the pseudorange and pseudorange rate calculated by using the current position and velocity information, respectively; 
${\textit{\textbf{F}}}_{i} =[ {0, \; {1}; \; 0, \; 0} ] $; 
$\begin{smallmatrix} {[ e_{x}^{( i) } } & {e_{y}^{( i) } } & {e_{z}^{( i) }] } \end{smallmatrix}$ is the unit line of sight (LOS) vector of the navigation satellite relative to the high-orbit spacecraft, which is expressed in the Earth-centred Earth-fixed (ECEF) frame; 
$[ w_{x} , \; w_{y} , \; w_{z} ] $ represents the process noise caused by the dynamics of the high-orbit spacecraft; w b and w d are the process noise caused by clock bias and clock drift, respectively.
Regarding the measurement model, the outputs of the discriminator in DLLs and FLLs are utilised as the measurements. The measurement equation can then be mathematically written as
$$\begin{align}{\textit{\textbf{Z}}}={\textit{{\textbf{h}}}}( {\textit{{\textbf{X}}}} ) +{\textit{{\textbf{V}}}} \end{align}$$
where 
${{\textbf{Z}}}=[ \delta \tau _{1} , \; \; \delta f_{1} , \; \; \; \cdots \, , \; \; \delta \tau _{n} \; , \; \; \delta f_{n} ] $ is measurement vector; 
$\delta \tau _{i} $ and δ f i are the residuals of the code phase and carrier Doppler shift output by the discriminator in the DLL and FLL, respectively; 
${\textit{\textbf{h}}} ( \bullet ) $ is the nonlinear measurement function; and 
${\textit{{\textbf{V}}}}$ is the measurement noise.
3.2. Factors affecting measurement noise
With the aim of numerically evaluating the impact factors of measurements accuracy, the statistical characteristics of code phase error and Doppler frequency error are derived and analysed. The correlation results of the in-phase branch and quadrature branch at epoch k can be expressed as
$$\begin{align} \left\{{\begin{array}{@{}l} I_k =\sqrt {2C \mathord{\left/ {\vphantom {C {N_0 \cdot }}} \right. \kern-\nulldelimiterspace} {N_0 \cdot }T} D\Re \left( {\Delta \tau _k } \right)\sin c\left( {\pi \Delta f_k T} \right)\cos \left( {\Delta \varphi _k +\pi \Delta f_k T} \right)+\eta _I \\ Q_k =\sqrt {2C \mathord{\left/ {\vphantom {C {N_0 \cdot }}} \right. \kern-\nulldelimiterspace} {N_0 \cdot }T} D\Re \left( {\Delta \tau _k } \right)\sin c\left( {\pi \Delta f_k T} \right)\sin \left( {\Delta \varphi _k +\pi \Delta f_k T} \right)+\eta _Q \\ \end{array}} \right. \end{align}$$
where T is the coherent integration time; 
$\Re ( \bullet ) $ is the code correlation result; 
$\Delta \tau _{k} $ is the code phase error; D is the navigation data; Δ f k is Doppler frequency error; 
$\Delta \varphi _{k} $ is the carrier phase error at the beginning of coherent integration; and ηI and ηQ are both the standard Gaussian noise and independent of each other.
3.2.1. Modelling for code phase error
The code phase error is produced by the non-coherent early minus late amplitude discriminator, for which the mathematical model can be written as
$$\begin{align}\Delta \tau _k =\frac{1}{2}\frac{\sqrt {I_{k,E}^2 +Q_{k,E}^2 } -\sqrt {I_{k,L}^2 +Q_{k,L}^2 } }{\sqrt {I_{k,E}^2 +Q_{k,E}^2 } +\sqrt {I_{k,L}^2 +Q_{k,L}^2 } } \end{align}$$
where the subscripts E and L represent the early code and late code, respectively. It can be seen from Equation (5) that, since the noise ηI and ηQ are both the standard Gaussian noise and independent of each other, 
$I_{k}^{2} +Q_{k}^{2} $ obeys the non-central χ2 distribution, and 
$\sqrt {I_{k}^{2} +Q_{k}^{2} } $ obeys the non-central χ distribution accordingly. The mean and variance are given as
$$\begin{align} \left\{ {\begin{array}{@{}l} \mu =\sqrt {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-\nulldelimiterspace} 2} \cdot L_n^\alpha \left( {{-\lambda ^2} \mathord{\left/ {\vphantom {{-\lambda ^2} 2}} \right. \kern-\nulldelimiterspace} 2} \right) \\[5pt] \sigma ^2=k+\lambda ^2-\mu ^2 \\ \end{array}} \right. \end{align}$$
where 
$L_{n}^{\alpha} ( x ) $ is n-order generalised Laguerre polynomial; α is the polynomial parameter; λ is the non-central parameter; and k is the degree of freedom.
In order to simplify the model, the scale constant 0.5 in Equation (6) is omitted to avoid tediousness, and the numerator and denominator in Equation (6) are recorded as M and N, respectively. Therefore, the two-order statistics of 
$M/N$ are as follows:
$$\begin{align} \left\{ {\begin{array}{@{}l} \mu _{M \mathord{\left/ {\vphantom {M N}} \right. \kern-\nulldelimiterspace} N} =\dfrac{\mu _M }{\mu _N }-\dfrac{{\rm cov}\left( {M,N} \right)}{\mu _N^2 }+\dfrac{\mu _M }{\mu _N^3 }\sigma _N^2 \\[15pt] \sigma _{M \mathord{\left/ {\vphantom {M N}} \right. \kern-\nulldelimiterspace} N}^2 =\left( {\dfrac{\mu _M }{\mu _N }} \right)^2\left( {\dfrac{\sigma _M^2 }{\mu _M^2 }+\dfrac{\sigma _N^2 }{\mu _N^2 }-\frac{2{\rm cov}\left( {M,N} \right)}{\mu _M \mu _N }} \right) \\ \end{array}} \right. \end{align}$$
where μM and μN are the means of M and N, respectively; 
$\sigma _{M}^{2} $ and 
$\sigma _{N}^{2} $ are the variances of M and N, respectively; and cov
$( {M, \; N} ) $ is the covariance matrix of M and N.
By combining Equations (6) and (8), the second-order statistics of 
$\Delta \tau _{k} $ are as follows:
$$\begin{align} \left\{ {\begin{array}{@{}l} \mu _{\Delta \tau } =\dfrac{{1}}{{2}}\mu _{M \mathord{\left/{\vphantom {M N}} \right. \kern-\nulldelimiterspace} N} \\[10pt] \sigma _{\Delta \tau }^2 =\dfrac{{1}}{{4}}\sigma _{M \mathord{\left/ {\vphantom {M N}} \right. \kern-\nulldelimiterspace} N}^2 \\ \end{array}} \right. \end{align}$$ According to Equations (5) and (9), the relationship between the measurement noise of VDLL and C/N0, the coherent integration time, the code phase error can be obtained. According to the simulation result, the code phase error is set to 0.1 chip, and the C/N0 increases uniformly from 11 dB-Hz to 45 dB-Hz in 40 s. The T are 1 ms, 2 ms, 5 ms and 10 ms, respectively. The varying curves of 
$\sigma _{\Delta \tau }^{2} $ are shown in Figure 2. These results indicate that the noise variance of VDLL will decrease monotonically with the increase of C/N0. If C/N0 is lower than 15 dB-Hz, the variance tends to be stable. Meanwhile, when C/N0 is constant, the noise variance would reduce gradually with the increase of T.
Varying curves of VDLL noise variance.
3.2.2. Modelling for Doppler frequency error
Compared with the phase locked loop, the FLL has a wider noise bandwidth and better dynamic tolerance performance (Shen and Zhang, Reference Shen, Zhang, Wong and Ma2013), according with the requirements of a weak and dynamic signal environment. The four-quadrant cross product discriminator is employed to estimate the carrier frequency error, for which the mathematical model can be written as
$$\begin{align} \Delta f_k =\frac{{\rm sign}( {P_{\rm {dot}} } ) \cdot P_{\rm {cross}} }{2\pi \cdot T\cdot ( {I_{k,P}^2 +Q_{k,P}^2 } ) } \end{align}$$
where 
${\rm sign}( \bullet ) $ is the symbolic function; and P dot and P cross are the coherent energy of the dot product and cross product respectively. The specific expressions of P dot and P cross are respectively,
$$\begin{align} \left\{ {\begin{array}{@{}l} P_{\rm {dot}} =I_{k-1,P} I_{k,P} +Q_{k-1,P} Q_{k,P} \\ P_{\rm {cross}} =I_{k-1,P} Q_{k,P} -Q_{k-1,P} I_{k,P} \\ \end{array}} \right. \end{align}$$where the subscript P represents the prompt code.
In order to simplify the carrier frequency error model, note,
$$\begin{align} \left\{ {{\begin{array}{@{}l} {\alpha =\Delta \varphi_{k-1} +\pi \Delta f_{k-1} T} \\ {\beta =\Delta \varphi_k +\pi \Delta f_k T\quad } \\ \end{array} }} \right. \end{align}$$
$$\begin{align} A=\sqrt {2C \mathord{\left/ {\vphantom {C {N_0 \cdot }}} \right. \kern-\nulldelimiterspace} {N_0 \cdot }T} D\Re \left( {\Delta \tau _k } \right)\sin c\left( {\pi \Delta f_k T} \right) \end{align}$$where α and β represent carrier phase errors at epoch k − 1 and epoch k, respectively.
When the VTL keeps the stable tracking of the signal, combining the output results of the VDLL and the VFLL, the spread-spectrum code (i.e., C/A code) and the navigation message information in the signal can be demodulated. In this case, A can be expressed as
$$\begin{align} A=\sqrt {2C \mathord{\left/ {\vphantom {C {N_0 \cdot }}} \right. \kern-\nulldelimiterspace} {N_0 \cdot }T} \end{align}$$ The scale factor term 
${\rm {sign}( {P_{\rm\ {dot}}\ }\ ) } / {2\pi T}$ in Equation (10) is omitted to avoid tediousness, and the numerator P cross and the denominator 
$I_{P, k}^{2} +Q_{P, k}^{2} $ in Equation (10) are recorded as M′ and N′, respectively. The specific expressions of M′ and N′ are respectively,
$$\begin{align} \left\{ {\begin{array}{@{}l} {M}'=\left( {A\cos \alpha +\eta _{I_P ,k-1} } \right)\left( {A\sin \beta +\eta _{Q_P k} } \right)-\left( {A\cos \alpha +\eta _{Q_P ,k-1} } \right)\left( {A\sin \beta +\eta _{I_P ,k} } \right) \\[5pt] {N}'=\left( {A\cos \beta +\eta _{I_P ,k} } \right)^2+\left({A{sin}\beta +\eta _{Q_P ,k} } \right)^2 \\ \end{array}} \right.\end{align}$$ Therefore, the two-order statistics of 
${{M}'} / {{N}'}$ are as follows,
$$\begin{align} \left\{ {{\begin{array}{@{}c} {\mu _{{{M}'} \mathord{\left/ {\vphantom {{{M}'} {{N}'}}} \right. \kern-\nulldelimiterspace} {{N}'}} =\dfrac{\mu _{{M}'} }{\mu _{{N}'}}-\dfrac{{\rm cov}\left( {{M}',{N}'} \right)}{\mu _{{N}'}^2 }+\dfrac{\mu _{{M}'}}{\mu _{{N}'}^3 }\sigma _{{N}'}^2 \quad \quad \quad } \\[15pt] {\sigma _{{{M}'} \mathord{\left/ {\vphantom {{{M}'} {{N}'}}} \right.\kern-\nulldelimiterspace} {{N}'}}^2 =\left( {\dfrac{\mu _{{M}'} }{\mu_{{N}'} }} \right)^2\left( {\dfrac{\sigma _{{M}'}^2 }{\mu _{{M}'}^2}+\dfrac{\sigma _{{N}'}^2 }{\mu _{{N}'}^2 }-\dfrac{2{\rm cov}\left( {{M}',{N}'}\right)}{\mu _{{M}'} \mu _{{N}'} }} \right)} \\ \end{array} }} \right.\end{align}$$
where μM′ and μN′ are the means of M′ and N′, respectively; 
$\sigma _{{M}'}^{2} $ and 
$\sigma _{{N}'}^{2} $ are the variances of M′ and N′, respectively; and cov
$( {{M}', \; {N}'} ) $ is the covariance matrix of M′ and N′.
By combining Equations (10) and (16), the second-order statistics of Δ f k are as follows,
$$\begin{align} \left\{ {\begin{array}{@{}l} \mu _{\Delta \tau } =\dfrac{\rm {sign}\left( {P_{\rm {dot}} } \right)}{2\pi T}\mu _{{{M}'} \mathord{\left/ {\vphantom {{{M}'} {{N}'}}} \right. \kern-\nulldelimiterspace} {{N}'}} \\[10pt] \sigma _{\Delta \tau }^2 =\dfrac{{1}}{\left( {2\pi T} \right)^{2}}\sigma _{{{M}'} \mathord{\left/ {\vphantom {{{M}'} {{N}'}}} \right. \kern-\nulldelimiterspace} {{N}'}}^2 \end{array}} \right. \end{align}$$ Similarly, the varying curves of 
$\sigma _{\Delta f}^{2} $ can be gained and shown in Figure 3. Unlike VDLL, the variance of VFLL holds a rapid descent at low C/N0 then turns to steady increase. From the vertical direction, the smaller C/N0 is, the more obviously the noise variance is affected by T.
Varying curves of VFLL noise variance.
The statistical models in Equation (9) and (17) indicate that measurement noise is mainly affected by C/N0 and tracking algorithm, which includes the coherent integration time and the discriminator results at the previous moment. Therefore, the adaptive filter algorithm is designed subsequently based on the above factors.
4. DESIGN OF AN ADAPTIVE VTL SCHEME
This section will focus on a non-Gaussian distributed measurement noise and the interdependence between different channels. To improve the tracking performance, an adaptive VTL scheme is proposed and described in detail.
4.1. Probability models
The continuous system model represented by Equations (1) and (4) is discretised to obtain,
$$\begin{align} \left\{ {\begin{array}{@{}l} {\textit{\textbf{X}}}_k ={\boldsymbol{\varPhi }}_{k,k-1} {\textit{\textbf{X}}}_{k-1} +{\boldsymbol{\varGamma}}_{k-1} {\textit{\textbf{W}}}_{k{-1}} \\[5pt] {\textit{\textbf{Z}}}_k ={\textit{\textbf{h}}}( {{\textit{\textbf{X}}}_k } ) +{\textit{\textbf{V}}}_k \end{array}} \right. \end{align}$$
where 
${\boldsymbol{\varPhi }}_{k, k-1} $ is the one-step transfer matrix and 
${\boldsymbol{\varGamma}}_{k-1} $ is the discrete process noise driven matrix.
It is assumed that both 
${\textit{\textbf{W}}}_{k{-1}} $ and 
${\textit{\textbf{V}}}_{k} $ obey Gaussian distribution, and the mean of 
${\textit{\textbf{W}}}_{k{-1}} $ is zero. According to Equation (18), the probability models based on the generalised Bayesian framework can be written as,
$$\begin{align} \left\{ {\begin{array}{@{}l} p( { {{\textit{\textbf{X}}}_k } \vert {\textit{\textbf{X}}}_{k-1} } ) =N( {{\boldsymbol{\varPhi }}_{k,k-1} {\textit{\textbf{X}}}_{k-1} , {\boldsymbol{\varGamma}}_{k-1} {\textit{\textbf{Q}}}_{k-1} {\boldsymbol{\varGamma}}_{k-1}^T } ) \\[5pt] p( { {{\textit{\textbf{Z}}}_k } \vert {\textit{\textbf{X}}}_k } ) =N( {\textit{\textbf{h}}}( {{\textit{\textbf{X}}}_k } ) +{\boldsymbol{\varepsilon }}_k ,{\textit{\textbf{{R}}}_k } ) \end{array}} \right. \end{align}$$
where 
${\textit{\textbf{Q}}}_{k} $ is the process noise covariance, and 
${\boldsymbol{\varepsilon }}_{k} $ and 
${\textit{\textbf{R}}}_{k} $ are the mean and variance of the measurement noise 
${\textit{\textbf{V}}}_{k} $, respectively.
Though the noises of the discriminators in each channel are Gaussian distributions, the mean and variance can be time-varying with the changing of the actual received signal strength according to Equations (5), (6) and (10). Due to the propagation of information among channels, the measurement noise would lose its independence and the non-diagonal elements of the covariance matrix are not equal to zero. Considering the uncertainty of measurement noise, the mean and variance of measurement noise are regarded as random variables. In view of the time-varying characteristics of the measurement noise in the VTL and the relativity between the channels, we assume that the measurement noise obeys the multiple Gaussian distribution. Because the multivariate Gaussian distribution and the inverse Wishart (IW) distribution are usually used as the conjugate prior distribution in Bayesian estimation to describe the uncertainty of the mean and variance of the multivariate Gaussian distribution (Bishop, Reference Bishop2006). Therefore, we assume that the prior distribution of the mean and variance of the measurement noise is the multivariate Gaussian distribution and IW distribution, respectively. Thus, the measurement noise models can be further written as,
$$\begin{align} \left\{ {{\begin{array}{@{}c} {p( { {{\boldsymbol{\varepsilon }}_k } \vert {\textit{\textbf{R}}}_k } ) =N( { {{\boldsymbol{\varepsilon }}_k } \vert {\boldsymbol{\eta }}_k , \beta_k {\textit{\textbf{R}}}_k } ) } \\[5pt] {p( {{\textit{\textbf{R}}}_k } ) ={IW}( { {{\textit{\textbf{R}}}_k } \vert \vartheta_k ,{\boldsymbol{\varLambda }}_k } ) \quad \;} \end{array} }} \right. \end{align}$$
where 
${\boldsymbol{\eta }}_{k} $ and 
$\beta _{k} {\textit{\textbf{R}}}_{k} $ are the mean and covariance of Gauss distribution; βk is the scale factor; ϑk is the degree of freedom for IW distribution; and 
${\boldsymbol{\varLambda }}_{k} $ is the symmetric positive definite variance matrix. IW distribution possesses the ability to make a full noise covariance matrix estimation which is suitable for handling the correlated noise among different channels of VTL.
Assuming the dynamic models of the states 
$p( { {{\textit{\textbf{X}}}_{k} } \vert {\textit{\textbf{X}}}_{k-1} } ) $ and the models of measurement noise variance 
$p( { {{\textit{\textbf{R}}}_{k} } \vert {\textit{\textbf{R}}}_{k-1} } ) $ are independent, then (Särkkä and Nummenmaa, Reference Särkkä and Nummenmaa2009)
$$\begin{align} p( { {{\textit{\textbf{X}}}_k ,{\textit{\textbf{R}}}_k } \vert {\textit{\textbf{ X}}}_{k-1} ,{\textit{\textbf{R}}}_{k-1} } ) =p( { {{\textit{\textbf{X}}}_k } \vert {\textit{\textbf{X}}}_{k-1} } ) \cdot p( { {{\textit{\textbf{R}}}_k } \vert {\textit{\textbf{R}}}_{k-1} } ) \end{align}$$
where 
$p( { {{\textit{\textbf{X}}}_{k} , \; {\textit{\textbf{R}}}_{k} } \vert {\textit{\textbf{X}}}_{k-1} , \; {\textit{\textbf{R}}}_{k-1} } ) $ is the joint distribution of 
${\textit{\textbf{X}}}_{k} $ and 
${\textit{\textbf{R}}}_{k} $.
Before we obtain the measurement at epoch k, the predictive distribution of 
${\textit{\textbf{X}}}_{k} $ and 
${\textit{\textbf{R}}}_{k} $ is given by the Chapman-Kolmogorov equation (Särkkä, Reference Särkkä2013),
$$\begin{align} p\left( {\left. {{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right|{{\textbf{Z}}}_{k-1} } \right) & =\int {p\left( {\left. {{{\textbf{X}}}_k } \right|{{\textbf{X}}}_{k-1} } \right)\cdot p\left( {\left. {{{\textbf{R}}}_k } \right|{{\textbf{R}}}_{k-1} } \right)\cdot } \nonumber\\ & p\left( {\left. {{{\textbf{X}}}_{k-1} ,{{\textbf{R}}}_{k-1} } \right|{{\textbf{Z}}}_{k-1} } \right)d{{\textbf{X}}}_{k-1} d{{\textbf{R}}}_{k-1} \end{align}$$
where 
$p( {\textit{\textbf{X}}}_{k}, \; {\textit{\textbf{R}}}_{k}\vert {{\textbf{Z}}}_{k-1}) $ is the joint posterior distribution of 
${\textit{\textbf{X}}}_{k} $ and 
${\textit{\textbf{R}}}_{k} $ at epoch k − 1.
After obtaining the measurement at epoch k, according to Bayes' theorem, the joint posterior distribution of 
${\textit{\textbf{X}}}_{k} $ and 
${\textit{\textbf{R}}}_{k} $ can be expressed as (Bishop, Reference Bishop2006)
$$\begin{align} p\left( {\left. {{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right|{{\textbf{Z}}}_k } \right) & =\frac{p\left( {\left. {{{\textbf{Z}}}_k } \right|{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right)\cdot p\left( {\left. {{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right|{{\textbf{Z}}}_{k-1} } \right)}{p\left( {\left. {{{\textbf{Z}}}_k } \right|{{\textbf{Z}}}_{k-1} } \right)}\nonumber \\ & =\frac{p\left( {\left. {{{\textbf{Z}}}_k } \right|{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right)\cdot p\left( {\left. {{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right|{{\textbf{Z}}}_{k-1} } \right)}{\int {p\left( {\left. {{{\textbf{Z}}}_k } \right|{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right)\cdot p\left({\left. {{{\textbf{X}}}_k ,{{\textbf{R}}}_k } \right|{{\textbf{Z}}}_{k-1} }\right)d{{\textbf{X}}}_{k-1} d{{\textbf{R}}}_{k-1} } }\end{align}$$
where 
$p( { {{\textit{\textbf{Z}}}_{k} } \vert {\textit{\textbf{X}}}_{k} , \; {\textit{\textbf{R}}}_{k} } ) $ is the likelihood function and 
$p( { {{\textit{\textbf{Z}}}_{k} } \vert {\textit{\textbf{Z}}}_{k-1} } ) $ is the marginal distribution. Due to the complexity of the integral of the denominator in Equation (23), it is difficult to obtain the analytical solution directly, which makes the process of calculating the joint posterior distribution 
$p( { {{\textit{\textbf{X}}}_{k} , \; {\textit{\textbf{R}}}_{k} } \vert {\textit{\textbf{Z}}}_{k} } ) $ according to Equation (23) more complex.
Instead of directly integrating the multi-dimensional probability density function, the VB inference method uses the newly introduced distribution to approach the true posterior distribution through iteration (Beal, Reference Beal2003]{bib1a}; Mbalawata et al., Reference Mbalawata, Vihola and Haario2015). We record the estimated parameters set as 
${\boldsymbol{\varTheta}}=( {{\textit{\textbf{X}}}_{k} , \; {\textit{\textbf{R}}}_{k} } ) $. According to the theory of VB, the logarithmic form of the marginal distribution is expressed as (Bishop, Reference Bishop2006),
$$\begin{align} \ln p\left( {\left. {{{\textbf{Z}}}_k } \right|{{\textbf{Z}}}_{k-1} } \right) & =\int {g\left( {\boldsymbol{\varTheta}} \right)\cdot \ln \left[ {\frac{p\left({{{\textbf{Z}}}_k ,{\boldsymbol{\varTheta}}} \right)}{p\left( {\left. {\boldsymbol{\varTheta}} \right|{{\textbf{Z}}}_k } \right)}} \right]{d}{\boldsymbol{\varTheta}}}\nonumber \\ & =\int {g\left( {\boldsymbol{\varTheta}} \right)\cdot \ln \left[ {\frac{p\left({{{\textbf{Z}}}_k ,{\boldsymbol{\varTheta}}} \right)}{g\left( {\boldsymbol{\varTheta}}\right)}} \right]{d}{\boldsymbol{\varTheta}}} +\int {g\left( {\boldsymbol{\varTheta}} \right)\cdot \ln \left[ {\frac{g\left( {\boldsymbol{\varTheta}}\right)}{p\left( {\left. {\boldsymbol{\varTheta}} \right|{{\textbf{Z}}}_k }\right)}} \right]{d}{\boldsymbol{\varTheta}}}\nonumber \\ & =T\left[ {g\left( {\boldsymbol{\varTheta}} \right)} \right]+D_{KL} \left[{\left. {g\left( {\boldsymbol{\varTheta}} \right)} \right\|p\left( {\left. {\boldsymbol{\varTheta}} \right|{{\textbf{Z}}}_k } \right)} \right]\end{align}$$
$$\begin{align}T\left[ {g\left( {\boldsymbol{\varTheta}} \right)} \right]=\int {g\left( {\boldsymbol{\varTheta}} \right)\cdot \ln \left[ {\frac{p\left( {{\textit{\textbf{Z}}}_k ,{\boldsymbol{\varTheta}}} \right)}{g\left( {\boldsymbol{\varTheta}} \right)}}\right]{d}{\boldsymbol{\varTheta}}}\end{align}$$
$$\begin{align}D_{KL} \left[ {\left. {g\left( {\boldsymbol{\varTheta}} \right)} \right\|p\left({\left. {\boldsymbol{\varTheta}} \right|{\textit{\textbf{Z}}}_k } \right)} \right]=\int{g\left( {\boldsymbol{\varTheta}} \right)\cdot \ln \left[ {\frac{g\left( {\boldsymbol{\varTheta}} \right)}{p\left( {\left. {\boldsymbol{\varTheta}} \right|{\textit{\textbf{Z}}}_k } \right)}} \right]{d}{\boldsymbol{\varTheta}}}\end{align}$$
where 
$p( {{\textit{\textbf{Z}}}_{k} , \; {\boldsymbol{\varTheta}}} ) $ is the joint distribution of 
${\textit{\textbf{Z}}}_{k} $ and 
${\boldsymbol{\varTheta}}$; 
$g( {\boldsymbol{\varTheta}} ) $ is an introduced approximate distribution for the true posterior distribution 
$p( { {\boldsymbol{\varTheta}} \vert {\textit{\textbf{Z}}}_{k} } ) $; 
$T[ {g( {\boldsymbol{\varTheta}} ) } ] $ is the free energy; and 
$D_{KL} [ { {g( {\boldsymbol{\varTheta}} ) } \Vert p( { {\boldsymbol{\varTheta}} \vert {\textit{\textbf{Z}}}_{k} } ) } ] $ is the Kullback-Leibler (KL) divergence between the approximate distribution 
$g( {\boldsymbol{\varTheta}} ) $ and the true posterior distribution 
$p( { {\boldsymbol{\varTheta}} \vert {\textit{\textbf{Z}}}_{k} } ) $. KL divergence is nonnegative and can evaluate the dissymmetry between two different distributions. If the approximate distribution 
$g( {\boldsymbol{\varTheta}} ) $ is similar to the true distribution
$p( { {\boldsymbol{\varTheta}} \vert Z} ) $, then 
$D_{KL} [ { {g( {\boldsymbol{\varTheta}} ) } \Vert p( { {\boldsymbol{\varTheta}} \vert {\textit{\textbf{Z}}}_{k} } ) } ] $ will be minimised.
Since the marginal distribution 
$p( { {{\textit{\textbf{Z}}}_{k} } \vert {\textit{\textbf{Z}}}_{k-1} } ) $ is independent of the introduced approximate distribution 
$g( {\boldsymbol{\varTheta}} ) $, the variation of the marginal distribution 
$p( { {{\textit{\textbf{Z}}}_{k} } \vert {\textit{\textbf{Z}}}_{k-1} } ) $ with respect to the approximate distribution 
$g( {\boldsymbol{\varTheta}} ) $ is zero. Therefore, when we adjust the parameters set 
${\boldsymbol{\varTheta}}$ to make 
$g( {\boldsymbol{\varTheta}} ) $ approximate the true posterior distribution 
$p( { {\boldsymbol{\varTheta}} \vert {\textit{\textbf{Z}}}_{k} } ) $, the marginal distribution is always constant. Moreover, minimising KL divergence is equivalent to maximising free energy 
$T[ {g( {\boldsymbol{\varTheta}} ) } ] $ (Bishop, Reference Bishop2006).
It is assumed that the number of parameters in the parameter set 
${\boldsymbol{\varTheta}}$ is d, and these parameters are independent of each other. The approximate distribution 
$g( {\boldsymbol{\varTheta}} ) $ can be decomposed as follows (Beal, Reference Beal2003),
$$\begin{align} g( {\boldsymbol{\varTheta}} ) =\mathop \prod \limits_{i=1}^d g( {\varTheta_i } ) \end{align}$$ Where Θ i is the ith parameter in the parameter set and 
$g( {\varTheta_{i} } ) $ is the approximate distribution of the corresponding parameter.
In this case, the free energy 
$T[ {g( {\boldsymbol{\varTheta}} ) } ] $ can be expressed as,
$$\begin{align} T\left[ {g\left( {\boldsymbol{\Theta}} \right)} \right] & =\int {g\left( {\boldsymbol{\Theta}} \right)\cdot \ln \left[ {\frac{p\left( {{\textit{\textbf{Z}}}_k ,{\boldsymbol{\Theta}}} \right)}{g\left( {\boldsymbol{\Theta}} \right)}}\right]{d}{\boldsymbol{\Theta}}}\nonumber \\ & =\int {\prod\limits_{i=1}^d {g\left( {\varTheta _i } \right)} \cdot \left\{{\ln \left[ {p\left( {{\textit{\textbf{Z}}}_k ,{\boldsymbol{\Theta}}} \right)}\right]-\ln \left[ {\prod\limits_{i=1}^d {g\left( {\varTheta _i } \right)} }\right]} \right\}{d}\varTheta _1 } {d}\varTheta _{2} \cdots{d}\varTheta _d \end{align}$$ Find the maximum value of the free energy 
$T[ {g( {\boldsymbol{\varTheta}} ) } ] $ with respect to the approximate distribution of each parameter, then the general solution of the approximate distribution of each parameter can be expressed as (Sato, Reference Sato2001; Beal,Reference Beal2003; Särkkä and Nummenmaa, Reference Särkkä and Nummenmaa2009).
$$\begin{align} g\left( {\varTheta _i } \right)=\frac{{exp}\left\{ {E_{g\left( {\varTheta_{j\not= i} } \right)} \left\langle {\ln \left[ {p\left( {{\textit{\textbf{Z}}}_k,{\boldsymbol{\varTheta}}} \right)} \right]} \right\rangle } \right\}}{\int{{exp}\left\{ {E_{g\left( {\varTheta _{j\not= i} } \right)} \left\langle{\ln \left[ {p\left( {{\textit{\textbf{Z}}}_k ,{\boldsymbol{\varTheta}}} \right)}\right]} \right\rangle } \right\}} {d}\varTheta _i }\end{align}$$
where 
$E_{g( {\varTheta_{j\not= i} } ) } \langle{\ln [ {p( {{\textit{\textbf{Z}}}_{k} , \; {\boldsymbol{\varTheta}}} ) } ] } \rangle $ is the expectation of joint distribution for other d–1 approximate distributions. It can be found that the approximate distribution of each parameter depends on the approximate distribution of other parameters. Since the free energy 
$T[ {g( {\boldsymbol{\varTheta}} ) } ] $ is a convex function of each approximate distribution, any Θ i can be initialised and the other Θ j can be calculated iteratively. The final parameter set 
${\boldsymbol{\varTheta}}$ will converge and make the approximate distribution close to a posteriori distribution. At the moment, the parameters in set 
${\boldsymbol{\varTheta}}$ are considered as the optimal estimation. They are used to update the prior distribution of measurements noise and match for the actual noise.
4.2. Adaptive VTL filter algorithm implementation
Based on the probability model of UKF, the adaptive algorithm is employed to identify the time-varying noise parameters 
$( {{\boldsymbol{\eta }}, \; \beta , \; \vartheta , \; {\boldsymbol{\varLambda }}} ) $ in Equation (20). The adaptive algorithm is a combination of the VB algorithm and C/N0 online estimation. The flow diagram of the adaptive filter algorithm is shown in Figure 4.
Flow diagram of adaptive filter algorithm.
The concrete implementation steps of the proposed adaptive algorithm are as follows:
Step 1 Initialisation
The purposes of initialisation mainly include two aspects. The first one is to give the initial value of recursive parameters in the VTL filter, which can be used to determine the prior distribution of the noise parameters in Equation (20). The other crucial one is when the tracking channel switches the satellite, it is necessary to carry out the reassessment operation on the noise level in order to adjust the prior distribution accurately.
The initial value of measurement noise can be given by Equation (9) and (17), in which the C/N0 is estimated by the power ratio method. It is low-complexity and widespread application in software receivers (Falletti et al., Reference Falletti, Pini and Lo Presti2011). Thus, the [C/N0]est can be obtained as
$$\begin{align} ( C/N_0) _{est}& =\frac{1}{T}\cdot \frac{\mu_{\rm {MP}} -1}{M-\mu _{\rm {MP}} } \nonumber\\ {\rm WBP}&=\sum_{i=1}^M {( {I_{P,i}^2 +Q_{P,i}^2 } ) } \nonumber\\ {\rm NBP}&=\left( {\sum_{i=1}^M {I_{P,i} } } \right)^2+\left({\sum_{i=1}^M {Q_{P,i} } } \right)^2\label{eq30} \end{align}$$where M is the number of integration times and μMP is the mean of ratio sequence for narrowband signal power (NBP) and wideband signal power (WBP).
Step 2 One-step prediction
We assume that the mean and variance of measurement noise are time-varying, but in general, the dynamic model of mean and variance of noise is unknown. For the consideration of time-varying of unknown noise parameters, one-step prediction of the parameters in the prior distribution can be expressed as
$$\begin{align} \left\{ {{\begin{array}{@{}c} {{\boldsymbol{\eta }}_{k/k-1} ={\boldsymbol{\eta }}_{k-1}, \; \beta _{k/k-1}={\beta _{k{-1}} } \mathord{\left/ {\vphantom {{\beta _{k{-1}} }\rho }} \right. \kern-\nulldelimiterspace} \rho } \\[5pt] {\vartheta _{k/k-1} =\rho \vartheta _{k-1}, \; {\boldsymbol{\varLambda }}_{k/k-1} =\rho {\boldsymbol{\varLambda }}_{k-1} }\end{array} }} \right. \end{align}$$ 
$\rho ( \rho \in ( 0, \; 1] ) $ is the attenuation factor introduced to reflect the fluctuation characteristics of noise parameters. The closer the attenuation factor ρ is to zero, the more obvious the instability of noise parameters is. In practical application, it can be set according to the fluctuation characteristics of noise (Beal, Reference Beal2003; Särkkä and Nummenmaa, Reference Särkkä and Nummenmaa2009).
One-step prediction equations of states are as follows:
$$\begin{align} \left\{ {\begin{array}{@{}l} {\textit{\textbf{L}}}_{k/k-1} ={\boldsymbol{\varPhi }}_{k/k-1} {\textit{\textbf{L}}}_{k-1} \\[5pt] {\textit{\textbf{P}}}_{k/k-1} ={\boldsymbol{\varPhi }}_{k/k-1} {\textit{\textbf{P}}}_{k-1} {\boldsymbol{\varPhi }}_{k/k-1}^T +{\boldsymbol{\varGamma}}_{k-1} {\textit{\textbf{Q}}}_{k-1} {\boldsymbol{\varGamma}}_{k-1}^T \end{array}} \right. \end{align}$$
where 
${\textit{\textbf{L}}}_{k-1} $ and 
${\textit{\textbf{P}}}_{k-1} $ are the mean and covariance of states, respectively.
Step 3 States estimation
Update the posterior distribution of states by the measurements,
$$\begin{align} \left\{ {{\begin{array}{@{}l} {{\textit{\textbf{L}}}_k ={\textit{\textbf{L}}}_{k/k-1} +{\textit{\textbf{K}}}_k \cdot \left[ {{\textit{\textbf{Z}}}_k -{\textit{\textbf{h}}}\left( {{\textit{\textbf{L}}}_{k/k-1} } \right)-E\left({{\boldsymbol{\varepsilon }}_k } \right)} \right]} \\[5pt] {{\textit{\textbf{P}}}_k ={\textit{\textbf{P}}}_{k/k-1} -{\textit{\textbf{K}}}_k {\textit{\textbf{H}}}_k{\textit{\textbf{P}}}_{k/k-1} \quad \quad \quad \quad \quad \quad }\end{array} }} \right. \end{align}$$
where the filter gain matrix 
${\textit{\textbf{K}}}_{k} ={\textit{\textbf{P}}}_{( {XZ} ) _{k/k-1} } \cdot {\textit{\textbf{P}}}_{( {ZZ} ) _{k/k-1} }^{-1} $.
Step 4 Update the measurement noise parameters
In analogous manner to Equation (33), the updated measurement noise parameters are
$$\begin{align} \left\{ {{\begin{array}{@{}l} {\beta_k =( {1+( {\beta_{k/k-1} } )^{-1}} )^{-1}\quad \quad \quad \quad \quad \quad \quad } \\[5pt] {{\boldsymbol{\eta }}_k ={\boldsymbol{\eta }}_{k/k-1} +\beta_k ( {{\textit{\textbf{Z}}}_k -{\textit{\textbf{h}}}( {E( {{\textit{\textbf{X}}}_k } ) } ) -{\boldsymbol{\eta }}_{k/k-1} } ) } \end{array} }} \right. \end{align}$$ According to Equation (20), 
${\boldsymbol{\varepsilon }}_{k} $ and 
${\textit{\textbf{R}}}_{k} $ are conditional independence. Therefore, the joint distribution 
$p( {{\boldsymbol{\varepsilon }}_{k} , \; {\textit{\textbf{R}}}_{k} } ) $ and the conditional distribution 
$p( { {{\boldsymbol{\varepsilon }}_{k} } \vert {\textit{\textbf{R}}}_{k} } ) $ can be used to calculate the ϑk and 
${\boldsymbol{\varLambda }}_{k} $ in IW distribution,
$$\begin{align} \left\{\begin{array}{@{}l} \vartheta _k =\vartheta _{k/k-1} +1\\{{\boldsymbol{\Lambda}}}_k ={{\boldsymbol{\Lambda}}}_{k/k-1} +\left( {1-\beta _k} \right)\left[ {{{\boldsymbol{\eta}}}_{k/k-1} -\left( {{{\textbf Z}}_k -{{\textbf{h}}}\left( {E\left( {{{\textbf X}}_k } \right)} \right)} \right)}\right] \\[5pt]\qquad\times \left[ {{{\boldsymbol{\eta}}}_{k/k-1} -\left( {{{\textbf Z}}_k -{{\textbf{h}}}\left( {E\left( {{{\textbf X}}_k } \right)} \right)} \right)}\right]^T+{\rm {\bf HP}}_k {{\textbf{H}}}^T \end{array} \right. \end{align}$$Step 5 Calculate the mathematical expectations
Based on the properties of Gauss distribution and IW distribution, the expectations of states and measurement noise parameters are as follows,
$$\begin{align} \left\{ {{\begin{array}{@{}l} {E( {{\textit{\textbf{X}}}_k } ) ={\textit{\textbf{L}}}_k \quad \quad \quad } \\[5pt] {E( {{\textit{\textbf{R}}}_k } ) =\dfrac{{\boldsymbol{\varLambda }}_k}{\vartheta_k -2n-1}} \\[15pt] {E( {{\boldsymbol{\varepsilon }}_k } ) ={\boldsymbol{\eta }}_k \quad \quad \quad } \end{array} }} \right. \end{align}$$where 2n is the dimension of measurements.
From Equations (33) to (35) in step 3 and step 4, it can be seen that the system state vector and measurement noise parameters are coupled with each other, so step 3 and step 4 need to be solved iteratively until the optimal parameter estimation is obtained after convergence. It is for these reasons that the ‘inner loop’ formed by steps 3 and 4 is necessary. The convergence of ‘inner loop’ is guaranteed by the convergence properties of VB inference, so the ‘inner loop’ can converge through iteration even if no more new measurements are input.
4.3. An innovative adaptive VTL scheme
As addressed in the last section, the adaptive VTL assisted by VB algorithm possesses the ability to estimate the time-varying noise. The architecture of the proposed innovative adaptive VTL scheme is shown in Figure 5.
Architecture of the proposed innovative adaptive VTL scheme.
The proposed scheme includes the following three main parts:
(1) The baseband signal processing. The input IF signal is mixed with the local carrier to obtain the in-phase branch and the quadrature branch signal. Afterwards, the correlation operation is implemented with the local code in the correlators. According to the correlation results, the VDFLL can identify the code phase error and Doppler frequency error.
(2) The estimation of states and measurement noise. The outputs of VDFLL are considered as the measurements of the VTL filter. According to the proposed adaptive algorithm, the time-varying noise can be estimated wholly through the recursive and update equations in the VB algorithm. In case of signal outage, it can effectively restrain the channel error divergence, and the estimated C/N0 can be utilised to re-assess the prior distribution function in order to guarantee the algorithm accuracy.
(3) Update the tracking parameters. The pseudorange and pseudorange rate would be acquired according to the GNSS orbit parameters provided by ephemeris. If there are more than four available satellites, the PVT information of the receiver can be calculated by the least square method. The tracking parameters would be updated in order to modify the numeric control oscillators.
5. SIMULATION RESULTS AND ANALYSIS
5.1. Simulation conditions
To validate the adaptive performance of the proposed VTL scheme under high-orbit environment, simulation tests were conducted subsequently, with the following conditions. The well-known satellite named AMSAT-OSCAR 40 (Moreau et al., Reference Moreau, Bauer, Carpenter, Davis, Davis and Jackson2002) is referred by the design of high-orbit spacecraft orbit. The total number of GNSS satellites is 110, including 24 GPS satellites, 35 BDS satellites, 27 Galileo satellites and 24 GLONASS satellites. The GNSS signal is generated by a software simulator with the emission power of 14.28 dBW. The free space propagation loss of GNSS signal can be calculated by the power budget equation based on the distance between the receiver and GNSS satellites (Capuano et al., Reference Capuano, Botteron and Farine2013). The noise spectrum power density N0 is −203.98 dBW/Hz.
The receiver antenna is assumed to be isotropic (7 dB). The IF and sampling frequency are 9.548 MHz and 38.192 MHz, respectively. The coherence integration time is 10 ms. The spacing between early and late code is 1 chip. The channel number is 8. The received power threshold is 15 dB-Hz. The updated frequency is 5 Hz and the simulation time is 200 s. The initial parameters of the filter are given as: 
${{\textbf{P}}} _{0} ={\rm diag}\{ {\; ( {10{m}} )^{2}, \; ( {{1m/s}} )^{2}, \; \cdots \; } \} $, 
${{\textbf{Q}}}_{0} ={\rm diag}\left\{{\; \left({\sqrt {20} {m}} \right)^{2}, \; \left({\sqrt {0.05} {m/s}} \right)^{2}, \; \cdots \; } \right\}$, 
${\textit{\textbf{R}}}_{0} $ can be gained by the middle equation of Equation (25), in which the variance matrix 
${\boldsymbol{\varLambda }}_{k} $ is equal to 
${\rm diag}\{ {\; 0.18^{2}, \; 0.66^{2}, \; \cdots } \} $ and the initial ϑ0 is 18.
5.2. GNSS signal availability
The signal simulation results are provided in this section to evaluate the availability of signal to the GNSS receiver. High Earth orbit can be divided approximately into two parts based on the flight altitude. There is more attention to the altitude higher than GNSS constellations, so the received signals in this orbit segment are collected and analysed in detail.
Visibility is the basal factor for availability in navigation solutions, and the former is implicitly satisfied when the latter is satisfied. The number of visible satellites is calculated in cases that received power thresholds equal to 15 dB-Hz, 20 dB-Hz, 25 dB-Hz and 30 dB-Hz, respectively. As shown in Figure 6, when the threshold is reduced to 15 dB-Hz, the number of visible satellites is always more than four. However, the number drops rapidly with the increase of threshold. It is clear that the number of visible satellites is lower than four in 30% of the simulation time with a threshold of 25 dB-Hz. In particular, most of the available satellites are concentrated at the beginning and end of the orbital period, which is near the perigee.
Number of visible satellites with different received thresholds.
The Doppler effect caused by the relative motion of spacecraft is widespread in the GNSS signal. Figure 7 shows the mean of Doppler and Doppler rate for all visible satellites at each time. It indicates that there is a conspicuous Doppler in the received signal, which can reach 11 kHz. By comparison, the Doppler rate is merely 3 Hz/s. For this reason, the influence of frequency changing within a transient coherent integration time can be largely ignored.
Doppler and Doppler rate of received signal.
Signal continuity is another factor that affects navigation solutions. Here, the power threshold is set to 15 dB-Hz. Taking GPS as an example, Figure 8 shows how long each satellite is continuously available. It indicates obviously that the available times of most satellites are less than 50%. And there is a high probability that the signal outage will occur during the high Earth orbit flight period, which would reduce the tracking accuracy.
Continuity of each GPS satellite.
5.3. Performance analysis
In this section,the tracking performance between the traditional and proposed VTL is compared. The traditional one utilises the Kalman filter VTL architecture assisted by C/N0 estimation. Of course, the receiver exploits the whole functionalities of STL until the first PVT is available. From then on, it turns to VTL architecture. Since the Doppler rate is low, it is assumed that the Doppler remains constant between the updates of two consecutive navigation solutions.
5.3.1. In case of signal outage
In this simulation, we focus on the tracking parameters changing under signal outage conditions. The outage duration of the first channel signal is sustained from 80 s to 90 s. The code phase error and Doppler are shown in Figures 9 and 10, respectively. It appears clear that the traditional VTL architecture can merely track the signal without outage condition. From then on, the tracking loop will lose the ability to correct errors and gradually turn to divergence. In the case of the proposed VTL architecture, it can restrain the growth of tracking error effectively. Thus, the code phase and Doppler would be estimated continuously during the outage duration by prediction. By comparsion, it also confirms that the VFLL is more vulnerable than VDLL in the condition of signal outage.
Code phase errors of the traditional and proposed algorithms.
Doppler estimation of the traditional and proposed algorithms.
5.3.2. In case of C/N0 changing
In this simulation, the signal C/N0 changing is considered and implemented by reducing the C/N0 uniformly with the rate of 0.5 dB-Hz/s from 60 s to 80 s, while the C/N0 is maintained at 30 dB-Hz for the rest of the time. Figures 11 and 12 show the estimation results of code phase error and Doppler. It is clearly illustrated that, with the decline of C/N0, the tracking error of the VTL architecture without adaptivity has risen rapidly and reached a remarkable level in the case of low C/N0. By comparison, from the adapting VTL assisted by VB algorithm, we can appreciate how the code phase error and Doppler continue to be calculated accurately during the C/N0 changing. It has benefited from the accurate estimation of the time-varying noise by VB algorithm.
Estimation of code phase errors.
Estimation of Doppler values.
5.4. Orbit determination accuracy
The orbit errors 
$\left({\Delta_{pos} =\sqrt {\Delta_{x}^{2} +\Delta_{y}^{2} +\Delta_{z}^{2} } } \right)$ of the traditional and proposed algorithms are depicted in Figure 13. It indicates that both algorithms can converge quickly, but the error of the proposed algorithm is just 29% of the traditional one, and the steadiness of orbit error has improved significantly. It is thus demonstrated that the proposed algorithm possesses the ability to make the navigation performance promisingly.
Comparison of orbit errors.
6. CONCLUSION
The high-orbit GNSS receiver faces the problems of discontinuous and time-varying strength of the received signal. To improve the tracking performance, an adaptive VTL scheme assisted by the variational Bayesian algorithm is proposed. The mean and variance of measurement noises are considered as the random variables in the hierarchical VB learning network. Based on the given prior distribution functions, they can be estimated adaptively to match for the actual noise and maintain the convergence of tracking errors.
The signal characteristics have been analysed for visibility, dynamics and continuity. Results show that the receiver can gain the signal adequately with the threshold as low as 15 dB-Hz. The high-orbit GNSS signal also meets the challenges of large Doppler and discontinuity. Performance evaluation for the proposed VTL scheme has been carried out in cases of signal outage and changing C/N0. Several simulation results are also provided and discussed in detail. As expected, the traditional algorithm is unable to overcome the above aggravated conditions effectively, especially under a low C/N0 environment. On the contrary, the matrix R can be estimated wholly in the adaptive algorithm. Thus, the relative noise between different channels would be calculated accurately. In the condition of signal outage, the proposed scheme can suppress the error growth in the corresponding channel. It also possesses the ability to estimate the time-varying noise and restrain the discriminator error. In conclusion, the proposed algorithm has significant advantages in accuracy and reliability in a high-orbit degraded GNSS signal environment.
ACKNOWLEDGEMENTS
This work was supported by the Natural Science Foundation of China (NSFC) (grant number 61673040, 61233005 and 61074157), the State Key Laboratory Foundation of Astronautical Dynamics of China (grant number 2012ADL-DW0201), the Joint Projects of NSFC-CNRS (grant number 61111130198), the Aeronautical Science Foundation of China (grant number 2015ZC51038, 20160812004 and 20170151002) and the Project of the Experimentation and Technology Research (1700050405). The work is also supported by the Open Research Fund of the State Key Laboratory of Space-Ground Integrated Information Technology under grant number 2015-SGIIT-KFJJ-DH-01.
