*Article* **Performance of BDS B1 Frequency Standard Point Positioning during the Main Phase of Different Classified Geomagnetic Storms in China and the Surrounding Area**

**Junchen Xue 1,\*, Sreeja Vadakke Veettil 2, Marcio Aquino 2, Xiaogong Hu 1, Lin Quan 3, Dun Liu 4, Peng Guo <sup>1</sup> and Mengjie Wu <sup>1</sup>**


**Abstract:** Geomagnetic storms are one of the space weather events. The radio signals transmitted by modern navigation systems suffer from the effects of magnetic storms, which can degrade the performance of the whole system. In this study, the performance of the BeiDou Navigation Satellite System (BDS) B1 frequency standard point positioning (SPP) in China and the surrounding area during different classes of storm was investigated for the first time. The statistical analysis of the results revealed that the accuracy of the BDS-2 B1 frequency SPP deteriorated during the storms. The probability of the extrema of the positioning error statistics was largest during strong storms, followed by moderate and weak storms. The positioning accuracy for storms of a similar class was found not to be at the same level. The root mean square error in positioning for the different classes of storm could be at least tens of centimeters in the east, north and up directions. The findings in this study could contribute toward the error constraint of BDS positioning accuracy during different classes of geomagnetic storm and be beneficial to other systems, such as BDS-3, as well.

**Keywords:** BDS; SPP; main phase; geomagnetic storms

### **1. Introduction**

A geomagnetic storm is defined as a period when the ring current becomes intense enough to exceed the key threshold of the disturbance storm time (Dst) index. Geomagnetic storms are induced by the intense and continuing interplanetary convection electric field and energy injection into the magnetosphere–ionosphere system [1]. The enhanced interplanetary convection electric field is motivated by a constant southward interplanetary magnetic field (IMF) [2]. Solar wind carries the coronal magnetic field out into the entire heliosphere, thus forming the IMF [3]. Based on the signatures in the magnetic field, a geomagnetic storm can be divided into three phases: initial, main and recovery. The main phase is the principal characteristic of a geomagnetic storm [1,4].

The largest global atmospheric effects can be activated by geomagnetic storms [5]. Storms can generate disturbances in the ionosphere, which varies with the location of the region under consideration, the local time (LT) of the geomagnetic storm onset and other parameters [6]. The Equatorial Ionization Anomaly also responds to the effects of storms [7]. The disturbed condition of the ionosphere during geomagnetic storms is known as an ionospheric perturbation, which can have great effects on radio propagation-dependent applications, especially for Global Navigation Satellite System (GNSS) single frequency users. These effects can usually be corrected by ionospheric models. In the case of BeiDou

**Citation:** Xue, J.; Veettil, S.V.; Aquino, M.; Hu, X.; Quan, L.; Liu, D.; Guo, P.; Wu, M. Performance of BDS B1 Frequency Standard Point Positioning during the Main Phase of Different Classified Geomagnetic Storms in China and the Surrounding Area. *Remote Sens.* **2022**, *14*, 1240. https://doi.org/10.3390/rs14051240

Academic Editor: Yunbin Yuan

Received: 25 January 2022 Accepted: 1 March 2022 Published: 3 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Navigation Satellite System (BDS-2) single frequency users, a Klobuchar-style ionospheric navigation model is applied for the corrections [8]. The mean correction precision of the model is better than 65% and it performs better for middle latitudes than low latitudes [9]. For BDS-3 single frequency users, the BeiDou global broadcast ionospheric delay correction model (BDGIM) is proposed for the ionospheric delay correction. Ionospheric errors can be mitigated by 80.9% in the China region and 77.6% at the global scale [10].

There are case studies demonstrating the adverse effects of geomagnetic storms on GNSS positioning. During storms, carrier phase cycle slips in Global Positioning System (GPS) signals may occur, which result from loss of lock (LOL) events [11]. Astafyeva et al. [12] demonstrated that the density of GPS LOL events can increase to 0.25% on the L1 band and 3% on the L2 band during super storms (Dst ≤ −250 nT) and 0.15% on L1 and 1% on L2 during intense storms (−250 nT < Dst ≤ −100 nT). Specifically, the tracking performance of GPS receivers in the high latitudes was investigated for the 2015 St. Patrick's Day strong storm. The significant scintillation caused by the storm contributed to a LOL on the GPS L2 band but had little influence on the tracking of the GPS L1 signal [13]. Kinematic GPS positioning could also be degraded during ionospheric disturbances induced by geomagnetic storms. The repeatability of kinematic positioning, which was estimated using a two-step approach based on double difference L3 phase measurements, reached 12.8, 8.1 and 26.1 cm in the north (N), east (E) and up (U) directions, respectively, during the 2003 Halloween storm [14]. The accuracy of real-time kinematic (RTK) positioning could also be deteriorated during a strong geomagnetic storm [15] and even during a weak storm at high latitudes [16]. Furthermore, positioning errors in network RTK and precise point positioning (PPP) techniques increased rapidly during the 2015 St. Patrick's Day strong storm [17]. An investigation was also performed into the effects of moderate and weak geomagnetic storms on the performance of GNSS–SBAS in low-latitude African regions using an SBAS emulator to simulate specific EGNOS-like messages. The SBAS performance in equatorial African regions showed a non-linear relationship with the geomagnetic storm indices [18]. Additionally, GPS instrumental biases, including receiver and satellite biases, are routinely estimated using a dual frequency geometry-free combination. These computations can also be affected during geomagnetic storms [19].

Even though previous studies have revealed the effects of individual, or several, geomagnetic storms on the Earth's upper atmosphere and GNSS applications, few papers have studied the performance of BDS-based applications during geomagnetic storms, especially during the different classes of storm. GNSS single frequency users are supposed to be more obviously affected compared to other positioning modes, such as PPP, during those periods. In this study, the effects of different types of storm on the performance of BDS-2 single frequency standard point positioning (SPP) are investigated comprehensively, especially for the most commonly used B1 frequency. In addition, the differences between the effects of separate storms are studied.

#### **2. Methodology**

Dst index can be used as a criterion to classify geomagnetic storms [4]. In this study, Dst indices were extracted from combined OMNI files that were obtained from the NASA database (https://omniweb.gsfc.nasa.gov, accessed on 20 January 2022). All storms in solar cycle 24 were analyzed and divided into three classes: strong, moderate and weak. Table 1 states the threshold conditions that were applied in the classification of storms (see [1,20]).

**Table 1.** The thresholds implemented in the classification of geomagnetic storms. Type refers to the storm classes, Dst is in nT and ΔT is in hours.


The basic strategy for selecting storms was that the Dst should be as minimal as possible and the duration of each storm should be more than 12 h. To ensure that each storm was independent and not influenced by another storm, a condition was applied that the Dst index for the ten days before and after the main phase day must be greater than the minimum value for each individual class of storm. Finally, five cases with noticeable main phases were chosen for each class of storm from 2015 to 2018 [20]. The principal characteristic of a geomagnetic storm is the main phase [4]. The details of the main phases of the storms, including the related Dst values, the start and end epoch and the duration, are presented in Table 2 (see [20]); Dates with a suffix of "0" refer to the start epoch while a suffix of "1" represents the end epoch. The same meanings are applied to the suffix used for the Dst values.

The earliest period for collecting BDS-2 observations from the Multi-GNSS Experiment (MGEX) network [21] in this study was 2015. BDS-2 observations were obtained for the chosen stations on the dates listed in Table 2. Those stations from the MGEX network, namely DAEJ, GMSD, JFNG, LHAZ, HKWS and HKSL, are distributed throughout China and the surrounding area (see Figure 1). The sampling interval was 30s. The related information, such as geodetic coordinates, receiver and antenna version, is shown in Table 3. The last two columns show the dates that the hardware, such as the receivers or antennae, was changed or updated.

**Figure 1.** The distribution of the selected MGEX stations throughout China and the surrounding area: a red dot indicates the location a station; the black dotted line indicates the geomagnetic equator. The latitude and longitude are in degrees).


*Remote Sens.* **2022** , *14*, 1240

**Table 2.** The main phases of the different classes of

geomagnetic

 storm from 2015 to 2018: Type, storm classes; MJD, modified Julian date;


**Table 3.** Information regarding the longitude, latitude and receiver and antenna versions of the stations. Site refers to the station name and the latitude and longitude are in degrees.

The data were processed in the kinematic mode of SPP using BDS-2 single frequency pseudorange observations. Considering the dispersive nature of the ionosphere, only the B1 pseudorange was used here. The pseudorange observation equation is illustrated as follows:

$$B\_1 = \rho + dt\_I - dt^s + T + I\_1 + db\_{r1} - db^{s1} + \varepsilon \tag{1}$$

where *B*<sup>1</sup> is the BDS-2 B1 pseudorange observation, *ρ* is the geometric range, *dtr* is the receiver clock error, *dt<sup>s</sup>* is the satellite clock error, *T* is the tropospheric delay, *I*<sup>1</sup> is the ionospheric delay, *dbr*<sup>1</sup> is the receiver differential code bias (DCB), *dbs*<sup>1</sup> is the satellite DCB and *ε* is the noise error.

A conventional option was set for the SPP program. The satellite orbit and clock were computed from IGS navigation data. The tropospheric delays were derived using the Saastamoinen model. The ionospheric delays were calculated using the broadcasted BDS-2 navigation ionospheric model. Time group delays in the broadcast ephemeris were extracted, converted and utilized to compute the satellite DCB. For each epoch, the station coordinates and receiver clock error were estimated using the Gauss–Newton least square method. The weight was set with the satellite elevation angle. The elevation mask angle was set to 10◦.

As a result, the station coordinates in the Cartesian coordinate system were compared to the precise solutions from SINEX files obtained from MGEX products. Positioning errors were obtained in the east (E), north (N) and up (U) directions of the local station coordinate framework. The related statistics were performed for the main phase periods using indices such as minimum (MIN), maximum (MAX), bias and root mean square error (RMSE). The MIN and MAX represent the minimum and maximum of the positioning errors for the three directions. The bias and RMSE were computed from the positioning errors for each component as well. The formulas are demonstrated as follows:

$$\begin{aligned} \text{MIN} &= \text{minimum} \{ \Delta POS\_i \} \\ \text{MAX} &= \text{maximum} \{ \Delta POS\_i \} \\ \text{BIAS} &= < \Delta POS\_i > \\ \text{RMSE} &= \sqrt{< \Delta POS\_i^2 >} \\ \text{\Delta POS\_i &= POS\_{ref,i} - POS\_{est,i}, i = 1, n} \end{aligned} \tag{2}$$

where <> is the average of the variable, *POSref* ,*<sup>i</sup>* is the precise solution from the SINEX files, *POSest*,*<sup>i</sup>* is the solution obtained from this processing and *n* is the total number of samples.

#### **3. Results and Discussion**

The accuracy of BDS-2 B1 frequency SPP positioning during the main phase of different classes of storm is analysed in this section. First, the positioning errors in the three directions for all stations during a period of 3 days' representative of each class of storm are presented in Figures 2–4. The periods shown in each of the figures cover the three days before and after the main phases of the individual storms. In each of the figures, the Dst time series is applied to indicate the activities of the storms. The slant total electron content (TEC) of each station was computed and converted into slant time delays (DION) in the B1 frequency. The positioning errors of the selected GNSS stations in the E, N and U directions are shown below. The main phases of the geomagnetic storms are indicated with the red dash–dot lines. The recovery phase is the period after the epoch indicated by the right-hand dash–dot line. From the figures, we can see that the positioning errors during the main phases were clearly different from those in other periods. This implies that the positioning errors were influenced during the main phases of the storms. In addition, the degradation of the positioning errors could occur at random epochs of the main phase with the decrease in Dst values. This could last until the recovery phase, as shown in Figure 2. Furthermore, the variations of positioning errors in the U direction were more significant than in the other directions.

For the strong storm shown in Figure 2, the positioning errors for the selected GNSS stations in the E, N and U directions represent the fluctuations during the main phase of the storm, especially for the U direction in which errors could reach approximately 10 m. The fluctuations of the positioning errors for the different stations varied. In general, the positioning errors of the stations could reach about 2 m in the E direction, about 5 m in the N direction and about 11 m in the U direction. The characteristics of the ionospheric delays were changed during the main phase. The maximum delays and the related epochs were different from other periods. Moreover, in this event, a larger degradation of the positioning errors occurred at the beginning of the recovery phase. The reason for this can be attributed to the fact that the Dst values were lower, the geomagnetic disturbance remained intense throughout and the ionospheric error correction was insufficient.

It is worth noting that there were some sharp increases in the positioning errors of the LHAZ station in the E, N and U directions during the main phase of the storm. Nonetheless, there were no similar increases observed for the HKWS and HKSL stations, which are also located at lower latitudes. The DIONs shown in Figure 2 stayed at a low level, thus suggesting that the correction of ionospheric error could be sufficient. The increases were observed at nearly 1 LT and lasted for about 30 min, thus indicating that this occurred locally and temporarily. The time series of F10.7 cm radio flux was also checked for the same period to discover whether there could be sudden ionospheric disturbances. However, no sudden variations were observed and the values were under 75 sfu (1 sfu = 10−<sup>22</sup> w/m2/hz), thus implying that there were no abrupt changes in solar activity. Solar radio burst (SRB) events can also affect GNSS positioning [22] and, therefore, the data from the NOAA websites (see ftp://ftp.swpc.noaa.gov/pub/indices/events, accessed on 20 January 2022) were checked to detect any possible occurrences of SRB events. It was noticed that no SRB events occurred during that epoch.

Further, BKG Ntrip Client (BNC) software [23] was used for the data quality check at the LHAZ station. The number of satellites and the position dilution of precision (PDOP) were computed and are shown in Figure 5. It was found that there were no variations in the number of satellites during the periods that sharp increases in the PDOP values were observed. This indicates that the degradation was directly caused by the poor geometrical structure of the constellations. The orbits of most of the observed satellites were geostationary earth orbits (GEOs) and inclined geosynchronous orbits (IGSOs). There was only one medium earth orbit (MEO) satellite, named C14, which was at a high elevation during this time. These factors might have contributed to the sharp jumps in the PDOP. However, frequent jumps in the number of satellites were observed near the end of the whole jump period. On checking the PDOP values at other stations, no such jumps were observed in the same period. The reason for this could be attributed to receiver signal tracking issues during the storm. In addition, it can be seen that the PDOP values varied and most values were greater than 5 during the whole storm period, thus suggesting that the influence of this storm on PDOP was strong.

**Figure 2.** The time series of the positioning errors for BDS-2 B1 frequency during a strong storm around DOY 238, 2018. The times series from top to bottom are for Dst, DION and the positioning errors (POS Err) in E, N and U directions, respectively. The brown dotted line is for DAEJ, the green dotted line is for GMSD, the purple dotted line is for JFNG, the yellow dotted line is for LHAZ, the black dotted line is for HKWS and the cyan dotted line is for HKSL. The red dash–dot lines indicate the borders of the main phase, the left-hand line is the start epoch and the right-hand line is the end epoch. The gray dash–dot lines indicate the borders of the days. The *X*-axis is epoch (DOY) in GPST and the *Y*-axis is Dst in nT. DION and POS Err are in meters.

**Figure 3.** The time series of the positioning errors for BDS-2 B1 frequency during a moderate storm around DOY 086, 2017. The times series from top to bottom are for Dst, DION and the positioning errors (POS Err) in E, N and U directions, respectively. The brown dotted line is for DAEJ, the green dotted line is for GMSD, the purple dotted line is for JFNG, the yellow dotted line is for LHAZ, the black dotted line is for HKWS and the cyan dotted line is for HKSL. The red dash–dot lines indicate the borders of the main phase, the left-hand line is the start epoch and the right-hand line is the end epoch. The gray dash–dot lines indicate the borders of the days. The *X*-axis is epoch (DOY) in GPST and the *Y*-axis is Dst in nT. DION and POS Err are in meters.

**Figure 4.** The time series of the positioning errors for BDS-2 B1 frequency during a weak storm around DOY 032, 2017. The times series from top to bottom are for Dst, DION and the positioning errors (POS Err) in E, N and U directions, respectively. The brown dotted line is for DAEJ, the green dotted line is for GMSD, the purple dotted line is for JFNG, the yellow dotted line is for LHAZ, the black dotted line is for HKWS and the cyan dotted line is for HKSL. The red dash–dot lines indicate the borders of main phase, the left-hand line is the start epoch and the right-hand line is the end epoch. The gray dash–dot lines indicate the borders of the days. The *X*-axis is epoch (DOY) in GPST and the *Y*-axis is Dst in nT. DION and POS Err are in meters).

The moderate storm event indicated by the Dst time series is shown in Figure 3. The positioning errors of most of the stations varied in all directions during the main phase in comparison to other periods. The positioning errors were degraded to some extent during the main phase. The degradation in the U direction was more obvious than in the other two directions. The maxima of the positioning errors were about 2 m in the E direction, 4 m in the N direction and 10 m in the U direction. When comparing the maxima of the positioning errors from the selected moderate storm and strong storm, the former were lower than the latter. It was found that the characteristics of the positioning errors of HKWS and HKSL were different from those of the other stations. This could be attributed to the different version of receiver hardware used at those two stations (LEICA version, see Table 3). The ionospheric delays were slightly enhanced during the main phase compared to other periods. In addition, there were also changes in the positioning errors along the E, N and U directions during the recovery phase. This suggests the there were influences of the geomagnetic storm on the positioning errors when the Dst index was still at the lower level.

**Figure 5.** The number of satellites and the PDOP for the LHAZ station during a strong storm around DOY 238, 2018. The green dots indicate the number of satellites and the blue dots indicate the PDOP. The red dash–dot lines indicate the borders of the main phase, the left-hand line is the start epoch and the right-hand line is the end epoch. The gray dash–dot lines indicate the borders of the days. The *X*-axis is epoch (DOY) in GPST, the left *Y*-axis is the number of satellites and the right *Y*-axis is PDOP.

The time series of the positioning errors during the selected weak storm are shown in Figure 4. From the figure, we can see that the activity of this storm was at a low level in terms of the Dst time series. The ionospheric delays were at a reduced status as well. There were slight variations in the positioning errors in the three directions during the main phase compared to the other periods. The positioning errors could reach about 2 m in the E direction, 2 m in the N direction and 8 m in the U direction. The statistical results seemed lower than those from the strong and moderate storms. This suggests that the influence of the weak storm on the positioning errors was less than those of the strong and moderate storms. It can also be seen that the variation in the U direction was greater than in the other directions. The variations in the U direction were stronger near the end of the main phase when the Dst was at its minimum. Additionally, there were also fluctuations in the positioning errors along the U direction during the recovery phase.

Tables 4–6 show the statistics for the BDS-2 B1 frequency positioning errors during the main phases of the independent storms. From the statistical results, the probability of the extrema in the four statistical indices was largest during strong storms, overall followed by moderate and weak storms. The difference in positioning errors between different classes of storm was greatest in the U direction.

**Table 4.** The statistical indices for the BDS-2 B1 frequency positioning errors during the main phases of strong storms: MJD, modified Julian date; SITE, station name; MIN, MAX, BIAS and RMSE, statistical indices. The three columns for each index are for results in the E, N and U directions, respectively. The last two rows are the mean and median of the statistics in each column. All indices are in meters).


From Table 4, we can see that the maximum positioning error during the strong storms was up to 33 m, the minimum was nearly −12 m, the bias approached 5 m and the RMSE for the E, N and U directions reached 1.52, 3.25 and 7.92 m, respectively. The mean of all RMSEs was 0.95, 1.32 and 3.50 m in the E, N and U directions, respectively, whilst the median was 1.03, 1.12 and 2.96 m for the E, N and U directions, respectively. Furthermore, the positioning accuracy was not comparable between the different strong storms; the accuracy during some storms was better than during others. This indicates that not all strong storms had a similar influence on the positioning accuracy. The same feature was observed during the moderate and weak storms as well.

From Table 5, we can see that the maximum positioning error during moderate storms reached 12 m, the minimum was close to −11 m, the bias approximated 4 m and the RMSE could be up to 1.83, 1.87 and 5.40 m for the E, N and U directions, respectively. The mean and median for the RMSEs of the E, N and U directions were 0.94, 1.09, 3.21 m and 0.90, 1.07, 3.06 m, respectively.



From Table 6, we can see that the positioning errors during weak storms were generally lower than those observed during moderate and strong storms. However, it can be noticed from Table 6 that the LHAZ, HKWS and HKSL stations presented an irregular behavior on MJD 57717 (DOY 329, 2016). Among the three stations, the maximum positioning error reached 50.08 m in the N direction of HKSL, while the minimum reached −93.21 m in the N direction of HKWS. The corresponding bias and RMSE were also large. There were no such big changes for the GMSD and JFNG stations, although they lie in higher latitudes. The reasons for this irregular behavior in the positioning errors of the three stations are analyzed in details later. Except for the anomalous results, the maximum was up to 10 m, the minimum was −20 m, bias was approximately 3 m and the RMSE could reach 2.15, 1.95 and 4.87 m for the E, N and U directions, respectively.

The corresponding time series of the irregular positioning errors during the selected weak storm are showed in Figure 6. The ionospheric delays were slightly depressed during the main phase. The positioning errors of all stations showed large fluctuations at the beginning of the storm, which were more noticeable in the N and U directions. There were also minor fluctuations in the E direction. The LT corresponding to the start epoch of the storm was around 14 h, which is when the ionospheric activity is the most intense during a day. There were large jumps in the positioning errors, especially in the N and U directions. Furthermore, there were no positioning estimations in any direction for many epochs. The strongest jumps occurred at 20 LT and lasted until 8 LT the next morning. Moreover, there were other jumps at around 20 LT in the recovery phase. These jumps were clearly visible in all three directions.


**Table 6.** The statistical indices for the BDS-2 B1 frequency positioning errors during the main phases of weak storms: MJD, modified Julian date; SITE, station name; MIN, MAX, BIAS and RMSE, statistical indices. The three columns for each index are for results in the E, N and U directions, respectively. The last two rows are the mean and median of the statistics in each column. All indices are in meters.

It was initially supposed that the irregular behavior could be related to the fact that ionospheric activity at low latitudes is more intense than at high latitudes [24]. However, the influence of ionospheric activity could not reach such a level (tens of meters) after the correction of ionospheric model (see Figure 6 and reference [25]). This can also be proved by the positioning errors at other epochs of the main phase. Further, the number of satellites and the PDOP for the three stations are shown in Figure 7. It can be seen that there were large jumps in the PDOP. The maximum could be at the level of several hundred counts. Combined with the number of satellites, these jumps were directly caused by a loss of satellite tracking. A similar feature was observed on the next day (DOY 330, 2016) as well. There are many reasons for the failure of tracking satellites, such as issues with receiver hardware or software, signal strength, space weather, etc. In this study, the main reason for the jumps in the positioning errors for the three stations (LHAZ, HKWS and HKSL) could be attributed to a comprehensive effect. The issues with that specific receiver (LEICA version, see Table 3) might be the most possible reason, which caused similar jumps for these three stations. In addition, the time series of the related space weather indices were checked. The southward IMF Bz had high-frequency variations during the storm period. The solar wind speed increased before the end of the main phase and subsequently decreased. There were continual and consistent variations in Kp and Dst during that period. The AE had large jumps near the epoch of the minimum Dst. There were no sharp

changes in F10.7 or the solar wind pressure time series, nor any occurrence of SRBs events. The variations of the space weather indices imply that the storm activity became complex during this period.

**Figure 6.** The time series of the positioning errors for BDS-2 B1 frequency during the main phase around DOY 330, 2016. The times series from top to bottom are for Dst, DION and the positioning errors (POS Err) in E, N and U directions, respectively. The green dotted line is for GMSD, the purple dotted line is for JFNG, the yellow dotted line is for LHAZ, the black dotted line is for HKWS and the cyan dotted line is for HKSL. The red dash–dot lines indicate the borders of the main phase, the left-hand line is the start epoch and the right-hand line is the end epoch. The gray dash–dot lines indicate the borders of the days. The *X*-axis is epoch (DOY) in GPST and the *Y*-axis is Dst in nT. DION and POS Err are in meters.

After removing these three stations, the final statistics are illustrated in Table 7. As shown in the table, the mean of the RMSEs for the B1 positioning errors in the three directions was 0.92, 1.06 and 2.70 m and the median was 0.82, 0.96 and 2.60 m. The statistics were lower than those for strong and moderate storms.

Futhermore, the 3D RMSEs of each station were computed around the selected events. By comparing the 3D RMSEs of the pre-storm period and main phase, the decrease in positioning performance was derived. The mean relative percentages of the 3D RMSEs were 5.38%, 1.53% and 0.59% for the strong, moderate and weak storms in this study, respectively.

**Figure 7.** The number of satellites and the PDOP for the HKSL, HKWS and LHAZ stations (top to bottom) around DOY 330, 2016. The green dots indicate the number of satellites and the blue dots indicate the PDOP. The *X*-axis is GPST, the left *Y*-axis is the number of satellites and the right *Y*-axis is PDOP.

**Table 7.** The statistical indices for the BDS-2 B1 positioning errors during the main phases of weak storms without the anomalous stations: (MJD, modified Julian date; MIN, MAX, BIAS and RMSE, statistical indices. The three columns for each index are for results in the E, N and U directions, respectively. The last two rows are the mean and median of the statistics. All indices are in meters).


#### **4. Conclusions**

In this study, a statistical analysis of BDS-2 B1 frequency SPP positioning errors during the main phases of different classes of storm in China and the surrounding area was conducted. From the results, it was observed that the positioning accuracy was affected to different degrees during the storms. Some relevant conclusions can be drawn from the analyses. Firstly, the probability of the extrema of the positioning error statistics was greatest during strong storms, followed by moderate and weak storms. Secondly, during the same class of storm, the positioning accuracy could vary. Thirdly, the positioning accuracy could be influenced even during the recovery phase of a storm.

The findings of this study will hopefully contribute toward the error constraint of BDS positioning accuracy during different strengths of geomagnetic storms. Additionally, the influence of storms could be comparable to other GNSS systems; thus, the findings could also be beneficial to those systems. However, it should be noted that the analysis in this study could have been limited by the uncertainties of the error models of SPP, although the further study of those uncertainties would be beneficial to the improvement of SPP during these storms. In addition, since the study period was in the descending phase of solar cycle 24, the effects on the positioning accuracy might not be entirely apparent. Thus, the study needs to be extended to the beginning of solar cycle 25 and with the addition of more storm events. Investigations into BDS-3 applications, such as the BDGIM model [10], could be performed with more BDS observations as well. Further work could also be focused on case studies of the effects and analyses should be performed with combinations of geophysical indices, drivers and sources, etc.

**Author Contributions:** Conceptualization, J.X.; data curation, J.X.; formal analysis, J.X., S.V.V., M.A., X.H., L.Q. and D.L.; funding acquisition, J.X.; investigation, J.X., S.V.V. and M.A.; methodology, J.X.; project administration, J.X.; resources, J.X.; software, J.X.; supervision, S.V.V. and M.A.; validation, J.X., S.V.V. and M.A.; visualization, J.X.; writing—original draft, J.X.; writing—review & editing, J.X., S.V.V., M.A., X.H., L.Q., D.L., P.G. and M.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 11703066, and the China Scholarship Council, grant number 201704910002.

**Data Availability Statement:** The datasets analysed during this study are available in the IGS MGEX repository (ftp://cddis.gsfc.nasa.gov/pub, accessed on 20 January 2022), NASA OMNI (https://omniweb.gsfc.nasa.gov, accessed on 20 January 2022) and the NOAA database (ftp://ftp. swpc.noaa.gov/pub/indices/events, accessed on 20 January 2022).

**Acknowledgments:** The authors would like to thank the anonymous referees for their valuable suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

