GPS Receiver VFLL-Assisted PLL High Dynamic Weak Signal Tracking Based on a Maximum Likelihood Estimation

Abstract

This paper proposes a GPS receiver vector frequency-locked loop-assisted phase-locked loop (VFAPLL) structure based on the maximum likelihood estimation (MLE) method for highly dynamic weak-signal scenarios. In this structure, the loop structure does not include a frequency discriminator, and the signal is directly input to the navigation filter after down-conversion, coherent integration, and other processing to avoid nonlinear noise error. Due to the high dimension and nonlinearity of the cost function of the MLE algorithm, the Levenberg Marquardt (LM) algorithm is used to optimize it. The proposed VFAPLL is compared with the VFAPLL implemented based on the extended Kalman filter (EKF) algorithm and the frequency locked loop assisted phase locked loop (FAPLL) implemented based on MLE. Through simulation verification, it was shown that the VFAPLL (MLE) has higher tracking accuracy, lower loss-of-lock threshold, and better robustness to the input signal than the other two loops.

1. Introduction

The traditional frequency locked loop assisted phase locked loop (FAPLL) couples the frequency locked loop (FLL) and the phase locked loop (PLL), where the FLL is re sponsible for the coarse acquisition of the frequency, and the PLL is responsible for the accurate tracking of the low dynamic part of the signal. It can better solve the signal tracking problem under the high dynamic stress of the receiver. However, the tracking performance of the FAPLL is unsatisfactory in extremely challenging environments, such as urban canyons, extremely high dynamics in the military field, and automatic driving.

In fact, each signal channel in the receiver is related. If we can establish a connection between each channel through the position and speed of the receiver to realize information sharing, then the tracking performance of the receiver can be improved to a certain extent. Especially in highly dynamic scenes and under weak signal conditions with a low carrier-to-noise ratio, this information-sharing loop shows better signal locking ability and higher tracking accuracy than the scalar tracking loop; that is, a strong signal can assist the tracking of a weak signal, and it can maintain a certain precision of tracking even when the signal is blocked and can realize fast recapture when the signal is recovered, which is the standard approach to vector tracking. The concept of a vector tracking loop was first proposed by Spilker [1]; subsequently, vector tracking technology gradually attracted the attention of researchers in related fields who carried out research on this technology.

Pany et al. proposed a carrier-to-noise ratio estimation method suitable for vector tracking under weak signal conditions [2]. Dr Lashley studied the tracking performance of vector tracking in weak signal and highly dynamic environments and further demonstrated the advantages of vector tracking technology [3,4,5]. Domestic scholars have also contributed to the development of vector tracking technology. The more representative studies were produced from Tsinghua University and the National University of Defense Technology [6,7].

Today’s vector tracking loops include the vector frequency locked loop (VFLL), vector delay locked loop (VDLL) and vector phase locked loop (VPLL). The combination of a FAPLL and vector tracking technology is a very promising development direction. In this paper, VFLL is used to assist the PLL. The advantages of vector tracking and the FLL are that they are fully utilized in the auxiliary loop, which avoids time-consuming and energy-consuming reacquisition after signal interruption, can tolerate the high dynamic stress of users more robustly, and can track signals with lower signal-to-noise ratios.

Scholars have been working to improve the tracking performance and robustness of FAPLLs. This loop essentially aids PLLs by estimating the Doppler frequency shift, whereas Kalman loops or other circuit forms can also be used for assistance in estimating frequency. Lashley proposed several carrier tracking loop structures combining vectors and scalars in [8], including the VFLL-assisted PLL, local carrier phase Kalman filter cascaded VDFLL, and cooperative carrier tracking. The VFLL-assisted PLL structure combines the carrier frequency error of the navigation processor with the phase discriminator. The position and clock offset of the navigation processor are predicted and updated by the code discriminator in the VFLL-assisted PLL structure. The speed measurement value is obtained by combining the carrier frequency and discriminator of each channel. Roncagliolo proposed a new loop structure called the unambiguous frequency-assisted phase locked loop (UFA-PLL). In terms of tracking ability and antijamming, this loop has the same advantages as commonly used coupled loops, but its design and implementation are simpler. An optimal method for smoothing phase estimation has also been proposed [9]. Stefan et al. proposed a vector receiver combined with an FAPLL [10], where the VFLL receiver outputs the frequency error used in the FAPLL. The loop proposed in this paper is more robust to carrier phase tracking than conventional FAPLLs, even in harsh signal environments.

Yang Rong et al. in China proposed an iterative FAPLL designed in the state space based on the minimum mean square error algorithm through constructing a high-fidelity carrier signal model [11]. Another study combined the FLL and PLL through an ensemble filter. If a decentralized state estimation framework is used for the FAPLL, then the measurements can be distributed and filtered in the local state estimator and are thus able to combine the local estimates by a weighting matrix to update the global carrier state [12].

To address the problem that the VFLL is not sufficiently accurate and unable to lock the carrier phase and demodulate the received signal correctly, a phase-locked loop-assisted vector frequency-locked loop method was proposed [13]. The method adds the output of the loop filter in the PLL to the prediction of the carrier frequency, which improves the prediction accuracy of the VFLL carrier phase. After that, the accuracy of the demodulated value is improved by compensating for the residual error of the carrier phase.

The loop filter of the FAPLL in [14] is realized by bilinear transformation from the S domain to the Z domain, and a weight adjustment module using an automatic adjustment strategy is designed to calculate the weight factor to adjust the effects of the FLL and PLL. In [15], the authors proposed a tracking loop structure to achieve carrier signal phase tracking in the VFLL, in which a loop filter and reset module in the channel will assist in tracking carrier phase and the signal processing period is matched with a preprocessing filter. Based on fuzzy control theory, an algorithm for the FLL auxiliary PLL structure is proposed, the algorithm uses the output of discriminator as the measured residuals of the Kalman filter [16].

In the above examples from the literature, in-depth research has been carried out on the FAPLL control loop in the receiver, and the performance of the corresponding receiver was improved. However, in the extremely harsh environment of highly dynamic weak signals and some application fields that require extremely high navigation accuracy, further improvement of the loop performance is needed.

From the perspective of the algorithm, the performance improvement of the carrier loop generally occurs in one of two ways. One way is to replace the loop filter with a nonlinear filtering algorithm for the traditional discriminator and loop filter structure to adaptively adjust the loop bandwidth. This improves the tracking performance of the system while considering the dynamic performance and steady-state performance of the system [17]. In [18], the authors evaluated the tracking performance of an adaptive third-order look-up table direct state Kalman filter in the loop bandwidth control algorithm FAPLL. Firstly, the system measurement model, state space model, and the transfer function of FAPLL were explained, and the relationship between the direct state Kalman filter and the FAPLL was analyzed. Then, by solving the discrete algebraic Riccati equation, the convergence of Kalman gain was calculated, and the look-up table direct state Kalman filter was derived, thus reducing the complexity of direct state Kalman filter. However, as the discriminator contained in the loop has only a small linear interval, nonlinear noise errors were introduced, and the loop performance degrades. The other way is to design the carrier tracking loop based on parameter estimation theory, which removes the discriminator in the traditional loop and replaces it with a nonlinear filter and uses the coherent integration result of the signal as the input of the navigation filter, which can improve the loop performance, noise immunity, and sensitivity.

In the above literature, the FLL or VFLL auxiliary PLL typically uses a discriminator to discriminate the phase or frequency difference between the input carrier and the replica carrier. As the discriminator has only a small linear interval, when the carrier-to-noise ratio C/No falls to an error close to or below the level required for the Kramer–Roman lower bound, the discriminator’s measurement error increases rapidly. This leads to a degradation in the tracking performance of closed-loop structures [19].

Therefore, according to parameter estimation theory, this paper proposes a VFLL-assisted PLL based on the maximum likelihood estimation algorithm. The loop no longer uses the discriminator but directly uses the correlation result as the input of the navigation filter and then uses the maximum likelihood estimation to process the input. The nonlinear filtering method can increase the linear interval compared with the direct use of the discriminator, while the nonlinear operational noise is lower [20]. Additionally, unlike Kalman filters, ML estimators are not recursive, and each estimate uses enough signal data to achieve good weak signal performance [19].

However, if the cost function of MLE algorithm has high dimensionality and nonlinearity, this will lead to low computational efficiency and no closed form of solution. A gradient method, such as steepest descent method, Levenberg–Marquardt method, etc., can be used to minimize the cost function through an iterative process. In [21], a maximum likelihood estimation method based on the Levenberg–Marquardt algorithm and message tree transmission idea was proposed. This algorithm not only has high positioning accuracy but can also allow convergence to the same solution as the centralized algorithm. A maximum likelihood estimation method for cases where the parameters are matrices and vectors are not measurable was proposed in [22]. Compared with the multivariable recursive extended least squares algorithm, the multivariable maximum likelihood iterative least squares algorithm has higher estimation accuracy and lower computational complexity. The maximum likelihood method is very important for parameter estimation and system modeling. In [23], the maximum likelihood principle was combined with data filtering technology to estimate the parameters of a class of bilinear systems, with the proposal of a maximum likelihood iterative least squares algorithm based on filtering, which was compared with an iterative algorithm based on least squares. The simulation results showed that the proposed algorithm had higher accuracy and efficiency under different noise variances.

The content in this paper is arranged as follows: Section 2 gives the operation mechanism of the proposed VFLL-assisted PLL carrier tracking loop and uses the jerk model to model this loop under high dynamics. The state equation and observation equation of the loop are also provided. Section 3 represents the implementation process of the maximum likelihood estimation method and its corresponding optimization algorithm, the Levenberg-Marquardt algorithm, and the specific operation process of the VFAPLL based on the MLE-LM algorithm. The intermediate frequency data are repeatedly correlated until maximum likelihood convergence is achieved. Section 4 employs a software receiver in MATLAB to verify the proposed loop. Section 5 summarizes the paper.

2. VFLL-Assisted PLL Carrier Tracking Loop Design

Vector tracking technology is considered to be a promising next-generation receiver technology. It couples all signal tracking channels together through a certain method and directly estimates the user’s position and clock error. The estimated state information is further converted into a control quantity for feedback control of the tracking loop. By realizing information sharing among each signal channel, the vector tracking loop can provide better performance in low signal-to-noise ratios and highly dynamic environments without adding any external assistance. Although the traditional vector tracking loop couples the output of the discriminator in each channel, the measurement residuals of each signal channel are still estimated independently [19]. Based on this situation, the VFLL in this paper no longer uses the frequency discriminator [24], and the coherent integration result of the I/Q branch is directly used as the observation value of the filter so that the receiver will easily operate under the condition of a low signal-to-noise ratio due to the nonlinear noise introduced by the frequency discriminator, which promotes its anti-noise performance. In addition, this paper combines the maximum likelihood estimation method with vector tracking in the VFLL and uses this VFLL to assist the PLL to reduce the tracking threshold of the carrier tracking loop.

To adapt to the highly dynamic and weak-signal environment, a second order VFLL-assisted third order PLL was designed with reference to the FAPLL structure. Figure 1 is a block diagram of the VFLL-assisted PLL based on the maximum likelihood estimation algorithm, and Figure 2 is a block diagram of the FAPLL based on the maximum likelihood estimation algorithm.

The carrier tracking loop of the i1 channel is shown in Figure 2, the digital intermediate frequency signal is first multiplied and correlated by the local signal and stripped the carrier and doppler frequency. The integral clearer integrates input correlation signals i(n) and q(n) and outputs integral results I(n) and Q(n) as the observations of the navigation filter. The navigation filter uses the maximum likelihood estimation method to process the input information and outputs the estimated speed and receiver clock error. The line-of-sight projection converts the receiver speed and clock error estimates into frequency estimates using the projection relationship of the speed error and frequency error between the receiver and the satellite. Then, the frequency estimation value obtained by the VFLL and the phase estimation value output by the PLL are simultaneously input to loop filter 2, and the output result enters the NCO to obtain a local signal similar to the receiver input intermediate frequency signal. Figure 2 shows the structure of the traditional FAPLL. In the figure, the input signal in the FLL and PLL is multiplied by the local replica signal, and the output is a real baseband signal containing only data bits. To further increase the signal-to-noise ratio, the coherent integration results I(n) and Q(n) are output by the integration cleaner, whose coherent integration time is $T_{coh}$. The phase detector/frequency detector uses the coherent integration result to estimate the phase or frequency difference. The error result is filtered by the loop filter and then fused to adjust the numerical control oscillator so that the frequency or phase of the output signal is consistent with that of the received signal. Because the FLL in the combined loop can track the higher-order components corresponding to high dynamic stress in the input signal change, the tracking and locking of the signal by the phase-locked loop in the combined loop will become relatively easy, thus improving the tracking ability of the loop for the high dynamic signal [25].

Based on the theory of parameter estimation, the maximum likelihood estimation method is applied to VFAPLL (MLE) loop to improve the robustness of the loop in high dynamic environment. The loop is modeled using the Jerk model of high-speed and highly maneuverable targets. In the model, the state vector is composed of the Doppler frequency $f_0$, the first-order rate of change of the Doppler frequency $f_1$, and the second-order rate of change of the Doppler frequency $f_2$. The state equation is as follows [26]:

$${\mathit{X}}_k=\left[\begin{array}{ccc}1& T& \frac{1}{2}{T}^2\\ 0& 1& T\\ 0& 0& 1\end{array}\right]{\mathit{X}}_{k-1}+\left[\begin{array}{c}\frac{1}{2}{T}^2\\ T\\ 1\end{array}\right]{\mathit{W}}_{k-1}$$

(1)

where ${\mathit{X}}_k={\left[ f_0 f_1 f_2\right]}^T$, ${\mathit{W}}_{k-1}$ is the state disturbance noise, and its covariance matrix $\mathit{Q}\left(k\right)=E\left[{\mathit{W}}_k{\mathit{W}}_k^T\right]$. The expression of $\mathit{Q}\left(k\right)$ is shown in Formula (2). The jerk model adds another dimension to the acceleration model; that is, it can estimate the jerk, so a more accurate estimation of the acceleration parameters can be obtained. Consequently, the jerk model is very suitable for modeling highly dynamic carriers.

$$\mathit{Q}\left(k\right)=\left[\begin{array}{c} \frac{T^5}{\mathit{20}} \frac{T^4}{\mathit{8}} \frac{T^3}{\mathit{6}}\\ \frac{T^4}{\mathit{8}} \frac{T^3}{\mathit{3}} \frac{T^2}{\mathit{2}}\\ \frac{T^3}{\mathit{6}} \frac{T^2}{\mathit{2}} T\end{array}\right]{\mathit{N}}_{\mathit{y}}T$$

(2)

where ${\mathit{N}}_{\mathit{y}}$ represents the random jitter of the third derivative of the carrier frequency. After coherent integration processing by the integral clearer, the output integral signal with phase and frequency deviation is [17]:

$$g\left(k\right)=\left[\begin{array}{c}{g}_I\left(k\right)\\ g_Q\left(k\right)\end{array}\right]=\left[\begin{array}{c}AD\left(k\right)R\left(\Delta \tau \right)\frac{\sin \left(f_e{T}_{coh}\right)}{f_e{T}_{coh}}cos\theta \left(k\right)\\ AD\left(k\right)R\left(\Delta \tau \right)\frac{\sin \left(f_e{T}_{coh}\right)}{f_e{T}_{coh}}sin\theta \left(k\right)\end{array}\right]+\left[\begin{array}{c}{n}_I\left(k\right)\\ n_Q\left(k\right)\end{array}\right]$$

(3)

where A is the signal amplitude, $D\left(k\right)$ represents the data code, and $R\left(\Delta \tau \right)$ represents the C/A code autocorrelation function with a maximum value of 1. As the code delay of the instant branch is within one-half chip when the tracking is initialized, R(Δτ) can be normalized to 1. Variable $\Delta \tau$ is the phase difference between the copied C/A code and the received C/A code, $n_I\left(k\right)$, $n_Q\left(k\right)$ are zero-mean white Gaussian noise, and the $\sin \mathrm{c}\left(f_e{T}_{coh}\right)$ term is considered to be approximately equal to 1.

Therefore, the observation equation of the loop is:

$$Z\left(k\right)=\left[\begin{array}{c}Acos({\theta }_{me})\\ Asin({\theta }_{me})\end{array}\right]+B\left(k\right)$$

(4)

where ${\theta }_{me}\left(k\right)$ is the average residual phase over time $T$.

$$\begin{gathered} {\theta }_{me}\left(k\right)=\frac{1}{T}{\int }_{\left(k-1\right)T}^{kT}\left[{\theta }_c\left(k\right)+2\pi f_c\left(k\right)T+2\pi {f_c}^{\prime }\left(k\right)\frac{T^2}{2}+2\pi {f_c}^{″}\left(k\right)\frac{T^3}{6}\right]d\tau \\ -\frac{1}{T}{\int }_{\left(k-1\right)T}^{kT}\left[{\theta }_{nco}\left(k\right)+2\pi f_{nco}\left(k\right)T+\Delta {\theta }_{nco}\right]d\tau \\ ={\theta }_e\left(k\right)+\pi f_c\left(k\right)T+\pi {f_c}^{\prime }\left(k\right)\frac{T^2}{3}+\pi {f_c}^{″}\left(k\right)\frac{T^3}{12}-\pi f_{nco}\left(k\right)T \end{gathered}$$

(5)

The covariance of the observation noise matrix $V\left(k\right)$ is $R_k$, where $R_k=E\left[B_k B_k^T\right]$, such that

$${\mathit{R}}_{\mathit{k}}\mathit{=}\mathit{\left[}\begin{array}{c}{\mathit{\sigma }}_{\mathit{n}}^{\mathit{2}} \mathit{0}\\ \mathit{0} {\mathit{\sigma }}_{\mathit{n}}^{\mathit{2}}\end{array}\mathit{\right]}$$

(6)

3. VFAPLL Based on Maximum Likelihood Estimation

The IF signal model received by the receiver is expressed in Formula (7) [27]:

$$g\left(t\right)=h\left(t\right)+n\left(t\right)$$

(7)

where $g$(t) is the observed signal, $h$(t) is the real signal, and n(t) is the noise. Assuming the signal amplitude, code phase delay, Doppler shift, carrier phase, and other parameters change very slowly during the observation time interval so that they can be regarded as unknown constants during this time interval, the digital model of the signal is as follows:

$$g\left(k\right)=A·C\left(k-\tau \right)\cos \left(2\pi T\left(f_{IF}+f_d\right)k+\phi \right)+n\left(k\right)$$

(8)

where for the meaning of the parameters A, C(k), τ, T, f, φ, please see “List of symbols”. ($A$, $\tau$, $f_d$, $\phi$) are unknown signal parameters.

The maximum likelihood estimation method is to find a value of $\theta$ suitable for maximizing the possibility of obtaining the observed data, that is, if there is a sample of random vector $X$ with observation quantity $x$, find a solution of $\theta$ such that the corresponding joint probability density function or probability function $p\left(\mathit{x};\theta \right)$ reach their maximum.

For probability $\mathrm{P}\left(x|\theta \right)$, there are two cases: if x is a variable, $\theta$ is invariant, and this function is termed the probability function and describes the probability for different sample points x; on the contrary, if x is invariant, $\theta$ is a variable, and this function is termed the likelihood function and describes the probability of having sample point x for different model parameters $\theta$.

If corresponding to the paper, $\mathrm{P}\left(r|A,\tau ,f_d,\phi \right)$ represents that for a set of given sample values [$g\left(0\right),g\left(1\right)$,$\dots ,g\left(M\right)$], find a set of values of the signal parameters $\left(A,\tau ,f_d,\phi \right)$, so that the probability of the sample values is maximum, that is, find the optimal estimation of the signal parameters $\left(A,\tau ,f_d,\phi \right)$. If the elements of $g$ are all independent of each other, then the joint distribution

$$\mathrm{P}\left(r|A,\tau ,f_d,\phi \right)={\displaystyle \prod }_{i=1}^{M}p\left(r_i|A,\tau ,f_d, \phi \right)$$

(9)

is the probability density in Formula (9) in the paper.

The product operation entails large amounts of computation, so it is common to take logarithms of the likelihood function, as shown below:

$$\ln [{\displaystyle \prod }_{i=1}^{M}p\left(r_i|A, \tau , f_d, \phi \right)]={\displaystyle \sum }_{i=1}^M\ln [p\left(r_i|A, \tau , f_d, \phi \right)]$$

(10)

where $\ln [\prod _{i=1}^{M}p\left(r_i|A, \tau , f_d, \phi \right)]$ is the Log—likelihood function, corresponding to $U\left(A,\tau ,f_d,\phi |r\right)$ in this paper.

The joint probability density of the received signal at $N$ sampling points is [27]:

$$p\left(r|A,\tau ,f_d,\phi \right)=\frac{1}{{\left(\pi N_0\right)}^N}\exp (-\frac{1}{N_0}{\left(g-\widehat{g}\right)}^{T}W\left(g-\widehat{g}\right))$$

(11)

where $\widehat{g}$ denotes the estimated value of $g$ obtained from Formula (8) in the absence of noise, and $W$ denotes the weighting matrix of $g$. Assuming that the observations and residuals are causally reversible, and the residuals are independent of each other, the log-likelihood cost function of the above joint probability density function can be defined as follows [27]:

$$U\left(A,\tau ,f_d,\phi |g\right)=-\frac{1}{N_0}{\displaystyle \sum }_{k=0}^{N-1}{w}_k{\left|g\left(k\right)-\widehat{g}\left(k\right)\right|}^2$$

(12)

The maximum likelihood estimate of the signal parameters can be obtained by maximizing the log-likelihood cost function $U$ as follows [27]:

$$\frac{\partial U\left(\theta |g\right)}{\partial \theta }=0$$

(13)

$$\begin{gathered} L\left(\tau ,f_d|g\right)=\frac{1}{N_0}{\left\{{\displaystyle \sum }_{k=0}^{N-1}g\left(k\right)C\left(k-\tau \right)\cos \left(2\pi T\left(f_{IF}+f_d\right)k\right)\right\}}^2 \\ +\frac{1}{N_0}{\left\{{\displaystyle \sum }_{k=0}^{N-1}g\left(k\right)C\left(k-\tau \right)\sin \left(2\pi T\left(f_{IF}+f_d\right)k\right)\right\}}^2 \end{gathered}$$

(14)

For the derivation of formula (14), please see the Appendix A. The Doppler frequency $f_d$ is:

$$f_d\left(k\right)=f_0+f_{1}k+\frac{f_2}{2}{k}^2$$

(15)

where $f_0$, $f_1$, and $f_2$ represent the Doppler frequency, first-order rate of change of Doppler frequency, and second-order rate of change of Doppler frequency, respectively. Then, the maximum likelihood cost function is as follows:

$$\begin{gathered} L(f_{0 }, f_1, f_2|g)=\frac{1}{N_0}{\left\{{\displaystyle \sum }_{k=0}^{N-1}g\left(k\right)·C\left(k-\tau \right)\cos {\displaystyle (}2\pi T\left(f_{IF}k+f_{0}k+\frac{f_1}{2} k^2+\frac{f_2}{6}{k}^3\right){\displaystyle )}\right\}}^2 \\ +\frac{1}{N_0}{\left\{{\displaystyle \sum }_{k=0}^{N-1}g\left(k\right)·C\left(k-\tau \right)\sin {\displaystyle (}2\pi T\left(f_{IF}k+f_{0}k+\frac{f_1}{2} k^2+\frac{f_2}{6}{k}^3\right){\displaystyle )}\right\}}^2 \end{gathered}$$

(16)

The so-called maximum likelihood estimation of signal frequency performs a three-dimensional search on $f_0,f_1,f_2$ to find a set of frequency estimates when the maximum likelihood function reaches the peak value. However, the three-dimensional search of frequency requires considerable computation. As shown in expression (19), the high dimensionality and nonlinearity of the cost function will entail a large amount of computation for the MLE algorithm, which is also the main disadvantage of this algorithm. The peak value of the cost function can be obtained through iterative calculation using the gradient method.

For optimizing the maximum likelihood estimation method, we use the Levenberg-Marquardt algorithm. The LM algorithm belongs to optimization algorithm as well as an iterative maximum likelihood estimation method, and it is the most widely used nonlinear least squares algorithm and is an algorithm for finding the extreme value by using gradient, which belongs to the “mountain climbing” method. It has the advantages of both the gradient method and the Newton method. When λ is very small, the step size is equal to the Newton method step size, and when λ is large, the step size is approximately equal to the step size of the gradient descent method.

The algorithm requires calculating the gradient matrix and Hessian matrix of the log-maximum likelihood function and using the gradient to find the extreme value. The update equation of the algorithm is as follows [27]:

$${\widehat{\theta }}_{ML}^{i+1}={\widehat{\theta }}_{ML}^i-{\left(H_i+d_i\right)}^{-1}{G}_{i } i=0, 1, \dots$$

(17)

where ${\theta }_{ML}$ represents the state vector, which is equal to ${\left[g,f_d\right]}^T$. $G_i$ represents the gradient matrix, $H_i$ is the Hessian matrix, and $i$ represents the number of iterations.

$$G_i={\left[\frac{\partial L\left(\theta |g_N\right)}{\partial \theta }\right]}_{\theta ={\widehat{\theta }}^i}$$

(18)

$$H_i={\left[\frac{{\partial }^{2}L\left(\theta |g_N\right)}{\partial {\theta }^2}\right]}_{\theta ={\widehat{\theta }}^i}$$

(19)

where

$$\frac{\partial L}{\partial \theta }=\mathit{\left[}\begin{array}{c}\frac{\partial L}{\partial \tau }\\ \frac{\partial L}{\partial f_d}\end{array}\mathit{\right]} , \frac{{\partial }^{2}L}{\partial {\theta }^2}=\mathit{\left[}\begin{array}{cc}\frac{{\partial }^{2}L}{\partial {\tau }^2}& \frac{{\partial }^{2}L}{\partial \tau \partial f_d}\\ \frac{{\partial }^{2}L}{\partial f_d\partial \tau }& \frac{{\partial }^{2}L}{\partial f_d^2}\end{array}\mathit{\right]}$$

(20)

The value of the diagonal matrix $d_i$ must ensure that $H_i+d_i$ is always positive such that ${(H_i+d_i)}^{-1}$ has a solution, the ${(H_i+d_i)}^{-1}{G}_i$ term is the correction term for each iteration, and the function is equivalent to the discriminator in the traditional loop. In addition, the second-order derivative matrix—the Hessian matrix—is of great significance because it can determine whether the obtained solution is the minimum or maximum value. If the Hessian matrix is negative at the point corresponding to the solution, then this solution is the maximum. The Hessian matrix also defines the covariance matrix of the parameter estimates, and its relationship with the covariance matrix $V\left(\theta \right)$ is $V\left(\theta \right)={\left\{-E\left[H\left(\theta \right)\right]\right\}}^{-1}$.

Assuming that the C/A code is stripped from the received signal, the expression of the received signal is changed from (8) to the following form:

$$\tilde{g}\left(k\right)=A\cos \left(2\pi T\left(f_{IF}+f_d\right)k+\phi \right)+\tilde{n}\left(k\right)$$

(21)

where $\tilde{g}\left(k\right)=g\left(k\right)C\left(k-\widehat{\tau }\right) \mathrm{and} \tilde{n}\left(k\right)=n\left(k\right)C\left(k-\widehat{\tau }\right)$ represent the received signal and noise after stripping the C/A code, respectively, and the loop can be simplified into a carrier tracking loop. At this time, the gradient matrix and Hessian matrix in the LM algorithm are as follows [27]:

$$\frac{\partial L}{\partial f_d}=-\frac{4\pi T}{N_0}\left\{{\displaystyle \sum }_k\tilde{g}·\cos \right\}\times \left\{{\displaystyle \sum }_k\tilde{g}·k·\sin \right\}+\frac{4\pi T}{N_0}\left\{{\displaystyle \sum }_k\tilde{g}·\sin \right\}\times \left\{{\displaystyle \sum }_k\tilde{g}·k·\cos \right\}$$

(22)

$$\begin{gathered} \frac{{\partial }^{2}L}{\partial {f_d}^2}=\frac{8{\pi }^2{T}^2}{N_0}{\left\{{\displaystyle \sum }_k\tilde{g}·k·\sin \right\}}^2+\frac{8{\pi }^2{T}^2}{N_0}{\left\{{\displaystyle \sum }_k\tilde{g}·k·\cos \right\}}^2-\frac{8{\pi }^2{T}^2}{N_0}\left\{{\displaystyle \sum }_k\tilde{g}·\cos \right\}\times \left\{{\displaystyle \sum }_k\tilde{g}·k^2·\cos \right\} \\ -\frac{8{\pi }^2{T}^2}{N_0}\left\{{\displaystyle \sum }_k\tilde{g}·\sin \right\}\times \left\{{\displaystyle \sum }_k\tilde{g}·k^2·\sin \right\} \end{gathered}$$

(23)

where $\sin =\sin \left(2\pi T\left(f_{IF}+f_d\right)k\right)=\sin {\displaystyle (}2\pi T\left(f_{IF}k+f_{0}k+\frac{f_1}{2} k^2+\frac{f_2}{6}{k}^3\right){\displaystyle )},$ $\cos =\cos \left(2\pi T\left(f_{IF}+f_d\right)k\right)=\cos {\displaystyle (}2\pi T\left(f_{IF}k+f_{0}k+\frac{f_1}{2} k^2+\frac{f_2}{6}{k}^3\right){\displaystyle )}$, $\tilde{g}=\tilde{g}\left(k\right), \sum _k(·)=\sum _{k=0}^{N-1}(·)$. This process eliminates the related computation of τ in the original gradient matrix and Hessian matrix, thus reducing the amount of computation.

The purpose of the iterative process is to use an efficient method to gradually find the maximum value. This process can be divided into three steps. First, a set of initial values is set, and there is a distance between this set of values and the peak value. Second, it approaches the peak in some way from the starting value. Third, some method is needed to test whether the obtained value is a peak (i.e., the maximum likelihood solution). This corresponds to testing a new set of estimated values to assess whether the maximum value has been reached. If the result of the judgment has reached the maximum value, then the search can be stopped, and the set of values can be used as the maximum likelihood estimation solution. However, if the judgment result has not reached the maximum value, then the second and third steps need to be repeated until the maximum likelihood solution is found.

The following is the specific iterative process of the LM algorithm [27]:

(1)

First, set the initial values of $\widehat{\theta }$ and matrix $d$. When selecting initial values, we can consider the range of variables in the model and the assumed range of real parameters and then randomly select a set of initial values, or we can obtain them based on prior information. In this paper, we select a set of initial values based on previous research, where matrix $d$ should be initialized as a diagonal matrix.

(2)

After obtaining the initial value, the maximum likelihood solution is gradually found by iteration. In this paper, the gradient is used to find the maximum value, and this algorithm uses a climbing method. Calculate the gradient matrix and the Hessian matrix based on the initial value, and then use Formula (l) to update ${\widehat{\theta }}^i$.

(3)

Calculate the maximum likelihood estimation cost function $L({\widehat{\theta }}^i|r_N)$ and compare it with the previous cost function. If it is greater than or equal to the previous one, then accept this iteration and go to the next step. In contrast, if the cost function is less than the previous one, then increase $d_i$ and repeat the calculation. The iteration is terminated when the change in $L({\widehat{\theta }}^i|r_N)$ is less than the tolerance level $\delta$; in most cases, the tolerance level is set to approximately 0.01 or lower. The lower the tolerance level is, the more accurate the result, and the closer the obtained solution is to the true value.

(4)

Check for convergence. If convergence is not achieved, repeat step 2 to calculate ${\widehat{\theta }}^{i+1}$ until convergence.

4. Simulation and Verification

The receiver adopts a combination of vector and scalar tracking loops. To compare the VFAPLL (MLE), VFAPLL (EKF) and FAPLL (MLE) under the same conditions, parameters such as their loop filter bandwidth and the predetection integral were selected to be consistent; that is, the predetection integration time was 1 ms, $B_{PLL}=15$ Hz, and $B_{FLL}=25$ Hz. The MATLAB software was used to simulate and generate the GPS digital intermediate frequency signal using a signal frequency $f_{IF}=4.092$ MHz, sampling frequency $f_s=16.368$ MHz, ratio between carrier and noise CNR =35 dB/Hz, and a 1 ms C/A code period.

During 0–25s, the carrier moved in a uniform straight line. Then, the number of channels with signals in the system was 12. During 26–59s, the carrier continued to move in a uniform straight line, but during this period, the number of signal channels decreased to three or less, and this weak signal state lasted until 120s. The carrier periodically performed the high dynamic motion as shown in Figure 2 from 60s to 85s, during which the carrier was in an extremely harsh high dynamic weak signal environment. After 85s, the carrier returned to uniform linear motion, and after 120s, the number of channels with signal returned to 12 [28].

In this paper, the highly dynamic model of the Jet Propulsion Laboratory (JPL) in the United States was used to simulate the highly dynamic part in the motion state of the carrier, as shown in Figure 3a,b. During 0–3 s, the receiver performed uniform acceleration movement at −25 g, and the jerk was zero. At 3 s, the jerk increased rapidly to 100g/s and lasted for 0.5s, meanwhile, the acceleration increased from −25 to 25 g, and the acceleration motion was then made uniform at 25 g. When it reached 5.5 s, the jerk sharply changed from 0 to −100 g/s, and lasted for 0.5 s, whereas the acceleration decreased from 25 to −25g, and the acceleration motion was then made uniform at −25 g [28].

The simulation tests the tracking performance of the VFAPLL and FAPLL based on the maximum likelihood estimation algorithm and VFAPLL based on the extended Kalman filter. Figure 4 and Figure 5 show the speed and position errors generated by the three loops in the movement process described above.

Figure 4 and Figure 5 show that the tracking performance of the three loops was good during the initial period of 0–24 s, but during the period of 25–120 s, the loops were in a weak signal environment, and the position and speed errors of the three loops significantly increased: at the same time, there was a gap in tracking performance. The VFAPLL (MLE) had the best performance both in a weak signal environment and in a harsh, highly dynamic weak signal environment during 60–85 s; however, the FAPLL (MLE) could not adapt to the high dynamic stress due to the limited frequency discrimination range of the FLL and began to lose lock after 63 s. After 120 s, the number of channels with signals was restored to 12; accordingly, the tracking errors of the VFAPLL (MLE) and VFAPLL (EKF) were greatly reduced and gradually returned to the state during 0–24 s. During the operation of the above loops, the maximum speed error difference between the VFAPLL (MLE) and VFAPLL (EKF) was 2.69 m/s, and the maximum position error difference was 57 m.

The frequency tracking accuracy of the loop was analyzed through the relationship between the root mean square ($RMS)$ frequency error and the carrier-to-noise ratio (CNR) of the loop, and the frequency tracking accuracy of VFAPLL (MLE), VFAPLL (EKF) and FAPLL (MLE) was compared. The definition of $RMSE$ was as follows [29]:

$$RMSE=\sqrt{\frac{1}{M}{\displaystyle \sum }_{i=1}^M{( X_n\left(k\right)-{\widehat{X}}_n^i(k|k))}^2}$$

(24)

where $X_n\left(k\right)$ and ${\widehat{X}}_n^i(k|k)$ represent the real value and filtered estimated value at the same time, respectively. Figure 6 shows the relationship between the RMS frequency error and the carrier-to-noise ratio of the three loops of the VFAPLL (MLE), VFAPLL (EKF) and FAPLL (MLE). Figure 7 shows the loss-of-lock probability of these three loops to the signal Doppler frequency under different carrier-to-noise ratios.

Figure 6 shows that the VFAPLL (MLE) loop has significantly improved frequency tracking accuracy compared with the VFAPLL (EKF) and FAPLL (MLE), whether in the low CNR range of 20–25 dB/Hz or in the extremely low CNR range of 15–20 dB/Hz. When the carrier -to-noise ratio of the input signal is at a very low value of 15 dB/Hz, the tracking accuracy of VFAPLL (MLE) was improved by 8 and 18 dB compared with the other two loops respectively. When C/N0 = 21.5, the tracking accuracy of the VFAPLL (MLE) was 6 dB higher than that of the VFAPLL (EKF) and 15 dB higher than that of the FAPLL (MLE). When C/N0 = 22.5, the mean square error of the frequency estimation of the VFAPLL (MLE) was 28 Hz, which was 6 dB higher than that of the VFAPLL (EKF) and 15 dB higher than that of the FAPLL (MLE). However, with the increase in the carrier noise ratio C/N0, this improvement gradually decreases.

Figure 7 shows that the tracking threshold of the VFAPLL (MLE) was 22.5 dB/Hz, and the tracking thresholds of the VFAPLL (EKF) and FAPLL(MLE) were 24.5 and 25.5 dB/Hz, respectively, which indicates that compared with the other two loops, the VFAPLL (MLE) achieved 2 and 3 dB/Hz improvements, respectively, on the loss-of-lock threshold. In addition, it can be seen from the figure that when CNR ≤ 18.5dB/Hz, the probability of FAPLL (MLE) losing lock was approximately 1, while VFAPLL (MLE) still had the possibility of locking signal under extremely low CNR of 15 dB/Hz, and when CNR was 24.5 dB/Hz or more, the probability of VFAPLL (MLE) losing lock became very low.

In summary, the VFAPLL (MLE) has higher tracking accuracy than the VFAPLL (EKF) and FAPLL (MLE) in both low-SNR and very low-SNR conditions.

The following is an analysis of the robustness of the simulated loops. Based on the highly dynamic scenarios given in Figure 3a,b, we conducted a statistical comparison of the measurement accuracy of the Doppler frequency for the VFAPLL (MLE), VFAPLL (EKF), and FAPLL (MLE), as shown in Figure 8 and Figure 9.

Figure 8 and Figure 9 show a comparison of the average values of the root mean square (RMS) Doppler frequency errors of VFAPLL (MLE), VFAPLL (EKF), and FAPLL (MLE) under different acceleration and jerk scenarios, respectively. The results show that the VFAPLL (MLE) had a smaller average RMS Doppler frequency error than the other two loops in the abovementioned highly dynamic scenarios, indicating that the VFAPLL (MLE) had the highest robustness of the three loops. This was mainly due to the efficient use of iterations by the LM algorithm and the advantages of the vector tracking loop.

5. Convergence Analysis

The convergence of iteration mainly depends on three aspects: convergence standard, maximum number of iterations and initial value.

(1)

Convergence standard: The convergence criteria are defined as

$$\left|L\left({\widehat{\theta }}^{i+1}|r_N\right)-L\left({\widehat{\theta }}^i|r_N\right)\right|<0.01$$

(25)

That is, if the difference between the logarithmic likelihood function at time I + 1 and the logarithmic likelihood function at time i of the maximum likelihood estimation is less than the set value, it is considered unnecessary to iterate further, that is, the algorithm has converged.

(2)

Maximum number of iterations: The maximum number of times that the iteration process terminates, given in advance, automatically terminates if it does not converge. This is set to 10 times for this article.

(3)

Initial value: set with predictor.

In order to evaluate the convergence speed and performance of the algorithm, the following two kinds of square root errors are introduced in this paper [29],

$${\mathrm{RMSE}}_{time}\left(k\right)=\sqrt{\frac{1}{M}{\displaystyle \sum }_{j=1}^M{[{\widehat{\theta }}_k^j-{\theta }_k^j]}^2}$$

(26)

$${\mathrm{RMSE}}_{total}\left(j\right)=\sqrt{\frac{1}{N}{\displaystyle \sum }_{k=1}^N{[{\widehat{\theta }}_k^j-{\theta }_k^j]}^2}$$

(27)

where, ${\widehat{\theta }}_j\left(k\right)$ represents the estimated value of $\theta$ in the jth Monte Carlo simulation at time $k$, $N$ is the length of discrete time, and $M$ is the number of Monte Carlo simulations. $RMS{E}_{time}\left(k\right)$ represents the average convergence speed of algorithms in simulations, $RMS{E}_{total}$ indicates the performance of the algorithm in a simulation, the smaller the value, the better the performance of the algorithm.

The Figure 10 compares the convergence performance of VFAPLL (MLE) and VFAPLL (EKF) when CNR = 23 dB/Hz. It is known from the figure that the average convergence time of VFAPLL (MLE) loop is shorter than that of VFAPLL (EKF).

Table 1. The average consumption time and Average value of ${\mathrm{RMSE}}_{total}$ of VFAPLL (MLE) when CNR is 15, 19, 20, 21, 23, and 25 dB/Hz.

Carrier-to-Noise Ratio (C/N0, dB/Hz)	Average Consumption Time of VFAPLL (MLE) (s)	Average Value of $\mathbf{R}\mathbf{M}\mathbf{S}{\mathbf{E}}_{\mathit{t}\mathit{o}\mathit{t}\mathit{a}\mathit{l}}$ of VFAPLL (MLE) (s)
25	0.023	0.02
23	0.026	0.029
21	0.036	0.035
20	0.042	0.046
19	0.056	0.05
15	0.065	0.053

and

Table 2. The average consumption time and Average value of ${\mathrm{RMSE}}_{total}$ of VFAPLL (EKF) when CNR is 15, 19, 20, 21, 23, and 25 dB/Hz.

Carrier-to-Noise Ratio (C/N0, dB/Hz)	Average Consumption Time of VFAPLL (EKF) (s)	Average Value of $\mathbf{R}\mathbf{M}\mathbf{S}{\mathbf{E}}_{\mathit{t}\mathit{o}\mathit{t}\mathit{a}\mathit{l}}$ of VFAPLL (EKF) (s)
25	0.026	0.029
23	0.035	0.039
21	0.039	0.042
20	0.053	0.05
19	0.06	0.055
15	0.065	0.06

compare the performance of the VFAPLL (MLE) and VFAPLL (EKF) using ${\mathrm{RMSE}}_{total}$ and average consumption time as the index, as follows:

It can be seen from the table that the average time of running VFAPLL (MLE) was less than that of running VFAPLL (EKF), which proves that the overall complexity of VFAPLL (MLE) is lower than that of VFAPLL (EKF) (Including frequency discriminator). The comparison of ${\mathrm{RMSE}}_{total}$ shows that the loop tracking performance of VFAPLL (MLE) is better than that of VFAPLL (EKF).

6. Conclusions

In this paper, a VFLL-assisted PLL for global navigation satellite system (GNSS) receivers based on maximum likelihood estimation was proposed. The loop makes full use of the advantages of the VFLL and PLL and can realize a trade-off between dynamic robustness and tracking accuracy in the loop. Compared with the traditional VFAPLL, the loop performance is improved. The new architecture not only takes advantage of the vector tracking loop but also achieves a very low tracking threshold through ML estimation. In the process of ML estimation, the discriminator is removed, and the coherent integration result of the signal is directly used as the input of the navigation filter. The problem of the limited tracking range of the discriminator is avoided, and the tracking sensitivity of the loop is improved. The method proposed in this paper improves the performance of carrier tracking loop in high dynamic weak signal environments. The proposed loop structure improves the tracking accuracy and robustness of the loop.

Deficiencies in the Paper and Prospects for Future Work

The test signals used in the paper are all from the GNSS simulator or the high dynamic motion model of the JPL laboratory in the United States. If we can find in the future an external environment where the navigation signal strength and carrier motion meet the test conditions, it will be more helpful to combine theory with practice.

In addition, although the MLE algorithm shows better robustness and faster response speed in high dynamic scenarios, its computing load is greater than that of traditional algorithms. Because the computing load of the MLE (LM) algorithm increases almost linearly with the increase in the number of iterations, so we need to further optimize the MLE algorithm based on LM algorithm.