Multi-frequency GNSS robust carrier tracking for ionospheric scintillation mitigation

Ionospheric scintillation is the physical phenomena affecting radio waves propagating from the space through the ionosphere to earth. The signal distortion induced by scintillation can pose a major threat to some GNSS application. Scintillation is one of the more challenging propagation scenarios, particularly affecting high-precision GNSS receivers which require high quality carrier phase measurements; and safety critical applications which have strict accuracy, availability and integrity requirements. Under ionospheric scintillation conditions, GNSS signals are affected by fast amplitude and phase variations, which can compromise the receiver synchronization. To take into account the underlying correlation among different frequency bands, we propose a new multivariate autoregressive model (MAR) for the multi-frequency ionospheric scintillation process. Multi-frequency GNSS observations and the scintillation MAR are modeled in state-space, allowing independent tracking of both line-of-sight phase variations and complex gain scintillation components. The resulting joint synchronization and scintillation mitigation problem is solved using a robust nonlinear Kalman filter, validated using real multi-frequency scintillation data with encouraging results.


Introduction
Global Navigation Satellite Systems (GNSS) is the technology of choice for most position-related applications when it is available (Dardari et al., 2015).A GNSS receiver relies on a constellation of satellites to estimate a set of range measures from which to compute its position.The main challenges of GNSS technology arise when operating in complex, harsh propagation scenarios which are naturally impaired by multipath, shadowing, high dynamics, or ionospheric scintillation.In the last decade, ushered by an ever increasing demand for availability, accuracy, and reliability, the mitigation of these challenges has steered intense research on advanced receiver design (Amin et al., 2016).
Among these challenges, in this paper we focus on ionospheric scintillation.Because such disturbance is not related to the local environment, as in the case of multipath or shadowing, it can degrade receiver performance even under ideal open-sky conditions.Ionospheric scintillation is caused by a disturbance in the portion of the ionosphere through which the GNSS signals propagates, an effect particularly relevant in polar and equatorial regions.It can produce rapidly varying constructive and destructive interference between multiple scattered signals.This turbulent behavior can render the carrier phase difficult to track by the receiver.One of the major challenges of severe scintillation propagation conditions is known as canonical fading, which is characterized by a combination of strong fading and rapid phase changes in a simultaneous and random manner (Kintner et al., 2009;Humphreys et al., 2009).From a synchronization standpoint, measuring the signal phase during such a fade is very challenging because rate of change of phase changes is at its highest when the signal amplitude is near its lowest (Humphreys et al., 2005).This may lead the receiver to momentarily lose carrier synchronization, resulting in either cycle-slips or loss of lock.Unfortunately, the accuracy improvement of modern high-precision receivers is based on the use of carrier phase-based positioning techniques, and so it is extremely sensitive to scintillation perturbations (Jacobsen & Dähnn, 2014;Prikryl et al., 2014;Banville & Langley, 2013, Banville et al., 2010).
Traditional synchronization relies on well-known closedloop DLL/PLL architectures (Kaplan, 2006).However, receiver architectures based on DLL/PLL are not always able to provide continuous and reliable phase tracking (Humphreys et al., 2005).Most of state-of-the-art (SoTA) techniques propose methods to robustify DLL/PLL schemes via fine tuning and augmentation (Yu et al., 2006;López-Salcedo et al., 2014), either by using adaptive bandwidth schemes (Skone et al., 2005) or combined architectures (Chiou, 2010;Won, 2014).It has already been shown that Kalman filter (KF)-based solutions provide better performances and robustness in these situations (Humphreys et al., 2005;Macabiau et al., 2012), being nowadays the performance benchmark for the development of new methodologies (Vilà-Valls et al., 2017b).The natural extension is to use adaptive KF-based solutions (Zhang et al., 2010;Won et al., 2012), which aim at sequentially adapt the filter parameters.However, even these advanced techniques do not solve the estimation versus mitigation trade-off (Vilà-Valls et al., 2015b), because the filter is not able to decouple the different phase contributions (i.e., dynamics and scintillation), and so may provide poor scintillation mitigation capabilities.A scalar tracking approach to overcome these limitations was proposed in Vilà-Valls et al. (2015a,b), where the scintillation physical phenomena is statistically modeled using autoregressive (AR) processes and embedded in the filter formulation.This approach was shown to outperform the SoTA techniques using synthetic data, being effectively able to discriminate between dynamic and scintillation contributions to the estimated carrier-phase.The Cornell Scintillation Model (CSM) (Humphreys et al., 2009(Humphreys et al., , 2010)), a popular synthetic signal simulator that generates scintillation amplitude and phase traces, was used to test the AR-based filtering methodology.Recent studies of real equatorial and highlatitude data have shown that scintillation characteristics are frequency-dependent (Carrano et al., 2012;Sokolova et al., 2015;Jiao et al., 2016).While phase disturbances may be correlated, deep amplitude fades tend to occur at different times on different frequencies, thus motivating the potential performance gain of using inter-frequency aiding and multifrequency architectures.
This paper exploits the statistical relationships between scintillation at different frequency bands to propose a multifrequency carrier tracking method that virtually mitigates scintillation component.The paper extends our previous work on single-frequency carrier tracking loops under scintillation conditions (Vilà-Valls et al., 2015a,b), for which a state-space formulation and a Kalman filter (KF) was considered.Here, an augmented state-space formulation considering multiple frequencies is investigated and a tracking method proposed on the top of it.Remarkably, the resulting method is validated with real data measurements showing notable performance.
To summarize, in this article, we further explore the idea and preliminary results in Vilà-Valls et al. (2017a) with the following four main contributions: multivariate AR ionospheric amplitude and phase scintillation modeling, featuring a joint fitting of the multifrequency scintillation time series.An optimal model order selection based on the Bayesian Information Criterion is used.This is discussed in Section 3; the novel multi-frequency joint carrier phase and scintillation state-space model (SSM) formulation, that takes into account the potential correlation at different frequency bands is introduced in Section 4; the proposed robust filter that performs state estimation on the aforementioned SSM is discussed in Section 5.The solution is based on a nonlinear, adaptive KF formulation; as opposed to previous works where synthetic data was used, the proposed method is analyzed using real ionospheric scintillation data in Section 6.

Ionospheric scintillation
Scintillation is caused by local variations in the structure of the ionospheric electron density along line-of-sight between the satellite and the receiver.The receiver observes signals that are refracted and diffracted such that the phase and amplitude of the signal exhibits large stochastic perturbations.The severity of the scintillation is traditionally quantified by two indices: an amplitude scintillation index, denoted S 4 ; and a phase scintillation index, denoted s f .The indices are computed on a per-signal basis and indicate average intensity of the signal variations over the preceding minute.They have been used for some decades and as a result there exist rich databases of historical data for a wide range of observation points.
Unfortunately these indices offer no insight into the timecorrelatedness of the instantaneous phase and amplitude of scintillation observed on multiple frequencies.A number of recent studies have provided empirical evidence that the phase scintillation exhibits a high degree of correlation between different GNSS frequencies, including L1, L2 and L5.It has been suggested that this correlation persists except in the case of strong amplitude scintillation, when deep amplitude fades occur (Sokolova et al., 2015;Carrano et al., 2012;Jiao et al., 2016).This implies that even when strong phase scintillation is observed, which is not uncommon at high latitudes, these variations in the carrier will be correlated across frequencies.When a receiver observes both amplitude and phase scintillation, then it has been observed that the phase variations are correlated depending on the scintillation intensity, but become decorrelated once deep amplitude fades begin to occur (Sokolova et al., 2015).These deep amplitude fades, it has been shown, tend to be uncorrelated across frequencies (Carrano et al., 2012).
From a carrier tracking perspective, these two observations can be exploited to improve and robustify carrier phase tracking.Firstly, the fact that the phase variations are highly correlated implies that a centralized, or joint tracking of the carrier across multiple GNSS frequencies can improve tracking accuracy (Sokolova et al., 2015).Secondly, the fact that deep amplitude fades occur at different times on different frequencies, implies that it is highly unlikely that a receiver will loose lock on multiple frequencies simultaneously, suggesting that it should be possible to assist one frequency with another (Sokolova et al., 2015;Carrano et al., 2012).These concepts are explored in the work presented here.
The experimental work presented in this manuscript, was conducted on data obtained from the Joint Research Centre (JRC) Scintillation Repository a collection of data with more than 10 hours of scintillation events recorded over Hanoi in March andApril, 2015 (Curran et al., 2015b).These datasets were based on a reconfigurable quad-channel front-end, Fourtune, which had been programmed to collect triplefrequency data, L1, L2 and L5, using 1-bit complex sampling at rates of 5 MHz, 5 MHz, and 30 MHz, respectively (Curran et al., 2015b).The recording system was configured to continuously record 50-minute datasets and to post-process each using an L1 software defined receiver for the purposes of basic scintillation detection.Those datasets in which severe scintillation was identified were then archived for postprocessing, and the others discarded (Curran et al., 2015a).
The post-processing stage employed a multi-frequency open-loop software receiver which exploited precise knowledge of the receiver location, and the well-disciplined reference oscillator to generate accurate reference carrier and code local replicas (Curran et al., 2015b).Once demodulated to complex baseband, the correlator values corresponding to each observed GNSS signal were processed to estimate the phase and amplitude perturbations induced by ionospheric activity.Being an open-loop post-processing scheme, a batch estimation of the amplitude and phase was possible, providing highly accurate and reliable characterization (Curran et al., 2015b).From these 1 kHz estimates ofcarrier phase and amplitude, traditional measures of scintillation activity were computed, such as S 4 and s f .Figure 1 shows an example of real equatorial ionospheric scintillation amplitude data for a scintillation event at three different GPS frequency bands on one particular satellite.This is an example of strong scintillation, exhibiting a few fades below -30 dB, and many prolonged fades below -10 dB.

Multivariate AR scintillation signal model
Ionospheric scintillation can be modeled as a multiplicative channel j s (t), with the received base-band signal being xðtÞ ¼ j s ðtÞsðtÞ þ wðtÞ; j s ðtÞ ¼ r s ðtÞe ju s ðtÞ ; ð1Þ where r s (t) and u s (t) are the ionospheric scintillation envelope and phase components, respectively; s(t) is the baseband equivalent of the transmitted signal; and w(t) is a randomnoise term including both thermal noise and any other disturbance.
The amplitude scintillation strength is typically described by the scintillation index S 4 , and is usually considered within three main regions (Humphreys et al., 2009) It was shown in Vilà-Valls et al. (2013, 2015a,b) for the single-frequency case that both scintillation components can be well approximated using an autoregressive (AR) process.The scalar AR scintillation model is given by with k a constant value, and the mean of the amplitude scintillation process equal to m ¼ k / ð1 À P q i¼1 g i Þ; {g i , b i } are the set of AR coefficients, and {s 2 h r ,s 2 h u } the driving noise variances.In Vilà-Valls et al. (2017a), we considered this scalar AR approximation to build a multi-frequency SSM, where the model order selection and AR parameters were obtained from synthetic scintillation data.Using this scalar approximation in a multi-frequency scenario is a simplistic approach, because it does not take into account the possible correlations among frequencies, and therefore loses the characteristics of the physical phenomena (see Section 2).To overcome the limitations of the scalar AR processes model in this context, we propose a new multivariate AR (MAR) ionospheric scintillation formulation.Real ionospheric scintillation data was used, selecting the model order and parameters.
Definition 1 (MAR(p)) A multivariate AR model of order p is defined as where p is the model order;  value of The two main MAR modelling challenges are: (i) selectionthe correct model order selection p, and (ii) the time-series analysis to obtain an estimate of the MAR model parameters, fw; fA i g p i¼1 ; S e g.Both points are discussed in the remainder of this section.

MAR model order selection
In the statistical and signal processing literature, the Akaike Information Criterion (AIC) (Akaike, 1974), extensions of the AIC such as the corrected AIC (AICc) and the generalized information criterion (GIC), and the Bayesian Information Criterion (BIC) (Schwarz, 1978), are widely used metrics for model order selection.The latter has become the method of choice mainly because: i) while the AIC has a nonzero overfitting probability, this tends to zero in the BIC when the sample size grows, and ii) it is an asymptotic approximation of the optimal maximum a posteriori (MAP) rule (Stoica and Selen, 2004).In this contribution we use the BIC as the model order selection criterion.
In this subsection we define the quantities for a generic time-series composed of N samples, y ∈ R N , and assume u is the vector parameters of the model with order p, which we aim to select.

Definition 2
(BIC) The Bayesian Information Criterion is defined as (Stoica and Selen, 2004) with y being the available data, u the parameter vector (whose dimension can vary with the model order) and p corresponds to the model order.p p (y, u) is the likelihood function which depends on the parameter vector u and the model order p, and the maximum likelihood estimate of the parameters is given by û ¼ arg max u ln p p ðy; uÞ: The BIC rule selects the order p that maximizes (7).Definition 2 corresponds to the BIC when a single timeseries is fitted.We are, however, interested in jointly fitting multiple time-series (from the various frequencies).This concept is generalized in Definition 3.

Definition 3
(Multivariate BIC) The Bayesian Information Criterion for the multiavariate AR model MAR(p) in Definition 1 yields to a model order estimate according to Neumaier and Schneider (2001).
where L p ¼ ln det ðN À dp À 1Þ Ŝe ðpÞ and Ŝe ðpÞ is an estimate of the noise covariance matrix obtained from the N samples.
The estimates of the noise covariance matrix Ŝe ðpÞ are typically obtained for a range of model order values, p min p p max .An efficient iterative way to compute an estimate of L p using a QR factorization and a stepwise downdating strategy is shown in Neumaier and Schneider (2001).

MAR model parameters
Several MAR model parameters estimation methods based on the minimization of the prediction error are available in the literature.A thorough comparison of the most common estimators is available in Schlögl (2006).The main conclusion is that when enough samples are available (i.e., N ≥ 400), the performance of the different methods is similar.In this contribution we propose to use the stepwise least squares estimation method proposed in Neumaier and Schneider (2001), which is available as a Matlab TM package named ARfit (Schneider & Neumaier, 2001).This method first rewrites the MAR model as a regression model, and then obtains the parameters using a least squares method.

Ionospheric scintillation MAR formulation
The MAR multi-frequency ionospheric scintillation model considering M frequency bands (i.e., B 1 to B M ) is given by Notice that each of these matrices has dimension M Â M, given by the number of frequency bands.It is worth saying that the MAR model noise covariance matrices, S r and S u , provide a good indicator of the interfrequency correlation of the time-series.It has been found empirically from the estimated noise covariances that while S r is almost diagonal (i.e., no correlation in the ionospheric scintillation amplitude), significant values off the main diagonal appear in S u (i.e., higher inter-frequency correlation in the ionospheric scintillation phase).
A key point to exploit the MAR scintillation approximation within a filtering framework is the model order selection.In Vilà-Valls et al. (2013, 2015a,b), the selection was done inspecting the power spectral densities (PSD), and subsequently in Vilà-Valls et al. (2017a), we used a statistical test based on the partial autocorrelation functions (PAF).Both analyses were done using synthetic data.In this contribution we use the BIC in Definition 3 and perform the fitting on real multi-frequency scintillation data.
Figure 2 shows the C/N 0 (carrier-to-noise density ratio) and amplitude scintillation index estimations, S 4 , at three different GPS frequency bands (i.e., L1, L2, and L5).These results correspond to two ionospheric scintillation events recorded over Hanoi in March and April 2015 (cf.Sect.2).While the first event (top figure) corresponds to a moderate to low scintillation scenario, the second one (bottom figure) is clearly a strong scintillation case.Thus, the real data available (over 3 hours of equatorial scintillation) contains both low, moderate and severe scintillation scenarios.
To obtain the model order for both amplitude and phase scintillation processes, we processed over 4 hours of 3dimensional time-series to compute the BIC metrics in (8).BIC results are obtained from non-overlapped sequences of N = 500 samples (5 s).In contrast to the single-frequency fitting to synthetic data given in Vilà-Valls et al. (2013, 2015a,b), the conclusion from the MAR analysis is that a suitable model order to fit the majority of scintillation scenarios is given by q = 6 and p = 5, i.e. a MAR(6) and MAR(5) for the scintillation amplitude and phase, respectively.These results are in contrast to the analysis of previous contributions (based on synthetic data), where lowerorders were sufficient.To justify the choice of the model order, we plot in Figure 3 the histograms of the order selection obtained from real multi-frequency ionospheric scintillation data.(Spilker, 1996) where T s is the sampling rate, k stands for the discrete time t k = kT s , A k is the signal amplitude at the output of the correlators, d k is the data bit, R(⋅) is the code autocorrelation function and Dt k ; Df d;k ; Du k n o are, respectively, the code delay, Doppler shift and carrier phase errors.Focusing on carriertracking, we work under the assumption of perfect timing synchronization and data wipe-off, resulting in a signal model where amplitude a k and phase u k may include any non-nominal harsh propagation disturbance.For instance, including the ionospheric scintillation into the signal model leads to The dynamics phase term u d,k refers to the line-of-sight (LOS) phase evolution due to the relative movement between the satellite and the receiver.Standard carrier tracking architectures, being proportional plus integral controller in the form of a PLL, model u d,k using a Taylor approximation of the time-varying phase evolution according to the expected dynamics.Typically, a 3 rd order representation of the dynamics phase evolution is considered where u 0 (rad) is a random phase value, f d,k (Hz) the carrier Doppler frequency, and f r,k (Hz/s) the Doppler frequency rate.

Multi-frequency GNSS state-space model
In a multi-frequency architecture, a single satellite transmits simultaneously at different frequencies, therefore at the receiver side several measurements are available.Although the methodology in this paper is general for any GNSS signal, we focus on a modern GPS triple band satellite, which transmits on the L 1 , L 2 and L 5 frequencies.In this case, the measurements are given by which may be useful to rewrite in real form as where for each frequency, The different LOS carrier-phase evolutions, when referencing all of the dynamics to L 1 , are given by with f L 1 , f L 2 and f L 5 being the carrier frequencies at the different bands.It is important to notice that the same Doppler frequency terms, f d,k and f r,k , appear in both frequencies.The multi-frequency noise diagonal covariance matrix is In general, x k refers to the state to be tracked and is defined as The LOS dynamic carrier-phase state evolution stands for possible uncertainties or mismatches on the dynamic model.The covariance Q d,k is usually a priori fixed according to the problem under consideration.
From the MAR scintillation model in ( 9), the multifrequency scintillation amplitude state evolution is where the elements not specified are 0, I M refers to the identity matrix of dimension M Â M, v r;k $ Nð0; Q r;k Þ with block diagonal covariance matrix Q r,k = blkdiag (S r , 0 M(q-1)ÂM(q-1) ) and k r is a null vector except for the first M values given by k.
Equivalently, we define the scintillation phase state evolution from (10) as ).Finally, the complete state evolution equation reads, where the Gaussian process noise has a block diagonal covariance matrix ). Equations ( 15) and ( 23) define the new SSM formulation.Notice that in contrast to the SSM in Vilà-Valls et al. (2017a), in this case the possible correlation between scintillation amplitude and phase at different frequency bands is accounted for.

New robust filtering methodology
In this article we elaborate on the previous KF-based architecture proposed in Vilà-Valls et al. (2017a), and extended in Section 4, where the scintillation physical phenomena is mathematically modeled and embedded into the state-space formulation.A nonlinear extended KF (EKF) formulation allows to directly operate with the received signal samples, avoiding the problems associated to the use of a phase discriminator, that is, nonlinearities, loss of Gaussianity and possible saturation at low C / N 0 .Using the SSM in ( 15) and ( 23), and the MAR formulation in Section 3.3, recall that the main parameters of the filter are which mainly depend on the amplitude MAR(6) and phase MAR (5) set of parameters obtained from fitting to real multi-frequency data, C ¼ ½fG i g 6 i¼1 ; k T ; fD i g 5 i¼1 ; S r ; S u .The complete state x k (n x = 38) to be tracked is defined in equation ( 25), Q d,k is a priori fixed by the user from the expected dynamics, and the measurement noise variances (i.e., s 2 n L1 ;k , s 2 n L2 ;k and s 2 n L5 ;k ) can be sequentially adjusted from the nominal C / N 0 estimators at each frequency band, being available at the receiver.
The main idea behind the EKF is to linearize the nonlinear measurement function around the predicted state estimate, and then use the standard linear KF equations (Anderson and Moore 1979).In this case, the measurement function is given by with the different matrices given in equations ( 26)-( 28).Note that the matrix is evaluated at the predicted state, 2 k|k-1 refers to the predicted estimate at time k using measurements up to time k -1, and k|k refers to the estimated value at time k using measurements up to time k.
For the sake of completeness, the EKF formulation of the general carrier-phase tracking problem is sketched in Algorithm 1.Notice that the Kalman gain, K k , is computed from the possibly time-varying system noise matrices, Q k and R k , which implies that the filter is adapting its behavior to the actual working conditions and has an implicit adaptive bandwidth, in contrast to the traditional PLL with constant coefficients (i.e., fixed bandwidth) (Vilà-Valls et al. 2017b).

Algorithm 1 Carrier-phase tracking EKF formulation
In practice, we compute an estimate of the MAR parameters using a subset of the real multi-frequency ionospheric scintillation data available.To obtain meaningful results, in the simulations we use a small sample set to obtain the model parameters and process the rest of the scintillation time-series to assess the filter performance.In other words, we test the robustness of the method by decoupling the model fitting from the state tracking.Finally, the block diagram of the new multi-frequency nonlinear EKF-based architecture for efficient scintillation mitigation is sketched in Figure 4.It includes both the state estimation procedure and the update of the prediction/estimation error covariance matrices.

Results and discussion
The performance of the proposed method (named MAR-EKF in the simulations) is analyzed using the real equatorial multi-frequency ionospheric scintillation data described in Section 2. We processed 6 sets of data (i.e., over 4 hours of data), each corresponding to a different scintillation event, containing low, moderate and severe ionospheric scintillation propagation conditions (see Eq. ( 2)).The different events are described in Table 1, with the corresponding estimation of the scintillation index S 4 shown in Figure 5. From the state vector in (25), the parameters of interest are the phase components related to the dynamics of the receiver, i.e., obtaining a good estimate of these parameters directly implies good ionospheric scintillation mitigation capabilities.Therefore, this is the final goal of the new receiver architecture.Note that in some applicationsit may be interesting to obtain an estimate of the scintillation amplitude and phase components, r s,k and u s,k .
The results obtained with the MAR-EKF are compared to the current state-of-the-art techniques, and the latest contributions using scalar AR scintillation model approximations: -3 rd order PLL with equivalent bandwidth B w = 5 Hz.We consider a Doppler profile given by an initial random phase in [-p, p], initial Doppler f d,0 = 50 Hz and rate f r,0 = 100 Hz/s, and a relatively low nominal C / N 0 = 30 dB-Hz.The signal is sampled at T s = 10 ms.The measure of performance is the root mean square error (RMSE) on the carrier-phase of interest, i.e., u L 1 d;k , u L 2 d;k and u L 5 d;k .We define the RMSE (radians) over time for each frequency L j as which is obtained from long real data sets of N t samples.The proposed methodology is characterized via an in depth analysis, considering the steady-state performance for several representative scenarios: (i) a low scintillation scenario which is used as an architecture validation, (ii) scintillation mitigation capabilities in challenging moderate and severe scintillation scenarios, (iii) robustness to modeling mismatch, where we take into account possible scintillation under and overestimation, and (iv) the estimation performance for the ionospheric scintillation amplitude and phase components.

Low scintillation scenario: architecture validation
In general, it is held that low ionospheric scintillation conditions are given for S 4 < 0.4.From Figure 5, we easily identify these conditions in the scintillation events #1 (t = 1500 -2500s), #2 (t = 100 -1800s), #3 (t = 100 -1700s) and #5 (t = 100 -1900s).A desired behavior of any tracking method designed to mitigate moderate and severe scintillation conditions is that it performs at least as good as the standard tracking techniques when low or no scintillation is present.Therefore, analyzing the performance of the proposed MAR-EKF in such scenarios can be understood as a representative architecture validation procedure.
We show the RMSE t (L j ) in radians, obtained for two different low scintillation scenarios (500s of data), in Table 2.In this case, we can see that the performance obtained with the two EKF-based solutions, i.e., AEKF-AR and MAR-EKF, is similar but slightly better with the MAR-EKF.Moreover, these methods provide a lower RMSE than the one given by the PLL and the three discriminator-based KF solutions, i.e., KF, AKF and AKF-AR.Therefore, we have two main conclusions: i) the proposed architecture is a valid approach even if no scintillation is present in the signal, that is, the filter is robust and including the scintillation components into the state to be tracked does not induce an error on the estimation of the dynamics related phase, and ii) it may be preferable to avoid discriminator-based architectures and operate directly with the received signal samples.Moreover, in this case we can include the scintillation amplitude into the state formulation, thus being aware of possible fading and amplitude variations.

Assessing the scintillation mitigation capabilities in moderate and severe scintillation scenarios
The main goal of the proposed architecture is to effectively counteract the undesired ionospheric effects in moderate and severe scintillation scenarios, without sacrificing performance in nominal scenarios (cf.Section 6.1).In this section, we analyze the performance of the MAR-EKF under harsh propagation conditions as per Table 1.

Moderate scintillation
For the sake of completeness, we plot an example of the steady-state instantaneous dynamics phase estimation error for the moderate scintillation event #1 (t = 200-600s) at L5, obtained with the AKF-AR, AEKF-AR and the MAR-EKF.Notice that we do not plot the results for the other methods because their performance is significantly poorer (see Table 3 for RMSE results).In view of the results, the benefits of the new multi-frequency approach with respect to the other filters using scalar AR approximations are apparent.
To obtain statistically more representative results, we show the RMSE t (L j ), obtained for events #1 (t = 500 -1000s) and #4  (t = 1300 -1800s), in Table 3.These results verify the appropriate behavior of the new architecture, indicating that it outperforms the other methods under moderate scintillation conditions.Again, there is a significant performance improvement attained when avoiding the use of a discriminator.Specifically, the performance obtained with the AEKF-AR and MAR-EKF is improved with respect to the othermethods.This is due to the Gaussianity assumption,which is does not hold when working with a discriminator under fading conditions (Vilà-Valls et al. 2017b).

Severe scintillation
Finally, we assess the performance of the MAR-EKF in challenging, severe scintillation conditions.Figure 7 shows an example of the steady-state instantaneous dynamics phase estimation error for the severe scintillation event #6 (t = 1800-2200s) at L2, obtained with the AKF-AR, AEKF-AR and the MAR-EKF.The best estimation performance is given by the MAR-EKF.
The RMSE t (L j ), obtained for events #3 (t = 2000 -2500s), #4 (t = 2000 -2500s) and #6 (t = 1700 -2200s), is shown in Table 4.These results confirm that the new architecture is superior to the rest of tracking methods also under severe scintillation, providing a significant performance improvement with respect to standard approaches and the scalar EKF-based solution in Vilà-Valls et al. (2015a), which justifies the interest and importance of considering a multi-frequency approach.Overall, the MAR-EKF appears to be relatively robust against a wide range scintillation propagation conditions.

Robustness to modeling mismatch
To fully characterize the proposed methodology, we analyze the impact of an incorrect MAR model parameters estimation on the filterperformance, that is, fitting the scintillation amplitude MAR(6) and phase MAR(5) to a set of real data with a different scintillation behavior than the actual signal being processed.Two cases are considered: i) underestimation of the scintillation conditions, which implies that while the filter deals with a severe scintillation scenario, the MAR model parameters are fitted to moderate scintillation data; and ii) overestimation of the scintillation conditions, which refers to themoderate scintillation case with MAR model parameters fitted to severe scintillation data.The RMSE t (L j ), obtained for the severe scintillation event #6 (t = 1700 -2200s) considering an underestimated fitting, and for the moderate scintillation event #6 (t = 1300 -1700s) with an overestimated fitting, is shown in Table 5.It is interesting to note that the performance degradation when having a model mismatch is reasonably bounded.In practice, it is always desirable to overestimate, rather than underestimate, the scintillation intensity.

Scintillation amplitude and phase estimation
In this section we analyze the estimation performance for the scintillation amplitude and phase components, r s,k and u s,k .We define the corresponding RMSE over time for each frequency L j as which are obtained from N t data samples.Notice that standard methods track the complete phase of the signal, therefore we do not have an estimate of the scintillation components.The results obtained with the MAR-EKF are compared to the scalar AKF-AR (i.e,only scintillation phase) and AEKF-AR.First we plot an example of the scintillation amplitude and phase estimates delivered by the MAR-EKF, for the severe scintillation event #3 at L5, in Figure 8.We can observe a remarkable recovery of scintillation components by the MAR-EKF filter.Both RMSE t ðr s;L j Þ and RMSE t ðu s;L j Þ are shown in Table 6.Again, we confirm that the MAR-EKF provides the best estimation performance for both scintillation amplitude and phase components.The performance gain with respect to the other methods using scalar AR approximations is clear for the severe scintillation scenario.In the moderate scintillation case, the performance obtained with the AEKF-AR and the MAR-EKF are equivalent.

Cycle slip performance analysis
To conclude the characterization of the new methodology, we assess its robustness to cycle slips, which is particularly necessary in carrier-based positioning techniques such as realtime-kinematic (RTK) and precise-point-positioning (PPP) that rely on the integrity of carrier phase measurements.
In terms of cycle slips, the most challenging scenario is the severe scintillation event #6.We can see in Figure 5 that this event has the highest values for the scintillation index S 4 .We show the phase error for the different methods in Figure 9.While some cycle slips appear with the discriminator-based PLL, KF, AKF and AKF-AR, the methods including both scintillation components into the SSM formulation, AEKF-AR, MAR-EKF, are more robust (i.e., no cycle slips).These results confirm the enhanced filtering performance of the new approach.Notice that for this particular scenario, the number of cycle slips with the AKF is higher than with the standard PLL and KF.This is because the filter is not able to cope with the strong fadings during severe scintillation due to the fast changes on the estimated C/N 0 , which controls the noise variance at the output of the discriminator.This is also the case with the discriminator-based AKF-AR, which does not take into account the scintillation amplitude.The number of cycle slips for event #6 and the three GPS frequency bands is shown in Table 7, which again confirm the superiority in terms of robustness of the new MAR-EKF.

Conclusions
Carrier-phase measurements can be severely degraded due to ionospheric scintillation.This paper proposes a new tracking methodology that is able to extract the desired phase component, used for positioning, from the scintillation component.A key point is first to model both ionospheric scintillation components, amplitude and phase, using an AR model formulation, and then include them into the state to be tracked by the filter.This allows to decouple the phase contribution due to the LOS dynamics from the phase perturbations induced by the ionosphere.Particularly, we want to exploit the advantages of multi-frequency architectures, because the scintillation characteristics are frequency-dependent.To take into account possible correlations among frequencies we must resort to multivariateAR models.A new state-space model is introduced and a Kalman-type state estimator operating on the outputs of the correlators implemented.The new method is validated using real scintillation data recorded in Hanoi, March and April 2015.It has been shown that the new filter outperforms state-of-theart solutions in the tested scenarios, both in terms of mean square error and robustness to cycle slips.This approach provides a solution to the estimation/mitigation challenge when dealing with scintillation corrupted measurements.
One of the main challenges in the design of the proposed solution is the estimation of the multivariate AR model parameters.While in the analysis conducted in this article these parameters have been obtained by fitting the desired model to a small part of the real time-series available (i.e., then processing a different part of the data), in real-life conditions, where the scintillation conditions evolve over time, we should consider an iterative procedure to adapt the filter to the actual working conditions.Another important point is how to fix the LOS phase dynamics covariance matrix, which depends on the expected dynamics.To fully characterize the new method, it would be desirable to have a larger set of data, recorded at different points over the earth, at different latitudes and containing several events.Even if the set of data available is representative of different scintillation conditions, it would be interesting to obtain a statistically more representative analysis.

Fig. 1 .
Fig.1.An example of 100 s of real equatorial scintillation amplitude data for a severe scintillation event at three different GPS frequency bands.

Fig. 2 .
Fig. 2. C/N 0 and S 4 scintillation index estimation at three GPS frequency bands for two ionospheric scintillation events recorded over Hanoi in March and April 2015.

Fig. 3 .
Fig. 3. Empirical model order probability density functions obtained from real multi-frequency ionospheric scintillation data.

-
Standard discriminator-based KF-based tracking.-Adaptive discriminator-based KF (AKF), adjusting the measurement noise variance from the C/N 0 estimate.-The scalar AKF-AR presented in Vilà-Valls et al. (2015b), only tracking thescintillation phase, and using a C/N 0 estimate to tune the measurement noise variance at the output of the discriminator.-The scalar AEKF-AR presented in Vilà-Valls et al. (2015a), tracking both scintillation amplitude and phase, treating each frequency separately.
the kth sample of a d-dimensional time-series; A 1 to A p ∈ R dÂd are the coefficient matrices of the MAR model; the random vector term e k = [e k,1 , ... , e k,d ] ⊤ is a zero-mean Gaussian random variable with covariance matrix S e ∈ R dÂd .e k is assumed uncorrelated in time (i.e., E {e k e ℓ } = 0 for k ≠ ℓ), but possibly correlated among time-series (i.e., S e can containoff-diagonal elements); w ∈ R d is a parameter vector of intercept terms to allow a nonzero mean time-series, which implies an expected

Table 1 .
Real equatorial multi-frequency ionospheric scintillation data characterization.

Table 4 .
Root mean square phase tracking error considering real scintillation data at different GPS frequency bands for the strong

Table 7 .
Number of cycle slips for the real ionospheric scintillation event #6.