Analysis of the BDGIM Performance in BDS Single Point Positioning

: The broadcast ionospheric model is mainly used to correct the ionospheric delay error for single-frequency users. Since the BeiDou global ionospheric delay correction model (BDGIM) is a novel broadcast ionospheric model for BDS-3, its performance was analyzed through single point positioning (SPP) in this study. Twenty-two stations simultaneously receiving B1C, B2a, B1I and B3I signals were selected from the International GNSS Service (IGS) and the International GNSS Monitoring and Assessment System (iGMAS) tracking networks for the SPP experiments. The differential code bias (DCB) parameters were used to correct the hardware delays in the signals of B1C and B2a. The results showed that the BDGIM performs the best in high-latitude areas, and can effectively improve the positioning accuracy compared with the Klobuchar model. The average 3D positioning accuracy of the four civil signals can reach 3.58 m in high-latitude areas. The positioning accuracies with the BDGIM in the northern hemisphere are better than those in the southern hemisphere, and the global average 3D positioning accuracy of the four civil signals is 4.60 m. The performance of the BDGIM also shows some seasonal differences. The BDGIM performs better than the Klobuchar model on the days of spring equinox and winter solstice, while the opposite is true on the days of summer solstice and autumn equinox. On the day of winter solstice, the average 3D accuracies with the BDGIM on the signals of B1C, B2a, B1I and B3I are 4.13 m, 5.32 m, 4.40 m and 4.49 m, respectively. Although the SPP accuracies are to some extent affected by the geomagnetic storm, the BDGIM generally performs better and are more resistant to the geomagnetic storm than the Klobuchar model. with the BDGIM in the southern hemisphere are not as good as their counterparts in the northern hemisphere. The differences between the accuracies with the BDGIM and the GIM are 1.11 m, 2.16 m and 3.79 m in high-, mid- and low-latitude areas, respectively.


Introduction
Ionosphere is one of the main error sources affecting the signals of global navigation satellite system (GNSS), and the magnitude of ionospheric effects could reach tens of meters [1,2]. In GNSS data processing, the ionospheric error should not be negligible, and efforts have been made to avoid its adverse effects on positioning. Although the ionosphere-free combination is usually employed by dual-frequency users to eliminate the first-order ionospheric term, it always leads to significant increase in the noise level [3,4]. When only single-frequency measurements are available, the ionospheric error has to be corrected with proper ionospheric models [2]. The choice of ionospheric models is of great importance in obtaining precise positioning results. The American global positioning system (GPS) and the European Galileo navigation satellite system (Galileo) adopt the Klobuchar model and improved NeQuick-G model, respectively [1,5,6]. An approach was proposed to reduce the computational load of the NeQuick-G algorithm by increasing the update interval of the ionospheric corrections [7]. China's BeiDou navigation satellite system (BDS) adopts both the Klobuchar model and the BeiDou global ionospheric delay correction model (BDGIM) [2,6,8]. It is worth noting that the Klobuchar models with GPS and BDS are different in the term of coordinate system. The Klobuchar model with BDS is referenced to the geographic coordinate system, while that with GPS is to the geomagnetic coordinate system [9,10].
Due to the simplicity in calculation and the performance in mid-latitude areas, the Klobuchar model has been widely used and extensively studied. By improving the nighttime term and the amplitude of the cosine term, Bi et al. [11] proposed a modified Klobuchar model without any additional coefficients and the new model can reduce the ionospheric error by 60% over the polar region. Based on the k-means clustering of ionospheric daily variations, a correction model was proposed and was able to improve the performance of Klobuchar model by 36.24% in mid-latitude areas [12]. Liu et al. [13] adjusted the number of model parameters to adapt areas of different scales, and the regional correction rates of the Klobuchar-like models exceeded 90%. Considering the applicability and the issue of broadcasting, the traditional Klobuchar model with eight parameters are still the most popular.
The information of the eight-parameter Klobuchar model used by BDS-2 is contained in D1 and D2 navigation messages and is broadcast every two hours [14]. The D1 navigation message is broadcast by the B1I signals of the medium Earth orbit (MEO) and the inclined geo synchronous orbit (IGSO) satellites, while the D2 navigation message is broadcast by the B1I signals of the geostationary Earth orbit (GEO satellites) [14]. Adopted by BDS-3, the BDGIM contains nine parameters, which are broadcast through the navigation messages by the B1C and B2a signals, respectively [15,16].
To evaluate the performance of ionospheric parameters of BDS, Zhang et al. [17] took the global ionospheric map (GIM) model released by the Center for Orbit Determination in Europe (CODE) as a reference and found that more than 70% of the ionospheric error could be corrected by the BDS Klobuchar model, which performs better in the northern hemisphere and in the mid-latitude areas [10]. Over the polar area, the eight-parameter Klobuchar model of BDS is remarkably outperformed by that of the GPS [18]. The initial performance assessment by Yuan et al. [8] shows that the BDGIM is able to mitigate the ionospheric errors by 80.9% in China and by 77.6% on the global scale, outperforming both the GPS Klobuchar model and the NeQuick Galileo model. With precise ionospheric grid data released by the CODE and slant ionospheric delay derived from dual-frequency observations, Zhou et al. [10] conducted comparison among the performances of the BDS Klobuchar model, the GPS Klobuchar Model and the BDGIM. It was found that the BDGIM in China region is comparable to the BDS Klobuchar model, and that the BDGIM exhibit an overall advantage on a global scale. Similar results were achieved by Yu et al. [2], and the global accuracy of the BDGIM is about 3.6 total electron content units (TECU). Wang et al. re-estimated the BDGIM parameters and compared the vertical total electron content (VTEC) with the VTECs from the GIM and Jason-2/3 altimetry missions to find that the overall correction capability of the BDGIM was better than 75% for over 90% of all the observed samples [19].
Since the adoption of ionospheric model could lead to an improvement in the precision of positioning with single-frequency observables, the results of single point positioning (SPP) are used to evaluate the performance of different ionospheric models [9,[20][21][22]. The SPP experiments based on GPS L1 C/A code pseudoranges suggest that the positioning performance with the BDGIM is comparable to the NeQuick-C corrections and is better than GPS ionospheric correction algorithm in both quiet and disturbed conditions [19]. The SPP results of B1 and B2 code measurements with the BDS Klobuchar model are 26.3% and 28.2% better than those with the GPS Klobuchar model [23]. For the users in the China region, the accuracy of SPP with BDS Klobuchar model is better than that with GPS Klobuchar model [24][25][26]. The opposite results are also achieved due to the heavy multipath and poor orbit accuracy of the GEO satellites, and a regional satellite clock model and a regional ionospheric model are able to help improve the SPP accuracy Remote Sens. 2021, 13, 3888 3 of 19 by 50% [27]. For BDS-3 users, SPP accuracy with the BDGIM is overall better than that with the BDS Klobuchar model, especially in vertical component [10,28]. However, the performances of both the BDGIM and the BDS Klobuchar model are undermined in some extreme ionospheric environments, such as magnetic storm and polar region [6,18,29].
The existing studies with the BDGIM performance are mainly based on either the GPS code pseudoranges or the legacy BDS signals, while more comprehensive research on the BDGIM with BDS-3 satellites and the novel signals are still necessary. This work is to further investigate the performance of the BDGIM in term of SPP accuracy. From the users' perspective, we not only compare the different ionospheric models, but also pay intensive attention on the applicability of the BDGIM to all the four civil signals of BDS. Immediately after this introduction section, the algorithm of ionospheric correction with the BDGIM is introduced, and the data processing strategies are described in detail. In the following section, the performances of SPP with the four civil signals of BDS-3 are analyzed, and the BDGIM are compared with the BDS Klobuchar model and the global ionosphere maps (GIM). In addition, the performance of SPP with the BDGIM during magnetic storm is discussed. This is followed by the conclusions.

The Mathematical Model of SPP with the BDGIM
Only code measurements are used for SPP, and the ionospheric delay can be significantly reduced through the linear combination of dual-frequency measurements. We focus on single-frequency scenarios so as to investigate the performance of the ionospheric models.

Observation Equation of SPP
Generally, the code measurement of BDS can be modeled as [30] where, s, r and j indicate the numbers or names of satellite, receiver and frequency, respectively, ρ denotes the geometric distance between satellite and receiver, T denotes the tropospheric delay, dt r and dτ s denote receiver clock and satellite clock, respectively, I is the ionospheric delay, d and D are the hardware delays at receiver and satellite ends, respectively, and ε represents the sum of the unmodeled errors such as observing noise and multipath effects. The variable of epoch time is omitted, and all the terms are in the dimension of length. When no confusion is caused, the superscripts and subscripts might be omitted in the following discussion without additional instruction.
To form the observation equation, the geometric distance is linearized through Taylor expansion at the approximate coordinate, and the corrections of coordinates will be estimated. The parameter of receiver clock is also to be estimated, and the hardware delay at receiver end will be lumped into it. Therefore, the observation equation can be rewritten as where ρ s 0r is the computed geometric distance between the satellite and approximate coordinate of the receiver, and d τ s is the satellite clock computed from broadcast ephemeris. Since the BDS satellite clock in the broadcast ephemeris is referenced to the B3I signal, there is d τ s = dτ s + D s B3I and the differential code biases (DCBs) δ s j should be considered for the other signals. ι x v x , ι y v y and ι z v z are the corrections of coordinates multiplied by their coefficients. d t r is the sum of the receiver clock and corresponding hardware delay at the receiver end. In our experiment, the tropospheric delay is computed with empirical model, and so is the ionospheric delay which will be discussed below in details. The left-hand side is actually the observed-minus-computed (OMC) values.

Correction of Hardware Delay
As is mentioned above, the hardware delay of the B3I code measurement has already been included in the broadcast satellite clock, and therefore no additional correction is needed, meaning δ s B3I equals to zero. For the other signals, the code hardware delays at satellite end have to be corrected either by the timing group delay (TGD) parameters in broadcast ephemeris or post-time DCB products.
According to the specification of the format RINEX 3.04, the TGD parameter for B1I signal is recorded in the broadcast ephemeris of BDS, while the corrections for B1C and B2a signals are not. DCB products can be used as the substitute of the TGD parameters, since their equivalence has been verified through both theoretic induction and experimental and analysis [30][31][32]. Taking the pilot components of the B1C and B2a signals as examples, the DCB parameters δ s B1C and δ s B2a in Equation (2) are calculated as where ∆ s C1P,C6I is the DCB between the C1P and C6I signals for satellite s, and ∆ s C1P,C5P is the DCB between the C1P and C5P signals. These two DCB parameters are those which can be looked up in the DCB products.

Algorithm of the BDGIM
Similar to the tropospheric delay, the ionospheric delay for single-frequency SPP users is also computed with empirical model, and it is dependent on the total electron content (TEC) and the signal frequency. The spherical harmonic functions are able to effectively describe the subtle variations in TEC, and are widely used in the modeling and prediction of global TEC in solar geomagnetic frame [33].
Based on the modified spherical harmonics method, the algorithm of the BDGIM is detailed in the interface control documents (ICDs) for open service signals B1C and B2a, released by China Satellite Navigation Office (CSNO) in December 2017 [15,16]. The BDGIM coefficients are divided into two groups [2,10]. The nine low-order spherical harmonic function coefficients are broadcast through the navigation messages, while the other 17 coefficients, along with the calculation method, are provided in the ICDs. The nine broadcast coefficients, playing a dominant role in the prediction of TEC, are calculated with the global ionospheric background prediction model and the ground tracking data of BDS master control stations and A/B monitoring stations in China [34]. The 17 non-broadcast coefficients are supposed to describe the small-scale variations of TEC [2].
With the BDGIM, the line-of-sight ionospheric delay (in meters) of the signal with frequency f can be calculated as where α i (i = 1 ∼ 9) are the nine broadcast coefficients, and f is the carrier frequency of the signal. The M F is the ionospheric mapping function for the conversion from vertical to slant TEC. The values of A i are calculated as [15,16] In Equations (5) and (6), ϕ and λ are the geomagnetic latitude and longitude of ionospheric pierce point (IPP) in the solar-fixed reference frame. a k,j and b k,j are the nonbroadcast coefficients of BDGIM, and T k are the corresponding periods for prediction. t p is the odd hour of the modified Julian day. P |n i |,|m i | are the un-normalized Legendre functions with degree |n i | and order |m i |, and N n i ,m i are the normalization functions. According to the ICDs, the values of n i and m i are listed in Table 1 [15,16].
The mapping function in Equation (4) can be expressed as where, E is the elevation angle of the satellite, R e is the mean radius of the Earth, and H ion is the altitude of the ionospheric single-layer shell. Shown in Figure 1 is the flow chart for computing the ionospheric delay with the BDGIM.
In Equations (5) and (6), ′ and ′ are the geomagnetic latitude and longitude of ionospheric pierce point (IPP) in the solar-fixed reference frame.
, and , are the non-broadcast coefficients of BDGIM, and are the corresponding periods for prediction.
is the odd hour of the modified Julian day. | |,| | are the un-normalized Legendre functions with degree | | and order | |, and , are the normalization functions. According to the ICDs, the values of and are listed in Table 1 [15,16].
The mapping function in Equation (4) can be expressed as where, is the elevation angle of the satellite, is the mean radius of the Earth, and is the altitude of the ionospheric single-layer shell. Shown in Figure 1 is the flow chart for computing the ionospheric delay with the BDGIM.

Time t
Geographic longitude and latitude of station

Processing Strategies
At present, the available signals of BDS include B1I, B3I, B1C and B2a, and it is necessary to investigate the adaptability of BDGIM on different signals. Since the performance At present, the available signals of BDS include B1I, B3I, B1C and B2a, and it is necessary to investigate the adaptability of BDGIM on different signals. Since the performance of BDGIM can be reflected in single-frequency positioning results, experiments of SPP with different civil signals are conducted.

Data Selection
Data was collected during the days of year (DOYs) 23-46, 2020 (i.e., from 23 January to 15 February 2020) at 22 global tracking stations of the International GNSS Service (IGS) multi-GNSS experiment (MGEX) and the international GNSS Monitoring and Assessment System (iGMAS). During the 24 days, the Kp index stayed around 2.0, the 3-h forecast value never exceeded 4.0. The distribution of the stations is shown in Figure 2, and the detailed information of the 22 tracking stations is listed in Table 2. All of the four civil signals are available at each station. The broadcast ephemeris with RINEX format Version 3.04 are obtained from the website of the IGS, and the coefficients of both the BDS Klobuchar model and the BDGIM are provided by the Test and Assessment Research Center (TARC) of CSNO.    The numbers of visible BDS satellites and the position dilution of precision (PDOP) at each station are shown in Figure 3. The average numbers of BDS-3 and BDS-2/3 satellites

SPP Strategies
When code measurements are used for SPP, the satellite equipment group delay ought to be considered. The equipment group delay of B3I is regarded as the reference and is included in the clock parameter. Therefore, no additional correction is needed for B3I code measurement. The differential equipment group delay between B1I and B3I signals is given as the TGD parameters broadcast in the navigation message. The signals of B1C and B2a contain pilot components and data components. For the pilot components of B1C and B2a signals, the equipment group delays also need to be compensated by the TGD parameters. If data component is used, the group delay differential between the data component and the corresponding pilot component has to be additionally compensated by the inter-signal correction (ISC) parameters besides the TGD parameter. Both the TGD and the ISC parameters are broadcast in the navigation message.
Since the TGD and the ISC parameters for B1C and B2a signals are not recorded in the current ephemeris files, the differential code bias (DCB) parameters were used as the substitute to correct the equipment delays in B1C and B2a code measurements. In the study, we adopted the multi-GNSS DCB products released by China Academy of Sciences (CAS).
The accuracy of the GIMs provided by the IGS analysis centers could reach 2-5 TECU, equivalent to a zenith delay of 0.32-0.80 m in GPS code measurements of L1 frequency [35,36]. Therefore, the GIM released by CODE were employed as the control. The GIM and broadcast ionospheric models were used in SPP, respectively, and the positioning results were compared and analyzed. The data processing strategies are listed in Table  3, and the 95th percentile of the absolute values of positioning errors was adopted to measure the daily positioning accuracy of each station. According to the user range accuracy

SPP Strategies
When code measurements are used for SPP, the satellite equipment group delay ought to be considered. The equipment group delay of B3I is regarded as the reference and is included in the clock parameter. Therefore, no additional correction is needed for B3I code measurement. The differential equipment group delay between B1I and B3I signals is given as the TGD parameters broadcast in the navigation message. The signals of B1C and B2a contain pilot components and data components. For the pilot components of B1C and B2a signals, the equipment group delays also need to be compensated by the TGD parameters. If data component is used, the group delay differential between the data component and the corresponding pilot component has to be additionally compensated by the inter-signal correction (ISC) parameters besides the TGD parameter. Both the TGD and the ISC parameters are broadcast in the navigation message.
Since the TGD and the ISC parameters for B1C and B2a signals are not recorded in the current ephemeris files, the differential code bias (DCB) parameters were used as the substitute to correct the equipment delays in B1C and B2a code measurements. In the study, we adopted the multi-GNSS DCB products released by China Academy of Sciences (CAS).
The accuracy of the GIMs provided by the IGS analysis centers could reach 2-5 TECU, equivalent to a zenith delay of 0.32-0.80 m in GPS code measurements of L1 frequency [35,36]. Therefore, the GIM released by CODE were employed as the control. The GIM and broadcast ionospheric models were used in SPP, respectively, and the positioning results were compared and analyzed. The data processing strategies are listed in Table 3, and the 95th percentile of the absolute values of positioning errors was adopted to measure the daily positioning accuracy of each station. According to the user range accuracy (URA) information in the broadcast ephemeris and the specifications in the ICDs, the a priori precisions of BDS code measurements were set as 2.40 m and they were weighted with the elevation angles [14][15][16]. For the selected stations, the coordinates extracted from the officially issued software independent exchange (SINEX) files were taken as the reference true position for comparison.

Results
The global performances of SPP with the four BDS civil signals were investigated, and the SPP results with the BDGIM, the BDS Klobuchar model and the GIM were compared. To further study the adaptability of the BDGIM in different seasons, the positioning performances on four typical days were analyzed. In addition, the performance of the BDGIM during the geomagnetic storm was also analyzed.

Global Positioning Accuracies
Shown in Figures 4 and 5 are the accuracies in each direction (i.e., east, north and up) and the three-dimensional (3D) accuracies of SPP with different models and different signals. In each panel, arranged in the far left is the northernmost station, and the other stations are sorted rightwards by latitude from north to south. In Figure 4, the results of SPP with the BDGIM, the Klobuchar model and the GIM are displayed in the panels of left, middle and right columns, and each row of panels represent the results with each of the four BDS signals. The SPP accuracies in east, north and up directions are displayed as red, green and blue bars, respectively. In Figure 5, the 3D accuracies of SPP with the BDGIM, the Klobuchar model and the GIM are displayed as red circles, green stars and blue squares, respectively. Generally, the horizontal accuracies are significantly better than the vertical for all the stations except CANB, and the east components are usually the most accurate among all the three directions. The SPP accuracies of most of the stations in the northern hemisphere are better than their counterparts in the southern hemisphere, which might be explained by the relatively abundant monitoring stations in the northern hemisphere. With more stations in a certain area, a larger amount of data can be used to monitor the ionosphere, and the model parameters are theoretically able to describe the ionospheric variations much better. It is also worth noting that around the polar regions both the BDGIM and the GIM perform better than the Klobuchar model.
To     From the perspective of average SPP accuracies, the Klobuchar model performs better at the low-latitude stations, while the BDGIM performs significantly better at the high-latitude stations. At the mid-latitude stations, the BDGIM performs better than the Klobuchar model for all the signals except B1I. Among the four signals, the BDGIM almost always performs the best for the B1C signal, either regionally or globally. At the high-latitude stations, the highest average positioning accuracy, reaching 2.98 m, is obtained with the BDGIM on the B1C signal, and the average positioning accuracy of the four signals reaches 3.58 m. Compared with the Klobuchar model, the average SPP accuracies at the high-latitude stations with the BDGIM on the signals of B1C, B2a, B1I and B3I are improved by 78.2%, 83.2%, 70.5% and 79.4%, respectively. When the Klobuchar model is adopted, the average SPP accuracy at the high-latitude stations seems the worst, especially at the two arctic stations, CNYR and KIRU, whose accuracies exceed 10 m, as are shown in Figures 4 and 5.
At the mid-latitude stations, the BDGIM performs 5.5%, 19.8% and 7.3% better than the Klobuchar model for the signals of B1C, B2a and B3I, respectively, while is 5.7% outperformed by the Klobuchar model only for the signal of B1I. For all of the four civil signals, the Klobuchar model outperforms the BDGIM at the low-latitude stations by around 30%. For each signal, the best accuracy of SPP with the BDGIM is obtained at the high-latitude stations, and the SPP accuracy seems to increase with the latitude.
To investigate the difference between the two hemispheres, each of the abovementioned station groups is divided into the northern and the southern subgroups. Shown in Figure 6   The sharp deterioration of the positioning accuracy with the Klobuchar m high-latitude areas can be explained from the following two aspects. In the Klobuchar model, the value of the parameter α is relatively large, and the latitude point is not constrained. The combined effect of the two factors leads to signif crease in the calculated amplitude in high-latitude areas. To avoid these adverse the BDGIM is based on the method of modified spherical harmonics and calcu amplitude with normalized Legendre function. Figure 7 presents the typical time series of positioning errors at the station C 25 January 2020 (DOY 25), and the three models are used to correct the ionospher in all of the four signals. For the signals of B2a and B3I, the BDGIM performs sign better than the Klobuchar model, while for the signals of B1C and B1I the perfo of the two models are at the same level. Although the positioning accuracies with are the highest among the three models, the GIM cannot be obtained in real tim The sharp deterioration of the positioning accuracy with the Klobuchar model in highlatitude areas can be explained from the following two aspects. In the BeiDou Klobuchar model, the value of the parameter α is relatively large, and the latitude of pierce point is not constrained. The combined effect of the two factors leads to significant increase in the calculated amplitude in high-latitude areas. To avoid these adverse impacts, the BDGIM is based on the method of modified spherical harmonics and calculates the amplitude with normalized Legendre function. Figure 7 presents the typical time series of positioning errors at the station CLGY on 25 January 2020 (DOY 25), and the three models are used to correct the ionospheric delays in all of the four signals. For the signals of B2a and B3I, the BDGIM performs significantly better than the Klobuchar model, while for the signals of B1C and B1I the performances of the two models are at the same level. Although the positioning accuracies with the GIM are the highest among the three models, the GIM cannot be obtained in real time. Since the GIM is a post-time high-precision product with a certain time delay, it is not as convenient as the broadcast ionospheric models.

Positioning Accuracies in Different Seasons
To

Positioning Accuracies in Different Seasons
To investigate the performance of the BDGIM in different seasons, four typical days (i.e., 20 March 2020, 20 June 2020, 22 September 2020 and 22 December 2019) were selected to represent spring, summer, autumn and winter. The average SPP accuracies of the 22 stations with the four signals on the four typical days are shown in Figure 8. Note that the

Positioning Accuracies in Different Seasons
To investigate the performance of the BDGIM in different seasons, four typical days (i.e., 20 March 2020, 20 June 2020, 22 September 2020 and 22 December 2019) were selected to represent spring, summer, autumn and winter. The average SPP accuracies of the 22 stations with the four signals on the four typical days are shown in Figure 8. Note that the values larger than 10 m are not displayed. It can be seen that the BDGIM generally performs the best on the day of winter solstice and the worst on the day of autumn equinox. On the day of winter solstice, the average 3D accuracies of the signals of B1C, B2a and B3I are 4.13 m, 5.32 m and 4.49 m, respectively, and they are better than those of all the other three days. For the signal of B1I, the 3D accuracies on the three days other than autumn equinox are almost the same, i.e., around 4.40 m, and that on the day of autumn reaches 6.31 m.

Positioning Accuracies during the Geomagnetic Storm
A geomagnetic storm can not only cause strong disturbances in the geomagnetic field, but also influence the ionosphere. According to the alerts issued by the Space Environment Prediction Center of CAS, the 3-h forecast Kp index reached 5.0 twice on 20 April, 2020 (DOY 111), indicating an event of a geomagnetic storm. To investigate the performance of the BDGIM around the geomagnetic storm, we conducted SPP experiments with the BDS data collected on the DOYs 110-112, 2020.
Shown in Figure 9 are the average positioning accuracies of the four signals before, during and after the occurrence of the geomagnetic storm. Generally, the geomagnetic storm has perceptible but slight influence on SPP with the BDGIM, and the two novel signals and the zenith direction are affected more obviously. On the very day of the geomagnetic storm, the positioning accuracies of B1C and B2a signals decrease by 5.6% and 6.7% with respect to the previous day. In addition, the positioning accuracies of the two days before and after the occurrence of the geomagnetic storm are almost the same. For the signals of B1I and B3I, the 3D accuracies stay around 5 m all through the three days, and the horizontal accuracies are around 2 m. It seems that the positioning accuracies with the two legacy signals are resistant to the influence of the geomagnetic storm.

Positioning Accuracies during the Geomagnetic Storm
A geomagnetic storm can not only cause strong disturbances in the geomagnetic field, but also influence the ionosphere. According to the alerts issued by the Space Environment Prediction Center of CAS, the 3-h forecast Kp index reached 5.0 twice on 20 April, 2020 (DOY 111), indicating an event of a geomagnetic storm. To investigate the performance of the BDGIM around the geomagnetic storm, we conducted SPP experiments with the BDS data collected on the DOYs 110-112, 2020.
Shown in Figure 9 are the average positioning accuracies of the four signals before, during and after the occurrence of the geomagnetic storm. Generally, the geomagnetic storm has perceptible but slight influence on SPP with the BDGIM, and the two novel signals and the zenith direction are affected more obviously. On the very day of the geomagnetic storm, the positioning accuracies of B1C and B2a signals decrease by 5.6% and 6.7% with respect to the previous day. In addition, the positioning accuracies of the two days before and after the occurrence of the geomagnetic storm are almost the same. For the signals of B1I and B3I, the 3D accuracies stay around 5 m all through the three days, and the horizontal accuracies are around 2 m. It seems that the positioning accuracies with the two legacy signals are resistant to the influence of the geomagnetic storm.  The averages of accuracies with the four civil signals at each station on the day of the geomagnetic storm (DOY 111) are listed in Table 8, and the positioning accuracies with the BDGIM and the Klobuchar model are compared. The mean accuracies of the 22 stations are also calculated. It can be seen that the mean 3D accuracy with the BDGIM is 6.22 m, much better than that with the Klobuchar model. Considering the poor performance The averages of accuracies with the four civil signals at each station on the day of the geomagnetic storm (DOY 111) are listed in Table 8, and the positioning accuracies with the BDGIM and the Klobuchar model are compared. The mean accuracies of the 22 stations are also calculated. It can be seen that the mean 3D accuracy with the BDGIM is 6.22 m, much better than that with the Klobuchar model. Considering the poor performance The averages of accuracies with the four civil signals at each station on the day of the geomagnetic storm (DOY 111) are listed in Table 8, and the positioning accuracies with the BDGIM and the Klobuchar model are compared. The mean accuracies of the 22 stations are also calculated. It can be seen that the mean 3D accuracy with the BDGIM is 6.22 m, much better than that with the Klobuchar model. Considering the poor performance of the Klobuchar model at the three high-latitude stations (CNYR, KIRU and ZHON), the averages without the three stations are also calculated and compared. Even after the three high-latitude stations are excluded, the mean SPP accuracies with the BDGIM and the Klobuchar model are 6.44 m and 7.66 m, respectively. The results show that the BDGIM performs better than the Klobuchar model even during the geomagnetic storm.

Conclusions and Discussion
In the study, the performance of the BDGIM in SPP was assessed, and the DCB parameters were used to correct the hardware delays in the signals of B1C and B2a. To quantitatively analyze the SPP accuracies, 24 days of data collected at 22 global stations were processed, and the positioning accuracies were investigated and compared.
On the global scale, the average accuracy of SPP with the BDGIM can reach 4.60 m, which is close to that with the GIM. Compared with Klobuchar model, the positioning accuracies of the four signals are improved by 15.4%-37.0% when the BDGIM is adopted. In high-latitude areas, the BDGIM performs much higher than the Klobuchar model. With the BDGIM, the positioning accuracy of the B1C signal can reach 2.98 m in high-latitude areas, and the average positioning accuracy of the four civil signals is 3.58 m.
In mid-latitude areas, the BDGIM outperforms the Klobuchar model by 5.5%, 19.8% and 7.3% for the signals of B1C, B2a and B3I, respectively. In low-latitude areas, the Klobuchar model performs better than the BDGIM by around 30%. Seeming to increase with the latitude, the accuracies of SPP with BDGIM are similar to those with the GIM in the northern hemisphere, and their counterparts in the southern hemisphere are not as good.
The positioning performance of the BDGIM show seasonal differences and the general SPP accuracy is the best in winter and the worst in autumn. The BDGIM performs better than the Klobuchar model on the days of spring equinox and winter solstice, while the opposite is true on the days of summer solstice and autumn equinox. On the day of winter solstice, the average 3D accuracies with the BDGIM on the signals of B1C, B2a, B1I and B3I are 4.13 m, 5.32 m, 4.40 m and 4.49 m, respectively, i.e., 16.1%, 27.9%, 6.2% and 20.5% better than those with the Klobuchar model.
Although the geomagnetic storm can to some extent affect the accuracies of SPP with the BDGIM, the influences are more perceptible only for the two novel signals and the zenith direction. The accuracies of SPP with the two legacy signals are less affected. On the very day of the geomagnetic storm, the average 3D accuracy of SPP with the BDGIM is 6.22 m, and is much better than that with the Klobuchar model.
In general, the BDGIM performs better in the northern hemisphere and in high-latitude areas than in the southern hemisphere and in low-latitude areas. Seasonal differences are also seen in the performance of the BDGIM, and it outperforms the Klobuchar model on the days of spring equinox and winter solstice. Despite some slight influences, the BDGIM are more resistant to the geomagnetic storm than the Klobuchar model. The advantages of the BDGIM are probably attributed to the spherical harmonics and more favorable distribution of the monitoring stations used in the generation of model parameters. It is also worth noting that the above results and analysis are derived from a relatively limited data set during the time of solar minimum. To evaluate the long-term performance of the BDGIM, further investigation will be conducted in the future.