Weak and Dynamic GNSS Signal Tracking Strategies for Flight Missions in the Space Service Volume

Abstract

Weak-signal and high-dynamics are of two primary concerns of space navigation using GNSS (Global Navigation Satellite System) in the space service volume (SSV). The paper firstly defines a reference assumption third-order phase-locked loop (PLL) as the baseline of an onboard GNSS receiver, and proves the incompetence of this conventional architecture. Then an adaptive four-state Kalman filter (KF)-based algorithm is introduced to realize the optimization of loop noise bandwidth, which can adaptively regulate its filter gain according to the received signal power and line-of-sight (LOS) dynamics. To overcome the matter of losing lock in weak-signal and high-dynamic environments, an open loop tracking strategy aided by an inertial navigation system (INS) is recommended, and the traditional maximum likelihood estimation (MLE) method is modified in a non-coherent way by reconstructing the likelihood cost function. Furthermore, a typical mission with combined orbital maneuvering and non-maneuvering arcs is taken as a destination object to test the two proposed strategies. Finally, the experiment based on computer simulation identifies the effectiveness of an adaptive four-state KF-based strategy under non-maneuvering conditions and the virtue of INS-assisted methods under maneuvering conditions.

1. Introduction

The demand for reducing the burden of ground TT&C (Tracking, Telemetry and Command) stations and surveying vessels stimulates the development of precise orbit and attitude determination using GNSS. On the other hand, to take full advantage of the service capacity of GNSS, more flight missions from different space agencies equipped with onboard GNSS receiving terminals will utilize satellite navigation when operating in the space service volume (SSV) [1]. Nevertheless, the volume covered from a height of 3000 km to 36,000 km is difficult for GNSS applications. It should be noted that the altitude of 3000 km is defined as the boundary of the terrestrial service volume (TSV) and SSV, but the altitude (>20,000 km) above the GPS constellation is the most challenging environment for satellite navigation. A research satellite of the Max-Planck Institut für Extraterrestrische Physik, Equator-S, equipped with a dual antenna Motorola Viceroy receiver was launched in December 1997 into a geostationary transfer orbit (GTO) by the German Space Agency (DLR). Its in-orbit experiment proved the reception of GPS signals above the GPS orbit, at an altitude of 34,000 km, is possible [2], but the quality of the physical signals and the data contents were not good enough for spacecraft in-orbit navigation.

The main problems of navigation above the GPS constellation concentrate on two aspects: (1) insufficient signal availability and poor dilution of precision (DOP); and (2) weak signal processing with high dynamic stress and poor ranging accuracy. Both of them degrade the robustness and reliability of GNSS space service. The first problem is mainly related to signal visibility, while the second one is related to signal power. For the state of the art, the former can be improved by an interoperable SSV based on the development of multi-GNSSs interoperability [3], and we specified the GNSS SSV characterization and service performance in terms of four GNSS constellations in our previous research publication [4], but for the second problem, the technical matters are definitely more complicated, which will be the discussion topics in this work.

In normal conditions, the sensitivity of a GNSS receiver and its dynamic performance interact with each other. Fortunately, there are few terrestrial GNSS terminal manufactures that claim both high receiver sensitivity and desirable dynamic performance simultaneously. However, for SSV users, the situation is quite different; herein, the issue about how to process weak GNSS signals with high dynamic stress have to be dealt with. Comparatively, the carrier tracking loop is more likely to lose lock than the code tracking loop in weak-signal and high-dynamic environments. Thus, a third-order phase-locked loop (PLL) is taken as the receiver reference assumption model (RRAM) for space applications. According to automatic control principle [5], the third-order loop is only sensitive to jerk rather than velocity and acceleration. Although the jerk for non-maneuvering space vehicle is quite small, the jerk for orbital maneuvering mission can be as high as 4 g/s (g = 9.8 m/s² is the Earth’s acceleration of gravity) or more. It means that the performance of the SSV RRAM should be analyzed quantitatively.

In order to keep the third-order loop stable, the coherent integration time (CIT) and loop noise bandwidth (LNBW) are designed with upper limitation. By calculating the loop measurement errors caused by thermal noise and dynamic stress respectively, the maximal bearable jerk under different carrier-to-noise density ratio (C/N₀) can be determined. Obviously, the loop of SSV RRAM is likely to lose its lock in harsh environments, for instance, when the C/N₀ is lower than 30 dB-Hz and the jerk is over 4 m/s³ at the same time.

Owing to this reason, an adaptive four-state Kalman filter (KF)-based algorithm is presented with the intention of maximizing the signal processing performance of closed loop (CL) form. On the basis of estimating the state of code phase and the third-order state vector of the carrier phase, the principle about how to regulate the filter gain matrix to vary the equivalent LNBW is introduced in this work. Furthermore, the relationship between noise and adaptive bandwidth is also discussed.

Actually, the CL form uses a feedback path to make corrections of code phase and carrier frequency errors. In contrast, open loop (OL) tracking updates the local code phase and carrier frequency estimates flexibly, not entirely dependent upon the feedback path [6]. Aided by supplementary measurements, the OL scheme is free of loop instability. Thus, it is suggested that an OL tracking design could be an alternative to be adopted in SSV user receivers to overcome the dilemma of settling the issues of weak-signal and high-dynamics together.

Various OL tracking methods have been developed up to now. A US patent (Patent No: US 6633255B2) [7] proposed four kinds of carrier frequency measuring approaches, based on a frequency doubler, a frequency discriminator, a block phase estimator, and a channelized filter, respectively. Some scholars demonstrated quasi-open loop architecture to update the frequency of the local oscillator (LO) every several epochs, instead of updating frequency estimations at each epoch [8], which can be recognized as an interim scheme between closed loop and open loop. Other literatures released an idea of FFT-based frequency-domain tracking [9], which is a novel high-sensitivity tracking method developed to perform accurate frequency parameter estimation by processing the spectral peak line and adjacent lines [10]. Without doubt, an INS-assisted open loop tracking strategy is another effective method to find the true signals with high dynamics and high noise level [11]. However, all of these methods are derived on the basis of some assumption that disregarding the coupling effect between the carrier phase and Doppler frequency. Then the technical matter comes down to put forward an improved methodology to generate the carrier phase and Doppler frequency estimates decoupled with each other for the LO.

This paper incorporates a GNSS orbit propagator with INS-assisted open loop tracking together to obtain the aiding information about the Doppler frequency. With the purpose of updating the Doppler frequency estimation that is insensitive to the corresponding carrier phase estimation, a non-coherent maximum likelihood estimator is used to eliminate the coupling relationship in this paper. A new gradient function from the two-dimensional log-likelihood cost function for code delay and Doppler frequency is first established, and the new cost function is totally independent of the carrier phase tracking error. Then we can get the gradient and Hessian expressions of this new cost function with respect to the Doppler frequency. When the frequency correction is computed with the gradient divided by the Hessian, the equation does not contain the component of carrier phase and proves the cancellation of the coupling effect. This paper illustrates the implementation block diagram of the open loop tracking structure using the maximum likelihood estimation (MLE) method, and introduces its detailed operational procedure.

In the part of testing, a typical flight mission, with large variation ranges of signal power and dynamic stress, is taken as a destination object. The trajectory of the object is hybrid, which consists of a period of normal status (non-maneuvering) followed with a period of orbital maneuvering. The simulation results verify the fact that the performance of an adaptive KF-based strategy is superior to the classical CL structure, especially when the object operates in a non-maneuvering status. However, in orbital maneuvering conditions, an adaptive KF-based strategy is incompetent to perform well as what the INS-assisted strategy does. Additionally, the non-coherent MLE method in the INS-assisted strategy provides more robustness to track the Doppler frequency.

The remainder of the paper is organized as follows: Section 2 establishes a reference assumption PLL structure as the baseline of state-of-art; in Section 3, a four-state KF-based signal tracking method is presented which can adjust the equivalent LNBW adaptively; Section 4 depicts the architecture of the INS-assisted open loop design, and further proposes the non-coherent MLE algorithm applied in the architecture. Simulation experiments are carried out to test the effectiveness and feasibility of our strategies in Section 5; and finally, the paper is concluded by outlining the distinctive benefits in Section 6.

2. Reference Assumption PLL Features

Generally, GNSS applications suffer from some technical problems in harsh environments. For instance, indoor localization is trying to invent high-sensitivity receivers to receive weak signals, while missile guidance requires a better capacity to bear high-dynamics. As a result, the technologies aimed at solving high-sensitivity and high-dynamics are developed separately. However, space navigation with GNSS in the domain of SSV is influenced by the problem of realizing simultaneous high-sensitivity and high-dynamics, especially when the space vehicles are making orbital maneuvers. Theoretically, it is much more difficult to settle both of the two problems which interact with each other. To evaluate the performance of a current onboard GNSS receiver convincingly, we must establish a unified standard loop named SSV RRAM first.

2.1. The Definition of SSV RRAM

The so-called ‘dynamics’, in a generic sense, are concerned with velocity, acceleration, and jerk of the relative movement between the GNSS satellite and the user vehicle. Considering the motion complexity of SSV users, a space vehicle equipped with a third-order PLL is treated as the baseline of current technology. A third-order PLL is recommended for SSV RRAM resulting from its unique attributes. First, a third-order PLL can track the variations of velocity and acceleration without any bias, which is prior to first-order and second-order PLLs. Second, the jerk is so small for non-maneuvering space vehicles that the steady state dynamic stress error caused by jerk is also small. With regard to the parameter design of the SSV RRAM, two core parameters, the LNBW and CIT, should be emphatically analyzed.

The interaction effect of high-sensitivity and high-dynamics is reflected by these two parameters. Reducing the tracking loop bandwidth is the most effective way to achieve high-sensitivity, but a narrower LNBW is not in the best interest of high-dynamics. Secondly, increasing CIT is another way to achieve high-sensitivity, but the Doppler frequency changes dramatically in high-dynamic conditions, which results in an unacceptable integration loss caused by frequency misalignment. Apparently, neither reducing the LNBW nor increasing CIT is able to satisfy high-sensitivity and high-dynamic performance at the same time. The contradiction analyzed above brings about an issue about how to select proper LNBW and CIT for our SSV RRAM.

2.2. The Parameter Design of SSV RRAM

The carrier tracking loop measurement errors consist of two portions, phase jitter (${\mathsf{\sigma }}_j$), and line-of-sight (LOS) dynamic stress error ${\mathsf{\theta }}_e$. The former is subdivided into three parts: the thermal noise ${\mathsf{\sigma }}_{tPLL}$, the vibration-induced oscillator jitter, and the Allan variance-induced oscillator jitter. When a user is transferred from a normal condition into a weak-signal and high-dynamic environment, the thermal noise increases as the C/N₀ decreases, while the dynamic stress error increases as the dynamics become greater. Apparently, for SSV users, the thermal noise and dynamic stress error are the two dominant PLL error sources. Both of them are far greater than the two kinds of oscillator jitters, so the influence of oscillator jitters can be ignored herein. Acting as the major constituent of phase jitter, the thermal noise is usually expressed in the following equation according to [12]:

$${\mathsf{\sigma }}_{tPLL}=\frac{180}{\mathsf{\pi }}\sqrt{\frac{B_L}{C/N_0}(1+\frac{1}{2T_{coh}\cdot C/N_0})}$$

(1)

where $B_L$ represents the LNBW, and $T_{coh}$ represents the CIT. It must be stressed that $C/N_0$ remains unchanged before and after coherent integration.

For the SSV RRAM, the steady state error caused by LOS dynamic stress is written as:

$${\mathsf{\sigma }}_d=\frac{1}{{\mathsf{\omega }}_n^3}\frac{d^{3}R}{d{t}^3}=\frac{{0.7845}^3}{B_L^3}{J}_{LOS}(m)=\frac{{0.7845}^3}{B_L^3}{J}_{LOS}\cdot \frac{360f_{L1}}{c}(\mathrm{deg})$$

(2)

where c is the speed of light and $f_{L1}$ is the carrier frequency of L1 signal. $\frac{d^{3}R}{d{t}^3}$ is the third-order derivative of LOS range with respect to time, which is equivalent to $J_{LOS}$, i.e., the projection of jerk between navigation satellite and user vehicle in the direction of LOS. ${\mathsf{\omega }}_n$ represents the natural frequency, and has a definite relationship with $B_L$ in the third-order PLL:

$${\mathsf{\omega }}_n=\frac{B_L}{0.7845}$$

(3)

2.2.1. The Determination of CIT

The length of the CIT is restricted to two factors, the length of navigation data bit and integration loss caused by frequency error. In weak-signal and high-dynamic environments, which of the two factors is dominant should be identified with quantitative calculation as follows. In order to keep the integration loss caused by frequency error $f_e$ below 3 dB [13], we have:

$$\left|\sin c(f_e{T}_{coh})\right|\le \frac{1}{\sqrt{2}}$$

(4)

Under high-dynamic circumstances, $f_e$ mainly comes from the relative movements between the navigation satellite and the user vehicle. Take $v_{LOS}$ and $a_{LOS}$ as projections of relative velocity and acceleration from the satellite to the vehicle in the LOS direction, respectively; the Doppler frequency is:

$$f_d=-\frac{f_{L1}\cdot v_{LOS}}{c}$$

(5)

The rate of Doppler frequency ${\dot{f}}_d$ equals the time derivative of $f_d$:

$${\dot{f}}_d=\frac{d{f}_d}{dt}=-\frac{f_{L1}\cdot a_{LOS}}{c}$$

(6)

Furthermore, the rate of ${\dot{f}}_d$ can be expressed with $J_{LOS}$:

$${\ddot{f}}_d=\frac{d^2{f}_d}{d{t}^2}=-\frac{f_{L1}\cdot J_{LOS}}{c}$$

(7)

Then the frequency error during the integration processing can be computed as:

$$f_e=f_d+{\dot{f}}_d\cdot T_{coh}+\frac{1}{2}{\ddot{f}}_d{T}_{coh}^2=-\frac{f_{L1}\cdot v_{LOS}}{c}-\frac{f_{L1}\cdot a_{LOS}}{c}{T}_{coh}-\frac{f_{L1}\cdot J_{LOS}}{2c}{T}_{coh}^2$$

(8)

Assume that the first two terms of the right side of Equation (8), related to $v_{LOS}$ and $a_{LOS}$, respectively, is completely compensated by the third-order loop; then the equation is simplified as:

$$f_e=-\frac{f_{L1}\cdot J_{LOS}}{2c}{T}_{coh}^2$$

(9)

Substituting Equation (9) into Equation (4), we have:

$$\left|\sin c(-\frac{\mathsf{\pi }{f}_{L1}{J}_{LOS}{T}_{coh}^3}{2c})\right|\le \frac{1}{\sqrt{2}}$$

(10)

Thus, the upper limitation of $T_{coh}$ is determined:

$$T_{coh}=1000\left[\sqrt[3]{\frac{0.886c}{\mathsf{\pi }{f}_{L1}{J}_{LOS}}}\right](ms)$$

(11)

where the operator $\left[x\right]$ represents the maximum rounding operation below the value $x$, for the reason that $T_{coh}$ equals the integer multiples of C/A code period (1 ms). Figure 1 presents the bound of $T_{coh}$ over the variation of $J_{LOS}$, and proves that CIT restricted to the length of navigation data bit is dominant when $J_{LOS}$ does not exceed 100 m/s³. With regard to the mission in Section 5, whose peak jerk is 40 m/s³, a CIT of 20 ms would not lead to an unacceptable power loss.

2.2.2. The Determination of LNBW

As discussed above, the errors from oscillator jitters are ignored in SSV applications, and the rule-of-thumb expression of the loop measurement error is expressed as the following equation [13]:

$${\mathsf{\sigma }}_{PLL}={\mathsf{\sigma }}_{tPLL}+\frac{1}{3}{\mathsf{\sigma }}_d=\frac{180}{\mathsf{\pi }}\sqrt{\frac{B_L}{C/N_0}(1+\frac{1}{2T_{coh}C/N_0})}+\frac{57.94f_{L1}{J}_{LOS}}{c{B}_L^3}$$

(12)

The SSV RRAM chooses a two-quadrant arctangent discriminator, which is the most accurate one among all kinds of phase detectors. Referring to prior discoveries [12], the one-sigma rule threshold for the two-quadrant arctangent discriminator is 15°. Then we get the maximum bearable jerk of the third-order PLL over the variation of $B_L$:

$$J_{LOS}=\frac{\left[{15}^{\circ }-\frac{180}{\mathsf{\pi }}\sqrt{\frac{B_L}{C/N_0}(1+\frac{1}{2T_{coh}C/N_0})}\right]c{B}_L^3}{57.94f_{L1}}$$

(13)

In the case that the details of data bit transition are unknown, we set the CIT at 20 ms. Considering the stable condition of the third-order loop is $B_L\cdot T_{coh}<1.0$ [14], $B_L$ changes from 1 Hz to 18 Hz. According to [4], the minimum received C/N₀ is about 14.59 dB-Hz for SSV users. Meanwhile, the normal C/N₀ is about 45 dB-Hz for ground users [13]. Thus, we choose the C/N₀ values ranging from 15 dB-Hz to 45 dB-Hz for weak signal analysis. Then we draw the fluctuation of the maximum bearable jerk in different power levels in Figure 2, where seven different C/N₀ values are selected from 15 dB-Hz to 45 dB-Hz with a uniform spacing of 5 dB-Hz. Without a doubt, a negative jerk is physically meaningless, but the negative jerks in Figure 2 represent the PLL is certainly out of lock under these conditions. The figure shows two important facts:

• The third-order PLL is vulnerable to bear any dynamics when the power becomes very weak (20 dB-Hz or lower);
• The threshold of maximum bearable jerk increases with the growth of LNBW when the received C/N₀ is greater than 30 dB-Hz, but vibrates under a weak signal environment. Thus, there is an optimal LNBW that minimizes the total loop measurement error.

From the two facts, we can conclude that our SSV RRAM is incompetent to realize both high-sensitivity and high-dynamics simultaneously.

3. Adaptive Bandwidth Tracking Strategy

Based on the second fact in the end of last section, this section introduces the optimization of $B_L$. To minimize the thermal noise error, a narrower $B_L$ is desirable. To minimize the LOS dynamic stress error, a wider $B_L$ is necessary. Obviously, $B_L$ should be optimized to bring down the weighted sum of thermal noise and dynamic stress error. The theoretical optimum LNBW and the method to get the optimal solution are introduced in this section.

3.1. Theoretical Optimum LNBW

Computing the partial derivative of ${\mathsf{\sigma }}_{PLL}$ to $B_L$, we get $\frac{\partial {\mathsf{\sigma }}_{PLL}}{\partial B_L}$. When $\frac{\partial {\mathsf{\sigma }}_{PLL}}{\partial B_L}=0$, the value of $B_L$ that make the ${\mathsf{\sigma }}_{PLL}$ minimal can be determined:

$$B_{L\_optimal}=\sqrt[7]{\frac{{173.8}^2{\mathsf{\pi }}^2{J}_{LOS}^2{f}_{L1}^2}{{90}^2{c}^2}\cdot \frac{2T_{coh}{(C/N_0)}^2}{2T_{coh}\cdot C/N_0+1}}$$

(14)

In practical SSV flight missions, both the C/N₀ and $J_{LOS}$ are changing continuously. To achieve an optimal $B_L$, a self-adaptive method should be utilized, and the adaptive KF-based algorithm is exactly such a method.

3.2. Adaptive Four-State KF-Based Algorithm

A three-state (phase error, Doppler frequency, and Doppler frequency rate) adaptive KF-based method was adopted in [15]. To improve the loop performance in weak-signal and high-dynamic environments, we introduce the state of code phase error and extend the three-state to four-state. Replacing DLL (Delay Lock Loop), FLL (Frequency Lock Loop), and PLL filters of conventional tracking structure with a combined adaptive Kalman filter, the block diagram of the four-state KF-based strategy is depicted in Figure 3. The error-state vector consists of code phase error $\mathsf{\delta }{d}_k$, carrier phase error $\mathsf{\delta }{\mathsf{\varphi }}_k$, Doppler frequency error $\mathsf{\delta }{f}_{d,k}$, and its rate error $\mathsf{\delta }{\dot{f}}_{d,k}$ at time $k$:

$${\mathbf{X}}_k={[\mathsf{\delta }{d}_k,\mathsf{\delta }{\mathsf{\varphi }}_k,\mathsf{\delta }{f}_{d,k},\mathsf{\delta }{\dot{f}}_{d,k}]}^T$$

(15)

$\mathsf{\delta }{f}_{d,k}$ equals the difference between the actual Doppler frequency and the frequency feedback to numerically controlled oscillator(NCO). For L1 C/A signal, the ratio of code chip rate $f_{code}=1.023 \mathrm{Mcps}$ and nominal radio frequency $f_{RF}=1575.42 \mathrm{MHz}$ is constant, so the state forward prediction is expressed as:

$${\widehat{\mathbf{X}}}_{k,k-1}={\mathsf{\Phi }}_{k,k-1}{\widehat{\mathbf{X}}}_{k-1}=\left[\begin{array}{cccc}1& 0& \frac{f_{code}{T}_{coh}}{f_{RF}}& \frac{f_{code}{T}_{coh}^2}{2f_{RF}}\\ & 1& 2\mathsf{\pi }{T}_{coh}& \frac{f_{code}{T}_{coh}^2}{2f_{RF}}\\ & & 2\mathsf{\pi }& 2\mathsf{\pi }{T}_{coh}\\ & & & 2\mathsf{\pi }\end{array}\right]{\widehat{\mathbf{X}}}_{k-1}=\left[\begin{array}{cccc}1& 0& \frac{T_{coh}}{1540}& \frac{T_{coh}^2}{2\cdot 1540}\\ & 1& 2\mathsf{\pi }{T}_{coh}& \frac{T_{coh}^2}{2\cdot 1540}\\ & & 2\mathsf{\pi }& 2\mathsf{\pi }{T}_{coh}\\ & & & 2\mathsf{\pi }\end{array}\right]{\widehat{\mathbf{X}}}_{k-1}$$

(16)

Since coherent integration is used to process weak signals, the measurements obtained from discriminators are the average phase differences during CIT, rather than instant phase errors. Thus, we have:

$$\begin{gathered} \mathsf{\delta }{\mathsf{\varphi }}_k=\frac{1}{T_{coh}}{\displaystyle {\int }_0^{T_{coh}}\mathsf{\delta }{\mathsf{\varphi }}_{k-1}}+2\mathsf{\pi }\mathsf{\delta }{f}_{d,k-1}\mathsf{\tau }+\frac{1}{2}2\mathsf{\pi }\mathsf{\delta }{\dot{f}}_{d,k-1}{\mathsf{\tau }}^{2}d\mathsf{\tau }=\mathsf{\delta }{\mathsf{\varphi }}_{k-1}+\mathsf{\pi }\mathsf{\delta }{f}_{d,k-1}{T}_{coh}+\frac{1}{3}\mathsf{\pi }\mathsf{\delta }{\dot{f}}_{d,k-1}{T}_{coh}^2 \\ =\left[\begin{array}{cccc}0& 1& \mathsf{\pi }{T}_{coh}& \frac{1}{3}\mathsf{\pi }{T}_{coh}^2\end{array}\right]X_{k-1} \end{gathered}$$

(17)

Likely, the measurement equation of code phase error is:

$$\mathsf{\delta }{d}_k=\left[\begin{array}{cccc}1& 0& \frac{\mathsf{\pi }{T}_{coh}}{1540}& \frac{\mathsf{\pi }{T}_{coh}^2}{3\cdot 1540}\end{array}\right]X_{k-1}$$

(18)

Rewriting Equations (17) and (18) in the form of matrix and get the measurement matrix:

$${\mathbf{H}}_k=\left[\begin{array}{cccc}1& 0& \frac{\mathsf{\pi }{T}_{coh}}{1540}& \frac{\mathsf{\pi }{T}_{coh}^2}{3\cdot 1540}\\ 0& 1& \mathsf{\pi }{T}_{coh}& \frac{1}{3}\mathsf{\pi }{T}_{coh}^2\end{array}\right]$$

(19)

The formal measurement equation is expressed as:

$${\mathbf{Z}}_k={\mathbf{H}}_k{\mathsf{\Phi }}_{k,k-1}{\mathbf{X}}_{k-1}+v_k$$

(20)

where $v_k$ stands for the difference between actual measurement obtained from discriminator and the theoretical forward prediction. The process noise matrix is computed as follows:

$$\begin{gathered} {\mathbf{Q}}_k=E[({\widehat{\mathbf{X}}}_k-{\widehat{\mathbf{X}}}_{k-1}){({\widehat{\mathbf{X}}}_k-{\widehat{\mathbf{X}}}_{k-1})}^T] \\ =S_f{f}_{L1}^2\left[\begin{array}{cccc}\frac{T_{coh}}{{1540}^2}& \frac{T_{coh}}{1540}& 0& 0\\ \frac{T_{coh}}{1540}& T_{coh}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]+S_g{f}_{L1}^2\left[\begin{array}{cccc}\frac{T_{coh}^3}{3\cdot {1540}^2}& \frac{T_{coh}^3}{3\cdot 1540}& \frac{T_{coh}^2}{2\cdot 1540}& 0\\ \frac{T_{coh}^3}{3\cdot 1540}& \frac{T_{coh}^3}{3}& \frac{T_{coh}^2}{2}& 0\\ \frac{T_{coh}^2}{2\cdot 1540}& \frac{T_{coh}^2}{2}& T_{coh}& 0\\ 0& 0& 0& 0\end{array}\right] \\ +q_a{\left(\frac{f_{L1}}{c}\right)}^2\left[\begin{array}{cccc}\frac{T_{coh}^5}{20\cdot {1540}^2}& \frac{T_{coh}^5}{20\cdot 1540}& \frac{T_{coh}^4}{8\cdot 1540}& \frac{T_{coh}^3}{6\cdot 1540}\\ \frac{T_{coh}^5}{20\cdot 1540}& \frac{T_{coh}^5}{20}& \frac{T_{coh}^4}{8}& \frac{T_{coh}^3}{6}\\ \frac{T_{coh}^4}{8\cdot 1540}& \frac{T_{coh}^4}{8}& \frac{T_{coh}^3}{3}& \frac{T_{coh}^2}{2}\\ \frac{T_{coh}^3}{6\cdot 1540}& \frac{T_{coh}^3}{6}& \frac{T_{coh}^2}{2}& T_{coh}\end{array}\right]+q_d\left[\begin{array}{cccc}{T}_{coh}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right] \end{gathered}$$

(21)

where $q_a$ is the spectral density of LOS acceleration random walk model, $q_d$ is the code density of LOS acceleration random walk model, $S_f$ and $S_g$ are relevant to the phase and frequency of clock oscillator. Generally, $S_f=2h_0$ and $S_g=8{\mathsf{\pi }}^2{h}_{-2}$. The coefficients $h_0$ and $h_{-2}$ are determined by the type of clock oscillator [16,17], which can be referred in

Table 1. Coefficients for oven-controlled crystal oscillator (OCXO) and temperature compensated crystal oscillator (TCXO).

Oscillator Type	${\mathit{h}}_{\mathbf{0}}$ (s2/Hz)	${\mathit{h}}_{-\mathbf{2}}$ (1/Hz)
OCXO	1.0 × 10−21	1.0 × 10−20
TCXO	2.51 × 10−26	1.0 × 10−22

.

Additionally, the measurement noise covariance matrix is related to the output of the discriminator, and it is a function of real-time signal power:

$${\mathbf{R}}_k=\left[\begin{array}{cc}{\mathsf{\sigma }}_d^2& \\ & {\mathsf{\sigma }}_{\mathsf{\varphi }}^2\end{array}\right]=\left[\begin{array}{cc}d(2-d){[\frac{1}{2}+\frac{1}{T_{coh}C/N_0(2-d)}]}^2& \\ & \frac{1}{2T_{coh}C/N_0}(1+\frac{1}{2T_{coh}C/N_0})\mathbf{I}\end{array}\right]$$

(22)

where $d$ is the code phase, ${\mathsf{\sigma }}_d^2$ is noise variance of the DLL discriminator, ${\mathsf{\sigma }}_{\mathsf{\varphi }}^2$ is noise variance of the PLL discriminator, and $\mathbf{I}$ represents the unit diagonal matrix. The iterative process of KF contains the computation of forward prediction covariance matrix ${\mathbf{P}}_{k,k-1}$, adaptive gain matrix ${\mathbf{K}}_k$, new state covariance matrix ${\mathbf{P}}_k$ and new state vector ${\mathbf{X}}_k$:

$$\{\begin{array}{c}{\mathbf{P}}_{k,k-1}={\mathsf{\Phi }}_{k,k-1}{\mathbf{P}}_{k-1}{\mathsf{\Phi }}_{k,k-1}^T+{\mathbf{Q}}_k\\ {\mathbf{K}}_k={\mathbf{P}}_{k,k-1}{\mathbf{H}}_k^T{[{\mathbf{R}}_k+{\mathbf{H}}_k{\mathbf{P}}_{k,k-1}{\mathbf{H}}_k^T]}^{-1}\\ {\mathbf{P}}_k=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_{k,k-1}\\ {\widehat{\mathbf{X}}}_k={\widehat{\mathbf{X}}}_{k,k-1}+{\mathbf{K}}_k{v}_k\end{array}$$

(23)

Compared to conventional tracking strategy, the advantage of an adaptive KF is that the equivalent LNBW can be easily adjusted by tuning the coefficients in the noise matrix. In-depth analysis was given by [18] about how the adaptive KF varies the LNBW in a time-varying optimal manner as the signal power and dynamics change. The equivalent LNBW in steady state is:

$${(B_{L\_eq})}_k\cong \frac{{\mathbf{K}}_k({\mathbf{Q}}_k,{\mathbf{R}}_k)}{c_3{T}_{coh}}=\frac{({\mathsf{\Phi }}_{k,k-1}{\mathbf{P}}_{k-1}{\mathsf{\Phi }}_{k,k-1}^T+{\mathbf{Q}}_k){\mathbf{H}}_k^T{[{\mathbf{R}}_k+{\mathbf{H}}_k({\mathsf{\Phi }}_{k,k-1}{\mathbf{P}}_{k-1}{\mathsf{\Phi }}_{k,k-1}^T+{\mathbf{Q}}_k){\mathbf{H}}_k^T]}^{-1}}{c_3{T}_{coh}}$$

(24)

where the filter coefficient $c_3$ for the third-order PLL equals 3.048 and the gain matrix ${\mathbf{K}}_k$ is a function of ${\mathbf{Q}}_k$ and ${\mathbf{R}}_k$. It can be proved than ${\mathbf{K}}_k$ is positively correlated with ${\mathbf{Q}}_k$ and negatively correlated with ${\mathbf{R}}_k$, thus ${\mathbf{K}}_k\propto \frac{{\mathbf{Q}}_k}{{\mathbf{R}}_k}$. This indicates the realization of adaptive bandwidth in weak-signal and high dynamic environment:

• When the signal turns weaker, the C/N₀ becomes lower, so that the measurement noise ${\mathbf{R}}_k$ becomes greater, ${\mathbf{K}}_k$ decreases further and, finally, the LNBW is narrowed;
• When the LOS dynamic increases, the process noise ${\mathbf{Q}}_k$ becomes greater, ${\mathbf{K}}_k$ increases further and, finally, the LNBW is widened.

4. INS-Assisted Open Loop Tracking Strategy

As an opposite method against the closed loop, the idea of the open loop has been discussed for many years. In our common sense, open loop tracking is considered as block processing or batch processing. Contrary to closed loop sequential processing, the open loop batch processing has no local loop updates [19]. The use of open loop batch processing allows parallel computations on the batches of samples, and brings the benefits of improved signal observability and improved tracking robustness [20]. This section illustrates the structure of INS-assisted tracking strategy and the frequency estimation algorithm adopted in this strategy.

4.1. Structure of INS-Assisted Tracking Strategy

It must be admitted that OL batch processing is less favorable than the CL sequential method with respect to computation efficiency, resulting from the complicated FFT calculation. Additionally, the OL batch scheme does not improve the measurement accuracy as compared to the CL sequential architecture. Instead, the OL scheme mainly promotes the robustness and reliability of the tracking loop, especially in challenging environment, such as SSV applications. Thus, it is necessary to put forward a new OL structure to ease the computation burden. In Figure 4, a suggested OL architecture is proposed with this purpose.

The proposed architecture is built on the basis of Doppler prediction utilizing inertial measurement unit (IMU) and GNSS ephemeris. IMU is a cluster of sensors, including accelerometers and gyros, which is responsible for obtaining the motion and attitude of the SSV user receiver, whereas the GNSS satellite ephemeris is used to calculate the position and velocity of the satellite by accurate orbit propagation models. Referring to [21], the Doppler frequency shift caused by the motion of the SSV user receiver relative to the GNSS satellite can be easily determined. The estimation of the Doppler frequency is used to reduce the local frequency search space, and the estimation accuracy depends on the quality of IMU. If the space vehicles are equipped with high quality IMUs, the expected Doppler frequency error accumulates slowly over time, which is adequate for SSV navigation. The ephemeris is obtained during visible intervals and saved into onboard registers for orbital propagation. It has been proved that the outage time when no satellite is in view does not exceed 108 min for SSV flight missions [1], and the length of the orbit propagation is relatively short so that the dynamic aiding is accurate enough.

After A/D (Analog-to-Digital) conversion, the digitized incoming IF signal is generally expressed as:

$$s(k)=AC(k-\mathsf{\tau })D(k)\cos [2\mathsf{\pi }kT(f_{IF}+f_d)+\mathsf{\phi }]+n(k)$$

(25)

where $k$ is the time index of each sampling time interval, and $T$ is the sampling time interval; $A$ is the signal amplitude, $C(\cdot )$ is the spreading code sequence and $D(\cdot )$ is the navigation data bit; the three parameters to be estimated are code delay $\mathsf{\tau }$, Doppler frequency shift $f_d$, and carrier phase $\mathsf{\phi }$; $f_{IF}$ is IF in Hz, and $n(k)$ is additive band-limited white Gaussian noise (AWGN).

Using the dynamic aiding with IMU and ephemeris described in Figure 4, we get the estimated Doppler frequency ${\widehat{f}}_d$. If the estimated initial phase is ${\widehat{\mathsf{\varphi }}}_0$, the two orthogonal carriers representing the outputs of carrier NCO are:

$$\{\begin{array}{c}{i}_{LO}(k)=\cos [2\mathsf{\pi }(f_{IF}+{\widehat{f}}_d)kT+{\widehat{\mathsf{\varphi }}}_0]\\ q_{LO}(k)=\sin [2\mathsf{\pi }(f_{IF}+{\widehat{f}}_d)kT+{\widehat{\mathsf{\varphi }}}_0]\end{array}$$

(26)

The received signals in Equation (25) are correlated with the local replica signals in Equation (26) to generate the correlated results as:

$$\{\begin{array}{c}I(k)=AR(\mathsf{\delta }\mathsf{\tau })D(k)\cos (2\mathsf{\pi }kT\cdot \mathsf{\delta }f+\mathsf{\delta }{\mathsf{\varphi }}_0)+n_i(kT)\\ Q(k)=AR(\mathsf{\delta }\mathsf{\tau })D(k)\sin (2\mathsf{\pi }kT\cdot \mathsf{\delta }f+\mathsf{\delta }{\mathsf{\varphi }}_0)+n_q(kT)\end{array}$$

(27)

where $\mathsf{\delta }f$ is the Doppler frequency estimation error, $\mathsf{\delta }{\mathsf{\varphi }}_0$ is the initial phase estimation error, $\mathsf{\delta }\mathsf{\tau }$ is the code delay error, $R(\cdot )$ is the autocorrelation function (ACF) of C/A code, $n_i$ is the in-phase noise and $n_q$ is the quadrature noise. After coherent integration, we get the integration results of $I_P$ and $Q_P$ [22]:

$$\{\begin{array}{l}{I}_P(k)=AR(\mathsf{\delta }\mathsf{\tau })D(k)\sin c(kT\cdot \mathsf{\delta }f)\cos [2\mathsf{\pi }\mathsf{\delta }f(kT+\frac{T_{coh}}{2})+\mathsf{\delta }{\mathsf{\varphi }}_0]\\ Q_P(k)=AR(\mathsf{\delta }\mathsf{\tau })D(k)\sin c(kT\cdot \mathsf{\delta }f)\sin [2\mathsf{\pi }\mathsf{\delta }f(kT+\frac{T_{coh}}{2})+\mathsf{\delta }{\mathsf{\varphi }}_0]\end{array}$$

(28)

It is obvious that the estimation of $\mathsf{\delta }{\mathsf{\varphi }}_0$ and $\mathsf{\delta }f$ are tightly coupled. The PLL discriminator output $\mathsf{\delta }{\mathsf{\varphi }}_0$ is influenced by $\mathsf{\delta }f$, while the FLL discriminator output $\mathsf{\delta }f$ is influenced by $\mathsf{\delta }{\mathsf{\varphi }}_0$. The calculation results obtained directly from $I_P$ and $Q_P$ in Equation (28) are unreliable, so a non-coherent MLE-based FLL discriminator is applied in the architecture shown in Figure 4 to settle the problem.

4.2. Non-Coherent MLE Algorithm

Hereinafter, a non-coherent MLE algorithm to estimate Doppler frequency error that is insensitive to carrier phase error is introduced. By computing the sum of the squares of $I_P$ and $Q_P$, we have the value of non-coherent integration (NI):

$$\begin{gathered} V=I_P^2(k)+Q_P^2(k)=A^2{[R(\mathsf{\delta }\mathsf{\tau })]}^2{[D(k)]}^2{[\sin c(kT\cdot \mathsf{\delta }f)]}^2 \\ =A^2{[R(\mathsf{\delta }\mathsf{\tau })]}^2{[\sin c(kT\cdot \mathsf{\delta }f)]}^2 \end{gathered}$$

(29)

It is obvious that NI brings about two positive features: First, the value of $V$ is not affected by data bit $D(k)$ regardless of whether $D(k)$ equals +1 or −1, so the influence of data bit transition is eliminated; Second, the value of $V$ has nothing to do with carrier phase error $\mathsf{\delta }{\mathsf{\varphi }}_0$, and it changes with $\mathsf{\delta }f$. The two features induce us to construct a maximum likelihood cost function for Doppler estimation using the NI result, and the expression is as follows:

$$L\left(\mathsf{\delta }f\right)={\displaystyle \sum _{k=1}^N[I_P^2(k)+Q_P^2(k)]}={\displaystyle \sum _{k=1}^N{A}^2{[R(\mathsf{\delta }\mathsf{\tau })]}^2{[\sin c(kT\cdot \mathsf{\delta }f)]}^2}$$

(30)

where $N$ represents the total number of IF data samples within a given CIT. According to the normal MLE solving method, $\mathsf{\delta }f$ reaches its MLE where the partial derivative of the likelihood cost function or its log-likelihood function [23] with respect to $\mathsf{\delta }f$ is 0. However, obviously, it is too complicated to write the analytic expression of the MLE of $\mathsf{\delta }f$. Thus, we rewrite the cost Function (30) as:

$$L\left(\mathsf{\delta }f\right)={\displaystyle \sum _{k=1}^N[I_P^2(k)+Q_P^2(k)]}={\displaystyle \sum _{k=1}^N{[s(k)\cdot C\cdot \cos ]}^2}+{\displaystyle \sum _{k=1}^N{[s(k)\cdot C\cdot \sin ]}^2}$$

(31)

For the sake of simplicity, we use $C$ as $C(k-\mathsf{\tau })$, $\cos$ as $\cos [2\mathsf{\pi }(f_{IF}+\mathsf{\delta }f)kT+\mathsf{\delta }\mathsf{\varphi }]$, and $\sin$ as $\sin [2\mathsf{\pi }(f_{IF}+\mathsf{\delta }f)kT+\mathsf{\delta }\mathsf{\varphi }]$. Then the gradient of the new cost function for the Doppler frequency is given by:

$$\begin{gathered} Gradient=\frac{\partial L}{\partial \mathsf{\delta }f}=-4\mathsf{\pi }T{\displaystyle \sum _{k=1}^N[s(k)\cdot C\cdot \cos ]\times }{\displaystyle \sum _{k=1}^N[k\cdot s(k)\cdot C\cdot \sin ]} \\ +4\mathsf{\pi }T{\displaystyle \sum _{k=1}^N[s(k)\cdot C\cdot \sin ]\times }{\displaystyle \sum _{k=1}^N[k\cdot s(k)\cdot C\cdot \cos ]}=-4\mathsf{\pi }T(I_P{\overline{Q}}_P-{\overline{I}}_P{Q}_P) \end{gathered}$$

(32)

where:

$$\{\begin{array}{c}{I}_P={\displaystyle \sum _{k=1}^{N}s(k)\cdot C\cdot \cos }\\ Q_P={\displaystyle \sum _{k=1}^{N}s(k)\cdot C\cdot \sin }\\ {\overline{I}}_P={\displaystyle \sum _{k=1}^{N}k\cdot s(k)\cdot C\cdot \cos }\\ {\overline{Q}}_P={\displaystyle \sum _{k=1}^{N}k\cdot s(k)\cdot C\cdot s}in\end{array}$$

(33)

The Hessian [24] of the proposed cost function with regard to $\mathsf{\delta }f$ is given by:

$$Hessian=\frac{{\partial }^{2}L}{\partial \mathsf{\delta }{f}^2}=8{\mathsf{\pi }}^2{T}^2({\overline{I}}_P^2+{\overline{Q}}_P^2-I_P{\stackrel{=}{I}}_P-Q_P{\stackrel{=}{Q}}_P)$$

(34)

where:

$$\{\begin{array}{c}{\stackrel{=}{I}}_P={\displaystyle \sum _{k=1}^N{k}^2\cdot s(k)\cdot C\cdot \cos }\\ {\stackrel{=}{I}}_P={\displaystyle \sum _{k=1}^N{k}^2\cdot s(k)\cdot C\cdot s}in\end{array}$$

(35)

Referring to [25], the frequency error is estimated with the gradient divided by the Hessian:

$$\mathsf{\delta }f=-\frac{Gradient}{Hessian}=\frac{I_P{\overline{Q}}_P-{\overline{I}}_P{Q}_P}{2\mathsf{\pi }T({\overline{I}}_P^2+{\overline{Q}}_P^2-I_P{\stackrel{=}{I}}_P-Q_P{\stackrel{=}{Q}}_P)}$$

(36)

The calculation is completed in the block named MLE-based FLL discriminator in Figure 4. Finally, the coupling effect between Doppler frequency error and carrier phase error is cancelled.

In addition, the code delay $\mathsf{\tau }$ is still estimated with CL scheme using a carrier-aided DLL. There are two early-minus-late (EML) type correlators with outputs:

$$\{\begin{array}{l}{I}_{EML}={\displaystyle \sum _{k=1}^{N}s(k)\cdot \dot{C}\cdot \cos }\\ Q_{EML}={\displaystyle \sum _{k=1}^{N}s(k)\cdot \dot{C}\cdot s}in\end{array}$$

(37)

where the new symbol $\dot{C}$ is the derivative of $C$ with respect to $\mathsf{\tau }$. For a definite correlator spacing $d$, it is calculated with the following equation:

$$\dot{C}=\frac{\partial C(k-\mathsf{\tau })}{\partial \mathsf{\tau }}\approx \frac{C(k-\mathsf{\tau }-\frac{d}{2})-C(k-\mathsf{\tau }+\frac{d}{2})}{d}$$

(38)

5. Simulation and Experiment Results

Under normal conditions, the comparison between loop measurement error and its threshold is used to judge whether the loop is out of lock. Considering the balance of both measurement accuracy and rapid response capability, a threshold of 15° is regarded as a judgment indicator about whether the loop is out of lock or not in this section. Hereafter, a testing scenario is built up at first, then the experimental system and initial settings are introduced, and the Doppler frequency estimation error is used to weigh the performance of our tracking strategies.

5.1. Scenario Settings

In order to identify the tracking quality of our proposed strategies, we establish a scenario of a lunar upper stage on a HwaCreat™ GNSS signal simulator (produced by Beijing Hwa Creat Technology Corporation Ltd., Beijing, China). In the scenario, the lunar exploration probe operates in its Earth phasing orbit [4], whose trajectory traverses both TSV and SSV in a geostationary transfer orbit (GTO). One certain GPS satellite is chosen as the signal source, and the signal emitted from the satellite is used for analysis when it is visible to our object vehicle. The simulated scenario is drawn in Figure 5. We can see the lunar upper stage receives the signal emitted by the GPS satellite and transmitted over the limb of the Earth, which is consistent with the basic characterization of SSV navigation. The simulation time lasts for 320 s, which is divided into four phases. During the simulation time, the lunar probe operates at an altitude of about 35,000 km, and transfers from its original phasing orbit to another phasing orbit by means of an orbital maneuver. In the first 50 s, the GPS satellite is not visible to the lunar probe. From the 51st s to 300th s, the GPS satellite enters the visible area of the lunar probe. In this period, the probe operates under gravity, and completely free from its own drag force. In the third phase, from the 301st s to the 310th s, the onboard motor starts to work and generates a ramp LOS jerk. The value of the LOS jerk reaches its maximum at the 310th s and keeps constant in the fourth phase. The dynamic performance requirements for the GNSS receiver of the lunar upper stage can be seen in

Table 2. The requirements of bearable dynamics for onboard GNSS receiver of the lunar probe.

Dynamic Type	Sustainable Range
Velocity	0~12 km/s
Acceleration	0~15 g
Jerk	0~4 g/s

. If the Doppler shift caused by the LOS velocity and acceleration are completely compensated by the third-order PLL, the bearable jerk of 4 g/s becomes the most important requirement. Thus, we set the maximal LOS jerk at 40 m/s³ with a little excess. The mission control sequence (MCS) of the scenario is shown in

Table 3. The mission control sequence (MCS) of the lunar probe during the simulation time.

Phase	Time Span	Visibility of GPS Satellite	Orbital Maneuvering Status
I	0~50 s	Not in view	No maneuver is made
II	51~300 s	In view	No maneuver is made
III	301~310 s	In view	Make amaneuver with a ramp LOS jerk (the peak is 40 m/s3)
IV	311~320 s	In view	Make amaneuverwith a constant LOS jerk of 40 m/s3

.

Seen from Figure 6a, the received C/N₀ and LOS jerk drawn in different colors are both changing over time. In the second phase, the LOS jerks are about 2.9~3.6 × 10⁻⁵ m/s³, which are too minor to be plotted clearly in the global view. Thus, we repaint the close-up view of LOS jerk from the 50th to 300th s in Figure 6b.

5.2. Experimental System and Initialization

This subsection puts emphasis on the introduction of the experimental equipment being applied and their initialization settings. There are three instruments in the testing system: a HwaCreat™ GNSS signal simulator, an IF sampler, as well as a SDR (software-defined receiver). Their connection is shown in Figure 7. As aforementioned, we build up the scenario containing the lunar upper stage and a certain GPS satellite in the GNSS signal simulator. The radio frequency (RF) export of the simulator is connected to the IF sampler, so that the simulated signals can be collected and stored in the sampler. The sampler is capable of playback to ensure that the recorded data can be used for testing repeatedly, and make it possible that sufficient data can be obtained for statistical analysis. Finally, the digitized IF data is delivered into the SDR for processing. Note that the SDR is designed as three different structures in the following order: SSV RRAM, adaptive four-state KF-based structure illustrated in Figure 3, and INS-assisted structure illustrated in Figure 4.

For the conventional RRAM of CL form, the bandwidth of DLL, FLL, and PLL are 0.1 Hz, 2 Hz, and 18 Hz, respectively. The IF of the SDR is set at 4.092 MHz, whereas the sampling frequency is 16.368 MHz.

For the adaptive four-state KF-based strategy of CL form, the initial covariance matrix ${\mathbf{P}}_0=diag(P_{00},P_{11},P_{22},P_{33})$, where $P_{00}=1chi{p}^2,P_{11}=1,P_{22}={(50Hz)}^2,P_{33}={(1Hz/s)}^2$, and the initial state ${\mathbf{X}}_0={(0,0,0,0)}^T$. Note that the values of the initialization can only affect the convergence speed of the filter instead of the estimation accuracy.

For the INS-assisted strategy of the OL form, the local Doppler frequency is predicted with the help of the input IMU measurements and GPS ephemeris. There are four different grade IMUs being used in the experiments: MEMS (micro electro mechanical system) grade, civil grade, tactical grade, and navigation grade. Their parameters are provided in

Table 4. Input IMU parameters [27] for local frequency prediction in SDR.

IMU Parameters		MEMS IMU	Civil IMU	Tactical IMU	Navigation IMU
gyro errors	constant bias, °/h	300	50	1	0.01
accelerometer errors	constant bias, mg	10	1	0.1	0.01

. If the initial misalignment error is $\mathsf{\delta }{a}_0$, the cumulative acceleration error can be modeled as a simplified form:

$$\mathsf{\delta }{a}_{LOS}=(\mathsf{\delta }{a}_0+\nabla +g\cdot \mathsf{\epsilon }\cdot \Delta t)\cos \mathsf{\theta }$$

(39)

where $\nabla$ is the accelerometer error, $\mathsf{\epsilon }$ is the gyro error and $\Delta t$ is the drift time after the latest IMU calibration [26], and $\mathsf{\theta }$ is the real-time projection angle between the IMU measurement vector and the LOS direction. In the case when the IMU correction to prevent drift error is unavailable, the estimated frequency error over drift time is:

$$\mathsf{\delta }{f}_d=-\frac{f_{L1}}{c}{\displaystyle {\int }_0^{\Delta t}(\mathsf{\delta }{a}_0+\nabla +g\cdot \mathsf{\epsilon }\cdot \mathsf{\tau })}\cos \mathsf{\theta }(\mathsf{\tau })d\mathsf{\tau }$$

(40)

5.3. Result Comparison

Seen from Figure 8, during the second phase, the loop measurement error obtained by an adaptive KF-based structure is roughly the same as that obtained by the conventional structure, namely SSV RRAM. However, in the third and fourth phases, the adaptive KF makes some improvement to reduce the measurement error, compared to SSV RRAM. It is clearly shown that the measurement errors of both the two structures exceed the tracking threshold of 15° after the 300th s. The result reveals several conclusions. First, the two schemes of CL form can work properly in non-maneuvering status when the C/N₀ is over 16.85 dB-Hz. Second, adaptive KF performs better than SSV RRAM due to its flexibility of regulating LNBW in harsh environment. Third, both the two structures lose their locks in orbital maneuvering status and they are incapable of accomplishing navigation in the orbital transfer period of the lunar upper stage.

Figure 9a shows a comparison of the Doppler frequency estimations of different tracking schemes in the same experiment scenario. The results are obtained through repeated trials for each scheme. In Figure 9b, no maneuver is made, all of the six schemes keep track of the target signal but the estimation accuracies diverge from each other. In Figure 9c, the probe makes an orbital maneuver, then the two CL schemes, SSV RRAM, as well as the adaptive KF-based method, are unavailable because the loop is out of lock. To demonstrate the performances of all these schemes clearly, we zoom in the plot around the 300th s with higher resolution in Figure 9d, and it is obvious that the estimation accuracy of INS-assisted strategy is determined by the quality of the IMUs. The IMU of the navigation grade is the best one, the tactical grade is secondary, the civil grade is tertiary, and the MEMS grade is the worst.

Figure 10 provides the comparison of the estimation error of the six schemes. Through statistical computing of the estimation errors in the second phase, the root mean square (RMS) is 2.3188 Hz for SSV RRAM and 0.4324 Hz for the adaptive KF-based method. This leads to the conclusion that the adaptive KF-based method can get more precise estimates of the Doppler frequency than the conventional SSV RRAM when the space vehicle operates under normal or non-maneuvering condition. For the INS-assisted schemes of OL form, the frequency errors accumulate continuously. The maximum frequency errors are shown in

Table 5. The Doppler estimation errors of four OL tracking schemes with a drift time of 270 s.

OL Tracking Schemes	Aided by MEMS IMU	Aided by Civil IMU	Aided by Tactical IMU	Aided by Navigation IMU
Maximum frequencyError (Hz)	1564.6000	260.6740	5.2246	0.0535

after a drift time of 270 s. It is obvious that the cumulative errors of the schemes aided by MEMS IMU and civil IMU grow quickly, which does not meet the performance requirements of SSV navigation. If the drift time without IMU calibration is short enough, e.g., within 77 s, the scheme aided by tactical IMU performs a little better than the adaptive KF-based method. However, once the drift time is over 78 s, the estimation error using the tactical IMU is above 0.4324 Hz which is inferior to the adaptive KF-based method. Although the scheme aided by navigation IMU is more accurate, the estimation error would exceed 0.4324 Hz after a cumulative time of 775 s. Therefore, if the drift time is too long or the dynamic is not too great, the adaptive KF-based method has its specific advantage compared to INS-assisted schemes. Therefore, what the OL tracking strategy actually improves is the loop robustness of GNSS signals tracking function, owing to the fact that OL tracking with high-quality IMU can work properly under orbital maneuvering conditions.

6. Conclusions

This section summarizes the contributions of this paper as follows. We first interpret why it is difficult to achieve high-sensitivity and high-dynamic GNSS signal tracking at the same time. Then a conventional design of a third-order PLL is put forward, and taken as the baseline, i.e., SSV RRAM. Taking both the low C/N₀ and high LOS jerk into account, an optimized strategy of CL form is recommended, and this strategy can be implemented by an adaptive four-state KF-based algorithm, which can adjust its LNBW according to the received signal power and LOS dynamics. The paper also describes the relationship between the filter gain matrix and equivalent LNBW. In order to prevent the PLL from losing lock, an INS-assisted strategy is adopted, and the decoupling of the Doppler frequency estimate and carrier phase estimate should be attributed to the non-coherent processing. Finally, a mission of a lunar upper stage in its Earth phasing orbit is used to test the validity of these strategies. The simulation results substantially prove the superiority of an adaptive four-state KF-based strategy under non-maneuvering conditions and the advantage of an OL tracking strategy aided by high-quality IMUs under maneuvering conditions.

It must be noted that the analyses in the paper are all based on the steady state performance of the tracking structure, but the transient response is another key point that requires a lot of in-depth studies. Therefore, future work will focus on the improvement of transient capability of the signal tracking strategy in the SSV.