Performance Analysis of GPS/BDS Broadcast Ionospheric Models in Standard Point Positioning during 2021 Strong Geomagnetic Storms

: The broadcast ionospheric model is one of the main methods for eliminating ionospheric delay errors for the Global Navigation Satellite Systems (GNSS) single-frequency users. GPS Klobuchar model (GPSK8) is the widely used broadcast ionospheric model for GPS, while BDS usually implements the BDS Klobuchar model (BDSK8) and BeiDou Global Broadcast Ionospheric Delay Correction Model (BDGIM). Geomagnetic storms may cause interference within the ionosphere and near-Earth space, compromising the accuracy of ionospheric models and adversely affecting the navigation satellite systems. This paper analyzes the static Standard Point Positioning (SPP) accuracy of GPS and BDS by implementing the broadcast ionospheric models and then investigates the impact of strong geomagnetic storms occurring in 2021 on positioning accuracy. The results show that the global 3D positioning accuracy (95%) of GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM are 3.92 m, 4.63 m, and 3.50 m respectively. BDS has a better positioning accuracy in the northern hemisphere than that of the southern hemisphere, while the opposite is valid for GPS. In the mid-latitude region of the northern hemisphere, BDS + BDSK8 and BDS + BDGIM have similar positioning accuracy and are both better than GPS + GPSK8. The positioning accuracy after applying those three broadcast ionospheric models shows the superior performances of winter and summer over spring and autumn (based on the northern hemisphere seasons). With the exception of during winter, nighttime accuracy is better than that of daytime. The strong geomagnetic storm that occurred on the day of year (DOY) 132, 2021 has an impact on the positioning accuracy for only a small number of stations; however, the global average positioning accuracy is not signiﬁcantly affected. The strong geomagnetic storms that occurred in DOY 307 and DOY 308 have a signiﬁcant impact on the positioning accuracy of dozens of stations, and the global average positioning accuracy is affected to a certain extent, with some stations experiencing a serious loss of accuracy. Decreased degrees in positioning accuracy is proportional to the intensity of the geomagnetic storm. Of the 33 IGS Multi-GNSS Experiment (MGEX) stations worldwide, those located in the low and mid-latitudes are more signiﬁcantly affected by the geomagnetic storms compared with higher latitudes. Evident ﬂuctuations of the positioning errors existed during the strong geomagnetic storms, with an increase in extreme values, particularly in the up direction


Introduction
The rapid development of GNSS has led to a wide range of applications in areas such as precise surveying, vehicle navigation, intelligent agriculture, machinery control, and maritime rescue.However, GNSS positioning accuracy might be compromised in an environment with poor observation conditions, such as urban canyons, geomagnetic storms, tunnels, heavy rainfall, etc. [1][2][3][4].The ionospheric delay is one of the significant error sources in GNSS navigation and positioning [5].The ionosphere is an essential part of the Earth's atmosphere, located between 60 and 1000 km above the Earth's surface.Under extreme ultraviolet, ultraviolet, X-rays, etc., the ionosphere contains a large number of free electrons and positively charged ions.An important source of these charged particles is the ionization of the corresponding thermosphere composition.During geomagnetic storms [6], enormous amounts of energy are injected into the upper atmosphere employing enhanced currents, electric fields, and energetic particle fallout.This is accompanied by a complex response of the ionosphere and thermosphere to the geomagnetic storm through the propagation of energy-momentum and thermospheric-ionospheric coupling [7,8].Due to the prevalent thermosphere-ionosphere coupling process, perturbations in the thermosphere may affect the behavior of the ionosphere (both positive and negative ionospheric responses) during geomagnetic storms through wind field transport and compositional changes [9,10].The ionospheric response can be severe enough to cause satellite signals to be interrupted [11][12][13].Therefore, as for the single-frequency navigation users, the ionospheric delay needs to be corrected using the broadcast ionospheric model, for instance, GPS Klobuchar Model (GPSK8) for GPS, BDS Klobuchar Model (BDSK8) and BDGIM for BDS, and NeQuickG model for Galileo.The Klobuchar model is widely applied for its computational simplicity.The reference frame used for the GPSK8 is a geomagnetic coordinate system, while the BDSK8 is a geographic coordinate system, and the ionospheric single layer heights for both models are 350 km and 375 km respectively.The BDGIM is a global ionospheric model based on a spherical harmonic function that uses non-broadcast and broadcast parameters to compute the ionospheric delay corrections.The NeQuickG is a three-dimensional global model capable of calculating the electron density at any given height in space and can also give total electron content (TEC) on the signal propagation path by numerical integration [14].The input parameters of the NeQuickG model include time, position, and solar activity index; meanwhile, the NeQuickG model replaces the solar activity index with an effective ionization level factor and uses a quadratic polynomial to fit three ionospheric parameters that are broadcast in the broadcast ephemeris [15].
Many scholars have assessed and analyzed the performances of the broadcast ionospheric models to improve the positioning accuracy of single-frequency users.By investigating the performance of the GPSK8, BDSK8, and NeQuickG models in the distance delay and position domains in China, Wang et al. [16] found that those three models could eliminate 64.8%, 65.4%, and 68.1% of the ionospheric delay, respectively.According to the evaluation of GPSK8 versus BDSK8 global standard point positioning accuracy by Wu et al. [17], in the northern hemisphere, including Asia, Europe, and North America, the 3D positioning accuracy was improved by 7.8-35.3%using BDSK8 compared to GPSK8, but the accuracy was reduced due to the lack of BDS monitoring stations in the southern hemisphere.Zhao et al. [18] performed a preliminary assessment of BDGIM using the International GNSS Monitoring and Assessment System (iGMAS) data from multiple stations and showed that BDGIM improved the performance of Klobuchar at night and in polar regions, with a global ionospheric correction rate of nearly 80%, providing the best global description of the ionosphere and performing optimally over most of the northern hemisphere.According to Yang et al. [19], the NTCM-BC, BDGIM, and GPSK8 can reduce 72.6%, 69.8%, and 61.7% of ionospheric delays globally, respectively.Wang et al. [20] carried out the experimental analysis of the performance of BDGIM in SPP.Compared with BDSK8, the positioning accuracy of the four signals could be improved by 15.4% to 37.0% when utilizing BDGIM, and the SPP accuracy was slightly affected by geomagnetic storms.Though the Klobuchar model has been improved and optimized by various scholars [21][22][23][24][25][26] and many new models have been proposed [27][28][29][30][31], the advantages of the original model in terms of simplicity of structure and easy implementation still make it widely used.In addition to improving the accuracy of the broadcast ionospheric model, the accuracy of the SPP can be improved by increasing the accuracy of satellite orbits and clock errors in the broadcast ephemeris or by adding additional augmentation information [32][33][34][35][36][37].
Geomagnetic storms, caused by the enhanced solar wind and the interactions with magnetosphere-ionosphere-thermosphere systems, are massive disturbances to the Earth's near-space environment [38].During geomagnetic storms, large amounts of energetic particles are fed into the Earth's magnetosphere.Some of these particles follow the Earth's magnetic field and then settle into the ionosphere at high latitudes, resulting in steep ionospheric density gradients and irregularities [39,40].Liu et al. [41] found that the geomagnetic storm that occurred on 22 February 2015 perturbed the ionosphere and coincided with ionospheric irregularities at mid-and high-latitudes, demonstrating that the perturbation of the rate of TEC index (ROTI) was caused by the geomagnetic storm.The Broadcast ionospheric model accuracy and GNSS positioning performance are significantly affected by the increased and fluctuating VTEC variations during geomagnetic storms [42,43].Xue et al. [44] investigated for the first time the accuracy of BDS B1 frequency standard point positioning during different levels of geomagnetic storms in and around China and showed that there was a decrease in positioning accuracy.Poniatowski and Nykiel [45] reported a significantly large root mean square (RMS) of position errors for several tens of centimeters for dynamic precision point positioning during the geomagnetic storm of 17 March 2015.Alcay [46] selected 12 intense geomagnetic storms for PPP experiments between 2000 and 2018 and showed that the degree of degradation of the dynamic PPP accuracy depends on the strength of TEC fluctuations.It has no significant effect on PPP positioning accuracy when TEC anomalies are at a relatively low level.The occurrence probability of positioning error extreme values is positively correlated with geomagnetic storm level and is generally greater at low latitudes than at high latitudes [47].In addition, satellite signals could be lost due to geomagnetic storms [48].
The Kp index is a widely used geomagnetic index.According to the National Oceanic and Atmospheric Administration (NOAA) space weather scale, a strong geomagnetic storm is considered to occur when the Kp index equals 7 [49].There were two strong geomagnetic storm events in 2021, on 12 May and 3 November, corresponding to the day of year (DOY) 132 and 307, respectively, and the DOY 307 geomagnetic storm lasted into DOY 308.Exploring the loss of standard point positioning accuracy during these three storms is useful to further understand the impact of geomagnetic storms on navigation satellite systems, and thus reduce such effects.In addition, it has been two years since BDS-3 was completely built, and thus assessing the SPP performance is of some practical importance to BDS users.According to Zhang et al. [50], the BDSK8 is capable of obtaining reliable SPP results on a global scale.This paper first evaluates the global standard point positioning accuracy of GPSK8, BDSK8, and BDGIM during the absence of strong geomagnetic storms in 2021, then analyses the impact on the SPP accuracy of broadcast ionospheric models during strong geomagnetic storms, gives experimental results and discussions, and finally draws conclusions.

Materials and Methods
The determination of the receiver position by distance intersection using the satellite orbits and satellite clock errors provided by the broadcast ephemeris, as well as the code pseudo-range observations, is called Standard Point Positioning (SPP).The observation equations for SPP are given below, along with three broadcast ionospheric models.

Observation Equation of SPP
Generally, the pseudorange observation equation is illustrated as follows: where P is the pseudorange measurements, ρ is the geometric distance between satellite and receiver, dt r and dt s denote receiver and satellite clock error in meters respectively, T is the tropospheric delay, I is the ionospheric delay, d r and d s denote the receiver and the satellite differential code bias (DCB), ε is the noise error and the multi-path error.Measurement noise is generally mitigated by improving hardware performance and extending the observation time.Multipath errors are related to the observation time and station environment.They are difficult to correct with the traditional GNSS mathematical model.In data processing, measurement noise and multipath effects are usually treated together as random measurement noise and attenuated using methods such as adjustment and filtering.Time group delays in the broadcast ephemeris were extracted, converted, and utilized to compute the satellite DCB.The geometric distance in the above equation is expanded by a Taylor series and truncated appropriately to give the following linearised observation equation: where l, m, n are the directional cosines of the approximate position of the station to the satellite, and v x , x y , v z are the correction values on the three position directions, respectively.The constant term is: where ρ 0 is the distance between the position of the satellite at the moment of transmission and the approximate position of the receiver at the moment of arrival of the signal.When four satellites are observed simultaneously, the system of equations can be solved directly to obtain the receiver coordinate increments and receiver clock error; when more than four are observed, the least squares adjustment is required.

Broadcast Ionosphere Models
Broadcast Ionospheric Model (BIM) is an ionospheric delay correction model for single-frequency GNSS users.The BDSK8 is an ionospheric regional model, mainly for the Asia-Pacific region, and GPSK8 and BDGIM are global models [51].Table 1 lists the key differences between the three models.The GPSK8 uses the assumption of a thin ionosphere.It is calculated in the geomagnetic coordinate system using an eight-parameter cosine function expression to reflect the daily ionospheric amplitude and phase variations.It treats the nighttime ionospheric TEC as a constant and specifies that the total number of electrons in the ionosphere peaks at 14:00 daily.The GPSK8 is given by [52] where T iono refers to the ionospheric time delay at L1 frequency.The equations for the other symbols in Equation ( 4) are as follows where F refers to the obliquity factor, α n and β n are the satellite transmitted data words with n = 0, 1, 2 and 3. φ m is the geomagnetic latitude of the Earth projection of the ionospheric intersection point (mean ionospheric height assumed 350 km).

BDS Klobuchar Model (BDSK8)
The BDSK8 is an improved model derived from the Klobuchar model, which is mainly applicable to the Asia Pacific region.The model parameters are obtained by solving the measured dual frequency data from the Chinese regional monitoring network.The user calculates the vertical ionospheric delay correction with BDSK8 as follows [53]: where I (t) is the vertical ionospheric delay in seconds for B1I, t is the local time (range 0~86,400 s) for the place under the intersection point (M) of the ionosphere and the direction from receiver to satellite.A 2 is the amplitude of the Klobuchar cosine curve in the daytime computed from the α n .
A 4 is the period of the cosine curve in seconds, which is computed from the β n .
where, φ M is the geographic latitude of Earth projection of the ionosphere intersection point in radians.I (t) can be converted to the ionospheric delay along the B1I propagation path I B1I (t) through the equation as follows and the unit is seconds.
where, R is the mean radius of the Earth (6378 km), E is the satellite elevation from the user's location in radians, h is the height of the ionosphere (375 km).

BeiDou Global Broadcast Ionospheric Delay Correction Model (BDGIM)
The BDGIM is based on the modified spherical harmonics method.The model consists of broadcast and non-broadcast parameters, nine broadcast parameters are transmitted in the broadcast ephemeris and the non-broadcast parameters are solidified at the user end.In this case, the number of parameters is reduced while retaining the contribution of the higher order spherical harmonic coefficients as much as possible, ensuring the global TEC descriptive power of the spherical harmonic function.According to the BDGIM, the user shall compute the ionospheric delay correction by using the equation as follows [54]: where, T ion is the line-of-sight (LOS) ionospheric delay along the signal propagation path from the satellite to the receiver (in meters).M F is the ionospheric mapping function for the conversion between vertical and slant TEC, f is the carrier frequency of the current signal (in Hertz); α i (i = 1 ∼ 9) are the BDGIM parameters (in TECU); A i is calculated by the geomagnetic latitude and longitude of the ionospheric puncture point (IPP) in the solar-fixed reference frame, A 0 is the predictive ionospheric delay (in TECU).Refer to the BDS interface control documents for detailed calculation processes.

Experiment Process 2.3.1. Data Selection
To analyze the SPP performance by applying the GPSK8, BDSK8, and BDGIM, 33 MGEX stations globally distributed were selected for the experiment, as shown in Figure 1.The following factors were taken into account when selecting the stations: (1) they were evenly distributed across the globe at different latitudes; (2) the station receivers were capable of receiving satellite signals from GPS, BDS-2, and BDS-3 simultaneously; and (3) the quality of the observed data could satisfy the requirements of this experiment.A total of 123 days of Observation data from January, April, July, and October 2021 were collected.According to the sun storm alert issued by the Space Environment Prediction Center of CAS, the Kp index was below 7 during the selected period, and no strong geomagnetic storm occurred.The average number of visible satellites and the average Position Dilution of Precision (PDOP) values for GPS and BDS at an elevation cut-off angle of 10 • were obtained from the January 2021 observation data, as is shown in Figure 2.  The number of visible GPS satellites at each station is mostly between 7 and 13, the PDOP values are below 2.5, and the variation between stations in different regions is not obvious.The overall number of visible BDS satellites is larger than that of GPS.This is because BDS currently has 7 GEO satellites and 10 IGSO satellites distributed in the Asia Pacific region, and therefore more BDS satellites can be observed in the eastern hemisphere than that of the western hemisphere.Although the number of visual BDS satellites is less for stations in South America, the PDOP values do not exceed 4. Observation data and GPS broadcast ephemeris used in the experiment were obtained from the IGS Data Centre of Wuhan University and the BDS broadcast ephemeris from the Test and Assessment Research Center of China Satellite Navigation Office.

SPP Strategies
All stations are almost static in the Earth-Centered-Earth-Fixed (ECEF) system, and thus it is not a dynamic scene.Static single-frequency SPP experiments use the L1 signal from GPS and B1I signal from BDS, respectively.The GPSK8 was used only for GPS, while the BDSK8 and BDGIM were used for BDS.The Saastamoinen model was utilized consistently for the tropospheric delay correction.The hardware delay is corrected using the TGD parameters in the broadcast ephemeris, and the detailed positioning strategy is presented in Table 2. To assess the accuracy of the SPP, the station coordinates from the IGS Software Independent EXchange (SINEX) products are extracted as the reference values, and the 95th percentile of the absolute value of the positioning error is used to measure the daily positioning accuracy of each station, and the average value over the study period is considered as the final positioning accuracy.

Results
The SPP performance of GPSK8, BDSK8, and BDGIM was investigated, comparing the results for different latitudes and different seasons and periods.In addition, the positioning performance of three broadcast ionospheric models during strong geomagnetic storms was analyzed.

Solar and Geomagnetic Conditions
In order to analyze the positioning performance during geomagnetic storms, January + April + July + October 2021 is referred to as the reference period.The Kp and F10.7 indices for the reference period are shown in Figure 3.The F10.7 is overwhelmingly below 110 sfu and the Kp index is overwhelmingly below 6, which shows the low level of solar and geomagnetic activity during the reference period.The Kp index during geomagnetic storms is shown in Figure 4.There are two "7" in DOY 132 and one "7" in DOY 307.In DOY 308, except for two "7", there are two "6" and one "5".This indicates that the geomagnetic storm activity was more intense on DOY 308 than on the other two days.The F10.7 indices for DOY 132, 307, and 308 are 75 sfu, 89 sfu, and 94 sfu, respectively, showing an increasing trend, but still at low activity levels.

Global Positioning Accuracies
Standard point positioning is currently at the meter level of accuracy, subject to the accuracy of the pseudo-range measurements, the accuracy of the satellite orbits and clock, and atmospheric delay errors.Figure 5 shows the average positioning accuracy for each station over the study period, from top to bottom for east (E), north (N), up (U) direction, and 3D positioning accuracy, with all stations sorted by latitude from north to south.The E and N directions of the three positioning solutions are generally more accurate than the U direction, most of the E and N directions are below 4 m and most of the U directions are below 6 m.
To analyze the differences in SPP accuracy at different latitudes, this paper considered 0-30 • as low latitude, 30-60 • as mid-latitude, and 60-75 • as high latitude, and then calculated the average horizontal and vertical accuracy of stations at different latitudes in the northern and southern hemispheres, as shown in Table 3. BDS + BDGIM performs best both globally and in China/Asia-Pacific region, and GPS performs best in the southern hemisphere, which is related to broadcast ionospheric model performance and modeling data sources, satellite visual conditions, etc. BDS has better positioning accuracy in the northern hemisphere than that of the southern hemisphere, while GPS has the opposite.The global 3D accuracy of GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM is 3.92 m, 4.63 m, and 3.50 m respectively.BDS + BDSK8 and BDS + BDGIM have the best positioning accuracy at mid-latitudes in the northern hemisphere and GPS + GPSK8 at mid-latitudes in the southern hemisphere.
Although the BDSK8 is an ionosphere regional model, it still has a usable positioning accuracy on a global scale, similar to BDS + BDGIM in the Chinese region with 3.84 m and 3.58 m respectively, which is better than GPS + GPSK8 with 5.23 m.  3. Statistical results of the average positioning accuracy of stations in different latitudes (unit: m).N stands for Northern Hemisphere, S for Southern Hemisphere, H for High Latitude (60~75 • ), M for Mid Latitude (30~60 • ), L for Low Latitude (0~30 • ), H for Horizontal, V for Vertical, and 3D for three dimensions.

Positioning Accuracies in Different Seasons and Periods
The positioning performance was further analyzed for different seasons and periods in 2021, with January, April, July, and October representing winter, spring, summer, and autumn respectively; with 08:00-20:00 local time each day as the daytime and the rest as the nighttime, the statistical analysis of the average 3D accuracy of all stations was carried out and the results are shown in Figure 6.Except for winter, all positioning solutions showed better nighttime accuracy than that the of daytime because of the absence of solar illumination, lower ionospheric electron content, and quieter ionospheric activity at night.
For GPS + GPSK8, the 3D positioning accuracy was 4.46 m, 3.80 m, 4.72 m, and 4.19 m for spring, summer, autumn, and winter, respectively; 5.10 m, 3.87 m, 5.26 m, and 4.34 m for BDS + BDSK8; and 3.89 m, 3.20 m, 3.96 m and 3.70 m for BDS + BDGIM.The positioning accuracy is better in winter and summer than that in spring and autumn.The global positioning performance of BDS + BDGIM is the best in all seasons.

Positioning Accuracies during Strong Geomagnetic Storms
Strong geomagnetic storms occurred in DOY 132, 307, and 308 in 2021, and SPP experiments were conducted using observation data from these three days.The three geomagnetic storms are DOY 308, DOY 307, and DOY 132 in descending order of degree.The average positioning accuracy for the four months of 2021 obtained in the above paper is regarded as the positioning accuracy for the reference period.
Daily positioning accuracy was counted and the statistics of GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM are shown in Tables 4-6 respectively, with the stations in the tables sorted by latitude from north to south and the average results for all stations in the last row of each table.The improvement in positioning accuracy during strong geomagnetic storms compared with the reference periods was calculated, and the results are presented in Figure 7.It can be seen that strong geomagnetic storms have a significant degradation of positioning accuracy for most stations and that the three strong storms produce different effects.Compared with the reference period, the strong geomagnetic storm of DOY 132 caused a maximum of 1.34 m (37.9%), 0.93 m (25.4%), and 1.84 m (43.8%) decrease in 3D accuracy for GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM respectively, which occurred at ABPO, ORID and OHI3 respectively.However, according to the average results of all stations, the geomagnetic storms on this day did not have a significant impact on the overall positioning results of the global stations.The strong geomagnetic storms occurring in DOY 307 and DOY 308 had a relatively more pronounced impact, with a greater number of stations affected and an increased level of impact.Compared with the reference period, the geomagnetic storms of DOY 307 caused a maximum of 2.29 m (46.0%), 4.54 m (102.3%), and 4.53 m (107.6%)decrease in 3D accuracy for the three positioning solutions, which occurred at LMMF, HKWS, and HKWS, respectively.According to the average results of the global stations, the 3D positioning accuracy of the three positioning solutions was reduced by 11.1%, 8.7%, and 12.9% respectively.
The strong geomagnetic storm of DOY 308 caused a maximum 3.05 m (61.2%), 4.43 m (90.8%), and 6.42 m (136.8%)decrease in 3D accuracy for the three positioning solutions, respectively, compared with the reference period, all at the LMMF station.According to the average results of the global stations, the 3D positioning accuracy of the three positioning solutions decreased by 10.8%, 12.6%, and 33.2% respectively.
To investigate the distribution of stations with significant accuracy degradation, the stations where the 3D positioning accuracy reduced by more than 30% in DOY 132, 307, and 308 are shown in Figure 8 for the three positioning solutions.The impact of the strong geomagnetic storm in DOY 132 was not significant, the vast majority of stations were not significantly affected, and the accuracy of BDS + BDSK8 was not degraded by more than 30% at all stations, while for GPS + GPSK8 and BDS + BDGIM, there were one and two stations respectively where the accuracy of the position was reduced by 30~50%.Strong geomagnetic storms of DOY 307 and DOY 308 had a significant impact on some stations.For DOY 307, there were 5, 6, and 6 stations with an accuracy decline of more than 30% for the three positioning solutions respectively, mainly in the middle and low latitudes of Asia and South America; for DOY 308, there were 3, 11, and 15 stations with accuracy reduction of more than 30% for the three positioning solutions respectively, mainly in the middle and low latitudes outside North America.To further analyze the positioning errors during strong geomagnetic storms, the positioning error time series of DOY 100, DOY 132, DOY 307, and DOY 308 are plotted in Figure 9, for which DOY 100 is a representative of the reference period.The positioning accuracy in each direction is marked in Figure 9.It can be found that the positioning error is relatively smooth when no geomagnetic storms occur, at DOY 132 the positioning error is slightly affected, and at DOY 307 and 308 the extreme values of the positioning error increase, and the error series shows significant fluctuations, more pronounced in the U-direction.
In summary, the impact of the strong geomagnetic storm in DOY 132 on SPP was local and controllable, with a few stations showing some degradation in accuracy and no significant degradation in global average accuracy; the impact of the geomagnetic storms in DOY 307 and DOY 308 was obvious and severe, with more stations showing significant degradation in accuracy at mid-and low-latitudes and some degradation in global average accuracy.In addition, there were significant fluctuations in the positioning error during the storms, and the probability of extreme values increased significantly, especially in the U direction.

Discussion
Through the experiments, we have found that the SPP accuracy of the same positioning solution is variable at different latitudes, which is mainly related to the ionospheric differences at different latitudes.The VTEC tends to increase roughly from high latitudes to low latitudes, where direct solar radiation makes the ionospheric activity more intense and the Equatorial Ionization Anomaly (EIA) occurs during the day, and broadcast ionospheric models are less accurate in this region, resulting in reduced positioning accuracy.At high latitudes, the mathematical shortcomings of the Klobuchar model have a detrimental effect, the SPP accuracy of GPSK8 and BDSK8 at high latitudes is not ideal.For the BDSK8, the accuracy is higher only at mid-latitudes in the northern hemisphere thanks to the limitations of the modeling data source.Although, BDGIM is limited by the modeling data source, the superiority of the BDGIM is rooted in the ability of the spherical harmonic function in describing the physical field of the spherical shell, and thus the BDGIM performs better than the BDSK8 at high latitudes.Temporal variation of global TEC causes variations in positioning accuracy across seasons and periods.With an 11-year solar cycle and 2021 being a low solar year, global TEC is typically lower in winter and summer than that in spring and autumn, a phenomenon known as an ionospheric semiannual anomaly, thus the positioning accuracy is higher in winter and summer than that in spring and autumn.Generally, nighttime positioning accuracy is better than that of daytime.The day and night positioning accuracies are similar in winter, due to the relatively high levels of ionospheric TEC at night.In summary, the SPP accuracy is closely related to ionospheric variability and broadcast ionospheric model performance.
It was found that strong geomagnetic storms have a significant impact on SPP.We checked the auroral electrojet index during geomagnetic storms.According to the World Data Center (WDC) for Geomagnetism, the real-time AE index reached a maximum of 1500 nT in DOY 132 and stayed at 500~1500 nT from 12:00 to 16:00.In DOY 307, it took a maximum of 2000 nT and remained in 500~2000 nT from 21:00 to 23:00.For DOY 308, it exceeded 2000 nT several times and kept in 500~2000 nT for most of the time between 01:00 and 15:00.It is clear that the geomagnetic storms from DOY 307 to 308 produced more severe and prolonged effects.Therefore, DOY 307 and DOY 308 show a more severe degradation of SPP accuracy than DOY 132.In addition, the 95th percentile of the absolute value of the positioning error is chosen as the positioning accuracy evaluation index, which is more sensitive to the extremes of the positioning error, and the increase of extremes during geomagnetic storms will be reflected in the final accuracy results.
Typically, the ionosphere is affected by upper-level solar activity, geomagnetic activity as well as thermosphere coupling, and the lower neutral atmosphere.During geomagnetic storms, the magnetosphere electric field penetrates the ionosphere, changing the electron density and causing a response in the ionosphere.The level of perturbation of the geomagnetic field is a useful proxy for the level of perturbation of other geospatial phenomena.As the most significant space weather event, geomagnetic storms can also cause serious hazards to communication and power systems, near-Earth spacecraft, etc.Therefore, it is critical to monitor geomagnetic disturbances, explore their variability and provide accurate alerts.

Conclusions
In the study, the performance of GPS and BDS broadcast ionospheric models in SPP was evaluated.33 MGEX stations globally distributed and a total of four months of observation data in 2021 were selected to analyze the positioning performance at different latitudes, seasons, and periods around the globe.
The global 3D positioning accuracy (95%) of GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM are 3.92 m, 4.63 m, and 3.50 m respectively.BDS + BDSK8 and BDS + BDGIM have the best positioning accuracy at mid-latitudes in the northern hemisphere, while GPS + GPSK8 has the best positioning accuracy at mid-latitudes in the southern hemisphere.Although the BDSK8 is a regional model, it still has usable positioning accuracy on a global scale.BDS + BDSK8 has similar positioning accuracy with BDS + BDGIM at mid-latitudes in the northern hemisphere and is preferable to GPS + GPSK8.All three positioning solutions, except winter, show better nighttime positioning accuracy than daytime.In addition, the winter and summer positioning accuracy are better than the spring and autumn.For GPS + GPSK8, the positioning accuracy was To analyze the positioning accuracy during strong geomagnetic storms, the SPP experiments of DOY 132, DOY 307, and DOY 308 of 2021 were carried out and the results showed that the degree of decrease in positioning accuracy during geomagnetic storms is proportional to the degree of geomagnetic storms.
Compared with the reference period, the strong geomagnetic storm of DOY 132 caused a maximum of 1.34 m (37.9%), 0.93 m (25.4%), and 1.84 m (43.8%) decrease in 3D accuracy for GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM respectively, and no degradation in global average positioning accuracy.
The strong geomagnetic storms occurring in DOY 307 and DOY 308 had a relatively more pronounced impact, with a greater number of stations affected and an increased

Figure 1 .
Figure 1.Distribution of the 33 MGEX stations used in this study.

Figure 2 .
Figure 2. Average number of visible satellites and average PDOP values for GPS and BDS in January 2021 with a cut-off elevation of 10 • .(a) Average numbers of Visible GPS satellites; (b) Average PDOPs of GPS satellites; (c) Average numbers of Visible BDS satellites; (d) Average PDOPs of BDS satellites.

Figure 5 .
Figure 5. Average east (E), north (N), up (U) and 3D positioning accuracy for each station during the reference period, i.e., January + April + July + October 2021.The stations in the figure are sorted by latitude from north to south.

Figure 7 .
Figure 7. Improvement in 3D positioning accuracy for all stations during strong geomagnetic storms compared to reference periods for GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM.When the value is negative it means a decrease in accuracy.

Figure 8 .
Figure 8. 3D positioning accuracy decreases in percentage during strong geomagnetic storms (only stations with more than 30% decrease in accuracy are plotted), from top to bottom three rows are GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM, and from left to right three columns are the accuracy decreases of DOY 132, DOY 307 and DOY 308, it should be noted that BDS + BDSK8 does not have stations with more than 30% decrease in accuracy in DOY 132.

Figure 9 .
Figure 9.Time series of positioning errors at LMMF for GPS + GPSK8, BDS + BDSK8, and BDS + BDGIM with DOY 100 as a reference period and strong geomagnetic storms at DOY 132, 307, and 308.In each panel, the time series and statistical results correspond to the E, N, and U directions, denoted in red, blue, and green respectively.

Table 2 .
The data processing strategies of static single-frequency SPP.