1 Introduction

When radio frequency (RF) signals pass through the plasma density irregularities in the ionosphere, their intensity and phase may suffer from rapid and random fluctuations, a phenomenon known as ionospheric scintillation. Scintillation has a higher probability of occurrence in the auroral to polar region and in the equatorial to low latitude region (Basu et al. 1988; Aarons 1982), roughly modulated by the solar and geomagnetic activities (Aarons et al. 1980; Aarons 1997). The characteristics of the ionospheric irregularities causing scintillation are different in these two regions. At high latitudes, scintillation is associated with the fast movement of large-scale ionospheric plasma structures along geomagnetic field lines (Basu et al. 2002; Forte and Radicella 2002) and is usually dominated by phase fluctuations (Jiao and Morton 2015). On the other hand, equatorial scintillation is associated with the post-sunset small-scale F-region irregularities and can be severe in both signal intensity and phase (Basu et al. 2002; Hysell and Kudeki 2004; Jiao and Morton 2015).

Signal quality and receiver performance of Global Navigation Satellite System (GNSS), such as Global Positioning System (GPS), GLONASS, Galileo and BeiDou Navigation Satellite System (BDS), can be adversely affected under scintillation (Skone et al. 2001; Chen et al. 2008; Sreeja et al. 2012), thereby leading to increased positioning errors. This usually leads to a worse accuracy of GNSS global broadcast ionospheric delay correction model for single-frequency users, which results in greater positioning error (Klobuchar 1987; Prieto-Cerdeira et al. 2014; Yuan et al. 2019). The effects of scintillation on GNSS positioning have been extensively studied. Based on the scintillation data recorded in the European Arctic region during the November 2004 storms, Meggs et al. (2008) showed that the rapid phase fluctuations caused by scintillation are highly correlated with losses of lock in GPS receivers. Jacobsen and Andalsvik (2016) analyzed the impact of ionospheric disturbances on the network real-time kinematic (RTK) positioning and Precise Point Positioning (PPP) techniques using the data collected in the Norway region during the 2015 St. Patrick’s day storm. Results show that the positioning errors of both RTK and PPP increase rapidly under ionospheric disturbances and that PPP can achieve a higher accuracy compared with the RTK positioning. Pi et al. (2017) calculated the kinematic PPP using the data collected by a scintillation monitoring receiver located in Alaska in 2013. It was pointed out that the dual-frequency PPP errors could be one to two orders of magnitude worse during high latitude scintillation. Both cycle slip and/or loss of carrier phase lock caused by scintillation can increase the positioning errors. Juan et al. (2018) investigated the scintillation effects on GNSS positioning by processing the data collected from 35 receivers at low and high latitudes. Results show that due to the difference in the characteristics of scintillation at low and high latitudes, the impacts of scintillation on GNSS positioning also differ. At high latitudes, cycle slips caused by scintillation are less frequent, thus it is possible to achieve precise positioning even during strong scintillation. By contrast, achieving precise positioning is more challenging in low latitude regions, as scintillation is more severe and produces more cycle slips in these regions.

Scintillation mitigation on GNSS positioning has been widely studied in recent years. Zhang et al. (2014) pointed out that the failure of cycle slip detection and abnormal blunders could lead to a sudden accuracy degradation under scintillation conditions. Thus, an improved cycle slip detection approach was proposed, which can help to avoid unnecessary reinitialization of phase ambiguity and to achieve an accuracy of 0.2–0.3 m in vertical direction under equatorial scintillation. A similar approach is proposed in Luo et al. (2019), where an improved cycle slip threshold model was developed aiming to decrease the frequent false alarm of cycle slip under scintillation. With the new model, the accuracy was improved by up to 47.9% in the up direction of BDS PPP. A new scintillation index is proposed in Juan et al. (2017) using the ionospheric-free (IF) combination of the phase measurements. An approach was developed to isolate the measurements suffering from cycle slips, thus the positioning errors under ionospheric disturbances are only due to the IF combination measurement noise.

The development of new GNSS constellations and signals provides the opportunity to mitigate scintillation effects on GNSS positioning. Marques et al. (2016) investigated the performance of using GPS L2C measurements in kinematic PPP calculation and compared with the solution using L2P under equatorial scintillation. Results show that using L2C against L2P can improve the positioning accuracy by up to 59% under weak scintillation, while no improvement is observed in the presence of strong scintillation. Furthermore, Marques et al. (2018) analyzed the integrated GPS and GLONASS kinematic PPP under different scintillation levels in the equatorial region. They suggested that compared with GPS-only PPP, a multi-GNSS solution improved the positioning accuracy by up to 60% under strong scintillation. Dabove et al. (2019) investigated the benefits of multi-GNSS positioning under high latitude ionospheric activities. It was found that adding GLONASS and Galileo satellites in positioning could improve the positioning accuracy and reduce the convergence time in PPP calculation.

To detect and remove the measurements affected by scintillation in Brazil, Bougard et al. (2013) applied the receiver autonomous integrity monitoring (RAIM) method in the PPP calculation. Results show that the method can help to mitigate the scintillation effects and achieve an accuracy level comparable to scintillation quiet conditions. Instead of removing the measurements affected by scintillation, a strategy to modify the least square stochastic models used in GNSS positioning was introduced in Aquino et al. (2009). The stochastic model was modified using the receiver phase-locked loop (PLL) and the delay-locked loop (DLL) tracking error variances, estimated using the scintillation sensitive tracking error models (Conker et al. 2003) with scintillation indices calculated at a rate of 1-min. In the modified stochastic models, the measurement noise levels are represented by the tracking error variances, thus measurements affected by scintillation will have higher tracking error variances and lower weights in the stochastic model, which results in a decreased contribution of the scintillating measurements to the position estimation. This mitigation approach has been applied successfully in many studies (da Silva et al. 2010; Strangeways et al. 2011; Park et al. 2017; Vani et al. 2019; Luo et al. 2020; Vadakke Veettil et al. 2020). Aquino et al. (2009), da Silva et al. (2010) and Strangeways et al. (2011) validated the mitigation approach by processing the scintillation data collected in northern Europe. Either relative positioning of double difference with code and phase measurements or point positioning with only code measurements was performed. PPP, which relies on the more precise carrier phase measurements, however, has not been performed. In the later studies of Vani et al. (2019), Luo et al. (2020) and Vadakke Veettil et al. (2020), PPP was performed following the mitigation approach under low latitude scintillation. But due to fact that high and low latitude scintillation have different characteristics, their effects on GNSS positioning are also different. Therefore, it is essential to validate the performance of the mitigation approach using PPP under high latitude scintillation. Another issue is that in the approach of Aquino et al. (2009), 1-min tracking error variances are used to represent the noise level of GNSS instantaneous measurements. The analysis in this study indicates that the signal fluctuations caused by scintillation vary significantly and presents large second-to-second variability within 1 min. Therefore, a 1-min tracking error variance may be a bad representative of the instantaneous measurement noise. 1-s scintillation indices are necessary to be developed to better describe the signal fluctuation under scintillation. A 1-s index was presented in Park et al. (2017) to describe the signal intensity fluctuation under low latitude scintillation. However, the estimation of the proposed 1-s index still requires 1-min scintillation data. 1-s indices for phase fluctuation are also not attempted, which is particularly essential at high latitudes.

To address the issues above, a novel approach is proposed to calculate 1-s scintillation indices based on 1-s scintillation data. Following the scintillation mitigation approach introduced in Aquino et al. (2009), the high latitude scintillation effects on GPS PPP are mitigated by exploiting the estimated 1-s scintillation indices and tracking error variances. The performance of mitigation approaches, exploiting the newly proposed 1-s scintillation indices and the commonly used 1-min indices, are also compared. The calculation of the commonly used 1-min scintillation indices and the datasets analyzed in this study are given next, followed by the details of the proposed approach for the estimation of the 1-s scintillation indices. The receiver PLL and DLL tracking error variance calculation accounting for scintillation are introduced thereafter. Scintillation effects on GPS PPP calculation are shown subsequently, followed by the results of scintillation mitigation on GPS positioning. Conclusions and remarks are presented finally.

2 Scintillation indices

Scintillation is classified as amplitude and phase scintillation, referring to the fluctuations in signal power and phase, respectively. The intensity of amplitude scintillation is characterized by S4, defined as the standard deviation of the detrended signal power normalized by its mean over 60 s (Van Dierendonck et al. 1993; Van Dierendonck 1999). Phase scintillation is characterized by \({{Phi60}}\), which is the standard deviation of the detrended carrier phase calculated in 60 s (Van Dierendonck 1999). One of the aims of the signal detrending is to minimize the low-frequency variations in signal power and phase measurements caused by non-scintillation effects, such as satellite motion, antenna patterns and multipath (Van Dierendonck et al. 1993). It is interesting to note that due to the changes in receiver ambient temperature, short-term temporal variations in code and phase biases in GNSS measurements are observed, which is discussed in Zhang et al. (2017, 2019). However, due to the fact that this variability is a low-frequency variation which can last over a few hours, it can be removed in the data detrending process and has minor effects on the analysis in this work. A detailed description of the signal power and phase measurement detrending can be found in Van Dierendonck et al. (1993).

Spectral parameters are also calculated to describe the characteristics of phase scintillation. The power spectral density (PSD) of phase scintillation is modeled as (Rino 1979)

$$ S_{\phi } \left( f \right) = \frac{T}{{\left[ {f^{2} + f_{0}^{2} } \right]^{p/2} }}, $$
(1)

where \(f_{0}\) is the frequency corresponding to the outer scale size of irregularities. \(f\) is the frequency component composing the detrended phase measurements. \(p\) is the spectral index corresponding to the opposite of the PSD slope in log–log axis. \(T\) is the spectral strength, which is the value of phase PSD at 1 Hz. In practical cases, \(f \gg f_{0}\), Eq. (1) is simplified as \(S_{\phi } \left( f \right) \approx Tf^{ - p}\).

Phase variance is related to \(p\) and \(T\) by (Rino 1979; Conker et al. 2003; Aquino et al. 2007)

$$ \sigma_{\phi }^{2} = 2\mathop \int \limits_{{f_{{{\text{lower}}}} }}^{{f_{{{\text{upper}}}} }} Tf^{ - p} {\text{d}}f $$
(2)

If 50 Hz data are used in scintillation index estimation within an interval of 1 min, the \(f_{{{\text{upper}}}}\) and \(f_{{{\text{lower}}}}\) are given as 25 and 0.1 Hz, respectively. In that case, \(\sigma_{\phi }\) is also equivalent to \(Phi60\). Performing the integration of Eq. (2), the following equation is obtained (Aquino et al. 2007):

$$ \sigma_{\phi }^{2} = 2T\left[ {\frac{{25^{r} - 0.1^{r} }}{r}} \right], $$
(3)

where \(r = 1 - p\). Figure 1 presents the relationship between \(\sigma_{\phi }\), \(p\) and \(T\) according to Eq. (3). As the figure shows, for fixed values of \(p\), dramatical increases in \(T\) values can be observed with the increase in \(\sigma_{\phi }\), while the increase in \(p\) will not cause a significant increase in \(T\) when \(\sigma_{\phi }\) is unchanged, indicating that \(T\) is less sensitive to the variation of \(p\). This conclusion is important as it motivates the concept of 1-s scintillation index estimation introduced hereafter.

Fig. 1
figure 1

Relationship among phase scintillation indices \(\sigma_{\phi }\), \(p\) and \(T\)

Full size image

3 Datasets

The scintillation data analyzed in this study were collected at Longyearbyen station (Lat. 78.35°N, Long. 15.63°E), denoted as LYB0, in the Arctic and Sachs Harbour station (Lat. 71.99°N, Long. 125.26°W), denoted as SAC, in northern Canada (Jayachandran et al. 2009). Each station is equipped with a Septentrio PolaRxS Pro receiver, which is a specialized ionospheric scintillation monitoring receiver (ISMR) that can output 50 Hz signal intensity and carrier phase measurements, 1-min scintillation indices S4, \(Phi60\), \(p\), \(T\), as well as carrier to noise density ratio (\(C/N0\)) and total electron content (TEC). The data were recorded during the geomagnetic storm that took place from 29 August, Day of year (DOY) 241, to 3 September, DOY 246, in 2019. Figure 2 presents the variation of \(Kp\) index on those days. It can be seen that the \(Kp\) increases from around 1 on DOY 241, to a maximum of 6- on DOY 243, followed by a gradual decrease to less than 2 on DOY 246.

Fig. 2
figure 2

Variation of \(Kp\) index from DOY 241 to 246 in 2019

Full size image

The amplitude scintillation with \(S4 \ge 0.3\) and phase scintillation with \(Phi60 \ge 0.3\) rad observed on GPS L1C/A signals at LYB0 and SAC stations from DOY 241 to 246 are counted and presented in Table 1. Compared with amplitude scintillation, phase scintillation is more frequently observed at both the stations. Although the amplitude scintillation at LYB0 station reaches an occurrence of 59, the \(S4\) values are all less than 0.4, indicating very weak scintillation intensities. It should be noted that satellites with elevation lower than 30° are not considered in Table 1 in order to remove the contamination of non-scintillation related effects, like multipath.

Table 1 Occurrence of amplitude and phase scintillation observed on GPS L1C/A signals at LYB0 and SAC stations from DOY 241 to 246 in 2019
Full size table

Figure 3 further shows the phase scintillation occurrence of different levels at LYB0 and SAC stations. The daily phase scintillation occurrence in the top panels peaks simultaneously on DOY 243 at both stations. Strong phase scintillation with \(Phi60 > 1.2\) rad are observed on DOY 242, 243 and 244 at LYB0 station and on DOY 243 and 246 at SAC station. The bottom panels show the overall occurrence of phase scintillation as a function of \(Phi60\). With the increase in \(Phi60\), the occurrences decrease dramatically at both stations, while apparent increases are seen when \(Phi60\) increases beyond 1.2 rad.

Fig. 3
figure 3

Daily phase scintillation occurrence of different levels (top) and overall phase scintillation occurrence in relation to \(Phi60\) (bottom) observed on GPS L1C/A signals from DOY 241 to 246 in 2019 at LYB0 and SAC stations

Full size image

4 Toward 1-s scintillation occurrence

The approach to estimate 1-s amplitude and phase scintillation indices using the 50 Hz signal intensity and phase measurements logged by the ISMRs at LYB0 and SAC stations are introduced in this section. The 1-s amplitude scintillation index, denoted as \(S4^{ - }\), is calculated following the same procedure of \(S4\) calculation, except that the standard deviation and normalization is performed within 1 s instead of 1 min. Figure 4 shows an example of the variation of detrended signal intensity \(P_{det}\), \(C/N0\) and the estimated \(S4^{ - }\) under a very weak scintillation with \(S4 = 0.15\). It can be seen that the detrended intensity in the top panel varies slightly in the first 30 s, while from around the 37th to the 46th s, apparent fluctuations are observed due to scintillation. A similar variation is seen in \(C/N0\), which presents obvious fluctuations in the same period. In the bottom panel, \(S4^{ - }\) presents large second-to-second variability. Obvious increases in \(S4^{ - }\) are seen from around the 37th s, following the signal intensity fluctuations in the top panel.

Fig. 4
figure 4

Variation of the detrended signal intensity, \(C/N0\) (top) and 1-s amplitude scintillation index \(S4^{ - }\) (bottom) on GPS L1C/A signal of PRN 16 at 18:45 UTC on DOY 243 at LYB0 station

Full size image

The 1-s phase scintillation spectral parameters \(p\) and \(T\) cannot be computed directly using the 50 Hz phase measurements, as there are not enough samples to calculate the PSD curves within an interval of 1 s. To address this limitation, this study exploits the relationship between 1-min \(Phi60\) and \(p\) indices, as shown in Fig. 5. It can be seen that with the increase in \(Phi60\), the \(p\) value first increases roughly and then maintains between a range of 2–3.5 at both stations.

Fig. 5
figure 5

\(p\) values in relation to \(Phi60\) recorded on GPS L1C/A signals at LYB0 (left) and SAC (right) stations from DOY 241 to 246 in 2019

Full size image

A power function given by

$$ p = a*Phi60^{b} + c $$
(4)

is used to fit the trend between the \(p\) and \(Phi60\) indices, as the red lines shown in Fig. 5. The fitted coefficients, \(a, b\) and \(c\) at LYB0 and SAC stations are listed in Table 2. It is assumed that the trend represented by Eq. (4) is still true between the corresponding 1-s scintillation indices. Thus, Eq. (4) can be modified as

$$ p = a*\left( {\sigma_{\phi } } \right)^{b} + c $$
(5)
Table 2 Coefficients of the power function for Eq. (4)
Full size table

Based on Eq. (5), 1-s \(p\) can be approximated by 1-s \(\sigma_{\phi }\), which is measurable using 50 Hz phase measurements, given by:

$$ \sigma_{\phi } = \sqrt {\langle{\phi^2_{det}\rangle_{1s} - \langle\phi_{\det}\rangle^{2}_{1s}}}, $$
(6)

where \(\phi_{det}\) is the detrended phase measurement. With Eqs. (3), (5) and (6), 1-s \(p\) and \(T\) can be estimated.

It is worth mentioning that although the \(p\) values are roughly estimated using Eq. (5), it can be used to calculate the 1-s \(T\) and the receiver tracking error variance, because (1) the phase scintillation-induced tracking error variance estimated using \(p \) and \(T\) maintains a relatively stable value when \(p\) varies from 2 to 3, which is shown in the next section; (2) according to the discussion of Fig. 1, \(T\) is not sensitive to the variation of \(p\). As a result, the rough estimation of \(p\) is sufficient to be used for the analysis in this study. Additionally, it should be noted that the coefficients given in Table 2 are only applicable to the scintillation datasets analyzed in this study. However, Eq. (5) can still be generally applied to other scintillation datasets at other stations, but the corresponding coefficients need to be re-estimated based on their observed scintillation indices.

Figure 6 shows an example of the estimated 1-s phase scintillation indices along with the variation of the detrended phase measurements. In this minute, although \(Phi60\) reaches a very high value of 1.49 rad, the phase measurement in the top panel remains relatively steady in the first 20 s, after which strong fluctuations are seen, leading to the increases in 1-s \(\sigma_{\phi }\) and T values shown in the bottom panel. Large second-to-second variabilities are also observed in both panels. Thus, compared with 1-min scintillation indices, 1-s indices can better describe the signal fluctuations under scintillation in this minute.

Fig. 6
figure 6

Variation of the detrended phase measurements (top) and 1-s scintillation indices (bottom) captured on GPS L1C/A signal of PRN 16 at 18:45 UTC on DOY 243 at LYB0 station

Full size image

The 1-s phase scintillation indices are estimated for all visible satellites observed at LYB0 and SAC stations from DOY 241 to 246 in 2019. Figure 7 shows the percentage of the estimated 1-s phase scintillation index \(\sigma_{\phi }\) of different levels as a function of \(Phi60\) at both stations. As the figure shows, the 1-s \(\sigma_{\phi }\) lower than 0.3 rad accounts for significant parts of all the \(Phi60\) levels. Even when \(Phi60\) is higher than 1.2 rad, around 60% of the 1-s \(\sigma_{\phi }\) are less than 0.4 rad at both stations. This indicates that most of the phase fluctuations are not severe within an interval of 1 s in the presence of the high latitude scintillation analyzed in this work.

Fig. 7
figure 7

Percentages of 1-s phase scintillation index \(\sigma_{\phi }\) of different levels as a function of 1-min index \(Phi60\) levels estimated on GPS L1C/A signals observed at LYB0 (left) and SAC (right) stations from DOY 241 to 246 in 2019

Full size image

To demonstrate the difference between the proposed 1-s phase scintillation indices and the 1-min indices output by the ISMRs, the 1-s indices are decimated to a rate of 60 s and aligned with the timestamp of the 1-min indices. The differences between the decimated 1-s \(\sigma_{\phi }\), \(p\) and \(T\) and the 1-min \(Phi60\), \(p\) and \(T\), denoted as \(\Delta \sigma_{\phi }\), \(\Delta p\) and \(\Delta T\), are calculated. The mean and standard deviation of \(\Delta \sigma_{\phi }\), \(\Delta p\) and \(\Delta T\) are estimated on each day, as shown by the error bar plots in Fig. 8. It can be seen that the standard deviation of \(\Delta \sigma_{\phi }\) and \(\Delta T\) increase obviously on DOY 242, 243 and 244 at LYB0 station and on DOY 243, 244 and 246 at SAC station. By contrast, on DOY 241 and 246 at LYB0 station and DOY 241, 242 and 245 at SAC station when less scintillation occurs, the differences \(\Delta \sigma_{\phi }\) and \(\Delta T\) tend to be smaller. This can be explained as due to the fact that when there is weak or no phase scintillation, the phase fluctuations are generally stable, thus the 1-min scintillation indices are comparable to the 1-s indices. However, under strong phase scintillation, large second-to-second variabilities in phase may appear, which results in the obvious differences between the 1-min and 1-s indices. On the other hand, the standard deviation of \(\Delta p\) is around 0.16 to 0.26 over all the days at both stations. This is reasonable as a rough estimation of 1-s \(p\) is used in the comparison. However, it has very limited effects on the estimation of phase scintillation-induced tracking error variance, which is introduced in detail in the next section.

Fig. 8
figure 8

Daily means of the differences between the estimated 1-s phase scintillation indices \(\sigma_{\phi }\), \(p\) and \(T\) and the ISMR output 1-min indices on DOY 241 to 246 at LYB0 (left) and SAC (right) stations. The standard deviations of the differences are also shown as error bars in each panel

Full size image

It is worth mentioning that since the 50 Hz signal intensity and carrier phase measurements for GPS L2P signals are not logged by the ISMRs used here, the scintillation indices for L2P signals are approximated using the following equations (Rino 1979; Hegarty et al. 2001; Conker et al. 2003)

$$ S4^{ - } \left( {L2} \right) = \left( {\frac{{f_{L1} }}{{f_{L2} }}} \right)^{1.5} S4^{ - } \left( {L1} \right) $$
(7)
$$ \sigma_{\phi } \left( {L2} \right) = \left( {\frac{{f_{L1} }}{{f_{L2} }}} \right)\sigma_{\phi } \left( {L1} \right) $$
(8)
$$ p\left( {L2} \right) = p\left( {L1} \right) $$
(9)

with the estimated 1-s scintillation indices, the tracking error variances on GPS L1C/A and L2P signals are calculated using the tracking error models described in the next section.

5 Scintillation sensitive tracking error variance

According to Knight and Finn (1998) and Conker et al. (2003), the tracking error variances at the output of PLL and DLL for GPS L1C/A signals are given by

$$ \sigma_{{{\text{PLL}}\_L1C/A}}^{2} = \frac{{B_{n} \left[ {1 + \frac{1}{{2\eta_{{{\text{PLL}}}} c/n_{0} \left( {1 - 2*S4^{2} } \right)}}} \right]}}{{c/n_{0} \left( {1 - S4^{2} } \right)}} + \frac{\pi T}{{kf_{n}^{p - 1} \sin \left( {\frac{{\left[ {2k + 1 - p} \right]\pi }}{2k}} \right)}} + \theta_{A}^{2} , \,\,\,S4 < \frac{\sqrt 2 }{2}, 1 < p < 2k $$
(10)
$$ \sigma_{{{\text{DLL}}\_L1C/A}}^{2} = \frac{{B_{L} d\left[ {1 + \frac{1}{{\eta_{{{\text{DLL}}}} c/n_{0} \left( {1 - 2*S4^{2} } \right)}}} \right]}}{{2c/n_{0} \left( {1 - S4^{2} } \right)}}\left( chip^{2} \right), $$
(11)

where \(B_{n}\) is the PLL one-side noise-equivalent loop bandwidth and equals to 15 Hz. \(\eta_{{{\text{PLL}}}}\) denotes the PLL integration time, equaling to 10 ms. The \(c/n_{0}\) is fractional form of \(C/N0\) measured on GPS L1C/A, given by \(c/n_{0} = 10^{{0.1{*}C/N0}}\). \(k\) is the PLL loop order, equaling to 3. \(f_{n}\) is the loop natural frequency, given by \(f_{n} = \frac{1.2}{{2\pi }}B_{n}\). \(\theta_{A}^{2}\) is the receiver oscillator noise-induced tracking error variance and equals to \(9.2 \times 10^{ - 6} \;{\text{rad}}^{2}\) (Irsigler and Eissfeller 2002). \(B_{L}\) is the DLL one-side bandwidth and equals to 0.25 Hz. \(d\) is the correlator spacing in C/A chips. \(\eta_{{{\text{DLL}}}}\) is the DLL integration time, equal to 100 ms. In Eqs. (10) and (11), the scintillation indices \(S4\), \(p \) and \(T \) are all measured on GPS L1C/A signal.

The second component in the right part of Eq. (10) is the phase scintillation-induced tracking error variance, denoted as \(\sigma_{pha}^{2}\). Variation of \(\sigma_{pha}^{2}\) in relation to \(p\) and \(\sigma_{\phi }\) is presented in Fig. 9. It can be seen that when \(\sigma_{\phi }\) is fixed, \(\sigma_{pha}^{2}\) decreases gradually with an increase of \(p\) from 1 to around 2. When \(p\) keeps increasing, \(\sigma_{pha}^{2}\) is seen to vary slightly, indicating the phase scintillation-induced tracking error variance \(\sigma_{pha}^{2}\) is not sensitive to \(p\) higher than 2.0.

Fig. 9
figure 9

Variation of phase scintillation-induced tracking error variance \(\sigma_{pha}^{2}\) in relation to phase scintillation indices \(p\) and \(\sigma_{\phi }\)

Full size image

Considering the case where the phase measurements on L1C/A and P are uncorrelated, the PLL and DLL tracking error variances for L2P signal are given by (Conker et al. (2003):

$$ \sigma_{{{\text{PLL}}\_L2P}}^{2} = \frac{{B_{n} \left[ {1 + \frac{1}{{2\eta_{Y} \left( {c/n_0} \right)_{ L1P} \left( {1 - S4\left( {L1} \right)^{2} } \right)}}} \right]}}{{\left( {c/n_0} \right)_{ L2P} \left( {1 - S4\left( {L2} \right)^{2} } \right) }} + \frac{\pi T}{{kf_{n}^{p - 1} \sin \left( {\frac{{\left[ {2k + 1 - p} \right]\pi }}{2k}} \right)}}\left( {\alpha^{2} + {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {\alpha^{2} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\alpha^{2} }$}}} \right) + \theta_{A}^{2} , S4\left( {L1} \right) < 0.687 $$
(12)
$$ \sigma_{{{\text{DLL}}\_L2P}}^{2} = \frac{{B_{L} \left[ {1 + \frac{1}{{2\eta_{Y} \left( {c/n_0} \right)_{{{ }L1P}} \left( {1 - S4\left( {L1} \right)^{2} } \right)}}} \right]}}{{2\left( {c/n_0} \right)_{{{ }L2P}} \left( {1 - S4\left( {L2} \right)^{2} } \right){ }}}, $$
(13)

where the PLL and DLL parameters \(B_{n}\), \(k\), \(f_{n}\), \(B_{L}\) are the same as in Eqs. (10) and (11) but for the GPS L2P signal tracking. In this case, \(B_{n} = 0.25{\text{ Hz}}\), \(k = 2\), \(B_{L} = 0.3{\text{ Hz}}\). \(\eta_{Y}\) is the integration for GPS Y code, equal to \(1.96 \times 10^{ - 6}\) s. \(\left( {c/n_0} \right)_{{{ }L1P}}\) and \(\left( {c/n_0} \right)_{{{ }L2P}}\) are the fractional form of \(C/N0\) values on L1P and L2P signals, respectively. \(\alpha\) is the ratio of the carrier frequencies of L2P and L1C/A. The \(C/N0\) on L2P in this analysis is obtained from the ISMR output, while the \(C/N0\) on L1P is approximately by \(\left( {C/N0} \right)_{L1P} \cong \left( {C/N0} \right)_{L1C/A} - 3{\text{dB}}\). Except for the \(S4\)(\(L1\)) which is measured on L1C/A signal, the \(S4\left( {L2} \right)\), \(p\) and \(T\) in Eqs. (12) and (13) are all for L2P signal and estimated using Eqs. (7), (8) and (9), respectively. It should be noted that all the parameters related to the receiver PLLs and DLLs are known by the ISMR configurations in the scintillation data collection. Additionally, it can be seen that the value of \(S4\) needs be less than \(\frac{\sqrt 2 }{2}\) in Eq. (10) and less than 0.687 in Eq. (12). In this analysis, a value of 0.70 and 0.65 is, respectively, used to replace the \(S4\) value when it is over \(\frac{\sqrt 2 }{2}\) when using Eq. (10) and over 0.687 when using Eq. (12) in the tracking error variance calculation.

By substituting the proposed 1-s scintillation index into these tracking error models, 1-s tracking error variances can be estimated. Figure 10 presents the 1-s DLL and PLL tracking error variances estimated with the 1-s scintillation indices shown in Figs. 4 and 6. As the figure shows, when there is no amplitude and phase scintillation during the first 20 s, the DLL and the PLL tracking error variances are relatively stable. Under scintillation occurrence between the 37th and the 43rd s, they increase to around 0.013 \({\mathrm{m}}^{2}\) and 0.018 \({\mathrm{rad}}^{2}\), respectively, in the top and bottom panels, which are in agreement with the signal intensity and phase fluctuations caused by scintillation shown in Figs. 4 and 6. Thus, these 1-s estimates of the tracking error variance can successfully describe the signal intensity and phase noise levels under the scintillation in this minute. The investigation of scintillation mitigation on GPS positioning with the 1-s tracking error variances is given hereafter.

Fig. 10
figure 10

Variation of 1-s DLL and PLL tracking error variances for GPS L1C/A signal on PRN 16 at 18:45 UTC on DOY 243 at LYB0 station

Full size image

6 PPP under high latitude scintillation

This section presents the high latitude scintillation effects on GNSS PPP, which is performed by the University of Nottingham in-house POINT software (Mohammed 2017). Code and carrier phase measurements on GPS L1C/A and L2P signals are used to form the IF linear combination. International GNSS Service (IGS) precise orbits and clock products estimated by the Centre for Orbit Determination in Europe (CODE) are used to constrain the PPP model. CODE Differential Code Biases (DCB) and satellite and receiver antenna corrections from IGS are used. The Melbourne-Wübbena Wide-Lane (MWWL) combination, involving code and carrier phase measurements on L1C/A and L2P signals, and ionospheric total electron contents rates (TECRs), involving the geometry-free combination of L1 and L2 carrier phase measurements, are both calculated to jointly detect the cycle slips (Liu 2011). When a cycle slip is detected, the biased carrier phase ambiguity will be reset and not used in positioning estimation. Details of the PPP calculation are summarized in Table 3.

Table 3 Strategies for PPP calculation
Full size table

Following the strategies in Table 3, PPP is calculated with the observations sampled at a rate of 1 s from DOY 241 to 246 at LYB0 and SAC stations. As no obvious scintillation is observed on DOY 241 by both stations, the precise coordinates obtained by static PPP based on 24 h data on this day are set as the reference coordinates for each station. The positioning errors obtained using kinematic PPP, which solves for the station coordinates on an epoch per epoch basis, are estimated against the reference coordinates using a satellite elevation-based weighting strategy. Figure 11 shows the kinematic PPP errors in the east (\(\Delta E\)), north (\(\Delta N\)) and up (\(\Delta U\)) directions along with the 1-s phase scintillation index \(\sigma_{\phi }\) on DOY 243 at LYB0 and SAC stations. The scintillation for different satellites in the top panels is distinguished by colors. When there is no obvious scintillation, the positioning errors in all three directions in the bottom panels are generally less than 0.2 m after convergence. However, under frequent strong phase scintillation from 18:00 to 20:00 UTC at LYB0 station and from 06:00 to 07:00 and 12:00 to 13:00 UTC at SAC station, positioning errors increase significantly to more than 0.5 m. Additionally, in the presence of moderate phase scintillation occurred from 08:00 to 14:00 UTC at LYB0 station and from 16:00 to 00:00 UTC at SAC station, only slight increases are observed in the positioning errors, which means that moderate and weak scintillation with \(\sigma_{\phi } < 0.5\) rad has less impacts on the positioning accuracy. Furthermore, the positioning errors in the up direction show more prominent fluctuations at both stations, indicating that the positioning accuracy in this direction is more susceptible to scintillation effects.

Fig. 11
figure 11

1-s phase scintillation index \(\sigma_{\phi }\) estimated on GPS L1C/A signals (top) and positioning errors in the east, north and up directions (bottom) when calculating kinematic PPP with an elevation weighting strategy on DOY 243 at LYB0 and SAC stations. The positioning errors are the deviations of the kinematic PPP results with respect to the reference coordinates

Full size image

7 Mitigating scintillation effects on PPP

The mitigation of scintillation effects on PPP is presented in this section. Kinematic PPP with a tracking error variance-based weighting strategy using code and phase measurements with a sampling rate of 1 s is first performed from DOY 241 to 246 at LYB0 and SAC stations. The 1-s tracking error variances, estimated by the newly proposed 1-s scintillation indices, are used to improve the least square stochastic models in positioning following the mitigation approach described in Aquino et al. (2009). Positioning errors are compared with respect to those of the satellite elevation weighting strategy. Figure 12 shows the kinematic PPP errors in the east, north and up directions using the tracking error variance weighting on DOY 243 at LYB0 and SAC stations. The hourly root mean squares (RMSs) of 3D positioning errors are also included. It can be seen that although the positioning errors in the top panels still obviously increase from 18:00 to 20:00 UTC at LYB0 station and from 06:00 to 08:00, 12:00 to 13:00 UTC at SAC station, the 3D positioning errors in the bottom panels are greatly reduced compared with those when an elevation weighting is used. The RMS of 3D errors in the 6th hour is mitigated from nearly 1 m, when using elevation weighting, to less than 0.2 m at SAC station, indicating a significant improvement in positioning accuracy by using the 1-s tracking error variance weighting approach in this hour.

Fig. 12
figure 12

Positioning errors in the east, north and up directions when calculating kinematic PPP using a tracking error variance weighting strategy (top) and hourly RMSs of the 3D positioning errors when using elevation and tracking error variance weighting strategies (bottom) on DOY 243 at LYB0 and SAC stations

Full size image

The daily 3D positioning errors of kinematic PPP when using elevation and tracking error variance weighting strategies from DOY 241 to 246 at LYB0 and SAC stations are calculated and shown in Fig. 13. Due to the PPP convergence process, the first hour of the positioning error time series on each day is not considered. With an elevation weighting strategy, the daily RMSs of 3D PPP errors are less than 0.05 m when there is only weak or no scintillation on DOY 241 and 246 at LYB0 station and on DOY 241 and 245 at SAC station. By contrast, when strong phase scintillation is frequently observed, the 3D positioning errors increase to 0.208 m on DOY 244 at LYB0 station and 0.873 m on DOY 246 at SAC station, which are more than 5 times higher than those on weak or no scintillation days. Additionally, it is interesting to note that the PPP errors peak on DOY 246 at SAC station, although the geomagnetic activity is quiet according to the \(Kp\) index shown in Fig. 2. This may be due to the different responses of the ionosphere at different high latitude regions to the geomagnetic storm and is worth analyzing further in future studies. On the other hand, the daily 3D positioning errors with tracking error variance weighting decrease at both LYB0 and SAC stations. In the left panel, the 3D positioning errors on DOY 244 at LYB station are reduced by 77.55% to 0.047 m. Meanwhile, a decrease of 90.30% is seen on DOY 246 at SAC station. Therefore, using the 1-s tracking error variance weighting can successfully improve the GPS positioning accuracy under the scintillation analyzed in this work. Although the kinematic PPP errors with the mitigation approach are slightly higher on DOY 241 at SAC station, the differences compared with the elevation weighting are less than 0.01 m, which is within the noise level and acceptable.

Fig. 13
figure 13

RMSs of 3D positioning errors of kinematic PPP using elevation and tracking error variance weighting strategies from DOY 241 to 246 at LYB0 (left) and SAC (right) stations. The kinematic PPP is estimated based on code and carrier phase measurements with a sampling interval of 1 s

Full size image

In order to compare the performance of tracking error variance weighting in scintillation mitigation by, respectively, exploiting the 1-min and 1-s scintillation index, kinematic PPP using the settings summarized in Table 3 is recalculated from DOY 241 to 246 at LYB0 and SAC stations, however, with a sampling interval of 60 s for code and carrier phase measurements. The 1-min tracking error variance is directly estimated using the 1-min scintillation indices and averaged \(C/N0\) values output by the ISMRs, while the 1-s tracking error variance is estimated using the proposed 1-s scintillation index but downsampled and aligned with the timestamp of code and carrier phase measurements. Figure 14 shows the up direction positioning errors for kinematic PPP using elevation, 1-min and 1-s tracking error variance weighting strategies from DOY 243 to 246 at LYB0 and SAC stations. The positioning errors are the deviations of kinematic PPP results against the reference coordinates. It is obviously seen that compared with the elevation weighting strategy, the 1-min and 1-s tracking error variance weighting can generally achieve better positioning accuracies in the up direction, especially on DOY 243 and 244 at LYB0 station and DOY 243, 244 and 246 at SAC station when phase scintillation is frequently observed.

Fig. 14
figure 14

Variation of kinematic PPP errors in the up directions from DOY 243 to 246 at LYB0 (top) and SAC (bottom) stations. The positioning error is obtained for PPP estimation using elevation, 1-min and 1-s tracking error variance weighting strategies

Full size image

The daily RMSs of 3D kinematic PPP errors with the three different weighting strategies at LYB0 and SAC stations are shown in Fig. 15. As the figure shows, on days when scintillation was frequently observed, the elevation weighting strategy has the largest 3D positioning errors reaching values of 0.200 m on DOY 244 at LYB0 station and 1.196 m on DOY 246 at SAC station compared with the tracking error variance weighting strategy. On the other hand, when tracking error variance weighting strategies are applied, the PPP errors are greatly reduced. On DOY 244 at LYB0 station, the positioning error is reduced, respectively, to 0.078 m and 0.043 m by 1-min and 1-s tracking error variance weighting strategies, manifesting great improvements with respect to the elevation weighting strategy. Additionally, although the 3D positioning error on DOY 246 at SAC station is only reduced to 0.400 m by the 1-min tracking error variance weighting, the 1-s one reduces the error to 0.210 m with an improvement of 88.26% in the positioning accuracy.

Fig. 15
figure 15

RMSs of 3D positioning errors of kinematic PPP using elevation, 1-min and 1-s tracking error variance weighting strategies from DOY 241 to 246 at LYB0 (left) and SAC (right) stations. The kinematic PPP is estimated based on code and carrier phase measurements with a sampling interval of 60 s

Full size image

Table 4 summarizes the improvements in 3D PPP errors using tracking error variance weighting compared with an elevation weighting corresponding to the results shown in Fig. 15. As it can be seen, both the 1-min and 1-s tracking error variance weighting strategies help to improve the PPP accuracy under scintillation, however, the greatest improvements mostly occur when the latter is used. An improvement of 78% and 83% in the 3D errors are, respectively, achieved with 1-s tracking error variance weighting on DOY 244 at LYB0 station and on DOY 246 at SAC station, when frequent and strong phase scintillation occurs. Consequently, it can be concluded that the 1-s tracking error variance weighting strategy performs better in GPS positioning calculation compared with the 1-min one under the high latitude phase scintillation analyzed in this study.

Table 4 Improvements in 3D PPP errors using tracking error variance weighting strategies with respect to using an elevation weighting strategy
Full size table

It is noticed that on DOY 241 and DOY 246 at LYB0 station, although the 3D positioning error for PPP using 1-min tracking error variance weighting is slightly worse than using elevation weighting, it only has a difference less than 0.005 m. Additionally, when there is only weak to moderate or no scintillation occurring on DOY 241, 245 and 246 at LYB0 station and DOY 241 and 245 at SAC station, comparable performance is observed when 1-s and 1-min tracking error variance weighting strategies are used. This may be due to the fact that under weak or moderate phase scintillation, the phase fluctuations are relatively stable and 1-min scintillation indices are comparable to the 1-s indices, which is also shown in Fig. 8, thus the 1-min tracking error variance achieves a better representative of the measurement noise level over the interval, which results in comparable results for the scintillation mitigation. Furthermore, although the 1-s tracking error variance weighting strategy significantly mitigates the strong phase scintillation effects on positioning, it can barely achieve a PPP accuracy comparable to non-scintillation days. As a result, scintillation mitigation tools can be further explored to mitigate the adverse scintillation effects on positioning. This will be the focus of the follow-on work.

8 Conclusion and remarks

This study analyzed the GPS scintillation data recorded by ISMRs operational at LYB0 station in the Arctic region and at SAC station in northern Canada during a geomagnetic storm in 2019. A new approach is proposed to estimate 1-s scintillation indices based on the 50 Hz signal intensity and carrier phase measurements logged by the ISMRs. It is found that the 1-s scintillation indices can better describe the detailed signal fluctuations and distortion under scintillation compared with the commonly used 1-min scintillation indices output by ISMRs. By exploiting the 1-s scintillation indices, the PLL and DLL tracking error variances considering scintillation effects were calculated. Those estimated tracking error variances are shown to successfully describe the signal intensity and phase noise levels under scintillation.

In the presence of the high latitude scintillation, the kinematic PPP is first calculated based on the code and phase measurements with a sampling rate of 1 s. Results show that when using an elevation weighting strategy, the RMS of 3D positioning error can be more than 5 times higher on days with strong phase scintillation than those days with only weak or no scintillation. The scintillation effects on PPP were then mitigated following the approach described in Aquino et al. (2009), that is, to modify the least square stochastic model using the 1-s tracking error variances. Results show that the 3D positioning errors can be mitigated by up to 77.55% at LYB0 station and up to 90.30% at SAC stations, indicating that the 1-s tracking error variance weighting can successfully mitigate GPS PPP errors under the scintillation analyzed in this study.

In order to compare the performance of tracking error variance weighting approaches in scintillation mitigation by, respectively, exploiting the 1-min scintillation indices output by ISMRs and the newly proposed 1-s scintillation indices, the kinematic PPP is recalculated based on the code and phase measurements with a sampling rate of 60 s. Results show that both the 1-min and 1-s tracking error variance weighting strategies help to improve the PPP accuracy under scintillation. However, the latter performs even better in GPS positioning calculation compared with the 1-min one under strong phase scintillation. Additionally, it is found that even when the 1-s tracking error variance weighting is applied, the PPP under scintillation can barely achieve an accuracy comparable to non-scintillation days. Therefore, the adverse scintillation effects on PPP can be further reduced by exploring new scintillation mitigation tools, which will be the focus of the follow-on work.