3.1. Inertial-aided Kalman filter-based tracking loop

A schematic diagram of one channel tracking loop aided by the inertial information is shown in Figure 1.

Figure 1. Basic structure diagram of inertial-aided Kalman filter-based tracking loop

The inertial-aided Kalman filter-based tracking loop consists of the following parts:

1. 1. The information of the inertial navigation system and the satellite ephemeris information are used to calculate the Doppler frequency parameters caused by the relative motion between the receiver and the satellites.

2. 2. The tracking loops are filtered using Kalman filters and estimate the error of control.

3. 3. The input of the local replica code NCO and local replica carrier NCO are calculated and generated.

3.2. Doppler frequency estimation

The inertial navigation system continuously outputs position and velocity information. At the same time, satellite position information can be obtained by using real-time demodulated satellite ephemeris. Based on Equation (1), the carrier Doppler frequency estimation can be obtained as:

(3)\begin{equation}{\tilde{f}_d} = \frac{{({\textbf{v}_I}\textbf{- }{\textbf{v}_s})}}{\lambda }\frac{{\boldsymbol{r}_I^{} - \boldsymbol{r}_s^{}}}{{|{\boldsymbol{r}_I^{} - \boldsymbol{r}_s^{}} |}}\end{equation}

where $\boldsymbol{r}_I^{}$ and ${\textbf{v}_I}$ are the position and velocity vector of the inertial navigation system; $\boldsymbol{r}_s^{}$ and ${\textbf{v}_s}$ are the position and velocity vector of the satellite; $\lambda$ is the carrier wavelength transmitted by the satellite.

3.3. Kalman filter

According to the characteristics of loop tracking, the state variates of Kalman filter are defined as:

(4)\begin{equation}\boldsymbol{x} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\Delta \tau }& {\Delta \phi } \end{array}}& {\Delta \omega }& {\Delta \dot{\omega }} \end{array}} \right]^T}\end{equation}

where $\Delta \tau$ is the phase difference between the local replica code and the received code, $\Delta \phi$, $\Delta \omega$ and $\Delta \dot{\omega }$ are the phase difference, angular frequency difference and the change rate of angular frequency difference between the local replica carrier and the received carrier, and T is the correlation integration time. The following relationship exists between the state vector ${\boldsymbol{x}_k}$ and ${\boldsymbol{x}_{k - 1}}$:

(5)\begin{equation}\left\{ \begin{array}{l} \Delta {\tau_k} = \Delta {\tau_{k - 1}} + \dfrac{1}{{763}}\Delta {\omega_{k - 1}}T - {u_{k - 1}}T\\ \Delta {\phi_k} = \Delta {\phi_{k - 1}} + \Delta {\omega_{k - 1}}T + \dfrac{1}{2}\Delta {{\dot{\omega }}_{k - 1}}{T^2} - {u_{k - 1}}T\\ \Delta {\omega_k} = \Delta {\omega_{k - 1}} + \Delta {{\dot{\omega }}_{k - 1}}T\\ \Delta {{\dot{\omega }}_k} = \Delta {{\dot{\omega }}_{k - 1}} \end{array} \right.\end{equation}

where $\boldsymbol{u}$ is the system input variable, that is, the control frequency of the local replica carrier and code. According to the relations above, the system state equation can be expressed as:

(6)\begin{equation}\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \Delta {\tau_k}\\ \Delta {\phi_k} \end{array}\\ {\Delta {\omega_k}}\\ {\Delta {{\dot{\omega }}_k}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1& 0& {\dfrac{1}{{763}}T}& 0\\ 0& 1& T& {\dfrac{1}{2}{T^2}}\\ 0& 0& 1& T\\ 0& 0& 0& 1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \Delta {\tau_{k - 1}}\\ \Delta {\phi_{k - 1}} \end{array}\\ {\Delta {\omega_{k - 1}}}\\ {\Delta {{\dot{\omega }}_{k - 1}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} \begin{array}{l} - T\\ - T \end{array}\\ 0\\ 0 \end{array}} \right]{\boldsymbol{u}_{k - 1}} + \left[ {\begin{array}{*{20}{c}} \begin{array}{l} 1\\ 0 \end{array}& \begin{array}{l} 0\\ 1 \end{array}& \begin{array}{l} \begin{array}{*{20}{c}} 0& 0 \end{array}\\ \begin{array}{*{20}{c}} 0& 0 \end{array} \end{array}\\ 0& 0& {\begin{array}{*{20}{c}} 1& 0 \end{array}}\\ 0& 0& {\begin{array}{*{20}{c}} 0& 1 \end{array}} \end{array}} \right]\boldsymbol{W}\end{equation}

where $\boldsymbol{W}$ stands for system noise.

The output of the phase discriminator is selected as the observation variable, and the phase discrimination result - currently represents the average of the phase difference from time to time, which can be expressed as:

(7)\begin{align}\Delta {\tau _{k,k - 1}} & = \frac{1}{T}\int_0^T {\left( {\Delta {\tau_{k - 1}} + \frac{1}{{763}}\Delta {\omega_{k - 1}}t} \right)dt} - \frac{1}{T}\int_0^T {{u_{k - 1}}tdt} \end{align}

(8)\begin{align}\Delta {\phi _{k,k - 1}} & = \frac{1}{T}\int_0^T {\left( {\Delta {\phi_{k - 1}} + \Delta {\omega_{k - 1}}t + \frac{1}{2}\Delta {{\dot{\omega }}_{k - 1}}{t^2}} \right)dt} - \frac{1}{T}\int_0^T {{u_{k - 1}}tdt} \end{align}

Therefore, the relationship between the observed variables $\Delta {\tau _{k,k - 1}}$ and $\Delta {\phi _{k,k - 1}}$ and the state variables can be obtained. By adjusting Equation (4), we can obtain:

(9)\begin{align} \Delta {\tau _{k,k - 1}} & =\dfrac{1}{T}\int_0^T {\left( {\left( {\Delta {\tau_k} - \dfrac{1}{{763}}\Delta {\omega_k}T + {u_{k - 1}}T} \right) + \dfrac{1}{{763}}({\Delta {\omega_k} - \Delta {{\dot{\omega }}_k}T} )t} \right)dt} - \dfrac{1}{T}\int_0^T {{u_{k - 1}}tdt} \notag\\ & =\dfrac{1}{T}\int_0^T {\left( {\Delta {\tau_k} + \dfrac{1}{{763}}\Delta {\omega_k}({t - T} )- \dfrac{1}{{763}}\Delta {{\dot{\omega }}_k}Tt} \right)dt} + \dfrac{1}{T}\int_0^T {{u_{k - 1}}({T - t} )dt} \notag\\ & =\Delta {\tau _k} - \dfrac{1}{{1526}}\Delta {\omega _k}T - \dfrac{1}{{1526}}\Delta {{\dot{\omega }}_k}{T^2} + \dfrac{1}{2}{u_{k - 1}}T \end{align}

(10)\begin{align} \Delta {\phi _{k,k - 1}} & =\dfrac{1}{T}\int_0^T {\left( {\left( {\Delta {\phi_k} - \Delta {\omega_k}T + \dfrac{1}{2}\Delta {{\dot{\omega }}_k}{T^2} + {u_{k - 1}}T} \right) + ({\Delta {\omega_k} - \Delta {{\dot{\omega }}_k}T} )t + \dfrac{1}{2}\Delta {{\dot{\omega }}_k}{t^2}} \right)dt}\notag\\ & \quad - \dfrac{1}{T}\int_0^T {{u_{k - 1}}tdt} \notag\\ & =\dfrac{1}{T}\int_0^T {\left( {\Delta {\phi_k} + \Delta {\omega_k}({t - T} )+ \Delta {{\dot{\omega }}_k}\left( {\dfrac{1}{2}{t^2} - Tt + \dfrac{1}{2}{T^2}} \right)} \right)dt} + \dfrac{1}{T}\int_0^T {{u_{k - 1}}({T - t} )dt} \notag\\ & =\Delta {\phi _k} - \dfrac{1}{2}\Delta {\omega _k}T + \dfrac{1}{6}\Delta {{\dot{\omega }}_k}{T^2} + \dfrac{1}{2}{u_{k - 1}}T \end{align}

The system observation equation is obtained as follows:

(11)\begin{equation}\left[ \begin{array}{l} \Delta {\tau_{k - 1,k}}\\ \Delta {\phi_{k - 1,k}} \end{array} \right] = \left[ {\begin{matrix} 1& 0& { - \dfrac{1}{{1526}}T}& { - \dfrac{1}{{1526}}{T^2}}\\ 0& 1& { - \dfrac{1}{2}T}& {\dfrac{1}{6}{T^2}} \end{matrix}} \right]\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \Delta {\tau_k}\\ \Delta {\phi_k} \end{array}\\ {\Delta {\omega_k}}\\ {\Delta {{\dot{\omega }}_k}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\dfrac{1}{2}T}\\ {\dfrac{1}{2}T} \end{array}} \right]{\boldsymbol{u}_{k - 1}} + \left[ {\begin{array}{*{20}{c}} 1& 0\\ 0& 1 \end{array}} \right]\boldsymbol{V}\end{equation}

where $\boldsymbol{V}$ represents measurement noise.

3.4. NCO generation

After the Kalman filter model is obtained, the state variables can be further used to calculate the local replica signal NCO generation to carry out feedback control, and the Kalman filter can replace the traditional loop filters.

Through the establishment of the Kalman filter model above, we can obtain:

(12)\begin{equation}\begin{array}{l} \Delta {\tau _{k + 1}} = \Delta {\tau _k} + \dfrac{1}{{763}}\Delta {\omega _k}T - {u_{k - 1}}T\\ \Delta {\phi _{k + 1}} = \Delta {\phi _k} + \Delta {\omega _k}T + \dfrac{1}{2}\Delta {{\dot{\omega }}_k}{T^2} - {u_{k - 1}}T \end{array}\end{equation}

The ultimate purpose of the tracking loop is to make the phase difference between the local replica signal and the received signal be zero. Therefore, set $\Delta {\tau _{k + 1}}$ and $\Delta {\phi _{k + 1}}$ to be zero in the equation, and the system input variable ${\boldsymbol{u}_k}$ can be obtained as follows:

(13)\begin{equation}{\boldsymbol{u}_k} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\dfrac{1}{T}}& 0& {\dfrac{1}{{763}}}& 0 \end{array}}\\ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0& {\dfrac{1}{T}} \end{array}}& 1& {\dfrac{1}{2}T} \end{array}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \Delta \tau \\ \Delta {\phi_k} \end{array}\\ {\Delta {\omega_k}}\\ {\Delta {{\dot{\omega }}_k}} \end{array}} \right]\end{equation}

The input of the code NCO and carrier NCO can be obtained as:

(14)\begin{equation}{\boldsymbol{f}_{NCO}}={\boldsymbol{u}_k} + {\tilde{f}_d}\end{equation}

So far, the tracking loop has realised the closed loop control.

4.1. Algorithm simulation

To achieve a comprehensive assessment of the proposed scheme, a simulation platform is arranged, including inertial navigation simulation system, intermediate frequency signal simulation system and satellite navigation receiver software system. The simulation platform structure is shown in Figure 2, including flight track generator, inertial IMU (inertial measurement unit) signal generator, inertial navigation calculation, intermediate frequency signal simulators and receiver signal processing module.

Figure 2. Diagram of the simulation platform

The IMU signal simulation includes the output data of the gyroscopes and the accelerometers. In the simulation process, the drift error of the gyroscope with time and the random error of the accelerometer are respectively considered.

The drift error of the gyroscope is defined as:

(15)\begin{equation}\varepsilon ={\varepsilon _b} + {\varepsilon _r}\end{equation}

where ${\varepsilon _b}$ and ${\varepsilon _r}$ represent random constant drift and first order Markov process.

(16)\begin{align}{\dot{\varepsilon }_b}& =\textrm{0}\end{align}

(17)\begin{align}{\dot{\varepsilon }_r}& = - \beta {\varepsilon _r} + {\omega _1}\end{align}

where $\beta$ is the time constant of inverse correlation and ${\omega _1}$ is white noise.

The random error of the accelerometer is a first order Markov process and is defined as:

(18)\begin{equation}\dot{\nabla }= - \mu \nabla + {\omega _2}\end{equation}

where $\mu$ is the time constant of inverse correlation and ${\omega _2}$ is white noise.

Simulation parameters are set as shown in Table 1.

Table 1. Inertial device simulation parameter settings

The inertial navigation solution includes the attitude matrix solution and the navigation solution, and the attitude matrix solution adopts the quaternion method. In order to verify and analyse the loop tracking effect of the designed Kalman filter structure, a comparative test is designed to test the filtering effect of the traditional loop filter and the inertial-aided Kalman filter. The user receiver is designed to move in a straight line with a variable acceleration of 1 g/s to the east. Two receiver channels are designed to track the same satellite at the same time. One channel adopts the traditional loop filter structure for tracking. The carrier loop filter bandwidth is set at 20 Hz. The other channel adopts the designed inertial-aided Kalman filter structure. The output of the carrier phase detector of the two channels is shown in Figure 3.

Figure 3. Phase discrimination results of traditional loop and inertial-aided Kalman filter-based tracking loop

As can be seen from Figure 3, in the first channel using the traditional loop filter, due to the constant loop parameters, the change rate of signal Doppler frequency increases with the increase of carrier acceleration, and the loop gradually loses the accurate tracking of signal frequency. The output value of the phase detector becomes larger, and the loop is close to losing the lock.

In the second channel using the inertial-aided Kalman filter structure, in about 12 s or so, with the completion of bit synchronisation and frame synchronisation signal, the pseudorange and ephemeris are obtained. With the condition of inertial information assistance, the frequency tracking error is reduced after inertial information assistance, and the phase jitter caused by thermal noise in the loop is reduced by Kalman filter. The output value of phase detector is significantly reduced, which can be equivalent to the decrease of loop bandwidth, and better tracking effect is obtained.

In addition, with the accumulation of time, the output value of the phase detector also increases with the adoption of the inertial-aided Kalman filter structure. This is because the errors of the inertial navigation system will gradually accumulate in the process of the navigation solution, resulting in an increasing error in the estimation of Doppler frequency. The position error and velocity error are shown in Figures 4 and 5.

Figure 4. Position error of the proposed algorithm

Figure 5. Velocity error of the proposed algorithm

As can be seen from these figures, the position error and velocity error gradually diverge, as does the carrier phase value in Figure 3.

4.2. Real data tests

To confirm the feasibility of the proposed tracking loop filter algorithm, real data tests on a Beidou/inertial combination navigation platform are conducted. The test platform is shown in Figure 6; it consists of a MEMS (micro-electromechanical system) inertial unit, a DSP-TMS320C6713B and a FPGA-EP4CE115F23, and an AD8347 used as the quadrature down-conversion mixer. The receiver sampling frequency is 62 MHz and IF (intermediate frequency) carrier frequency is 4⋅098 MHz, the carrier loop filter bandwidth is 20 Hz. The proposed tracking loop filter algorithm assisted by inertial information is realised in the receiver, meanwhile, the combined navigation information is used to correct the inertial navigation.

Figure 6. Beidou/inertial combination navigation test platform

The platform is held stationary for a long time. At the same time, the traditional receiver is used as a comparison. The positioning results of the two systems are shown in Figure 7. For a period of time, the positioning results of the traditional receiver fluctuate greatly, which is related to the change in the number of satellites involved in positioning calculation. The position results are statistically analysed as shown in Table 2.

Figure 7. Positioning results of the traditional algorithm and the proposed algorithm

Table 2. Mean square errors of positioning results of the traditional algorithm and the proposed algorithm

It can be seen that the positioning accuracy and stability of the proposed algorithm are better than that of the traditional algorithm, indicating that the tracking loop filter algorithm can effectively reduce the loop tracking error, improve the quality of measurement information, and further improve the positioning accuracy. Based on the static test results above, the feasibility of the proposed algorithm has been verified.