A Method to Mitigate the Effects of Strong Geomagnetic Storm on GNSS Precise Point Positioning

Geomagnetic storm can affect the performance of Global Navigation Satellite System (GNSS) precise positioning services. To mitigate the adverse effects of strong geomagnetic storms, we propose to establish the geometry‐free (GF) cycle slip threshold model based on ionospheric disturbance index rate of total electron content index to reduce the false detection rate of cycle slip in GNSS precise point positioning (PPP) during strong storm periods, thus improving the accuracy and reliability of GNSS PPP. The performance of our proposed model is validated by using 171 International GNSS Service (IGS) tracking stations data on 8 September 2017. The analysis indicates that compared with conventional PPP scheme, the proposed model can improve the positioning accuracy by approximately 14.0% (36.8%) and 23.1% (51.5%) in the horizonal and vertical components for global (high latitudes) stations. Furthermore, the availability of our proposed model is also validated by PPP experiments using 379 IGS tracking stations data during another strong storm occurred on 26 August 2018.

• Strong geomagnetic storm can significantly degrade the performance of Global Navigation Satellite System (GNSS) precise point positioning (PPP) • Geometry-free cycle slip detection observation shows close relationship to the ionospheric disturbance index rate of total electron content index • Positioning accuracy of GNSS PPP based on cycle slip threshold model is generally better than that of conventional PPP Correspondence to: X. Luo, luoxiaomin@cug.edu.cn
10.1029/2021SW002908 2 of 13 Berdermann et al., 2018;Zakharenkova & Cherniak, 2021). In the early stage, Odijk (2001) investigated the effects of geomagnetic storms on global positioning system (GPS) relative positioning. It is reported that the relative ionospheric delays cannot be neglected during strong storm periods even for rather short baselines (4 km). Statistical results in Odijk (2001) indicated that, without (with) ionospheric correction, the ambiguity success rates of instantaneous GPS precise positioning are 66% (95%) and 36% (77%) under minor and intense storms, respectively. The study shown in Rama Rao et al. (2009) demonstrated that during storm periods, the number of cycle slips increased significantly compared with those during quiet days. Astafyeva et al. (2014) further reported that most GNSS cycle slips occurred in geomagnetic storms should be associated with ionospheric irregularities. It is known that cycle slips caused by the rapid changes of TEC during storm periods can result in loss of locks of GNSS receivers, thus reducing the ambiguity success rates as well as the positioning performance of GNSS real-time kinematic (RTK) positioning (Odolinski & Teunissen, 2019).
Except GNSS RTK technique, many studies also investigated the geomagnetic storm effects on precise point positioning (PPP). Focusing on different levels of geomagnetic storm, Luo et al. (2018) Yasyukevich et al. (2020) reported that the time range of positioning degradation of PPP is longer than that during the St. Patrick's Day storm. It is also found that except high latitudes, the ionospheric irregularities of middle latitudes expanded from the auroral oval can also result in PPP degradation. Recently, Zakharenkova and Cherniak (2021) assessed the performance of ground-and space-based GPS PPP during the strong storm on 7-8 September 2017. Under the strong storm, the 3D error of GPS PPP for ground-based stations rose to several meters, while that on normal condition is only centimeter or decimeter level. In addition, the storm-induced irregularities can result in 5 times errors in kinematic orbit for Swarm satellites.
Although many studies have reported the adverse effects of geomagnetic storm on GNSS, few papers sought ways to mitigate these effects on GNSS precise positioning (Lu et al., 2020;Rodríguez-Bilbao et al., 2015). Since the number of cycle slip increases in storm conditions, correctly detecting or even repairing cycle slips should be the effective solution to improve GNSS PPP or RTK performance during storm periods. However, under complex ionospheric environment caused by geomagnetic storm, it is difficult to detect cycle slips correctly. It is known that the classic method as TurboEdit method, combined with the geometry-free (GF) and Melbourne-Wübbena (MW) wide-lane combinations, is usually applied in cycle slips detection (Blewitt, 1990). The GF combination is sensitive to cycle slips under normal ionospheric conditions, but it will be degraded during active ionospheric environments (Luo et al., 2020). That is why some researchers chose to set an empirical loose threshold for GF method (Zhang et al., 2014) or directly deactivated the GF method (Rodríguez-Bilbao et al., 2015) to improve GNSS PPP performance in geomagnetic storm. In this study, we propose to establish the GF cycle slip detection threshold model through considering ionospheric disturbance index ROTI (rate of total electron content index) to reduce the false detection rate of cycle slip, thus improving GNSS PPP accuracy under strong geomagnetic storm conditions.

Data and Methodology
This section first introduces the strong geomagnetic storm event occurred on 8 September 2017. Then, the geographical distribution of 417 IGS stations is presented. The ionospheric-free PPP model and cycle slip detection method of TurboEdit are introduced following. Finally, we give the cycle slip threshold model in detail.

Strong Geomagnetic Storm on 8 September 2017
Since two solar flares X2.2 and X9.3 erupted over the active region 2673 on 6 September 2017, a strong geomagnetic storm with double main phases was generated on 8 September 2017 (Yamauchi et al., 2018). Based on the observations derived from GNSS receivers, radars, and ionosondes, numerous studies have investigated LUO ET AL.

Experiment Data
The experiment is conducted by using GPS dual-frequency data collected at 417 IGS tracking stations on 8 September 2017. Figure 2 presents the geographical distribution of 417 IGS stations. We randomly select 246 stations (black dots) data for cycle slip threshold modeling, while the remaining 171 stations (red dots) data are used for PPP verification.

Ionospheric-Free PPP Model
Using precise satellite orbit and clock products from IGS as well as pseudorange and carrier phase measurements from single GNSS receiver, PPP technique can achieve decimeter or even millimeter levels of positioning accuracy  10.1029/2021SW002908 4 of 13 (Zumberge et al., 1997). The GNSS pseudorange and carrier phase measurements are expressed as below (Kouba & Héroux, 2001): where " ," " ," and " " represent the identification of GNSS receiver, satellite, and observation frequency, respectively; and Φ represent the pseudorange and carrier phase measurements, respectively; is the geometric distance (m); and are the receiver and satellite clock offset (s); trop is the tropospheric delay (m); ion, is the ionospheric delay at the frequency (m); and are the frequency-dependent signal delay for the receiver and satellite (m); is the wavelength of the satellite signal at the frequency (m); is the float ambiguity (cycle); is the phase windup error (cycle); and is the measurement noise (m).
In this study, GNSS PPP based on ionospheric-free linear combination model are performed to eliminate the first-order effect of ionospheric refraction. The model can be expressed as:

TurboEdit Method
The TurboEdit method is widely used in cycle slip detection and correction for dual-frequency GNSS measurements (Blewitt, 1990). It covers MW wide-lane combination as MW and GF combination as ΦGF for cycle slip detection: where MW is the MW wide-line combination of carrier phase measurements (m); MW = ∕( 1 − 2) ; and = 2 1 ∕ 2 2 .
For MW wide-lane combination, the cycle slip detection observation Δ MW and corresponding standard deviation can be expressed as: where is the current epoch in the data arc; MW is the mean value of MW , which can be calculated as: Generally, the cycle slip is considered to be detected when the following two conditions are satisfied: LUO ET AL. 10.1029/2021SW002908

of 13
For GF combination, the residual between GF and GF , which is the value of the polynomial fit of GF , is used to detect cycle slip. Since the GF is affected by the polynomial type, fitting order, and truncation error, we use Δ GF , which is the GF difference between two adjacent epochs, to detect cycle slips. The feasibility of this operation in cycle slip detection has been demonstrated by several previous studies (Chen et al., 2016;Zhang et al., 2014) as well as GNSS software gLAB (GNSS-Lab Tool; Ibáñez et al., 2018). The threshold of GF as ThGF is closely related to the ionospheric environment, data sampling interval, and satellite elevation angle. The empirical model of ThGF associated with data sampling interval has been applied in gLAB software as: where the max = 0.08, min = 0.034, 0 = 60, and step is the data sampling interval. When the sampling interval is 30 s, the values of ThGF is 0.05 m.
Under geomagnetic storm conditions, the empirical value ThGF as 0.05 m may result in the misjudgment of cycle slip detection. In the following, we provide a typical example to analyze the limitations of the empirical value. Figure 3 shows the time series of Δ MW , ΔΦGF , and ROTI for different GPS satellites observed at SCOR station (70.49°N, 21.95°W; MLAT: 71.43°N) on 8 September 2017. The threshold of MW wide-lane combination as 4 for different satellites are shown in the first panel. The ROTI value in the figure is calculated by using carrier phase measurements on GPS L1 and L2 with the sampling interval of 30 s, which is the same as in Pi et al. (1997). To guarantee that each data epoch has ROTI value, the calculation of ROTI is based on a sliding window for each Note that the ROTI will be overestimated if the cycle slip is not detected and corrected (Juan et al., 2017). In this study, we correct the cycle slip in L1 and L2 measurements based on TurboEdit method. Specifically, according to the analysis in Luo et al. (2020) and Zhang et al. (2014), under geomagnetic storm conditions, the relatively loose threshold as Th = 0.5 m are adopted in the cycle slip detection. Although the threshold of GF as 0.5 m is relative loose, the 4 threshold of MW guarantees that most cycle slips can be detected. After cycle slip detection, the candidates of cycle slip can be calculated by using the MW wide-lane and GF combinations. Through the residual test method, we can finally determine the proper cycle slip value in general. Note that under active ionospheric conditions, there are some complex cycle slips, such as the small cycle slips, that are difficult to detect and repair. For this situation, we also adopt the L1 and L2 measurements to calculate the ROTI since the small cycle slips cannot result in largely overestimated ROTI (Juan et al., 2017).
In Figure 3, it is seen that during 8:10-11:50 UT, no ionospheric disturbance is observed at SCOR station, and the time series of Δ MW and ΔΦGF fluctuate within ±4 cycle and ±0.05 m, respectively. For the other time periods on 8 September 2017, obvious disturbances occurred over SCOR station and the value of ionospheric disturbance index ROTI can reach around 6 TECU/min. Under this condition, however, the time series of Δ MW fluctuate within 4 and they do not show obvious variations with ROTI. Note that this situation is also found in many other stations during storm periods. That means Δ MW is not sensitive to ionospheric disturbance. For the ΔΦGF , significant fluctuations can be seen in the figure and many epochs of ΔΦGF are larger than the common threshold of 0.05 m (see the black lines), which means that the empirical cycle slip threshold (ThGF = 0.05 m) is not suitable under geomagnetic storm condition. For the figure, we may also note that the variations of ΔΦGF shows close relationship with ROTI. In view of this, we next establish the GF cycle slip threshold model based on abundant data of ΔΦGF and ROTI derived from 246 IGS stations. Figure 4 gives the distribution of ΔΦGF and ROTI derived from 246 IGS tracking stations data on 8 September 2017. Under non-ionospheric disturbance (ROTI ≤ 0.5 TECU/min), it is seen that most |ΔΦGF| are less than 0.1 m. With the increase of ROTI, the value of |ΔΦ | gradually increases. Note that the |ΔΦGF| value no longer increases significantly when the ROTI increases to around 3 TECU/min.

GF Cycle Slip Threshold Model
According to the distribution of ΔΦGF and ROTI, we choose the piecewise fitting strategy to model the GF cycle slip threshold. From Figure 4, ΔΦGF values are generally distributed within ±0.1 m when ROTI ≤ 0.5 TECU/ min. Therefore, a constant of threshold can be used in this ROTI range. For 0.5 < ROTI ≤ 3 TECU/min, ΔΦGF appears to increase with the increase of ROTI, so a second-order polynomial model of ThGF can be established based on ROTI values. Specifically, the selection of independent variable as ROTI is based on each 0.05 TECU/min interval in the range of 0.5-3 TECU/min, while the corresponding dependent variable as ThGF is based on the principle of 3σ (99.7%). For instance, for the interval of 0.625 < ROTI ≤ 0.675 TECU/ min (see Figure 5), the independent variable as ROTI is 0.650 TECU/min, and the dependent variable as ThGF is 0.141 m (3 × 0.047). According to this strategy, a series of independent and dependent variables can be obtained. Finally, we can get the coefficients of the second-order polynomial model. For ROTI > 3 TECU/min, ΔΦGF does not show obvious increase, so we set ThGF to a constant for this range.  From the above analysis, the GF cycle slip threshold model based on ionospheric disturbance index ROTI can be established as: Note that for higher ROTI values as ROTI > 3 TECU/min, ΔΦGF should also have a tendency of increase. Statistics indicate that in Figure 4, the number of ROTI ≤ 0.5 TECU/min, 0.5 < ROTI ≤ 3 TECU/min, and ROTI > 3 TECU/min accounts for 94.6%, 5.2%, and 0.2%, respectively. For such little data of ROTI > 3 TECU/min, we chose to set ThGF to a constant instead of modeling it. That means the performance of our model may be reduced at very high ROTI values.

Experimental Validations
To evaluate the performance of GF cycle slip threshold model, GPS data derived from 171 non-modeling IGS stations, are performed in two PPP schemes. Specifically, scheme 1 uses the common cycle slip threshold as ThGF = 0.05 m, while scheme 2 is based on our proposed GF cycle slip threshold model. Since Δ MW is not sensitive to ionospheric disturbances, the threshold of MW wide-line combination still takes the 4 in two schemes. In the data processing, the satellite elevation angle is 15° and the data sampling interval is 30 s. The precise satellite orbit and clock products are the final precise products of IGS. Table 1 presents the detailed information of models and methods used in GPS PPP. The experimental data is processed by FUSING software developed by Wuhan University (Shi et al., 2019;Yang et al., 2019).  Under ionospheric disturbances, the number of cycle slip detected by scheme 2 is obviously less than that of scheme 1, which reduces the frequently unnecessary ambiguity resets in GPS PPP, thus improving the PPP accuracy. To a single satellite system like GPS, if four to six satellites encounter cycle slip at the same epoch, the PPP will be reinitialized at this epoch in general, especially the case when less than four available satellites used in PPP positioning (Zhang & Li, 2012). From the figure, it is clearly seen that several reinitializations in PPP results of scheme 1 corresponding to lots of cycle slips very well. With the cycle slip threshold model, the number of misjudgment cycle slips is obvious reduced. Although two or three satellites suffer cycle slips simultaneously during strong storm periods, the PPP solution would not be deteriorated since most of ambiguities are not reset. That is why, compared with scheme 1, the time series of positioning results of scheme 2 show smoother variations. It also should be mentioned that when real cycle slips occurred in three or more satellites carrier phase measurements at the same epoch, our proposed model is unable to improve the PPP performance. For instance, at the epoch of 12:29:30 UT for YEL3 station, a reinitialization of positioning of scheme 2 PPP can be seen in the figure. For this situation, combining other GNSS measurements in PPP may be an optional strategy (Marques et al., 2018).

Single Station PPP
In addition, in the figure, we can see that the numerous points of positioning errors are in the convergence time period of PPP. Since the reinitialization of PPP in this study are generally associated with ionospheric disturbances during geomagnetic storm, the statistics of RMS are calculated by using all points of positioning errors for the two PPP schemes in this study. For the normal day, as shown in Figure 7, the positioning results of two PPP schemes are generally comparable. It is seen that scheme 2 based on cycle slip threshold model is effective for several stations located at high latitudes. Table 2 indicates that positioning accuracy of scheme 2 is generally better than that of scheme 1, which means our proposed model is valid under normally ionospheric environment.
For the abnormal day, as shown in Figure 8, we can see that the number of stations with darker colors in the bottom panels is less than that in the upper panels. Statistical results indicate that after applying cycle slip threshold model in GPS PPP, the number of stations with improved accuracy is 105 (61.4%), same accuracy is 39 (22.8%), and degraded accuracy is 27 (15.8%), respectively. Note that the positioning results of scheme 2 for the 27 degraded stations are comparable compared with scheme 1. Compared with scheme 1, as shown in Table 2, global GPS PPP based on cycle slip threshold model can improve the positioning accuracy by 14.0% and 23.1% in the horizontal and vertical directions, respectively. Compared with low and middle latitudes, the high latitudinal regions are generally close to the open magnetic field lines, so they are more easily affected by the geomagnetic storm. We further present statistical results of two PPP schemes for high latitudes stations in Table 2. It is seen that GPS PPP based on cycle slip threshold model at high latitudes can improve the positioning accuracy by approximately 36.8% and 51.5% in the horizontal and vertical components, over the conventional PPP solutions.
The above experiments mainly focus on the availability analysis of cycle slip threshold model for the strong storm on 8 September 2017. Next, we further present positioning results of two PPP schemes for another strong storm  (Bolaji et al., 2021;Yang et al., 2020). Using 379 IGS stations data on 26 August 2018, Figure 9 presents RMS values of two PPP schemes in the horizontal and vertical components, respectively. It is clearly seen that the positioning accuracies of scheme 2 are better than these of scheme 1. According to the averaged RMS values shown in Table 3, our proposed model can improve the positioning accuracy by 9.7% (19.8%) and 14.6% (39.7%) in the horizontal and vertical components for global (high latitudes) stations compared with conventional PPP solution.

Summary
As an important space weather, geomagnetic storm can affect the accuracy and reliability of GNSS PPP. To mitigate the effects of geomagnetic storm on GNSS precise positioning, we propose to establish the cycle slip threshold model based on ionospheric disturbance index ROTI to reduce the rate of false detection of cycle slip in GNSS PPP under the conditions of strong geomagnetic storm.

of 13
This study collects 246 IGS stations data during the strong geomagnetic storm on 8 September 2017 to analyze the relationship between GF cycle slip detection observation ΔΦGF and the ionospheric disturbance index ROTI. A piecewise function model of ThGF based on ROTI is established by using the polynomial fitting method. GPS PPP results of three IGS stations as PRDS, YEL3, and SCOR indicate that the positioning accuracy based on our proposed model is obviously better than that based on common cycle slip threshold. The PPP experiments using the data from global 171 IGS stations also demonstrate that compared with conventional PPP solution, our proposed model can improve the positioning accuracy by 14.0% and 23.1% in the horizontal and vertical directions, respectively; meanwhile, the PPP results of 105 (61.4%) stations are improved, especially those located at high latitudes. Furthermore, the availability of our proposed model is also verified by PPP experiments using 379 IGS stations data during another strong geomagnetic storm occurred on 26 August 2018. The analysis indicates that compared with conventional PPP scheme, the proposed model can improve the positioning accuracy by approximately 9.7% (19.8%) and 14.6% (39.7%) in the horizonal and vertical components for global (high latitudes) stations on 26 August 2018.