GPS Solutions (2022) 26:17 https://doi.org/10.1007/s10291-021-01193-0 ORIGINAL ARTICLE Performance assessment of a low‑complexity autoregressive Kalman filter for GNSS carrier tracking using real scintillation time series Sergi Locubiche‑Serra1 · Gonzalo Seco‑Granados1 · José A. López‑Salcedo1 Received: 28 January 2021 / Accepted: 11 October 2021 / Published online: 26 November 2021 © The Author(s) 2021 Abstract 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 experimenta- tion 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. Keywords Autoregressive modeling · Carrier tracking · GNSS · Ionospheric scintillation · Kalman filter Introduction entails moving into the arena of high-precision, safety–criti- cal applications such as aviation, maritime navigation and It is well known that global navigation satellite systems autonomous vehicle driving, just to mention a few. One key (GNSS) have become the key technology for positioning and enabling feature is that GNSS carrier phase measurements navigation purposes using satellite ranging signals. Their do provide ultra-precise positioning information, and thus, maturity, widespread deployment and accuracy in open-sky their exploitation becomes essential. However, the use of environments are leading to a growing demand for extend- GNSS in environments with many propagation impairments, ing GNSS beyond the limits of its original designs. This 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 next- * Sergi Locubiche‑Serra sergi.locubiche@uab.cat generation GNSS receivers. Some examples are weak signal reception and multipath. But in particular, one of the most Gonzalo Seco‑Granados gonzalo.seco@uab.cat challenging sources of error to deal with in GNSS is the so-called ionospheric scintillation effect (Lee et al. 2017). José A. López‑Salcedo jose.salcedo@uab.cat The ionosphere is the upper earth’s atmosphere ionized by solar radiation, and it has a significant influence on tran- 1 Department of Telecommunications and Systems sionospheric radio wave propagation. Ionospheric scintilla- Engineering, IEEC‑CERES, Universitat Autònoma de tion is a known effect of space weather whereby ionospheric Barcelona (UAB), 08193 Bellaterra, Barcelona, Spain 13 Vol.:(0123456789) 17 Page 2 of 14 GPS Solutions (2022) 26:17 electron density irregularities introduce rapid fluctuations From the above observations, the use of a KF with a onto GNSS signals when crossing the ionosphere in form kinematic model augmented with a linear Gaussian model of random amplitude fades and abrupt phase jumps (Kintner such as the AR becomes then the natural leap forward for et al. 2007). These have a detrimental effect on the perfor- optimally dealing with both phase dynamics and scintilla- mance and stability of GNSS receivers and the proper deter- tion experienced by received GNSS signals. This gives rise mination of the synchronization parameters, particularly the to a hybrid autoregressive Kalman filter (KF-AR) tech- carrier phase dynamics of interest. Ionospheric scintillations nique, an innovative approach pioneered by the authors can cause cycle slips, severe carrier phase jitter, loss of sig- of the present work in Locubiche-Serra et al. (2016) and nal even in line-of-sight (LoS) conditions, signal integrity Vilà-Valls et al. (2013), whereby the above-mentioned issues and degrade navigation performance. Furthermore, kinematic and AR models are hybridized into one single although ionospheric scintillation can occur anywhere on state-space formulation. The underlying idea is to discern earth during intense solar activity, it is more frequent in the between phase dynamics and scintillation and thus provide earth’s high- and low-latitude/equatorial regions (Kintner decoupled estimates of both while making explicit that et al. 2009) and it certainly poses serious concerns to the they must be determined in the presence of one another. widespread deployment of GNSS in emerging countries. The latter means, indeed, that the KF-AR can be employed Therefore, the design and validation of novel signal pro- to obtain either dynamics-robust phase scintillation esti- cessing techniques for providing carrier phase estimates mates or scintillation-robust phase dynamics estimates. that are robust to ionospheric scintillation have become of Even though remaining out of the scope of this work, the paramount importance. former does enable GNSS for ionosphere-related scien- In order to exploit carrier phase measurements, GNSS tific applications such as ionospheric scintillation monitor- receivers usually implement a carrier (along with code) ing (ISM) itself. Our focus is placed on the latter, where tracking stage that keeps them synchronized with the sat- robustness to scintillation is sought for positioning and ellites so that it accurately monitors any variation in the navigation purposes. The final goal is to locally recon- observed satellite-to-user dynamics. In most existing receiv- struct the phase scintillation disturbance so that it can be ers, this is implemented using the well-known phase-lock removed from the input carrier phase and result in a clean loop (PLL) technique, a simple closed-loop architecture that phase depending on the user’s dynamics only. compares the input carrier phase values to a local replica and The existing contributions on this topic are focused on drives the resulting error to zero by properly adjusting the the design and characterization of the above techniques. phase of the local oscillator (Gupta 1975). Notwithstanding, However, an analysis of their efficiency based on real data this technique is known to experience serious troubles in experimentation is still lacking. In that sense, our objec- the presence of unexpected propagation impairments, and it tive is to validate the applicability and suitability of KF-AR becomes the bottleneck in GNSS receivers under scintilla- techniques in realistic environments for scintillation-robust tion conditions (Lee et al. 2017). As a countermeasure, the GNSS carrier tracking. To this end, an extensive experimen- use of advanced tracking techniques such as those based on tation campaign is carried out using real scintillation data the Kalman filter (KF) is gradually being introduced from captures at high and low latitudes to introduce real-world- the general framework of optimal minimum mean square representative scintillation disturbances in the analysis, in estimation (MMSE). Some examples can be found in Maca- contrast to many existing contributions that are restricted biau et al. (2012), Won et al. (2012), Jiang et al. (2017), Susi to synthetic data. In this sense, the work leads to a three- et al. (2017) and Yang et al. (2017), where carrier phase fold contribution. First, we confirm the suitability of AR dynamics is estimated using KFs that encompass kinematic processes for modeling phase scintillation. Second, we pro- models in their state-space formulation. As a matter of fact, pose a novel adaptive, low-complexity KF-AR technique the application of KF-based GNSS carrier tracking turns into designed with the practical implementation of the approach an interesting approach for dealing with phase scintillation in the spotlight, aiming at facilitating its applicability in when noting that the random nature of the latter can be mod- practice. This implementation is intended to circumvent the eled through, for instance, a correlated Gaussian distribution performance and computational limitations of other fixed (Humphreys et al. 2010). This opens the door to phase scin- and adaptive implementations previously proposed in the tillation being accommodated by the KF in a natural man- literature and attain the theoretical performance given by ner, thus taking advantage of its optimality properties under the well-known Bayesian Cramér–Rao bound (BCRB) irre- Gaussian disturbances for linear models. A case of particu- spective of the input working conditions. Third, we provide lar interest is the class of autoregressive (AR) processes, some guidelines on the behavior of the proposed techniques which are gaining relevance for the problem at hand and are as a function of the environmental working conditions. We claimed to provide potentially promising results (Nunes and also evaluate the performance limits at which the proposed Sousa 2014; Fohlmeister et al. 2018; Morton et al. 2020). 13 GPS Solutions (2022) 26:17 Page 3 of 14 17 techniques can operate in practice and outperform conven- tional 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 con- version. Then, at its output, the digital processing of GNSS signals is usually performed through two sequen- tial stages: the acquisition and tracking stages, followed by the position, velocity and time (PVT) module, which ultimately determines the user’s position using the pro- cessed information. The acquisition stage is in charge of detecting the sat- ellites that are present in the received signal and finding a tentative estimate of their code delays (i.e., pseudor- Fig. 1  Block diagram of conventional GNSS code and carrier track- anges). As some Doppler shift is present owing to receiver ing architectures. Code tracking is used for fine code delay estima- oscillator inaccuracies and the relative motion between tion. Carrier tracking is used for fine carrier phase estimation 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 coherent in the PLL case). The result is fed to the dis- function. Therefore, when satellites are detected, the criminators, which provide a function that is proportional acquisition stage also provides coarse estimates of the to the local replica errors to be corrected. The loop filters code delay and Doppler shift at which the signal of interest (LF) aim at removing the noise at the discriminator output is located. These outputs are then provided to the tracking and providing a smoothed version of the corrections to stage. 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 Tracking stage work, we will focus on the PLL architecture, as detailed next. The tracking stage is a closed-loop architecture with a two- For a given satellite, consider the signal x(n) at the fold objective. First, to refine the synchronization param- input of the tracking stage with the following discrete- eters estimated in the acquisition stage. Second, to accu- time complex baseband equivalent signal (Seco-Granados rately follow any possible variation of these parameters et al. 2012), over time. The underlying idea is to find the code delay √ and carrier phase values to generate a local replica of the x(n) = Pd(n − 𝜏)r(n) + e(n) (1) signal of interest that aligns with the received one. This is where n is the discrete-time variable, P is the signal power, done by recursively comparing the local replica with the d(n) is the navigation message, 𝜏 is the code delay, and e(n) received signal and driving the resulting error to zero. To is the contribution due to receiver noise. The term r(n) is this end, a closed-loop system such as the one shown in the transmitted signal encompassing the receiver-to-satel- Fig. 1 is employed, formed by two coupled architectures, lite synchronization parameters, which can be expressed as one for code tracking and one for carrier tracking, usually (Seco-Granados et al. 2012), referred to as delay-lock loop (DLL) and PLL architec- tures, respectively. Both follow the same basis. First, the r(n) = c(n − 𝜏)ej𝜃d (n) (2) input signal is correlated with the local replica using a set of correlators, namely early (E), prompt (P) and late (L), that is, 𝜏 delays the pseudorandom sequence c(n) linked to formed by integrate-and-dump (I&D) blocks performing the satellite under study, and 𝜃d (n) is the carrier phase of a coherent ( Nc ) and non-coherent ( Ni ) integrations (only complex exponential that contains the Doppler shift due to 13 17 Page 4 of 14 GPS Solutions (2022) 26:17 signal dynamics. The PLL aims at providing an estimate Signal models of the latter, denoted as 𝜃̂d (n), and using it for generating the local replica in a way to compensate for the complex The term 𝜃(n) is the phase to be dealt with at the carrier track- exponential in (2), ing level. For this purpose, the corresponding signal models ̂ must first be defined. The models adopted in this work are r(n)e−j𝜃d (n) = c(n − 𝜏)ej𝜃𝜀 (n) (3) described next. where the phase error 𝜃𝜀 (n) ≐ 𝜃d (n) − 𝜃̂d (n) is driven to zero by iterating the loop. To that effect, the input signal is cor- Signal model for carrier phase dynamics related with the local replica c(n − 𝜏) ̂ through the P correla- tor in Fig. 1, denoted as yP (𝜏) ̂ , which is centered at 𝜏̂ , the For the problem of GNSS carrier tracking, we assume that delay estimated by the DLL. Note that the DLL employs the discrete-time evolution of the carrier phase due to signal two additional correlators, the E–L correlators, centered at dynamics can be approximated by a third-order kinematic 𝜏̂ plus some shift 𝛿 at each side, which facilitate the detection model, which has been found to be a self-standing approach of the main correlation peak. In contrast, the PLL employs of interest in many practical cases (Bar-Shalom et al. 2004), only the P one. 1 1 𝜃d (n) ≈ 𝜃d (n − 1) + Ts 𝜃̇ d (n − 1) + Ts2 𝜃̈d (n − 1) + Ts3 𝜃⃛d (n − 1) When feeding the result to the PLL discriminator, the 2 3! objective is to extract the phase error 𝜃𝜀 (n) of the complex (6) exponential in (3). For this purpose, several discriminators where Ts is the sampling time, and the terms 𝜃̇ d (n) , 𝜃̈d (n) can be employed, highlighting the PLL-, Costas PLL- and and 𝜃⃛d (n) are the Doppler shift, Doppler rate and Doppler frequency-lock loop (FLL)-type ones. Since we are interested acceleration (i.e., first, second and third derivatives of the in phase measurements, we focus on the PLL-type discrimi- carrier phase), whose evolution is given by, nators, specifically the well-known four-quadrant arctangent (ATAN2) discriminator as it duplicates the dynamic range of 1 𝜃̇ d (n) ≈ 𝜃̇ d (n − 1) + Ts 𝜃̈d (n − 1) + Ts2 𝜃⃛d (n − 1) the Costas PLL and is robust to phase variations. In addition, 2! 𝜃̈d (n) ≈ 𝜃̈d (n − 1) + Ts 𝜃⃛d (n − 1) (7) it is the optimal maximum likelihood (ML) phase extractor at high carrier-to-noise ratio (C/N0) (Proakis 2001), while it also 𝜃⃛d (n) ≈ 𝜃⃛d (n − 1) fits very well into GPS signals synchronized at bit level, as well as into the data-less nature of pilot signals which follows the structure in (7). This is then the signal ( )growingly used model to be used later to formulate the problem of GNSS nowadays. Its output, denoted as Δ𝜃ATAN2 𝜃𝜀 , is given by, carrier tracking. ( [ ]) ( ) Im yP (𝜏)̂ Δ𝜃ATAN2 𝜃𝜀 = arctan [ ] (4) Signal model for phase scintillation Re yP (𝜏)̂ [ ] [ ] When approximating the discrete-time evolution of phase where the terms Im yP (𝜏) ̂ and Re yP (𝜏) ̂ correspond to the imaginary and real parts of the prompt correlator, scintillation by a pth-order AR process, namely AR(p), the respectively. following signal model is adopted (Kay 1993), ∑ p Effect of ionospheric scintillation 𝜃s (n) = 𝛼m 𝜃s (n − m) + sp (n) (8) m=1 When the signal received from a given satellite is affected by { }p where 𝛼m m=1 is the set of p AR coefficients and sp (n) is the ionospheric scintillation, the effect can be modeled as a multi- so-called AR driving noise or prediction error ( of an)AR(p) plicative channel introducing signal amplitude As (n) and phase process, usually modeled as sp (n) ∼ N 0, 𝜎s2 with 𝜃s (n) variations to the transmitted signal r(n) in (Vilà-Valls [ ] p et al. 2013), E sp (i)s∗p (j) = 0 for i ≠ j . This is then the signal model to be used later to formulate the problem of phase scintillation As (n)ej𝜃s (n) r(n) = As (n)c(n − 𝜏)ej(𝜃d (n)+𝜃s (n)) = As (n)c(n − 𝜏)ej𝜃(n) tracking. Additionally, the power spectral density (PSD) of (5) an AR(p) process is given by (Kay 1993), in such a way that the carrier phase perceived by the PLL becomes the combination of the signal dynamics plus the � � 𝜎s2 p scintillation phase, that is, 𝜃(n) ≐ 𝜃d (n) + 𝜃s (n). S𝜃s ej𝜔 = (9) � ∑p �2 �1 − 𝛼k e−j𝜔k � � k=1 � 13 GPS Solutions (2022) 26:17 Page 5 of 14 17 with 𝜔 the discrete-time frequency variable. Expression (9) Detrended phase is used next to evaluate the feasibility of modeling phase 2 scintillation through an AR process. (rad) 0 -2 AR model fitting of phase scintillation 1.66 1.67 1.68 1.69 10 5 Once the signal models are defined, this section is intended Detrended power to show the suitability of modeling scintillation phase 5 0 through AR processes, a required step before formulating (dB) -5 the KF-AR. -10 -15 Scintillation time series 1.66 1.67 1.68 1.69 TOW (s) 5 10 Scintillation is reported to be frequent in low- and high- latitude regions and in the dark times of the day and year Fig. 2  Detrended phase and power of DAK.2014.335.22.GPS24 cap- (Aarons 1982). Following these baselines, we analyze two ture, representative of low-latitude scintillation. The focus is placed on the segment delimited within vertical stripped lines different scintillation data captures gathered at monitoring stations of the Monitor network from the European Space Agency (ESA). First, a capture from the station in Dakar, Detrended phase Senegal (lat. 14.765, long. 352.624°) as representative of 2 low-latitude scintillation, taken on December 1, 2014 from (rad) 22:00 to 23:00 h local time coinciding with Coordinated 0 Universal Time (UTC + 0) for GPS satellite vehicle (SV) -2 #24. This is henceforth denoted as DAK.2014.335.22. 3.53 3.54 3.55 3.56 GPS24, with 335 referring to the day of the year. Second, 10 5 a capture from the station in Kiruna, Sweden (lat. 67.743, Detrended power long. 21.060°) as representative of high-latitude scintillation, 5 taken on April 16, 2015 from 04:00 to 05:00 h local time 0 (dB) (UTC + 2) for GPS SV #22. This is henceforth denoted as -5 KIR.2015.106.04.GPS22. In both cases, the setup for scin- -10 tillation monitoring is formed by a high-grade antenna fol- -15 lowed by a Septentrio PolaRxS receiver providing binary 3.53 3.54 3.55 3.56 files with raw data and scintillation parameters at a rate of TOW (s) 10 5 50 Hz. These files are converted to raw American Standard Code for Information Interchange (ASCII) files contain- Fig. 3  Detrended phase and power of KIR.2015.106.04.GPS22 cap- ing the in-phase and quadrature (IQ) correlation samples, ture, representative of high-latitude scintillation. The focus is placed on the segment delimited within vertical stripped lines from which the scintillation power and phase time series are extracted following the data-detrending algorithm in Desh- pande et al. (2012). We also use the algorithm to eliminate index usually employed in the literature, defined as (Hum- residual low-frequency systematic effects initially affecting phreys et al. 2010), the captures, such as those from the troposphere, satellite √ √ [ geometry and receiver oscillator instabilities, which intro- √ E A2 (n)] − E[A (n)]2 √ S4 = √ s duce undesired fluctuations that add to the scintillation data s [ ]2 (10) and distort the measurements. E As (n) The look of the detrended phase and power is shown in Figs. 2 and 3 for the Dakar and Kiruna captures, respectively. with E[⋅] the expectation operator. On the other hand, phase We will focus on the segments delimited within the vertical scintillation is characterized through the standard devia- stripped lines as scintillation is more prominent there, with tion index 𝜎𝜃s of the phase series. At low latitudes, phase satellite elevation angles above 15° to avoid any propagation scintillation is correlated with amplitude scintillation, degradation effect owing to low satellite elevation. On the one with 𝜎𝜃s = 0.35 rad and power drops around 5–10 dB with hand, in order to characterize amplitude scintillation, S4 is the S4 = 0.3, being overall representative of mild scintillation. 13 17 Page 6 of 14 GPS Solutions (2022) 26:17 Table 1  Optimal AR parameters DAK.2014.335.22.GPS24 KIR.2015.106.04.GPS22 for phase scintillation time series 𝛼̂ p (unitless) 𝜎̂ s2 ­(rad2) 𝛼̂ p (unitless) 𝜎̂ s2 ­(rad2) p p AR(1) 0.9501 1.8658 × ­10−3 0.9606 3.0462 × ­10−3 AR(2) [0.7210, 0.2316] 1.7372 × ­10−3 [0.7515, 0.2085] 2.9008 × ­10−3 AR(3) [0.7024, 0.1637, 0.0874] 1.7148 × ­10−3 [0.7504, 0.1588, 0.0511] 2.8904 × ­10−3 10 10 0 0 -10 -10 -20 -20 -30 -30 -40 -40 (dB) (dB) -50 -50 -60 -60 -70 Periodogram of Dakar capture -70 Periodogram of Kiruna capture -80 AR(1) estimated PSD -80 AR(1) estimated PSD -90 AR(2) estimated PSD -90 AR(2) estimated PSD -100 AR(3) estimated PSD -100 AR(3) estimated PSD -110 -110 10 -2 10 -1 10 0 10 -2 10 -1 10 0 (rad) (rad) Fig. 4  Periodogram of DAK.2014.335.22.GPS24 versus estimated Fig. 5  Periodogram of KIR.2015.106.04.GPS22 versus estimated PSD of optimal AR(1, 2, 3) models. Low-latitude scintillation PSD of optimal AR(1, 2, 3) models. High-latitude scintillation However, a different situation is observed at high lati- observed at low frequencies that are, though, mainly due tudes, where even larger phase scintillations are perceived, to the low-pass filtering steps involved during the detrend- with 𝜎𝜃s = 0.65 rad, but no major power drops occur, with ing process. Note that the interest lies in the high-fre- S4 = 0.05, thus being representative of strong scintillation quency variations introduced by scintillation (Deshpande events in terms of phase. The behavior in both cases is in et al. 2012). In that sense, an AR process with any order p line with what is reported in the literature for low- and high- can be used in both captures, thus suggesting that a good latitude scintillation (Jiao et al. 2013; Lee et al. 2017). 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 Suitability of phase scintillation AR fitting errors are very similar for all orders. Therefore, these results show the suitability of modeling phase scintillation In order to evaluate how well the above live captures of events using an AR(1) process. phase scintillation can be modeled through an AR(p) pro- cess, the AR parameters that best fit the {time} series must p be first found, namely the coefficients 𝛼m m=1 and the prediction error power 𝜎s . By applying either the 2 Kalman filter equations p Yule–Walker equations or the least squares method (Kay 1993; Stoica and Moses 2005), the obtained results are As previously explained, there is a growing trend to replace shown in Table 1 for the Dakar and Kiruna data captures conventional PLL architectures (in particular the loop filter) considering AR(p) models with p up to three. These values by a KF-based one. Therefore, the objective of this section are used to determine the corresponding PSDs using (9), is to briefly present the generic KF equations, also termed which are then compared to the periodogram of the meas- state-space and observation models, which will be used later urements as a nonparametric estimator of the true PSDs. to formulate the KF-AR with the kinematic and AR signal The results are shown in Figs. 4 and 5 for each capture. A models introduced earlier to track carrier phase dynamics very tight match is obtained, despite some discrepancies and phase scintillation. 13 GPS Solutions (2022) 26:17 Page 7 of 14 17 The sequential estimation of the parameters of interest is [ ]T 𝐱d (n) ≐ 𝜃d (n) Ts 𝜃̇ d (n) Ts2 𝜃̈d (n) (14) governed by the following state transition equation (Yang et al. 2017), ⎡1 1 1∕2 ⎤ 𝐱(n + 1) = 𝐅𝐱(n) + 𝐆𝐯(n) (11) 𝐅d ≐ ⎢ 0 1 1 ⎥ (15) ⎢ ⎥ where 𝐱(n) is the state vector containing such parameters of ⎣0 0 1 ⎦ interest at time instant n. It propagates toward the next time ( ) The process noise is given by u(n) ≐ Ts3 𝜃⃛d (n) ∼ N 0, 𝜎u2 , instant n + 1 through the transition matrix F, plus the process which is intended to cover the higher-order terms that are noise 𝐯(n) that accounts for the possible modeling discrep- missing in the truncated model[ in (6). It propagates ]T to the ancy of the Kalman state-space model to the true one charac- ( 2) Kalman states through 𝐆d ≐ 1∕6 1∕2 1 that leads to terizing the[input samples. ] It usually follows 𝐯(n) ∼ N 0, 𝛔𝐯 𝐐d ≐ 𝜎u2 𝐆d 𝐆H . such that E 𝐯 (i)𝐯(j) = 0 for i ≠ j . As observed in (11), the H d On the other hand, the recursion in (8) for an AR(1) weight of 𝐯(n) onto the different Kalman states is adjusted process becomes readily the Kalman state transition equa- through a linear transformation given by matrix G that tion for the problem of phase scintillation tracking, where results in a zero-mean [ Gaussian ] process with the covari- the only state 𝜃s (n) propagates through 𝛼 , and the process ance matrix 𝐐 ≐ E 𝐆𝐯𝐯H 𝐆H . noise is now given by the AR driving noise sp (n), which The KF aims at providing clean estimates of 𝐱(n) based propagates entirely to the state. Then, by merging the kin- on the available noisy input measurements, which are mod- ematic and AR models, the KF-AR(1) state-space system eled through the following measurement equation (Yang model results into the following augmented state transition et al. 2017), equation, z(n) = 𝐇T 𝐱(n) + w(n) (12) [ ] [ ][ ] [ ] 𝐱d (n) 𝐅d 𝟎3×1 𝐱d (n − 1) 𝐆d 𝟎3×1 = + 𝐯(n) where the scalar measurements z(n) can be understood 𝜃s (n) 𝟎1×3 𝛼 𝜃s (n − 1) 0 1 as a linear transformation of the system states through (16) the observation matrix H. The measurements are cor- where the state vector and transition matrix are identified as, rupted by the presence of an [uncorrelated ] Gaussian noise [ ]T w(n) ∼ N(0, R(n)) such that E w(i)w (j) = 0 for i ≠ j. ∗ 𝐱(n) ≐ 𝐱dT (n) 𝜃s (n) (17) [ ] 𝐅d 𝟎3×1 Hybrid autoregressive Kalman filter (KF‑AR) 𝐅≐ 𝟎1×3 𝛼 (18) This section is intended to show the integration of the kin- The process noise becomes, ematic and AR signal models introduced earlier into the {[ ]T Kalman equations, leading to the hybrid autoregressive u(n) sp (n) if p = 1 𝐯(n) ≐ (19) Kalman filter (KF-AR), and present low-complexity imple- u(n) if p = 0 mentation proposed in this work. with p = 0 referring to a dynamics-only, non-AR KF that will be used by the adaptive technique in the next section. Baseline technique: the KF‑AR(1): state‑space The noise in (19) propagates to the Kalman [ states ]through and observation models the augmented process noise matrix 𝐆 ≐ 𝐆d 𝟎3×1 , result- 0 1 Recalling the kinematic signal model in (6)–(7), the states of ing in[a zero-mean Gaussian ] noise with covariance matrix interest for a three-dimensional KF are 𝜃d (n), 𝜃̇ d (n) and 𝜃̈d (n). 𝜎u2 𝐆d 𝐆H 𝟎 𝜎s2 . d 3×1 𝐐≐ Hence, the state-space model for carrier dynamics tracking 𝟎1×3 p can be written in normalized matrix notation as, When dealing with both phase dynamics and phase scin- tillation, we assume to have the following measurements ⎡ 𝜃d (n) ⎤ ⎡ 1 1 1∕2 ⎤⎡ 𝜃d (n − 1) ⎤ ⎡ 1∕6 ⎤ ⎢ Ts 𝜃̇ d (n) ⎥ = ⎢ 0 1 1 ⎥⎢ Ts 𝜃̇ d (n − 1) ⎥ + ⎢ 1∕2 ⎥Ts3 𝜃⃛d (n − 1) available, ⎢ 2̈ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎣ Ts 𝜃d (n) ⎦ ⎣ 0 0 1 ⎦⎣ Ts2 𝜃̈d (n − 1) ⎦ ⎣ 1 ⎦ z(n) = 𝐇T 𝐱 + w(n) = 𝜃d (n) + 𝜃s (n) + w(n) (20) (13) where by direct comparison with (11), the (3 × 1) state vec- where[ the augmented ]T observation matrix is identified as tor of the carrier phase dynamics and the transition matrix 𝐇 ≐ 1 0 0 1 . The measurement noise w(n) corresponds are identified, to the phase noise at the discriminator output. Assuming 13 17 Page 8 of 14 GPS Solutions (2022) 26:17 that the output of the prompt correlator is normalized to unit mean power, the variance of w(n) for the ATAN2 discrimina- tor is a scalar given by (Del Peral-Rosado et al. 2010), [ ] 1 1 R(n) = / 1+ / (21) 2Ts C N0 (n) 2Ts C N0 (n) 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/N0 (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/N0 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 bur- den. 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 eventu- ally 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 com- puting 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 imple- Fig. 6  Block diagram of AHL-KF-AR(0, 1) carrier tracking tech- mentation that uses a fixed set of AR model parameters, nique. The scintillation detector is used to set the AR model order to while partially preserving the other adaptive modules of the 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 AHL-KF-A2R(p): the AHL implementation remains as is, implementation 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 is compensated with the estimated scintillation phase (i.e., KF-AR(1) when scintillation is present. output of I&D block), the architecture implements a second The proposed technique corresponds to the carrier track- carrier discriminator that, since the carrier dynamics have ing architecture depicted in Fig. 6. As can be observed, the already been compensated, outputs only the input scintil- carrier discriminator output is the magnitude to be processed lation phase. This is used to detect the presence of scintil- by the KF, which corresponds to the innovation sequence lation to consequently enable or disable the AR module of i(n) that contains the error between the incoming phase the KF and update the Kalman gains K(n) and the transition and the Kalman predicted phase. Before the local replica matrix accordingly. The state correction is carried out, and 13 GPS Solutions (2022) 26:17 Page 9 of 14 17 the prediction for the next time instant is computed. The lat- 10 0 ter is then employed to obtain the measurement prediction No scint. and update the local replica, which is correlated with the 10 -1 =0.250 input signal, thus closing the loop. =0.500 =0.750 10 -2 =0.925 (rad 2) Feature 1: optimal fixed AR parameters Target BCRB 10 -3 2 Strong scint. s In order to determine the optimal AR parameters for the 10 -4 AHL-KF-AR(0, 1), the best achievable performance for each of the considered datasets must first be known. This 10 -5 information is readily provided by the BCRB, which can be recursively computed as (Van Trees and Bell 2007), 10 -6 [( ] 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 )−1 2 2 𝐉−1 B (n + 1) = 𝐐(n) + 𝐅(n)𝐉−1 T B (n)𝐅 (n) + R−1 (n)𝐇T (n)𝐇−1 (n) (rad ) s (22) with 𝐉B (n) the Bayesian information matrix. For the optimal Fig. 7  BCRB of KF-AR(1) as a function of AR model parameters. AR parameters in Table 1, and considering a high C/N0 of For 𝜎s2 fixed to track strong scintillation, the objective is to find the 𝛼 p 45 dB Hz, the BCRB for the phase dynamics of interest is that reaches the target BCRB located at 3.6 × ­10−3 and 1 × ­10−2 ­rad2 for the Dakar and Kiruna captures, respectively. We will take the smallest measurements. In that sense, determining the optimal order value to be our best possible target performance. That is to of a statistical model is a well-known problem in the field of say, we take 3.6 × ­10−3 ­rad2 as the MSE carrier phase perfor- signal processing often referred to as model-order selection mance to be attained by the AHL-KF-AR(0, 1). (Stoica and Selen 2004). Several criteria can be found in the In order to make the AHL-KF-AR(0, 1) robust, the opti- literature in this regard, such as Akaike’s information crite- mal AR parameter determination is carried out with strong rion (AIC) (Akaike 1974), the modified AIC, the generalized phase scintillation in the spotlight, such as the one from the information criterion (GIC) and the minimum description Kiruna capture analyzed herein. That is, the objective is to length (MDL) (Djuric et al. 1999). The proposed technique design the parameters such that the resulting tracking implements the latter, as it is reported to be a consistent scheme is able to cope with strong scintillation while attain- criterion. It computes the optimal p as, ing the aforementioned target MSE of 3.6 × ­10−3 ­rad2. For ( ( ) ) this purpose, note that the KF-AR is a closed-loop architec- p̂ MDL = arg min N log 𝜎̂ s2 + p log (N) (23) ture that can thus be characterized by its equivalent loop p p bandwidth, which must be large enough for adequate signal with N the number of samples in the data segment under tracking. As a matter of fact, such bandwidth is driven by analysis. In the proposed AHL-KF-AR(0, 1), this estimator the KF-AR process noise variance (Jwo 2001), given by 𝜎s2 p is implemented in form of a sliding window that evaluates for the case of phase scintillation tracking. Therefore, as a segments with 5 s of data as a function of time, and detects first step we can simply take the value of 𝜎s2 in Table 1 whether scintillation is present or not by determining 𝜎̂ s2 for p p required to track the scintillation in Kiruna as reference, p = 0 and p = 1 using 𝛼 = 0.925 and applying (23). namely 𝜎s2 = 3 × ­10−3 ­rad2. Then, it remains to determine the p coefficient 𝛼 , which is sought as the value for which, with such 𝜎s2 , the KF-AR attains the above target performance. Feature 3: adaptive hard‑limited measurement noise p variance By evaluating the BCRB for different candidate values of 𝛼 as shown in Fig. 7 it is found that, for 𝜎s2 = 3 × ­10−3 ­rad2, the p The adaptive hard-limited (AHL) implementation is thought performance of 3.6 × ­10−3 ­rad2 is attained for 𝛼 = 0.925. to deal with the nonlinear amplitude fades introduced by Therefore, this is the set of AR parameters to be used by the scintillation, whereby the KF-AR is affected by abnormal AHL-KF-AR(0, 1). measurements that may compromise the linearity of the phase discriminator, thus degrading its performance and Feature 2: real‑time scintillation detection even driving the technique to lose lock. The AHL exploits the fact that the amplitude fades make themselves explicit at When designing an AR(p) process, the model order p plays C/N0 level. Interestingly, this information can be provided to an important role for optimal AR fitting to a given set of the KF-AR through the measurement noise variance R(n) as 13 17 Page 10 of 14 GPS Solutions (2022) 26:17 seen in (21), so that the filter be dynamically self-adjusted KF-AR(1) PLL (expected) to match the actual working conditions. Therefore, the AHL adapts R(n) based on actual C/N0 measurements, subject to 10 0 AHL-KF-A2R(p) Target BCRB AHL-KF-AR(0,1) a hard-limiting threshold such that, { / ̂ � N0 (n) ≥ 𝛾 R(n) = R(n) if C (24) 10 -1 (rad) ∞ otherwise with ̂ ̂ / R(n) the measurement noise variance in (21) using C N0 (n), and where setting R(n) = ∞ drives the Kalman 10 -2 gains to zero whenever the C/N0 drops below a given thresh- old 𝛾 . 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 10 -3 0 100 200 300 400 500 600 the technique from abnormal measurements until the C/N0 [No - Yes] recovers to nominal conditions. 1 To estimate the C/N0, the new AHL-KF-AR(0, 1) imple- >99.99% 99.65% 93.68% 0 ments the narrow-wideband power ratio (NWPR) estimator, as it is a low-complexity algorithm of easy implementation 0 100 200 300 400 500 600 that provides reliable results (Van Dierendonck 1996). On Time (s) the other hand, taking into account that the C/N0 tracking threshold is located around 20 dB-Hz (Kaplan and Hegarty Fig. 8  Phase RMSE performance for DAK.2014.335.22.GPS24 (top), 2006), the threshold 𝛾 is set to 25 dB-Hz so as to leave some and success rate on detecting the presence or absence of scintillation margin for KF-AR operation at low C/N0. (bottom) Experimentation results and performance that sense, two sets of results are provided, namely for low evaluation and high latitudes. The metrics employed for performance evaluation are the root mean square error (RMSE) of the In this section, we assess the performance of fixed and estimated carrier phase dynamics, the number of cycle slips adaptive KF-AR techniques to validate the suitability and and the probability of loss of the tracking lock (LoL). The efficiency of the proposed approach for tracking carrier three metrics are computed by running the executions with phase dynamics under scintillation conditions. We evalu- 100 Monte Carlo realizations. ate the adaptability features of the techniques under time- varying scintillation working conditions and the closeness Experiment 1: KF‑AR performance for low‑ to the expected BCRB in comparison with the performance and high‑latitude scintillation provided by conventional PLLs. To this end, an experimen- tation campaign has been carried out using a MATLAB Figures 8 and 9 show the RMSE performance of the previ- software GNSS carrier tracking simulator for GPS L1 C/A ous KF-AR techniques and the new one proposed herein signals with a sampling time of Ts = 20 ms. The campaign for the Dakar and Kiruna data captures, respectively. Addi- simulates the environmental conditions in terms of signal tionally, Fig. 10 illustrates the error in cycles between the dynamics and receiver noise and adds the scintillation time estimated and true carrier phase of interest for some of the series captured in Dakar and Kiruna on top of these. We Monte Carlo realizations. consider a C/N0 of 45 dB Hz in a static receiver. The only As can be observed, all techniques outperform the RMSE experienced dynamics are due to the satellite motion, thus that would be obtained with a conventional PLL, as the lat- assuming a Doppler shift of 10 Hz, Doppler rate of 1 Hz/s ter estimates scintillation as part of the signal dynamics. and maximum Doppler acceleration of ± 2 × ­10−4 Hz/s2. The The PLL presents a LoL of about 20%, and it suffers from latter leads to 𝜎v2 = 3.4 × 10−17 ­rad2 as the variance of a uni- some tens of cycle slip occurrences over time, especially in form distribution delimited by such Doppler acceleration Kiruna, as shown in Fig. 10 top. More particularly, the mean values. The executions are 600 s long with a time-varying number of cycle slip occurrences in 300 s is found to be 45, presence of scintillation. Only AWGN is present during the meaning that a cycle slip is expected to occur approximately first 150 s, then scintillation appears between 150 and 450 s, every 6–7 s. On the contrary, KF-AR techniques are found to then scintillation disappears again during the last 150 s. In present no cycle slips, as shown in Fig. 10 bottom, and the 13 GPS Solutions (2022) 26:17 Page 11 of 14 17 10 1 through an AR process. Note that the ideal situation would KF-AR(1) PLL (expected) be to attain the BCRB obtained when scintillation is absent. AHL-KF-A2R(p) Target BCRB However, the KF-AR approach already outperforms the 0 AHL-KF-AR(0,1) conventional PLL by a factor of 6. The inconvenient of the 10 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 expe- (rad) 10 -1 rienced 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 scin- -2 10 tillation. 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 10 -3 the new condition. In that sense, the major advantage of 0 100 200 300 400 500 600 this technique is that it presents good agility and flexibility [No - Yes] 1 in front of time-varying conditions. Notwithstanding, the >99.99% 99.83% 91.88% main inconvenience of the AHL-KF-A2R(p) is that, besides 0 presenting a slightly larger LoL probability, it presents some RMSE gap that prevents attaining the BCRB by a factor 0 100 200 300 400 500 600 between 2 and 3. This is owing to the online AR param- Time (s) eter estimator providing biased estimates due to the pres- ence of AWGN. Furthermore, the computational burden of Fig. 9  Phase RMSE performance for KIR.2015.106.04.GPS22 (top), this technique is found to be 2.75 times that of the simplest and success rate on detecting the presence or absence of scintillation (bottom) 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) 10 when under AWGN only, and to the fixed KF-AR(1) when in PLL the presence of scintillation, thus attaining the correspond- (cycles) 5 0 ing BCRB values in both situations for both data captures. -5 The success rate on the detection of scintillation is shown -10 in the bottom plots of Figs. 8 and 9, with values above 90%, 0 100 200 300 400 500 600 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 compu- 10 KF-AR tational cost is reduced to 1.75 times that of the simplest (cycles) 5 0 KF-AR and the LoL probability is improved to below 1%, -5 thus potentially becoming one robust approach that can be -10 implemented in practice. 0 100 200 300 400 500 600 Time (s) Experiment 2: KF‑AR operational region The objective of this section is twofold. First, to extend the Fig. 10  PLL error between estimated and true carrier phase of inter- est (top), and KF-AR error (bottom). The figure reveals the presence analysis in the previous section and evaluate the KF-AR of cycle slips in the PLL performance as a function of the input C/N0 and the high- order dynamics that the KF can tolerate through 𝜎v2. Second, to quantify the limits of the proposed AHL-KF-AR(0, 1) LoL probability is below 1% as no Monte Carlo realization technique within which it provides a performance gain over has lost lock, except for the AHL-KF-A2R(p), which shows the conventional PLL and out of which there is no point in a LoL probability of 5% in Kiruna. using it over this latter. When scintillation is present, the fixed KF-AR(1) that By performing a dual-sweep analysis on both the above employs the optimal AR parameters in Table 1 can fully parameters, Fig. 11 shows a 3D plot of the AHL-KF-AR(0, reach the target performance given by the BCRB, thus 1) phase RMSE when considering low-latitude scintilla- confirming the validity of the modeling phase scintillation tion (Dakar). As can be observed, the RMSE increases 13 17 Page 12 of 14 GPS Solutions (2022) 26:17 Fig. 11  Phase RMSE of AHL-KF-AR(0, 1) as a function of C/N0 and Fig. 12  Phase RMSE of AHL-KF-AR(0, 1) as a function of C/N0 and 𝜎v2 for DAK.2014.335.22.GPS24. Comparison to expected PLL per- 𝜎v2 for KIR.2015.106.04.GPS22. Comparison to expected PLL perfor- formance given by gray plane mance given by gray plane when decreasing the C/N0 or when enhancing the dynam- latitudes in Fig. 11. The technique outperforms the PLL ics through the Kalman model discrepancy 𝜎v2 . The latter for a lower-bounded C/N0 of 25 dB Hz, that is, the value is explained by the fact that, for larger 𝜎v2 , the Kalman to which the AHL threshold has been set. At this point, must increase its dependability on the input measure- 𝜎v2 is limited to around ­10 −20 ­rad 2, similarly to the case ments rather than on its internal state-space model. As at low latitudes, corresponding to a maximum Doppler this occurs, it becomes more difficult for the Kalman to acceleration of ± 3 × ­10–6 Hz/s2. For C/N0 of 30 dB Hz and discern the scintillation phase from the dynamics phase of above, the tolerated higher-order dynamics is found to be interest, thus eventually leading to a larger RMSE. up to 𝜎v2 = 10−10 ­rad2, corresponding to a maximum Dop- On the other hand, the gray plane in Fig. 11 stands for pler acceleration of ± 3 × ­10−1 Hz/s2. Beyond this point, the best performance that could be expected from PLLs. the AHL-KF-AR(0, 1) provides similar performance to the The performance limits of the AHL-KF-AR(0, 1) are given PLL. The above results confirm the suitability of KF-AR by the boundary lines given where both KF-AR and PLL techniques for precise localization and cruise navigation RMSE plots intersect. The region where the proposed purposes, and autonomous vehicle driving applications technique outperforms the PLL is bounded by a C/N0 in with mild acceleration variations. the region of 30 dB Hz, with 𝜎v2 restricted to 10−20 ­rad2, whereas the technique supports up to 𝜎v2 = 10−12 ­rad2 for C/N0 above 35 dB Hz. For larger 𝜎v2 the technique performs Conclusions very similarly to the PLL, whereby phase scintillation is estimated as part of the carrier dynamics of interest. These We have tackled the problem of robust GNSS carrier track- values can be translated into more friendly magnitudes ing in the presence of ionospheric scintillation utilizing such as the maximum supported Doppler acceleration. For hybrid autoregressive Kalman filters. Two main contribu- uniformly distributed process noise, both are related as: tions have been presented. First, we proposed the AHL-KF- √ AR(0, 1) low-complexity technique to facilitate the imple- 12𝜎v2 mentation of the approach in practice. Second, we assessed ⃛ 𝜃dmax = ± (25) the performance of these techniques using representative 4𝜋Ts3 captures of real scintillation data at low and high latitudes, in Hz/s2. Therefore, in other words, the technique can toler- which are two different and complementary situations com- ate a Doppler acceleration of ± 3 × ­10−6 Hz/s2 at low C/N0, monly addressed in studies of scintillation disturbances. whereas it can support ± 3 × ­10−2 Hz/s2 at high C/N0. From the work performed, three main conclusions can be Figure 12 shows the 3D plot of the proposed AHL- drawn. First, the obtained results confirm the feasibility KF-AR(0, 1) phase RMSE when considering high-lati- of modeling the phase scintillations observed by the user tude scintillation (Kiruna). Here, it is remarkable to note through AR processes. Second, the proposed AHL-KF- that the outperformance region is larger than that at low AR(0,1) low-complexity technique has shown to clearly out- perform the conventional PLL architecture and other KF-AR 13 GPS Solutions (2022) 26:17 Page 13 of 14 17 implementations while conserving the optimality properties Deshpande K, Bust G, Clauer C, Kim H, Macon J, Humphreys T, of the Kalman filter. Third, the performance limits where the Bhatti J, Musko S, Crowley G, Weatherwax A (2012) Initial GPS scintillation results from CASES receiver at south pole, Antarc- technique provides a performance gain over the PLL have tica. Radio Sci 47(5):1–10 been quantified, thus becoming valuable guidelines during Djuric P, Kay S, Vijay K, Douglas B (1999) Spectrum estimation and the design of the tracking stage of GNSS receivers in real modeling. Digital signal processing handbook. CRC Press, Boca environments. Therefore, the work performed validates the Raton Fohlmeister F, Antreich F, Nossek J (2018) Dual Kalman filtering excellence of the proposed approach for robust GNSS carrier based GNSS phase tracking for scintillation mitigation. In: Pro- tracking under scintillation conditions. ceedings of the IEEE/ION position, location and navigation sym- posium (PLANS). pp 1151–1158 Acknowledgements This work has been supported in part by the Span- Gupta S (1975) Phase-locked loops. Proc IEEE 63(2):291–306 ish Ministry of Economy and Competitiveness under Grant TEC2017- Humphreys T, Psiaki M, Ledvina B, Cerruti A, Kintner P (2010) A 89925-R, in part by the ICREA Academia program and in part by the data-driven testbed for evaluating GPS carrier tracking loops FI-DGR PhD grant with file number 9015-376195/2014 from Gener- in ionospheric scintillation. IEEE Trans Aerosp Electron Syst alitat de Catalunya. The authors would also like to thank the European 46(4):1609–1623 Space Agency (ESA), and in particular Mr. Juan M. Parro, for pro- Jiang R, Wang K, Liu S, Li Y (2017) Performance analysis of a Kalman viding the Monitor scintillation time series employed to perform the filter carrier phase tracking loop. GPS Solut 21:551–559 analyses addressed in this work. Jiao Y, Morton Y, Taylor S, Pelgrum W (2013) Characterization of high-latitude ionospheric scintillation of GPS signals. Radio Sci Funding Open Access Funding provided by Universitat Autonoma de 48(6):698–708 Barcelona. Jwo D (2001) Optimization and sensitivity analysis of GPS receiver tracking loops in dynamic environments. Proc IEEE Radar Sonar Navig 148(4):241–250 Data availability The scintillation time series analyzed in this work Kaplan E, Hegarty C (2006) Understanding GPS: principles and appli- have been obtained from the Monitor network of the European Space cations, 2nd edn. Artech House, Boston Agency (ESA), a project aimed at monitoring ionospheric scintillation Kay S (1993) Fundamentals of statistical signal processing: estimation events from multiple ground stations located worldwide. Access to theory, vol 1. Prentice Hall, Hoboken these time series is publicly available and can be granted by ESA upon Kintner P, Levidna B, De Paula E (2007) GPS and ionospheric scintil- request through the Monitor website (http://​monit​or.​estec.​esa.​int/). lations. Space Weather. https://​doi.​org/​10.​1029/​2006S​W0002​60 Kintner P, Humphreys T, Hinks J (2009) GNSS and ionospheric scin- Open Access This article is licensed under a Creative Commons Attri- tillation. How to survive the next solar maximum. Inside GNSS bution 4.0 International License, which permits use, sharing, adapta- 4(4):22–30 tion, distribution and reproduction in any medium or format, as long Lee J, Morton Y, Lee J, Moon H, Seo J (2017) Monitoring and mitiga- as you give appropriate credit to the original author(s) and the source, tion of ionospheric anomalies for GNSS-based safety critical sys- provide a link to the Creative Commons licence, and indicate if changes tems: a review of up-to-date signal processing techniques. IEEE were made. The images or other third party material in this article are Signal Process Mag 34(5):96–110 included in the article's Creative Commons licence, unless indicated Locubiche-Serra S, Seco-Granados G, López-Salcedo J (2016) Doubly- otherwise in a credit line to the material. If material is not included in adaptive autoregressive Kalman filter for GNSS carrier tracking the article's Creative Commons licence and your intended use is not under scintillation conditions. In: Proceedings of the international permitted by statutory regulation or exceeds the permitted use, you will conference on localization and GNSS (ICL-GNSS). pp 1–6 need to obtain permission directly from the copyright holder. To view a Macabiau C, Deambrogio L, Barreau V, Vigneau W, Valette JJ, Artaud copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. G, Thevenon P, Ries L (2012) Kalman filter-based robust GNSS signal tracking algorithm in the presence of ionospheric scintilla- tions. In: Proceedings of the 25th ION ITM. pp 1–8 Morton Y, van Diggelen F, Spilker J, Parkinson B, Lo S, Gao G (2020) References Position, navigation, and timing technologies in the 21st century, volumes 1 and 2: integrated satellite navigation, sensor systems, and civil applications. Wiley, Hoboken Aarons J (1982) Global morphology of ionospheric scintillations. Proc Nunes F, Sousa F (2014) Practical simulation of GNSS signals in the IEEE 70(4):360–378 presence of ionospheric scintillation. In: IEEE/ION position, loca- Akaike H (1974) A new look at the statistical model identification. tion and navigation symposium (PLANS). pp 50–58 IEEE Trans Autom Control 19(6):716–723 Proakis J (2001) Digital communications, 4th edn. McGraw-Hill, New Bar-Shalom Y, Li X, Kirubarajan T (2004) Estimation with applica- York tions to tracking and navigation: theory algorithms and software. Seco-Granados G, López-Salcedo J, Jiménez-Baños D, López-Risueño Wiley, Hoboken G (2012) Challenges in indoor global navigation satellite systems: Curran J, Lachapelle G, Murphy C (2012) Improving the design of unveiling its core features in signal processing. IEEE Signal Pro- frequency lock loops for GNSS receivers. IEEE Trans Aerosp cess Mag 29(2):108–131 Electron Syst 48(1):850–868 Stoica P, Moses R (2005) Spectral analysis of signals. Prentice Hall, Del Peral-Rosado J, López-Salcedo J, Seco-Granados G, López- Hoboken Almansa J, Cosmen J (2010) Kalman filter-based architecture Stoica P, Selen Y (2004) Model-order selection. A review of informa- for robust and high-sensitivity tracking in GNSS receivers. In: tion criterion rules. IEEE Signal Process Mag 21(4):36–47 Proceedings of the 5th ESA workshop on satellite navigation tech- Susi M, Andreotti M, Aquino M, Dodson A (2017) Tuning a Kalman nologies and European workshop on GNSS signals and signal filter carrier tracking algorithm in the presence of ionospheric processing (NAVITEC). pp 1–8 scintillation. GPS Solut 21:1149–1160 13 17 Page 14 of 14 GPS Solutions (2022) 26:17 Van Dierendonck A (1996) GPS receivers. In: Bradford W, Spilker J working as a post-doctoral researcher and project engineer at UAB, (eds) Global positioning system theory and applications. Ameri- and he is also affiliated with the Institute of Space Studies of Catalonia. can Institute of Aeronautics and Astronautics, Reston He held a visiting appointment at the European Space Agency (ESA) Van Trees H, Bell K (2007) Bayesian bounds for parameter estimation in 2016. With more than 7 years’ experience in GNSS technologies, and nonlinear filtering/tracking. Wiley, Hoboken he has been involved in numerous projects funded by ESA and other Vilà-Valls J, López-Salcedo J, Seco-Granados G (2013) An interactive national and international space-related entities on advanced position- multiple model approach for robust GNSS carrier phase track- ing solutions for next-generation GNSS receivers. Since 2018 he has ing under scintillation conditions. In: IEEE international confer- been involved in the design of the signal steering algorithm for the ence on acoustics, speech and signal processing (ICASSP). pp EGNOS V3 system. 6392–6396 Won J, Pany T, Eissfeller B (2012) Characteristics of Kalman filters Gonzalo Seco‑Granados received the Ph.D. degree in Telecommu- for GNSS signal tracking loop. IEEE Trans Aerosp Electron Syst nications Engineering from the Universitat Politècnica de Catalunya, 48(4):3671–3681 in 2000, and the MBA degree from IESE Business School in 2002. Yang R, Ling K, Poh E, Morton Y (2017) Generalized GNSS sig- Until 2005, he was with the European Space Agency, involved in the nal carrier tracking: part I—modeling and analysis. IEEE Trans design of the Galileo system. Since 2006, he has been with Universitat Aerosp Electron Syst 53(4):1781–1797 Autònoma de Barcelona and affiliated with the Institute of Space Stud- ies of Catalonia. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. José A. López‑Salcedo received his M.Sc. and Ph.D. degrees in Tel- ecommunications Engineering from the Universitat Politècnica de Catalunya, in 2001 and 2007, respectively. Since 2006, he is with Uni- versitat Autònoma de Barcelona. He is also affiliated with the Institute Sergi Locubiche‑Serra received his M.Sc. and Ph.D. degrees in Tel- of Space Studies of Catalonia and held a visiting appointment at the ecommunications Engineering from the Universitat Autònoma de European Commission Joint Research Center in 2014. Barcelona (UAB), in 2014 and 2019, respectively. He is currently 13