A statistical approach to estimate Global Navigation Satellite Systems (GNSS) receiver signal tracking performance in the presence of ionospheric scintillation Application to Space Weather Analysis, Modelling,

– Ionospheric scintillation can seriously impair the Global Navigation Satellite Systems (GNSS) receiver signal tracking performance, thus affecting the required levels of availability, accuracy and integ-rity of positioning that supports modern day GNSS based applications. We present results from the research work carried out under the Horizon 2020 European Commission (EC) funded Ionospheric Prediction Service (IPS) project. The statistical models developed to estimate the standard deviation of the receiver Phase Locked Loop (PLL) tracking jitter on the Global Positioning System (GPS) L1 frequency as a function of scintillation levels are presented. The models were developed following the statistical approach of generalized linear modelling on data recorded by networks in operation at high and low latitudes during the years of 2012–2015. The developed models were validated using data from different stations over varying latitudes, which yielded promising results. In the case of mid-latitudes, as the occurrence of strong scintillation is absent, an attempt to develop a dedicated model proved fruitless and, therefore, the models developed for the high and low latitudes were tested for two mid-latitude stations. The developed statistical models can be used to generate receiver tracking jitter maps over a region, providing users with the expected tracking conditions. The approach followed for the development of these models for the GPS L1 frequency can be used as a blueprint for the development of similar models for other GNSS frequencies, which will be the subject of follow on research.


Introduction
The effect of the Earth's ionospheric environment currently represents the single largest contributor to the Global Navigation Satellite System (GNSS) error budget and abnormal ionospheric conditions can impose serious degradation on GNSS system functionality, including integrity, accuracy and availability. In particular, small-scale time varying plasma density irregularities in the ionosphere can introduce fluctuations in the amplitude and phase of the received GNSS signal, a phenomena known as scintillation (Kintner et al., 2001). Scintillation occurrence shows large day-to-day variability with dependence on local time, season, latitude, longitude, as well as solar and geomagnetic activity. The global morphology of ionospheric L-band scintillation occurrence during the solar maximum and minimum conditions is well known (Basu et al., 2002;Kintner et al., 2007;Alfonsi et al., 2011;Prikryl et al., 2015). It has been observed that the two regions where scintillation occurs more predominantly are the equatorial bands extending from about 20°N to 20°S geomagnetic latitudes and the high latitude regions extending from 65°to 90°geomagnetic latitudes. However, in these two regions, the processes that govern the generation of the irregularities causing scintillation are quite different, thereby leading to significant differences in the characteristics of the observed scintillation effects (Jiao & Morton, 2015).
Global Navigation Satellite System (GNSS) receivers need robust and continuous tracking of the satellite signals in order to compute the distance travelled by the incoming signal as accurately as possible. Demodulation of the navigation data from the incoming signal requires an exact replica of the signal, which can be achieved by tracking its phase through a Phase Locked Loop (PLL). However, ionospheric scintillation can degrade the GNSS receiver signal tracking performance. Rapid phase fluctuations associated with scintillation may cause the receiver PLL to lock onto a wrong phase while still tracking the signal, or even to lose lock completely (Humphreys et al., 2005). The receiver PLL performance is usually evaluated in terms of the phase tracking error variance, which is observed to correlate well with the scintillation levels at both high and low latitudes (Sreeja et al., 2012;Aquino & Sreeja, 2013). Although some high grade receivers, e.g. with an adaptive tracking loop, may be able to maintain track during moderate to strong scintillation conditions, the errors due to scintillation still propagate to the GNSS measurements leading to degradation in the positioning accuracy (Marques et al., 2018). Several efforts have been made to model the effects of scintillation on the Global Positioning System (GPS) receiver PLL performance. These can be summarized as: (i) development of scintillation sensitive PLL tracking models (Conker et al., 2003;Moraes et al., 2014), (ii) evaluation of PLL performance using a GPS signal simulator (Morrissey et al., 2004) and (iii) observation of receiver performance during actual events of scintillation in locations of interest (Groves et al., 2000).
With the growing reliance on GNSS for several modern life applications, including those demanding high positioning accuracy, precise ionospheric forecasts can contribute to the understanding and mitigation of the impact of significant ionospheric related events on our technology-based society. In this context, the research presented here was carried out under the ''Ionosphere Prediction Service'' (IPS) project funded by the European Commission (EC) under Horizon 2020, in the frame of the Galileo programme. The project is led by Telespazio (Italy) in collaboration with University of Nottingham (UoN, United Kingdom), Istituto Nazionale di Geofisica e Vulcanolgia (INGV, Italy), the University of Tor Vergata (UTOV, Italy), Nottingham Scientific Ltd (NSL, United Kingdom) and Telespazio VEGA (Germany). The main aim of this project is to translate the state of the ionosphere into GNSS user-devoted metrics through the design and development of an ionospheric prediction service prototype. The user-devoted products are tailored to meet the needs of the different communities, namely aviation, mass market, high accuracy and critical infrastructures monitoring, to which the service is targeted. The core scientific contribution of this project is represented by the research activities carried out by the project's research partners, namely UTOV, INGV and UoN, dealing with different topics of Space Weather phenomena from the Sun to the ionosphere affecting GNSS service providers and the user community.
In this framework, statistical models have been developed for the first time by UoN to estimate the standard deviation of the receiver PLL tracking jitter on the GPS L1 frequency separately at high and low latitudes as a function of the scintillation levels -these are presented and discussed. The research does not include the estimation of the scintillation levels, which is a topic being dealt with separately in the IPS project. Therefore the statistical models are developed based on real longterm data that is representative of historical scintillation levels actually observed in the different regions of the globe addressed in the project. The considerable size of the data set (12744760 (2.14 GB) and 11790943 (2.05 GB) data points respectively for the high and low latitudes) and the varied levels of scintillation that it covers provide a large and representative sample to ensure the statistical robustness of the models, which are validated in the research. The data and methodology are described in Section 2, followed by results and discussions in Section 3. Finally, the conclusions are presented in Section 4.

Data and methodology
The data collected by networks of specialised Ionospheric Scintillation Monitoring Receivers (ISMRs) that were in operation at the high and low latitudes during the years of 2012-2015 have been exploited to develop the statistical models to estimate the standard deviation of the receiver PLL tracking jitter, r PLL , on the GPS L1 frequency. In the development of these models, only measurements from satellites with an elevation angle greater than 20°were considered, in order to remove the contribution from non-scintillation related effects, such as multipath. Also, a minimum lock time threshold of 240 s was used to allow the convergence of the phase detrending filter. Two types of ISMRs have been used in this analysis, namely the Novatel GSV4004 and the Septentrio PolaRxS. The geographic latitude and longitude of the stations along with the GNSS receiver type used in the analyses are shown in Table 1. The data availability at the different stations are also shown in Table 1.
Statistical models have been developed to estimate r PLL as a function of scintillation indices and ROT rms , the latter defined as the root mean square of the Rate of change of slant Total Electron Content (STEC). Both types of ISMR used in this study use similar algorithms to provide the scintillation indices, namely the amplitude scintillation index, S4 and the phase scintillation index, r u . The analyses presented in Sreeja et al. (2011a) show that the scintillation indices recorded by these two types of ISMRs are comparable.
The S4 is defined as the standard deviation of the received signal power normalized by its mean value and r / is defined as the standard deviation of the detrended carrier phase using a high pass filter with 0.1 Hz cut off computed over 60s. In this study, the one minute r / and S4 indices recorded on the GPS L1 frequency per satellite-receiver link are used. The values of the scintillation indices values typically varies between 0 and 1. The uncalibrated STEC and the differential TEC (dTEC) values computed every 15 s per satellite-receiver link are also recorded by both types of ISMR (Van Dierendonck et al., 1993;Septentrio PolaRxS application manual, 2010). The dTEC values are computed from the carrier phase measurements and refer to the change in TEC over the four 15 s intervals during the last minute. The ROT rms per satellite-receiver link over a minute used in this study is computed by calculating the root mean square (rms) of the four dTEC values provided by the ISMR over each minute. The development of the models based on ROT rms will support the estimation of r PLL over regions where the scintillation indices or the ISMR data is unavailable, as ROT rms can be computed from uncalibrated STEC values that can be estimated using data from any dualfrequency GNSS receiver. The ROT estimation will remove any biases in the STEC measurements due to the differencing of the STEC values over time. The range of ROT rms values used in the model development is between 0 and 5.
As the range of validity of the developed models is driven by the model input data, the distribution of the one minute scintillation indices and ROT rms values covered by the four year long data used in the model development is shown in Table 2.
It can be observed from Table 2 that there are no major inherent limitations or bias in the input data, as all levels of scintillation are represented well. Therefore, the developed statistical models can be used under different scintillation conditions. Furthermore, the models are more likely to perform well during weak and moderate scintillation levels characterized by S4/r / < 0.7 and ROT rms < 2, where there is an abundance of input data.
The r PLL at high latitudes was computed using the scintillation sensitive PLL tracking model of Conker et al. (2003), hereafter referred to as the Conker model. However, this model is limited to weak-to-moderate levels of scintillation, S4 (L1) < 0.707, and hence cannot be applied for all levels of scintillation. This is particularly relevant for the equatorial and low latitudes, where very strong scintillation conditions are frequently encountered, with S4 (L1) reaching values greater than 0.8. To overcome this limitation and to enable the computation of r PLL over the equatorial and low latitudes, even under strong levels of scintillation, the PLL tracking model proposed in Moraes et al. (2014), hereafter referred to as the a-l model, was exploited. The parameters required to compute r PLL using either of these models are recorded by the ISMRs. The Conker model assumes a Nakagami-m distribution for the statistical characterization of amplitude scintillation, thus leading to the limitations in the model. The a-l model is an extension of the Conker model and is based on the a-l distribution, a more realistic distribution for amplitude scintillation (Moraes et al., 2013). The a and l coefficients can be computed either from the field recorded data, provided high rate 50 Hz data is available, or from an established empirical relationship between a and S4 (Moraes et al., 2014).
Statistical methods in the form of polynomial regression analysis have been applied previously to analyze the correlation of PLL tracking jitter with scintillation levels (Sreeja et al., 2012;Aquino & Sreeja, 2013). However, these correlation studies were limited to small data sets, usually from a limited number of individual days, not allowing for any attempt to develop robust statistical models that could be reliably used to estimate the receiver signal tracking performance as a function of scintillation levels. In this study, for the first time the statistical approach of Generalized Linear Modelling (GLM) was applied on a long-term, four years long data set, i.e. from 2012 to 2015, covering different seasons, solar and geomagnetic activity conditions, to develop models to estimate r PLL as a function of scintillation levels.
The GLM (McCullagh & Nelder, 1989) is a flexible generalization of linear regression models including response variables that do not necessarily have a Gaussian error distribution and uses the maximum likelihood estimation to estimate the model coefficients. The analysis of the error distribution curve for the response variable, r PLL in this case, showed that they follow a normal distribution, therefore this approach can be considered as linear regression, which is a special case of GLM. The standard errors of the statistical model coefficients are inversely proportional to the square root of the size of the dataset involved in the model development. In this study, using 4 years of data, the developed statistical models have smaller standard errors than previous models (Sreeja et al., 2012;Aquino & Sreeja, 2013) in the estimated model coefficients, thus resulting in an increased ability to estimate under various conditions. To determine the significance of each coefficient estimated in the model, a Matlab built in statistical test based on the t-statistic was used (Lehmann & Romano, 2005) at a 5% significance level, which is commonly used in this type of analysis.
For the validation of the developed models, a data sacrifice strategy was adopted, whereby selected data, covering scintillation and non-scintillation days during different seasons, has been left out of the development phase of the models. The criteria for categorizing a day as scintillation or non-scintillation was based on the temporal profiles of the phase and amplitude scintillation indices for the high and low latitudes respectively. The following criteria was applied to categorize the day as scintillation: (1) the threshold for amplitude and phase scintillation index is set as 0.3 as this study focuses on moderate to strong levels of scintillation; (2) the satellite elevation cut off is set at 20°; (3) the one minute scintillation index values remain above the threshold for more than 15 min; (4) for equatorial and low latitude stations, this 15 min interval refers to the post sunset period and for high latitudes, this 15 min interval refers to any time during the day. Days on which the scintillation index values were less than 0.3 for the whole day were categorized as non-scintillation days. In addition, manual visual inspections of the figures illustrating the temporal profiles of the scintillation indices were performed to confirm the selection of scintillation and non-scintillation days. The goodness of fit of the model estimations were evaluated by using the correlation coefficient, R, and the Root Mean

Results and discussions
The development of the statistical models following the approach of GLM to estimate r PLL at high and low latitudes is described in this section. The validation of the developed models is also described.

High latitude model
The data recorded at the high latitude stations, IQA and CBB in the Canadian High Arctic Ionospheric Network (Jayachandran et al., 2009) were used to develop the statistical models incorporating the dependence on the scintillation indices and the ROT rms . Scintillation at high latitudes frequently occurs in the cusp region with extension into dayside polar cap and in the nightside auroral oval (Prikryl et al., 2015). This was the rationale behind using data from stations CBB and IQA for model development, which are located in the cusp/polar cap and the auroral oval respectively.
To enable the computation of r PLL over the high latitudes, the Conker model has been exploited. Using the statistical test described in the previous section, the models developed for the high latitudes showed a dependence of r PLL on the phase scintillation index, r u as well as on the ROT rms . The models based on r u and ROT rms developed using data from IQA have the following forms, respectively: r PLL mm ð Þ ¼ 3:0941 þ 0:1452 Â ROT rms À 0:0226 The models based on r u and ROT rms developed using data from CBB have the following forms, respectively: r PLL mm ð Þ ¼ 3:092 þ 0:1463 Â ROT rms À 0:0111 It can be observed from the above equations that the coefficients of the developed models based on data from IQA and CBB are quite close to one another, thus indicating that either set of models (Eqs. (1) and (2) or Eqs. (3) and (4)) would be able to represent well the r PLL over varying latitudes. Therefore, the models developed using data from IQA were chosen as the high latitude models and used in the subsequent analysis. For illustration, a sample result showing the validation of the model represented in equations (1) and (2) on 03 March 2012, a day selected using the data sacrifice strategy described in the previous section, at station IQA is shown respectively in the top and bottom panels of Figure 1. In this figure, the black symbols indicate the observed values of r PLL , whereas the red symbols indicate the model estimations.
The RMSE of the residuals estimated for this validation day was 0.0664 mm and 0.0504 mm and the correlation coefficient R was 0.91 and 0.60 respectively for the models represented in equations (1) and (2), i.e. for the models based on r u and ROT rms respectively. This indicated that the models were able to correctly estimate both the values and the variability of the actual observations. The results from the days selected using the data sacrifice strategy described in the previous section for the validation of the models represented in (1) and (2) are summarized in Table 3. The table shows the average correlation coefficient R and the average RMSE in mm along with the number of days used for validation.
To further assess the suitability of the models represented in equations (1) and (2) to estimate r PLL at high latitudes, these models were also tested at the European high latitude stations of BRO and NYA, which did not contribute any data to the model development. The results from the validation days are shown in Table 4.
The fact that the choice of the cut off frequency used in detrending the phase can impact the estimation of the phase scintillation index, r u is known for a number of years (Forte & Radicella, 2002;Mushini et al., 2012). The 0.1 Hz cutoff frequency is widely implemented in ISMRs used for scintillation characterization. The results presented in this manuscript are based on data recorded by ISMRs. The impact of the results presented in Wang et al. (2018) could be on the developed high latitude models based on r u (Eqs. (1) and (3)). The assessment of the impact can only be carried out in a future study when a more definitive choice on the detrending cut off frequency is widely accepted by the community and therefore is now beyond the scope of this manuscript. The evidence presented in Wang et al. (2018) will not have an impact on the models based on ROT rms , as ROT rms is computed from uncalibrated STEC values estimated using data from any dual-frequency GNSS receiver.

Low latitude model
The data recorded at station PRU in Brazil was used to develop the statistical model for the low latitudes. This station was chosen because it lies close to the crest of the Equatorial Ionization Anomaly (EIA), where the occurrence of low latitude scintillation is frequent (Basu et al., 2002). To enable the computation of r PLL over the equatorial and low latitudes, the a-l model has been exploited. Based on the chosen statistical test described in the previous section, the models developed for the low latitudes showed a dependence on the amplitude scintillation index, S4 and on the ROT rms . These models have the following forms, respectively: The validation results of the models represented in equations (5) and (6) for the station PRU are summarized in Table 5.
As with the high latitude models, the suitability of the models represented in equations (5) and (6) were also tested at the Brazilian low latitude stations POA and PAL, which did not contribute any data to the model development. Table 6 summarizes the validation results for these two stations.
The validation results presented in Tables 3-6 indicate that both the developed high and low latitude models perform well under scintillation and non-scintillation conditions. The average values of R are mostly above 0.5, indicating correlation between the observed trend and the model estimated trend. The average RMSE values are all below 0.2 mm, which is one order of magnitude smaller than the precision of the GPS L1 carrier phase measurements of about 2 mm. This indicates that the values of the residuals are well below the L1 tracking noise level. Under these conditions, the value of the RMSE is the relevant criterion to evaluate the goodness of fit of the model estimations.
From the above tables, values of R observed at the low latitudes are on average lower than those observed at the high latitudes. This can be attributed to the fact that at low latitudes, during strong scintillation, there could also be fluctuations in the received signal to noise ratio (CN0), which have not been taken into account in the development of the above models. The statistical models for the estimation of r PLL have been developed exploiting their dependence only on the scintillation levels, as these are nowcasted by dedicated algorithms developed by INGV in the context of the IPS project. However, in reality, r PLL also has a dependence on the CN0, which is a receiver tracking dependent parameter. Clearly, the models can be further improved if the dependence on this parameter could be incorporated. However, incorporating the dependence on CN0 in the models is out of the scope of the IPS project, as nowcasted values of CN0 are not foreseen to be provided by the prototype.

Mid-latitude model
In the case of mid-latitude stations, the occurrence of strong scintillation is almost absent and hence r PLL values are very low. On the model development at these latitudes, the statistical test described in Section 2 indicated that none of the estimated coefficients were statistically significant, leading to conclude that any attempt using our approach to develop a model specific for the mid-latitudes would be meaningless. Nevertheless, during weak scintillation periods at the mid-latitude station of CYP, the variations in S4 are more pronounced than the variations in r u (Vadakke Veettil et al., 2017) Therefore, the models developed for the low latitudes were tested at this station. Concurrently, at the European mid-latitude station NOT, during weak scintillation periods the variations in r u are more pronounced than the variations in S4 (Aquino et al., 2005). Therefore, the models developed for the high latitudes were tested at this station. The results are summarized in Table 7.
The results shown in Table 7 indicate that the low and high latitude models can be used to estimate r PLL respectively at the mid-latitude stations CYP and NOT. It is to be noted that at both CYP and NOT stations, the occurrence of strong scintillation was observed to be almost absent, rendering impossible to categorize the days between scintillation and non-scintillation based on the criteria described in Section 2. Also, r PLL values are much smaller in these regions and therefore most of the time not likely to represent information that might be of relevance for users, as significant peaks are extremely rare. Additionally, it can be observed from Table 7 that the average RMSE values are less than 0.2 mm and the average R values are mostly above 0.5, thus indicating that the high and low latitude models may be used to estimate r PLL respectively at NOT and CYP. However, further analyses would be necessary to determine a latitude threshold to guide the decision as to what model to apply. If scintillation information can be nowcasted, the above developed models for the high and low latitudes can be used to nowcast r PLL on a regional or a global scale and thus used to generate tracking jitter maps (Sreeja et al., 2011b) over a region, which can provide GNSS users with the expected tracking conditions under scintillation. In the scope of the IPS project, partners INGV have developed algorithms to nowcast the scintillation and Proxy Scintillation Indices (PSI), which can then be used as inputs to the above developed models to generate r PLL maps. As an example, Figure 2 shows the nowcasted PSI values at a global scale on a regular grid of 2.5°· 5°respectively in latitude and longitude. The map is for 17 February 2018 at 05:15 UT. The PSI values in TECu/ min are defined in the color bar.
The nowcasted r PLL map, estimated on a global scale over the same grid size as that of Figure 2, using the models represented in equations (2) and (6), for the high and low latitudes respectively, is shown in Figure 3. r PLL values in mm are defined in the color bar.
It can be observed from Figures 2 and 3 that regions with higher values of PSI are in correspondence with the occurrence

Conclusion
The development and validation of statistical models to estimate the receiver r PLL for the GPS L1 signal developed by UoN under the EC funded IPS project have been presented. The GLM technique was applied to develop the statistical models to estimate r PLL as a function of the scintillation levels given by the conventional phase and amplitude scintillation indices and the ROT rms . Data collected by ISMR networks in operation at the high and low latitudes during the years of 2012-2015 was exploited to develop these models. The high latitude model showed a dependence on the phase scintillation index, whereas the low latitude model showed a dependence on the amplitude scintillation index. Also, statistical models were developed to incorporate the dependence of r PLL on ROT rms at both high and low latitudes. The validation of the developed models was carried out at different stations covering a range of latitudes, which yielded results that indicate the models' suitability to estimate r PLL over that range of latitudes. In the case of the mid-latitudes, as the occurrence of strong scintillation is almost absent, an attempt to develop a dedicated model proved fruitless. Nevertheless, to address the mid-latitudes, the models developed for the high and low latitudes were tested, respectively, at two stations in that region, namely Nottingham (53°N) and Cyprus (37°N). The rationale is that for these two stations the scintillation characteristics are somewhat similar to what is observed at high and low latitudes respectively. The results of the validation tests indicated that the developed high and low latitude models may be used to successfully estimate r PLL at mid-latitudes, but more detailed data analyses would be necessary to determine a latitude threshold to guide the decision as to what model to apply. The validation results for the high, mid and low latitude stations indicated average RMSE values below 0.2 mm and average R values above 0.5, thus suggesting that the developed models perform well under scintillation and non-scintillation conditions.