Performance assessment of a low-complexity autoregressive Kalman filter for GNSS carrier tracking using real scintillation time series

Ionospheric scintillation is one of the most challenging sources of errors in global navigation satellite systems (GNSS). It is an effect of space weather that introduces rapid amplitude and phase fluctuations to transionospheric signals and, as a result, it severely degrades the tracking performance of receivers, particularly carrier tracking. It can occur anywhere on the earth during intense solar activity, but the problem aggravates in equatorial and high-latitude regions, thus posing serious concerns to the widespread deployment of GNSS in those areas. One of the most promising approaches to address this problem is the use of Kalman filter-based techniques at the carrier tracking level, incorporating some a priori knowledge about the statistics of the scintillation to be dealt with. These techniques aim at dissociating the carrier phase dynamics of interest from phase scintillation by modeling the latter through some correlated Gaussian function, such as the case of autoregressive processes. However, besides the fact that the optimality of these techniques is still to be reached, their applicability for dealing with scintillation in real-world environments also remains to be confirmed. We carry out an extensive analysis and experimentation campaign on the suitability of these techniques by processing real data captures of scintillation at low and high latitudes. We first evaluate how well phase scintillation can be modeled through an autoregressive process. Then, we propose a novel adaptive, low-complexity autoregressive Kalman filter intended to facilitate the implementation of the approach in practice. Last, we provide an analysis of the operational region of the proposed technique and the limits at which a performance gain over conventional tracking architectures is obtained. The results validate the excellence of the proposed approach for GNSS carrier tracking under scintillation conditions.


Introduction
It is well known that global navigation satellite systems (GNSS) have become the key technology for positioning and navigation purposes using satellite ranging signals. Their maturity, widespread deployment and accuracy in open-sky environments are leading to a growing demand for extending GNSS beyond the limits of its original designs. This entails moving into the arena of high-precision, safety-critical applications such as aviation, maritime navigation and autonomous vehicle driving, just to mention a few. One key enabling feature is that GNSS carrier phase measurements do provide ultra-precise positioning information, and thus, their exploitation becomes essential. However, the use of GNSS in environments with many propagation impairments, different from those for which GNSS was initially conceived, hinders the determination of reliable carrier phase estimates and thus poses new technological challenges faced by nextgeneration GNSS receivers. Some examples are weak signal reception and multipath. But in particular, one of the most challenging sources of error to deal with in GNSS is the so-called ionospheric scintillation effect (Lee et al. 2017).
The ionosphere is the upper earth's atmosphere ionized by solar radiation, and it has a significant influence on transionospheric radio wave propagation. Ionospheric scintillation is a known effect of space weather whereby ionospheric electron density irregularities introduce rapid fluctuations onto GNSS signals when crossing the ionosphere in form of random amplitude fades and abrupt phase jumps (Kintner et al. 2007). These have a detrimental effect on the performance and stability of GNSS receivers and the proper determination of the synchronization parameters, particularly the carrier phase dynamics of interest. Ionospheric scintillations can cause cycle slips, severe carrier phase jitter, loss of signal even in line-of-sight (LoS) conditions, signal integrity issues and degrade navigation performance. Furthermore, although ionospheric scintillation can occur anywhere on earth during intense solar activity, it is more frequent in the earth's high-and low-latitude/equatorial regions (Kintner et al. 2009) and it certainly poses serious concerns to the widespread deployment of GNSS in emerging countries. Therefore, the design and validation of novel signal processing techniques for providing carrier phase estimates that are robust to ionospheric scintillation have become of paramount importance.
In order to exploit carrier phase measurements, GNSS receivers usually implement a carrier (along with code) tracking stage that keeps them synchronized with the satellites so that it accurately monitors any variation in the observed satellite-to-user dynamics. In most existing receivers, this is implemented using the well-known phase-lock loop (PLL) technique, a simple closed-loop architecture that compares the input carrier phase values to a local replica and drives the resulting error to zero by properly adjusting the phase of the local oscillator (Gupta 1975). Notwithstanding, this technique is known to experience serious troubles in the presence of unexpected propagation impairments, and it becomes the bottleneck in GNSS receivers under scintillation conditions (Lee et al. 2017). As a countermeasure, the use of advanced tracking techniques such as those based on the Kalman filter (KF) is gradually being introduced from the general framework of optimal minimum mean square estimation (MMSE). Some examples can be found in Macabiau et al. (2012), Won et al. (2012), Jiang et al. (2017), Susi et al. (2017) and Yang et al. (2017), where carrier phase dynamics is estimated using KFs that encompass kinematic models in their state-space formulation. As a matter of fact, the application of KF-based GNSS carrier tracking turns into an interesting approach for dealing with phase scintillation when noting that the random nature of the latter can be modeled through, for instance, a correlated Gaussian distribution (Humphreys et al. 2010). This opens the door to phase scintillation being accommodated by the KF in a natural manner, thus taking advantage of its optimality properties under Gaussian disturbances for linear models. A case of particular interest is the class of autoregressive (AR) processes, which are gaining relevance for the problem at hand and are claimed to provide potentially promising results (Nunes and Sousa 2014;Fohlmeister et al. 2018;Morton et al. 2020).
From the above observations, the use of a KF with a kinematic model augmented with a linear Gaussian model such as the AR becomes then the natural leap forward for optimally dealing with both phase dynamics and scintillation experienced by received GNSS signals. This gives rise to a hybrid autoregressive Kalman filter (KF-AR) technique, an innovative approach pioneered by the authors of the present work in Locubiche-Serra et al. (2016) and Vilà-Valls et al. (2013), whereby the above-mentioned kinematic and AR models are hybridized into one single state-space formulation. The underlying idea is to discern between phase dynamics and scintillation and thus provide decoupled estimates of both while making explicit that they must be determined in the presence of one another. The latter means, indeed, that the KF-AR can be employed to obtain either dynamics-robust phase scintillation estimates or scintillation-robust phase dynamics estimates. Even though remaining out of the scope of this work, the former does enable GNSS for ionosphere-related scientific applications such as ionospheric scintillation monitoring (ISM) itself. Our focus is placed on the latter, where robustness to scintillation is sought for positioning and navigation purposes. The final goal is to locally reconstruct the phase scintillation disturbance so that it can be removed from the input carrier phase and result in a clean phase depending on the user's dynamics only.
The existing contributions on this topic are focused on the design and characterization of the above techniques. However, an analysis of their efficiency based on real data experimentation is still lacking. In that sense, our objective is to validate the applicability and suitability of KF-AR techniques in realistic environments for scintillation-robust GNSS carrier tracking. To this end, an extensive experimentation campaign is carried out using real scintillation data captures at high and low latitudes to introduce real-worldrepresentative scintillation disturbances in the analysis, in contrast to many existing contributions that are restricted to synthetic data. In this sense, the work leads to a threefold contribution. First, we confirm the suitability of AR processes for modeling phase scintillation. Second, we propose a novel adaptive, low-complexity KF-AR technique designed with the practical implementation of the approach in the spotlight, aiming at facilitating its applicability in practice. This implementation is intended to circumvent the performance and computational limitations of other fixed and adaptive implementations previously proposed in the literature and attain the theoretical performance given by the well-known Bayesian Cramér-Rao bound (BCRB) irrespective of the input working conditions. Third, we provide some guidelines on the behavior of the proposed techniques as a function of the environmental working conditions. We also evaluate the performance limits at which the proposed techniques can operate in practice and outperform conventional PLLs.

Architecture of GNSS receivers
Apart from the antenna itself, the radiofrequency (RF) front-end is the first module of any GNSS receiver, where the analog processing and conditioning of the received signal for further digital processing is carried out. This involves low-noise amplification, filtering, frequency down-conversion to baseband and analog-to-digital conversion. Then, at its output, the digital processing of GNSS signals is usually performed through two sequential stages: the acquisition and tracking stages, followed by the position, velocity and time (PVT) module, which ultimately determines the user's position using the processed information.
The acquisition stage is in charge of detecting the satellites that are present in the received signal and finding a tentative estimate of their code delays (i.e., pseudoranges). As some Doppler shift is present owing to receiver oscillator inaccuracies and the relative motion between the satellite and the receiver, this search is performed in a two-dimensional grid, that is, in both time and frequency domains, by computing the so-called cross-ambiguity function. Therefore, when satellites are detected, the acquisition stage also provides coarse estimates of the code delay and Doppler shift at which the signal of interest is located. These outputs are then provided to the tracking stage.

Tracking stage
The tracking stage is a closed-loop architecture with a twofold objective. First, to refine the synchronization parameters estimated in the acquisition stage. Second, to accurately follow any possible variation of these parameters over time. The underlying idea is to find the code delay and carrier phase values to generate a local replica of the signal of interest that aligns with the received one. This is done by recursively comparing the local replica with the received signal and driving the resulting error to zero. To this end, a closed-loop system such as the one shown in Fig. 1 is employed, formed by two coupled architectures, one for code tracking and one for carrier tracking, usually referred to as delay-lock loop (DLL) and PLL architectures, respectively. Both follow the same basis. First, the input signal is correlated with the local replica using a set of correlators, namely early (E), prompt (P) and late (L), formed by integrate-and-dump (I&D) blocks performing coherent ( N c ) and non-coherent ( N i ) integrations (only coherent in the PLL case). The result is fed to the discriminators, which provide a function that is proportional to the local replica errors to be corrected. The loop filters (LF) aim at removing the noise at the discriminator output and providing a smoothed version of the corrections to be applied. This is then fed to the numerically controlled oscillators (NCO) that indicate such corrections to the code and carrier generators, thus closing the loop. In this work, we will focus on the PLL architecture, as detailed next.
For a given satellite, consider the signal x(n) at the input of the tracking stage with the following discretetime complex baseband equivalent signal (Seco-Granados et al. 2012), where n is the discrete-time variable, P is the signal power, d(n) is the navigation message, is the code delay, and e(n) is the contribution due to receiver noise. The term r(n) is the transmitted signal encompassing the receiver-to-satellite synchronization parameters, which can be expressed as (Seco-Granados et al. 2012), that is, delays the pseudorandom sequence c(n) linked to the satellite under study, and d (n) is the carrier phase of a complex exponential that contains the Doppler shift due to (2) r(n) = c(n − )e j d (n) Fig. 1 Block diagram of conventional GNSS code and carrier tracking architectures. Code tracking is used for fine code delay estimation. Carrier tracking is used for fine carrier phase estimation signal dynamics. The PLL aims at providing an estimate of the latter, denoted as ̂d (n) , and using it for generating the local replica in a way to compensate for the complex exponential in (2), where the phase error (n) ≐ d (n) −̂d(n) is driven to zero by iterating the loop. To that effect, the input signal is correlated with the local replica c(n −̂) through the P correlator in Fig. 1, denoted as y P (̂) , which is centered at ̂ , the delay estimated by the DLL. Note that the DLL employs two additional correlators, the E-L correlators, centered at ̂ plus some shift at each side, which facilitate the detection of the main correlation peak. In contrast, the PLL employs only the P one.
When feeding the result to the PLL discriminator, the objective is to extract the phase error (n) of the complex exponential in (3). For this purpose, several discriminators can be employed, highlighting the PLL-, Costas PLL-and frequency-lock loop (FLL)-type ones. Since we are interested in phase measurements, we focus on the PLL-type discriminators, specifically the well-known four-quadrant arctangent (ATAN2) discriminator as it duplicates the dynamic range of the Costas PLL and is robust to phase variations. In addition, it is the optimal maximum likelihood (ML) phase extractor at high carrier-to-noise ratio (C/N 0 ) (Proakis 2001), while it also fits very well into GPS signals synchronized at bit level, as well as into the data-less nature of pilot signals growingly used nowadays. Its output, denoted as Δ ATAN2 , is given by, where the terms Im y P (̂) and Re y P (̂) correspond to the imaginary and real parts of the prompt correlator, respectively.

Effect of ionospheric scintillation
When the signal received from a given satellite is affected by ionospheric scintillation, the effect can be modeled as a multiplicative channel introducing signal amplitude A s (n) and phase in such a way that the carrier phase perceived by the PLL becomes the combination of the signal dynamics plus the scintillation phase, that is, (n) ≐ d (n) + s (n).

Signal models
The term (n) is the phase to be dealt with at the carrier tracking level. For this purpose, the corresponding signal models must first be defined. The models adopted in this work are described next.

Signal model for carrier phase dynamics
For the problem of GNSS carrier tracking, we assume that the discrete-time evolution of the carrier phase due to signal dynamics can be approximated by a third-order kinematic model, which has been found to be a self-standing approach of interest in many practical cases (Bar-Shalom et al. 2004 where T s is the sampling time, and the terms ̇d (n) , ̈d (n) and ⃛ d (n) are the Doppler shift, Doppler rate and Doppler acceleration (i.e., first, second and third derivatives of the carrier phase), whose evolution is given by, which follows the structure in (7). This is then the signal model to be used later to formulate the problem of GNSS carrier tracking.

Signal model for phase scintillation
When approximating the discrete-time evolution of phase scintillation by a pth-order AR process, namely AR(p), the following signal model is adopted (Kay 1993), where m p m=1 is the set of p AR coefficients and s p (n) is the so-called AR driving noise or prediction error of an AR(p) process, usually modeled as s p (n) ∼ N 0, 2 This is then the signal model to be used later to formulate the problem of phase scintillation tracking. Additionally, the power spectral density (PSD) of an AR(p) process is given by (Kay 1993), with the discrete-time frequency variable. Expression (9) is used next to evaluate the feasibility of modeling phase scintillation through an AR process.

AR model fitting of phase scintillation
Once the signal models are defined, this section is intended to show the suitability of modeling scintillation phase through AR processes, a required step before formulating the KF-AR.

Scintillation time series
Scintillation is reported to be frequent in low-and highlatitude regions and in the dark times of the day and year (Aarons 1982 . This is henceforth denoted as KIR.2015.106.04.GPS22. In both cases, the setup for scintillation monitoring is formed by a high-grade antenna followed by a Septentrio PolaRxS receiver providing binary files with raw data and scintillation parameters at a rate of 50 Hz. These files are converted to raw American Standard Code for Information Interchange (ASCII) files containing the in-phase and quadrature (IQ) correlation samples, from which the scintillation power and phase time series are extracted following the data-detrending algorithm in Deshpande et al. (2012). We also use the algorithm to eliminate residual low-frequency systematic effects initially affecting the captures, such as those from the troposphere, satellite geometry and receiver oscillator instabilities, which introduce undesired fluctuations that add to the scintillation data and distort the measurements. The look of the detrended phase and power is shown in Figs. 2 and 3 for the Dakar and Kiruna captures, respectively. We will focus on the segments delimited within the vertical stripped lines as scintillation is more prominent there, with satellite elevation angles above 15° to avoid any propagation degradation effect owing to low satellite elevation. On the one hand, in order to characterize amplitude scintillation, S 4 is the index usually employed in the literature, defined as (Humphreys et al. 2010), with E[⋅] the expectation operator. On the other hand, phase scintillation is characterized through the standard deviation index s of the phase series. At low latitudes, phase scintillation is correlated with amplitude scintillation, with s = 0.35 rad and power drops around 5-10 dB with S 4 = 0.3 , being overall representative of mild scintillation. However, a different situation is observed at high latitudes, where even larger phase scintillations are perceived, with s = 0.65 rad, but no major power drops occur, with S 4 = 0.05 , thus being representative of strong scintillation events in terms of phase. The behavior in both cases is in line with what is reported in the literature for low-and highlatitude scintillation (Jiao et al. 2013;Lee et al. 2017).

Suitability of phase scintillation AR fitting
In order to evaluate how well the above live captures of phase scintillation can be modeled through an AR(p) process, the AR parameters that best fit the time series must be first found, namely the coefficients m p m=1 and the prediction error power 2 s p . By applying either the Yule-Walker equations or the least squares method (Kay 1993;Stoica and Moses 2005), the obtained results are shown in Table 1 for the Dakar and Kiruna data captures considering AR(p) models with p up to three. These values are used to determine the corresponding PSDs using (9), which are then compared to the periodogram of the measurements as a nonparametric estimator of the true PSDs. The results are shown in Figs. 4 and 5 for each capture. A very tight match is obtained, despite some discrepancies observed at low frequencies that are, though, mainly due to the low-pass filtering steps involved during the detrending process. Note that the interest lies in the high-frequency variations introduced by scintillation (Deshpande et al. 2012). In that sense, an AR process with any order p can be used in both captures, thus suggesting that a good fit can be obtained by just keeping the lowest order, namely AR(1). This result can already be anticipated by looking at the values in Table 1, where the prediction errors are very similar for all orders. Therefore, these results show the suitability of modeling phase scintillation events using an AR(1) process.

Kalman filter equations
As previously explained, there is a growing trend to replace conventional PLL architectures (in particular the loop filter) by a KF-based one. Therefore, the objective of this section is to briefly present the generic KF equations, also termed state-space and observation models, which will be used later to formulate the KF-AR with the kinematic and AR signal models introduced earlier to track carrier phase dynamics and phase scintillation.  The sequential estimation of the parameters of interest is governed by the following state transition equation (Yang et al. 2017), where (n) is the state vector containing such parameters of interest at time instant n. It propagates toward the next time instant n + 1 through the transition matrix F, plus the process noise (n) that accounts for the possible modeling discrepancy of the Kalman state-space model to the true one characterizing the input samples. It usually follows (n) ∼ N 0, 2 such that E H (i) (j) = 0 for i ≠ j . As observed in (11), the weight of (n) onto the different Kalman states is adjusted through a linear transformation given by matrix G that results in a zero-mean Gaussian process with the covariance matrix ≐ E H H . The KF aims at providing clean estimates of (n) based on the available noisy input measurements, which are modeled through the following measurement equation (Yang et al. 2017), where the scalar measurements z(n) can be understood as a linear transformation of the system states through the observation matrix H. The measurements are corrupted by the presence of an uncorrelated Gaussian noise

Hybrid autoregressive Kalman filter (KF-AR)
This section is intended to show the integration of the kinematic and AR signal models introduced earlier into the Kalman equations, leading to the hybrid autoregressive Kalman filter (KF-AR), and present low-complexity implementation proposed in this work.

Baseline technique: the KF-AR(1): state-space and observation models
Recalling the kinematic signal model in (6)-(7), the states of interest for a three-dimensional KF are d (n) , ̇d (n) and ̈d (n) . Hence, the state-space model for carrier dynamics tracking can be written in normalized matrix notation as, where by direct comparison with (11), the (3 × 1) state vector of the carrier phase dynamics and the transition matrix are identified, (11) (n + 1) = (n) + (n) The process noise is given by u(n) ≐ T 3 s ⃛ d (n) ∼ N 0, 2 u , which is intended to cover the higher-order terms that are missing in the truncated model in (6) (8) for an AR(1) process becomes readily the Kalman state transition equation for the problem of phase scintillation tracking, where the only state s (n) propagates through , and the process noise is now given by the AR driving noise s p (n) , which propagates entirely to the state. Then, by merging the kinematic and AR models, the KF-AR(1) state-space system model results into the following augmented state transition equation, where the state vector and transition matrix are identified as, The process noise becomes, with p = 0 referring to a dynamics-only, non-AR KF that will be used by the adaptive technique in the next section. The noise in (19)  When dealing with both phase dynamics and phase scintillation, we assume to have the following measurements available, where the augmented observation matrix is identified as ≐ 1 0 0 1 T . The measurement noise w(n) corresponds to the phase noise at the discriminator output. Assuming that the output of the prompt correlator is normalized to unit mean power, the variance of w(n) for the ATAN2 discriminator is a scalar given by (Del Peral-Rosado et al. 2010), where the second term inside the brackets at the right-hand side accounts for the squaring losses due to the nonlinear behavior of the discriminator at low C/N 0 (Curran et al. 2012).

Adaptive technique: the AHL-KF-AR(0, 1)
In order to deal with the time-varying nature of ionospheric scintillation, Locubiche-Serra et al. (2016) proposed a KF-AR implementation that consisted in three real-time adaptive algorithms on the baseline KF-AR state-space and observation models presented in the previous section. First, an estimator of the optimal AR model parameters. Second, an estimator of the optimal AR model order. Third, the use of adaptive hard-limited (AHL) Kalman gains to deal with the nonlinear C/N 0 drops introduced by scintillation. This KF-AR technique is referred to as AHL-KF-A2R(p). However, despite its agility and adaptability advantages, it is found to present two drawbacks. On the one hand, the simultaneous use of three adaptive mechanisms at each time instant incurs quite a high resource and computational burden. On the other hand, the available scintillation data are distorted by the presence of additive white Gaussian noise (AWGN) in the input measurements, which leads to a biased estimation of the optimal AR model parameters that eventually incurs some performance degradation, particularly when scintillation is not very strong. As a matter of fact, this latter adaptive implementation is also the major contributor to the computational cost enhancement, as it involves either the Yule-Walker or the least squares methods that require computing several autocorrelation functions and inverse matrices at each time instant. Therefore, in order to circumvent these limitations, we present the AHL-KF-AR(0, 1), a new implementation that uses a fixed set of AR model parameters, while partially preserving the other adaptive modules of the AHL-KF-A2R(p): the AHL implementation remains as is, whereas the AR model-order estimator will be employed as a scintillation detector switching between a kinematic-only KF when scintillation is absent, that is, KF-AR(0), and a KF-AR(1) when scintillation is present.
The proposed technique corresponds to the carrier tracking architecture depicted in Fig. 6. As can be observed, the carrier discriminator output is the magnitude to be processed by the KF, which corresponds to the innovation sequence i(n) that contains the error between the incoming phase and the Kalman predicted phase. Before the local replica is compensated with the estimated scintillation phase (i.e., output of I&D block), the architecture implements a second carrier discriminator that, since the carrier dynamics have already been compensated, outputs only the input scintillation phase. This is used to detect the presence of scintillation to consequently enable or disable the AR module of the KF and update the Kalman gains K(n) and the transition matrix accordingly. The state correction is carried out, and Fig. 6 Block diagram of AHL-KF-AR(0, 1) carrier tracking technique. The scintillation detector is used to set the AR model order to 0 or 1, thus disabling or enabling the AR module. The C/N0 estimator is used to deal with nonlinear signal power fades through the AHL implementation the prediction for the next time instant is computed. The latter is then employed to obtain the measurement prediction and update the local replica, which is correlated with the input signal, thus closing the loop.

Feature 1: optimal fixed AR parameters
In order to determine the optimal AR parameters for the AHL-KF-AR(0, 1), the best achievable performance for each of the considered datasets must first be known. This information is readily provided by the BCRB, which can be recursively computed as (Van Trees and Bell 2007), with B (n) the Bayesian information matrix. For the optimal AR parameters in Table 1, and considering a high C/N 0 of 45 dB Hz, the BCRB for the phase dynamics of interest is located at 3.6 × 10 −3 and 1 × 10 −2 rad 2 for the Dakar and Kiruna captures, respectively. We will take the smallest value to be our best possible target performance. That is to say, we take 3.6 × 10 −3 rad 2 as the MSE carrier phase performance to be attained by the AHL-KF-AR(0, 1).
In order to make the AHL-KF-AR(0, 1) robust, the optimal AR parameter determination is carried out with strong phase scintillation in the spotlight, such as the one from the Kiruna capture analyzed herein. That is, the objective is to design the parameters such that the resulting tracking scheme is able to cope with strong scintillation while attaining the aforementioned target MSE of 3.6 × 10 −3 rad 2 . For this purpose, note that the KF-AR is a closed-loop architecture that can thus be characterized by its equivalent loop bandwidth, which must be large enough for adequate signal tracking. As a matter of fact, such bandwidth is driven by the KF-AR process noise variance (Jwo 2001), given by 2 s p for the case of phase scintillation tracking. Therefore, as a first step we can simply take the value of 2 s p in Table 1 required to track the scintillation in Kiruna as reference, namely 2 s p = 3 × 10 −3 rad 2 . Then, it remains to determine the coefficient , which is sought as the value for which, with such 2 s p , the KF-AR attains the above target performance. By evaluating the BCRB for different candidate values of as shown in Fig. 7 it is found that, for 2 s p = 3 × 10 −3 rad 2 , the performance of 3.6 × 10 −3 rad 2 is attained for = 0.925 . Therefore, this is the set of AR parameters to be used by the AHL-KF-AR(0, 1).

Feature 2: real-time scintillation detection
When designing an AR(p) process, the model order p plays an important role for optimal AR fitting to a given set of measurements. In that sense, determining the optimal order of a statistical model is a well-known problem in the field of signal processing often referred to as model-order selection (Stoica and Selen 2004). Several criteria can be found in the literature in this regard, such as Akaike's information criterion (AIC) (Akaike 1974), the modified AIC, the generalized information criterion (GIC) and the minimum description length (MDL) (Djuric et al. 1999). The proposed technique implements the latter, as it is reported to be a consistent criterion. It computes the optimal p as, with N the number of samples in the data segment under analysis. In the proposed AHL-KF-AR(0, 1), this estimator is implemented in form of a sliding window that evaluates segments with 5 s of data as a function of time, and detects whether scintillation is present or not by determining ̂2 s p for p = 0 and p = 1 using = 0.925 and applying (23).

Feature 3: adaptive hard-limited measurement noise variance
The adaptive hard-limited (AHL) implementation is thought to deal with the nonlinear amplitude fades introduced by scintillation, whereby the KF-AR is affected by abnormal measurements that may compromise the linearity of the phase discriminator, thus degrading its performance and even driving the technique to lose lock. The AHL exploits the fact that the amplitude fades make themselves explicit at C/N 0 level. Interestingly, this information can be provided to the KF-AR through the measurement noise variance R(n) as (23) p MDL = arg min p N log ̂2 s p + p log (N) 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 s 2 (rad 2 ) 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1  Fig. 7 BCRB of KF-AR(1) as a function of AR model parameters. For 2 s p fixed to track strong scintillation, the objective is to find the that reaches the target BCRB seen in (21), so that the filter be dynamically self-adjusted to match the actual working conditions. Therefore, the AHL adapts R(n) based on actual C/N 0 measurements, subject to a hard-limiting threshold such that, with R (n) the measurement noise variance in (21) using Ĉ N 0 (n) , and where setting R(n) = ∞ drives the Kalman gains to zero whenever the C/N 0 drops below a given threshold . In this way, the KF-AR becomes totally isolated from the carrier discriminator output, and the loop relies only on the Kalman internal state-space model, thus protecting the technique from abnormal measurements until the C/N 0 recovers to nominal conditions.
To estimate the C/N 0 , the new AHL-KF-AR(0, 1) implements the narrow-wideband power ratio (NWPR) estimator, as it is a low-complexity algorithm of easy implementation that provides reliable results (Van Dierendonck 1996). On the other hand, taking into account that the C/N 0 tracking threshold is located around 20 dB-Hz (Kaplan and Hegarty 2006), the threshold is set to 25 dB-Hz so as to leave some margin for KF-AR operation at low C/N 0 .

Experimentation results and performance evaluation
In this section, we assess the performance of fixed and adaptive KF-AR techniques to validate the suitability and efficiency of the proposed approach for tracking carrier phase dynamics under scintillation conditions. We evaluate the adaptability features of the techniques under timevarying scintillation working conditions and the closeness to the expected BCRB in comparison with the performance provided by conventional PLLs. To this end, an experimentation campaign has been carried out using a MATLAB software GNSS carrier tracking simulator for GPS L1 C/A signals with a sampling time of T s = 20 ms. The campaign simulates the environmental conditions in terms of signal dynamics and receiver noise and adds the scintillation time series captured in Dakar and Kiruna on top of these. We consider a C/N 0 of 45 dB Hz in a static receiver. The only experienced dynamics are due to the satellite motion, thus assuming a Doppler shift of 10 Hz, Doppler rate of 1 Hz/s and maximum Doppler acceleration of ± 2 × 10 −4 Hz/s 2 . The latter leads to 2 v = 3.4 × 10 −17 rad 2 as the variance of a uniform distribution delimited by such Doppler acceleration values. The executions are 600 s long with a time-varying presence of scintillation. Only AWGN is present during the first 150 s, then scintillation appears between 150 and 450 s, then scintillation disappears again during the last 150 s. In that sense, two sets of results are provided, namely for low and high latitudes. The metrics employed for performance evaluation are the root mean square error (RMSE) of the estimated carrier phase dynamics, the number of cycle slips and the probability of loss of the tracking lock (LoL). The three metrics are computed by running the executions with 100 Monte Carlo realizations. Figures 8 and 9 show the RMSE performance of the previous KF-AR techniques and the new one proposed herein for the Dakar and Kiruna data captures, respectively. Additionally, Fig. 10 illustrates the error in cycles between the estimated and true carrier phase of interest for some of the Monte Carlo realizations. As can be observed, all techniques outperform the RMSE that would be obtained with a conventional PLL, as the latter estimates scintillation as part of the signal dynamics. The PLL presents a LoL of about 20%, and it suffers from some tens of cycle slip occurrences over time, especially in Kiruna, as shown in Fig. 10 top. More particularly, the mean number of cycle slip occurrences in 300 s is found to be 45, meaning that a cycle slip is expected to occur approximately every 6-7 s. On the contrary, KF-AR techniques are found to present no cycle slips, as shown in Fig. 10  LoL probability is below 1% as no Monte Carlo realization has lost lock, except for the AHL-KF-A2R(p), which shows a LoL probability of 5% in Kiruna. When scintillation is present, the fixed KF-AR(1) that employs the optimal AR parameters in Table 1 can fully reach the target performance given by the BCRB, thus confirming the validity of the modeling phase scintillation through an AR process. Note that the ideal situation would be to attain the BCRB obtained when scintillation is absent. However, the KF-AR approach already outperforms the conventional PLL by a factor of 6. The inconvenient of the fixed KF-AR(1) is, though, that it fails at providing optimal performance when pulled out of the test case for which it has been designed, particularly in terms of the transient experienced when scintillation disappears, requiring between 50 and 100 s to converge. This is the penalty incurred by forcing a fixed AR module in the KF in the absence of scintillation. This is overcome by the AHL-KF-A2R(p), which rapidly detects the absence of scintillation and disables the AR module, thus minimizing the transient period toward the new condition. In that sense, the major advantage of this technique is that it presents good agility and flexibility in front of time-varying conditions. Notwithstanding, the main inconvenience of the AHL-KF-A2R(p) is that, besides presenting a slightly larger LoL probability, it presents some RMSE gap that prevents attaining the BCRB by a factor between 2 and 3. This is owing to the online AR parameter estimator providing biased estimates due to the presence of AWGN. Furthermore, the computational burden of this technique is found to be 2.75 times that of the simplest KF-AR. This is where the new AHL-KF-AR(0, 1) technique proposed in this work comes into play, which shows very promising results. The technique commutes to the KF-AR(0) when under AWGN only, and to the fixed KF-AR(1) when in the presence of scintillation, thus attaining the corresponding BCRB values in both situations for both data captures. The success rate on the detection of scintillation is shown in the bottom plots of Figs. 8 and 9, with values above 90%, even though some transition time is seen to be particularly required when scintillation disappears. The agility of the AHL-KF-A2R(p) is preserved, while the required computational cost is reduced to 1.75 times that of the simplest KF-AR and the LoL probability is improved to below 1%, thus potentially becoming one robust approach that can be implemented in practice.

Experiment 2: KF-AR operational region
The objective of this section is twofold. First, to extend the analysis in the previous section and evaluate the KF-AR performance as a function of the input C/N 0 and the highorder dynamics that the KF can tolerate through 2 v . Second, to quantify the limits of the proposed AHL-KF-AR(0, 1) technique within which it provides a performance gain over the conventional PLL and out of which there is no point in using it over this latter.
By performing a dual-sweep analysis on both the above parameters, Fig. 11 shows a 3D plot of the AHL-KF-AR(0, 1) phase RMSE when considering low-latitude scintillation (Dakar). As can be observed, the RMSE increases when decreasing the C/N 0 or when enhancing the dynamics through the Kalman model discrepancy 2 v . The latter is explained by the fact that, for larger 2 v , the Kalman must increase its dependability on the input measurements rather than on its internal state-space model. As this occurs, it becomes more difficult for the Kalman to discern the scintillation phase from the dynamics phase of interest, thus eventually leading to a larger RMSE.
On the other hand, the gray plane in Fig. 11 stands for the best performance that could be expected from PLLs. The performance limits of the AHL-KF-AR(0, 1) are given by the boundary lines given where both KF-AR and PLL RMSE plots intersect. The region where the proposed technique outperforms the PLL is bounded by a C/N 0 in the region of 30 dB Hz, with 2 v restricted to 10 −20 rad 2 , whereas the technique supports up to 2 v = 10 −12 rad 2 for C/N 0 above 35 dB Hz. For larger 2 v the technique performs very similarly to the PLL, whereby phase scintillation is estimated as part of the carrier dynamics of interest. These values can be translated into more friendly magnitudes such as the maximum supported Doppler acceleration. For uniformly distributed process noise, both are related as: in Hz/s 2 . Therefore, in other words, the technique can tolerate a Doppler acceleration of ± 3 × 10 −6 Hz/s 2 at low C/N 0 , whereas it can support ± 3 × 10 −2 Hz/s 2 at high C/N 0 . Figure 12 shows the 3D plot of the proposed AHL-KF-AR(0, 1) phase RMSE when considering high-latitude scintillation (Kiruna). Here, it is remarkable to note that the outperformance region is larger than that at low Fig. 11. The technique outperforms the PLL for a lower-bounded C/N 0 of 25 dB Hz, that is, the value to which the AHL threshold has been set. At this point, 2 v is limited to around 10 −20 rad 2 , similarly to the case at low latitudes, corresponding to a maximum Doppler acceleration of ± 3 × 10 -6 Hz/s 2 . For C/N 0 of 30 dB Hz and above, the tolerated higher-order dynamics is found to be up to 2 v = 10 −10 rad 2 , corresponding to a maximum Doppler acceleration of ± 3 × 10 −1 Hz/s 2 . Beyond this point, the AHL-KF-AR(0, 1) provides similar performance to the PLL. The above results confirm the suitability of KF-AR techniques for precise localization and cruise navigation purposes, and autonomous vehicle driving applications with mild acceleration variations.

Conclusions
We have tackled the problem of robust GNSS carrier tracking in the presence of ionospheric scintillation utilizing hybrid autoregressive Kalman filters. Two main contributions have been presented. First, we proposed the AHL-KF-AR(0, 1) low-complexity technique to facilitate the implementation of the approach in practice. Second, we assessed the performance of these techniques using representative captures of real scintillation data at low and high latitudes, which are two different and complementary situations commonly addressed in studies of scintillation disturbances. From the work performed, three main conclusions can be drawn. First, the obtained results confirm the feasibility of modeling the phase scintillations observed by the user through AR processes. Second, the proposed AHL-KF-AR(0,1) low-complexity technique has shown to clearly outperform the conventional PLL architecture and other KF-AR  implementations while conserving the optimality properties of the Kalman filter. Third, the performance limits where the technique provides a performance gain over the PLL have been quantified, thus becoming valuable guidelines during the design of the tracking stage of GNSS receivers in real environments. Therefore, the work performed validates the excellence of the proposed approach for robust GNSS carrier tracking under scintillation conditions. Acknowledgements This work has been supported in part by the Spanish Ministry of Economy and Competitiveness under Grant TEC2017-89925-R, in part by the ICREA Academia program and in part by the FI-DGR PhD grant with file number 9015-376195/2014 from Generalitat de Catalunya. The authors would also like to thank the European Space Agency (ESA), and in particular Mr. Juan M. Parro, for providing the Monitor scintillation time series employed to perform the analyses addressed in this work.
Funding Open Access Funding provided by Universitat Autonoma de Barcelona.

Data availability
The scintillation time series analyzed in this work have been obtained from the Monitor network of the European Space Agency (ESA), a project aimed at monitoring ionospheric scintillation events from multiple ground stations located worldwide. Access to these time series is publicly available and can be granted by ESA upon request through the Monitor website (http:// monit or. estec. esa. int/).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.