Estimation and Analysis of BDS2 and BDS3 Differential Code Biases and Global Ionospheric Maps Using BDS Observations

: Following the continuous and stable regional service of BDS2, the BDS3 ofﬁcially announced its global service in July 2020. To fully take advantage of the new multi-frequency BDS3 signals in ionosphere sensing and positioning, it is essential to understand the characteristics of the differential code bias (DCB) of new BDS3 signals and BDS performance in global ionospheric maps (GIMs) estimation. This article presents an evaluation of the characteristics of 13 types of BDS DCBs and the accuracy of BDS-based GIM based on the data provided by the International GNSS Service (IGS) and International GNSS Monitoring and Assessment System (iGMAS) for the ﬁrst time. The GIMs and DCBs are estimated by the APM (Innovation Academy for Precision Measurement Science and Technology) method in a time efﬁcient manner, which can be divided into two main steps. The ﬁrst step is to produce GIMs based on BDS observations at the B1I, B2I and B3I signals, and the second step is to estimate DCBs among the other frequency bands by removing the ionospheric delay using the precomputed GIMs. Good agreement is found between the APM-based satellite DCB estimates and those from the Chinese Academy of Sciences (CAS) and the German Aerospace Center (DLR) at levels of 0.26 ns and 0.18 ns, respectively. The results, spanning one month, show that the stability of BDS DCB estimates among different frequency bands are related to the contributed observations, and the receiver DCB estimates represent larger STD values than the satellite DCB estimates. The differences in receiver DCB estimates between BDS2 and BDS3 are found to be related to the types of receivers and antennas and ﬁrmware version, and the bias of the JAVAD receivers reaches 1.03 ns. The results also indicate that the difference in the single-frequency standpoint positioning (SPP) accuracy using GPS-based and BDS-based GIMs for ionospheric delay corrections is less than 0.03 m in both the horizontal and vertical directions. between the APM-based BDS satellite DCB and daily DCBs by DLR and CAS validated in Section 3.2. To better understand the characteristics of the BDS3 satellite DCBs, mean values of the DCBs among the six different frequency bands covering all operational BDS3 satellites over a period of one month are shown in Figure 8. A close match of BDS3 DCBs between the same two frequency band signals can be noted from this ﬁgure. The DCBs between the B1I and B1C signals are conﬁned to a range of ± 5 ns, while the DCB values between B1I and other signals of the BDS3 satellites, with the exception of the PRN C33 satellite, are conﬁned to a range of − 31 to 55 ns.


Introduction
As the third global navigation satellite system (GNSS) following the great successes of the GPS and GLONASS systems, the BeiDou navigation satellite system (BDS) has evolved from a regional navigation satellite system (BDS2) to the third generation of the navigation satellite system (BDS3) with global coverage [1]. The BDS2 has provided services mainly in the Asia-Pacific area since 27 December 2012.The completion of BDS3 was officially announced on 31 July 2020, after which BDS3 started to provide services at its full operational capability, including positioning, timing and short message communications [1]. As shown in Table 1, BDS2 currently consists of three medium earth orbit (MEO) satellites, five geostationary earth orbit (GEO) satellites and seven inclined geosynchronous orbit (IGSO) satellites. As shown in Table 2, by June 2020, the BDS3 system consisted of three IGSO satellites, two GEO satellites and 24 MEO satellites [2]. The BDS2 satellites transmit triple frequency signals centered at B1I, B2I and B3I. In BDS3, signals with different modulation techniques and chip rates at a total of six frequencies are transmitted, including the new B1C, B2a, B2b and B2ab signals and the original B1I and B3I signals [3]. To track and analyze multifrequency observations and broadcast navigation-related messages from multi-GNSSs, several publicly accessible global networks of monitoring stations have been set up by various institutions, including the international GNSS service (IGS) multi-GNSS experiment (MGEX) network [4], the Chinese international GNSS monitoring and assessment system (iGMAS) tracking network [5] and the Australian continuously operating reference stations (AUSCORS) [6]. By June 2020, more than 240 receivers worldwide that support BDS2 signal tracking and 160 receivers that support BDS3 signal tracking were accessible. Multifrequency BDS signals tracked by globally distributed monitoring stations offer a basis for the independent monitoring of the global ionospheric total electron content (TEC).
Ionospheric effects need to be eliminated in GNSS data processing for precise navigation and positioning applications [7]. On the other hand, the computation of reliable global ionospheric maps (GIMs) of vertical TEC (VTEC) based on GNSS data provides valuable information in both science and technology fields, such as space weather monitoring and user navigation improvement [8,9]. Thus, the performance of the GIMs computed based on the GNSS data can be regarded as an important index that scales the whole performance of one specific navigation system. With the aim of continually monitoring the ionosphere, GIMs developed by different techniques have been routinely provided by the IGS since 1998 [8]. Although there are currently seven IGS-related ionospheric associate analysis centers (IAACs), the GIMs provided by these IAACs are almost exclusively produced using only GPS/GLONASS data [10]. In addition, although ionospheric modeling based on BDS data has been analyzed in previous studies [11,12], most of these studies used data combinations from BDS and other GNSS systems. That is, due to the limited number of BDS satellites and the uneven distribution of stations available to track BDS3 signals before the completion of the establishment of BDS3, BDS data have not been applied in computations of IGS-related GIM products, and the issue of the performances of GIMs calculated using only BDS observations has not yet been reported.
Considering the effects of differential code bias (DCB) on ionospheric delay modeling and satellite clock correction, satellite and receiver DCBs are thought to be important parameters in the pseudorange-based positioning, navigation and timing (PNT) of GNSSrelated and ionospheric research [13][14][15]. For the proper processing of linear combinations of code observations on multiple frequencies, two independent timing group delay (TGD) parameters, which are defined as the satellite DCBs between the antenna phase center of the B3I to B1I and B2I frequencies, namely TGD1 and TGD2, have been broadcast in navigation-related messages of the BDS [16][17][18][19]. In addition, two types of daily multi-GNSS DCB products have been routinely provided by the German Aerospace Center (DLR) and the Chinese Academy of Sciences (CAS) at the IGS CDDIS website [20]. The DLR-based DCBs are calculated by introducing an a priori ionospheric TEC [21], and the CAS-based DCBs are computed using combined estimations of DCBs and ionospheric parameters [22].
Due to the dispersive nature of ionospheric refractivity, the ionospheric observables are well known to be obtained from the geometry-free linear combination of dual-frequency GNSS measurements [23]. Thus, the accuracy of ionospheric TEC estimation can be improved by selecting two different frequencies with high-quality signals, and in this process, DCBs between different frequencies need to be precisely calibrated. To make full use of the new multifrequency BDS3 signals for positioning and ionospheric research, there is a strong need for precise handling of the BDS3 DCBs that are associated with the newly emerging frequencies. Although the estimation of BDS3 DCBs has been included in a few previous studies [5,[24][25][26][27], further research should be performed, as existing studies have mainly focused on experimental satellites and few BDS3 satellites and a limited number of stations were applied. On one hand, existing studies related to the estimation and analysis of BDS DCBs mainly concentrated on C2I-C7I and C2I-C6I DCBs. The characteristics of the DCBs related to the new BDS3 signals have not yet been systematically studied. On the other hand, proper knowledge of whether there are intersystem biases between BDS3 and BDS2 receiver DCBs is of interest for users who need to jointly process BDS2 and BD3 data [25]. However, different conclusions have been drawn from the existing studies. Li et al. [5] found that an obvious systematic bias existed in the receiver DCB differences between BDS3 and BDS2 when different networks were processed together. Wang et al. [24] found no evident systematic bias between BDS2 and BDS2+BDS3 receiver DCBs.
The main objective of this study is to present a comprehensive evaluation of the characteristics of satellite and receiver DCBs for the new BDS3 signals and the performance of the BDS in modeling global ionospheric TEC. To this end, we carried out experiments and analyses based on BDS2 and BDS3 multifrequency observations collected at 246 IGMAS and MGEX sites during two selected periods in 2020. This paper proceeds as follows. First, the estimation method, which is named APM (Innovation Academy for Precision Measurement Science and Technology), and was used to determine the GIMs and DCBs in the BDS satellites and receivers, is introduced in detail. Then, the characteristics of 13 types of BDS2 and BDS3 DCBs, which have been defined as DCB between the C1X, C1P, C1B, C1A, C5I, C5X, C5Q, C5P, C6I, C7I, C7Z, C7A and C8X signals with respect to the C2I signal, are analyzed in detail. Finally, the performance of the GIMs is validated and evaluated in terms of comparisons with IGS-related final GIMs and standard point positioning (SPP).

Methodology
The APM algorithm used to estimate the BDS-based GIMs and the satellite and receiver DCBs between different frequencies in a time-efficient manner includes three main processing steps. First, the ionospheric observables are extracted from the original BDS pseudorange and carrier-phase observations based on the pseudorange-leveled carrier-phase approach [28,29]. Second, the GIMs are calculated by establishing a global ionospheric VTEC model. Third, the satellite and receiver DCBs between different frequency bands, including those of the new BDS3 signals, are calculated using the GIMs computed in the second step.

Extraction of Ionospheric Observables
Ignoring the effects of multiple paths and noise, the code and carrier-phase measurements in units of length can be expressed as follows [30] where P s r,m (k) and L s r,m (k) are the code and carrier-phase measurements associated with receiver r, satellite s, carrier frequency band m and epoch k, respectively; ρ s r (k) = g s r (k) + dt r (k) − dt s (k) + T r (k) denotes the combination of all frequency-independent terms, including the geometric distance from the satellite to the receiver g s r (k), receiver clock errors dt r (k), satellite clock offset dt s (k) and tropospheric delay T r (k); I s r,1 (k) is the slant ionospheric delay on the first frequency; µ m = f 1 / f m 2 is the frequency-dependent factor, where f denotes the frequency of the carrier-phase; b r,m and d s ,m are the receiver and satellite code hardware delays, respectively; and N s r,m is the real-value ambiguity. Based on Equation (1), the geometry-free code and carrier-phase observation equation can be constructed as follows [31] P s r,mn (k) = P s r,m (k) − P s r,n (k) = (µ m − µ n ) · I s r,1 (k) + DCB r,mn − DCB s mn L s r,mn (k) = L s r,m (k) − L s r,n (k) = −(µ m − µ n ) · I s r,1 (k) + N s r,mn where P s r,mn (k) and L s r,mn (k) denote the geometry-free code and carrier-phase observables between frequency bands m and n, respectively; DCB r,mn = b r,m − b r,n and DCB s mn = b s ,m − b s ,n are the receiver and satellite DCBs, respectively, which are assumed to be constant over periods of 24 h; and N s r,mn denotes the geometry-free carrier-phase ambiguity.
As Equation (2) indicates, the non-dispersive effects ρ s r (k) have been removed and the ionospheric delay parameter I s r,1 (k) changes between epochs, while the remaining parameters are assumed to be constant over time. Thus, the so-called leveling constant, which amounts to the sum of the satellite and receiver DCBs and the ambiguities, can be determined by averaging P s r,mn (k) and L s r,mn (k) for a continuous arc that consists of a total of t epochs, as follows where · t is the operator that computes the average of the variables over t epochs (k = 1 · · · · · ·t). By combining Equations (2) and (3), the ionospheric observableÎ s r,1 (k) can be expressed as follows:Î As such, ionospheric observables, which are linear combinations of the original slant ionospheric delay and the satellite and receiver DCBs, are extracted.

Establishment of a Thin-Layer Global Ionosphere Model
As expressed in Equation (4), the estimable ionospheric observablesÎ s r,1 (k) are coupled with DCBs. To obtain clean slant ionospheric delays, one interrelated task is to separate the satellite and receiver DCBs from the ionospheric observables. It has been demonstrated that the spherical harmonic function model can efficiently model spatial and temporal variability in global ionospheric TEC [8]. In this study, the VTEC is represented by a maximum degree and order equal to 15 spherical harmonic expansions with a 2-h resolution based on the so-called single-layer model [10], which reads as follows x ∑ y=0 P xy (sin(ϕ)) · N xy · A xy cos(yλ) + B xy sin(yλ) (5) where I v (k) is the ionospheric VTEC at the ionospheric intersecting pierce point (IPP); ϕ and λ are the latitude and longitude of the IPP in the sun-fixed geomagnetism coordinate system; P xy (sin(ϕ)) is the classical unnormalized Legendre function; N xy = (x−y)!(2x+1)(2−δ 0y ) (x+y)! is the normalization function; δ is the Kronecker delta; and A xy and B xy are the unknown spherical harmonic coefficients to be estimated. With the reparameterized ionospheric delay I v replacing I s r,1 (i) in Equation (4), we can rewrite the ionospheric observables as follows is the shell mapping function used to convert the slant TEC (STEC) along a given ray path to the VTEC, with E, R e and H shell being the satellite elevation angle, mean radius of Earth and altitude of the ionospheric thin-layer shell (450 km), respectively, and P being the weight of the ionospheric observable, which is dependent on the elevation angle of the satellite.
In the parameter estimation process, piece-wise linear functions are used to guarantee the continuity of the TEC model between two neighboring sessions. Based on Equation (6), the spherical harmonic coefficients can be solved for based on least-squares adjustment, and can be further used to compute the VTEC values for the GIMs at the designated latitudes and longitudes. To improve the computational efficiency, only the ionospheric observables extracted based on the B1I, B2I and B3I signals are selected to estimate the spherical harmonic coefficients. An elevation cutoff angle of 15 • is used to mitigate the impacts of multipath and mapping function errors at low elevations.

Estimation of BDS DCBs between the Different Frequency Bands
By correcting the ionospheric delay based on the GIM calculated with Equation (5), the satellite and receiver DCBs between different frequency bands can be estimated using the following equation where I v,GI M in the second term is the VTEC at the IPP location calculated based on a bilinear interpolation in space and time on previously generated GIM. To eliminate the rank deficiency occurring between the satellite and receiver DCBs, the zero-mean constraint condition assuming that the sum of satellite DCBs of all satellites is equal to zero has been imposed in this study.

Results
In this chapter, we first describe the experimental datasets followed by the validation of the proposed APM method. Thereafter, analyses of the characteristics of satellite and receiver DCBs of the new BDS3 signals are conducted. Then, the performance of the GIMs estimated with only BDS data is analyzed.

Experimental Data
For our analysis, we selected two observation periods in 2020. The first period, covering six months, from day of year (DOY) 061 to 244, is used for validation of the APM method in terms of the accuracy of the APM-based satellite DCB estimates. The second period, covering 1 month, from DOY 214 to 244, is used for the analysis of DCBs for the BDS3 new signals and the BDS-based GIMs.
As shown in Figure 1, 246 globally distributed stations from the IGMAS and MGEX tracking networks were selected to perform the evaluation and analysis. The data sampling rate was chosen as 30 s. Table 3 summarizes the basic information of the tracking stations whose receiver types stayed unchanged during the period of DOY 214-244, 2020. The receiver hardware device information for stations whose receiver type or antenna type changed during the test period is listed in Table 4. Information on the BDS signals, including the signal frequencies, observation type and number of accessible stations, is listed in Table  5. Different colors used to denote stations in Figure 1 correspond to different receiver types in the first column, "Receiver type index", in Table 5. As shown in Table 5, more than 180 receivers can track triple-frequency legacy BDS signals. The C7I signal can be tracked by the majority of receivers except the CETC-54-GMR, BD070 and GNSS_GGR receivers. Considering that the C2I signal can be tracked by all receivers, all DCBs estimated in this study are formed with respect to the C2I signal.   Table 3.    Table 3.  As shown in Figure 2, 48 globally distributed MGEX stations were also selected to investigate the performance of the BDS-based GIM in the position domain. The station coordinates derived from the final IGS solutions were used as a reference to calculate the position errors. As shown in Figure 2, 48 globally distributed MGEX stations were also sel investigate the performance of the BDS-based GIM in the position domain. The coordinates derived from the final IGS solutions were used as a reference to calcu position errors.

Validation of the APM Method for BDS Satellite DCB Estimation
To test the effectiveness of the APM method, the overall performance of th based BDS DCB estimates was first validated, by comparison with daily DCBs p by DLR, CAS and TGD parameters collected from the BDS navigation message. bility, expressed as the standard deviation (STD) of satellite DCB estimates and t sistency among different DCB solutions for the period from DOY 061 to 244, 20 employed as the two performance indicators.
Since the two TGD parameters (TGD1 and TGD2) can be interpreted as C2I-C7I-C6I DCBs, the C2I-C7I DCB is, in theory, equal to the TGD1-TGD2 paramete sidering that the TGD1 parameter for BDS2 and BDS3 might refer to the respectiv ence datum [19], the BDS2 and BDS3 DCBs obtained from different solutions w normalized by a zero-mean constraint across the respective full sets of the BDS2 o satellites before comparison. Additionally, because the DCBs obtained from diffe lutions may be calculated based on different BDS satellites, this section presents o results for the common satellites. For example, since the DCBs for BDS3 GEO s

Validation of the APM Method for BDS Satellite DCB Estimation
To test the effectiveness of the APM method, the overall performance of the APMbased BDS DCB estimates was first validated, by comparison with daily DCBs provided by DLR, CAS and TGD parameters collected from the BDS navigation message. The stability, expressed as the standard deviation (STD) of satellite DCB estimates and the consistency among different DCB solutions for the period from DOY 061 to 244, 2020 were employed as the two performance indicators.
Since the two TGD parameters (TGD1 and TGD2) can be interpreted as C2I-C6I and C7I-C6I DCBs, the C2I-C7I DCB is, in theory, equal to the TGD1-TGD2 parameter. Considering that the TGD1 parameter for BDS2 and BDS3 might refer to the respective reference datum [19], the BDS2 and BDS3 DCBs obtained from different solutions were first normalized by a zero-mean constraint across the respective full sets of the BDS2 or BDS3 satellites before comparison. Additionally, because the DCBs obtained from different solutions may be calculated based on different BDS satellites, this section presents only the results for the common satellites. For example, since the DCBs for BDS3 GEO satellites (PRN C59 and C60) are not provided by the CAS and the DCBs for the PRN C60 satellite are not provided by the DLR, the statistical results do not cover these two satellites.
To understand the stability of APM-based BDS satellite DCB estimates, the time series of the weekly STDs of the C2I-C6I and C2I-C7I DCB estimates during the test period are first calculated. As shown in Figure 3, abnormal variations occurred in the BDS DCB estimates for the PRN C01 and C28 satellites, indicating that the satellite DCBs are not always stable enough. The abnormal weekly mean DCB estimates for the PRN C01 satellite during the periods from DOY 084 to 097 can be attributed to the decreased DCB values of satellite C01 on DOY 091. Obvious jumps in the weekly mean DCB time series shown at the bottom of Figure 3 can be attributed to the DCB variation on DOY 135 for the PRN C28 satellite. Apart from the two obvious jumps in the weekly mean DCB time series shown in Figure 3, DCB estimates for the two satellites appear to be rather stable. The weekly STD values for satellites except PRN C01 and C28 are plotted in Figure 4. The maximum values of the STDs of the DCB estimates obtained from APM, DLR and CAS are 0.17, 0.15 and 0.23 ns, which shows the consistency of the stability of DCB estimates from the APM and other solutions. The mean STDs for the three DCB products over the test period are all less than 0.1 ns, which shows good repeatability of the BDS satellite DCB estimates. weekly STD values for satellites except PRN C01 and C28 are plotted in Figure 4. maximum values of the STDs of the DCB estimates obtained from APM, DLR and are 0.17, 0.15 and 0.23 ns, which shows the consistency of the stability of DCB estim from the APM and other solutions. The mean STDs for the three DCB products ove test period are all less than 0.1 ns, which shows good repeatability of the BDS sat DCB estimates.   weekly STD values for satellites except PRN C01 and C28 are plotted in Figure 4. maximum values of the STDs of the DCB estimates obtained from APM, DLR and are 0.17, 0.15 and 0.23 ns, which shows the consistency of the stability of DCB estim from the APM and other solutions. The mean STDs for the three DCB products ove test period are all less than 0.1 ns, which shows good repeatability of the BDS sat DCB estimates.    Table 6. The statistics show that the TGD-based DCBs represent the largest differences with respect to the other solutions. The averaged RMS values of the differences between the APM-based DCB estimates and the DLR and CAS products are less than 0.23 and 0.32 ns, respectively, showing the encouraging consistency of these three DCB estimation methods. Due to the similar DCB estimation strategy, the APMand DLR-based DCBs have better agreement for both the BDS2 and BDS3 constellations than do the other DCB products. Table 6 shows the DCB estimates obtained from different solutions using BDS2 satellites show better consistency than those using BDS3 satellites. This finding is consistent with the study of Wang et al. [19], in which this phenomenon was attributed to the poor B1I-B3I dual-frequency coverage of the BDS3 satellites.
constellations than do the other DCB products. Table 6 shows the DCB estimate from different solutions using BDS2 satellites show better consistency than th BDS3 satellites. This finding is consistent with the study of Wang et al. [19], in phenomenon was attributed to the poor B1I-B3I dual-frequency coverage of the ellites.  As shown in Figure 5, the differences between the TGD-based DCBs and DCB solutions decreased after DOY 141, 2020, suggesting that an improved str used in the generation of broadcasted TGDs in BDS-based ground control seg addition, the RMS values of the differences of BDS2 DCBs between the TGDother solutions have an obvious jump during the period of DOY 091-103 in 202 the reason for the larger RMS for TGD-based DCBs, Figure 6 presents the RM BDS2 satellites. From Figure 6, it may be concluded that the observed abnormal values were largely driven by abnormal variations that occur in the RMS valu PRN C01 satellite.  Table 6. Statistics of the mean RMSs of the differences between the different satellite-based DCB solutions during the period of DOY 061-244 in 2020, expressed in nanoseconds.

DCB Type
System CAS-DLR As shown in Figure 5, the differences between the TGD-based DCBs and the other DCB solutions decreased after DOY 141, 2020, suggesting that an improved strategy was used in the generation of broadcasted TGDs in BDS-based ground control segments. In addition, the RMS values of the differences of BDS2 DCBs between the TGD-based and other solutions have an obvious jump during the period of DOY 091-103 in 2020. To find the reason for the larger RMS for TGD-based DCBs, Figure 6 presents the RMS for each BDS2 satellites. From Figure 6, it may be concluded that the observed abnormal daily RMS values were largely driven by abnormal variations that occur in the RMS values for the PRN C01 satellite.

CAS-TGD
The mean differences observed among the BDS-based DCBs of individual satellites during the period of DOY 061-244 in 2020 are illustrated in Figure 7, and the RMS values of the differences among varying DCB estimates for different orbit types are summarized in Table 7. As seen from the statistics, the BDS2-based DCB differences among the different products of the IGSO satellites are smaller than those of the GEO and MEO satellites; however, this conclusion is not suitable for BDS3-based DCBs. The mean differences observed among the BDS-based DCBs of individual satellites during the period of DOY 061-244 in 2020 are illustrated in Figure 7, and the RMS values of the differences among varying DCB estimates for different orbit types are summarized in Table 7. As seen from the statistics, the BDS2-based DCB differences among the different products of the IGSO satellites are smaller than those of the GEO and MEO satellites; however, this conclusion is not suitable for BDS3-based DCBs.    The mean differences observed among the BDS-based DCBs of individual satellites during the period of DOY 061-244 in 2020 are illustrated in Figure 7, and the RMS values of the differences among varying DCB estimates for different orbit types are summarized in Table 7. As seen from the statistics, the BDS2-based DCB differences among the different products of the IGSO satellites are smaller than those of the GEO and MEO satellites; however, this conclusion is not suitable for BDS3-based DCBs.

BDS Satellite-Based DCBs among Different Frequency Bands
The consistency between the APM-based BDS satellite DCB estimates and daily DCBs estimated by DLR and CAS has been validated in Section 3.2. To better understand the characteristics of the BDS3 satellite DCBs, mean values of the DCBs among the six different frequency bands covering all operational BDS3 satellites over a period of one month are shown in Figure 8. A close match of BDS3 DCBs between the same two frequency band signals can be noted from this figure. The DCBs between the B1I and B1C signals are confined to a range of ±5 ns, while the DCB values between B1I and other signals of the BDS3 satellites, with the exception of the PRN C33 satellite, are confined to a range of −31 to 55 ns.

BDS Satellite-Based DCBs among Different Frequency Bands
The consistency between the APM-based BDS satellite DCB estimates an DCBs estimated by DLR and CAS has been validated in Section 3.2. To better und the characteristics of the BDS3 satellite DCBs, mean values of the DCBs among different frequency bands covering all operational BDS3 satellites over a period month are shown in Figure 8. A close match of BDS3 DCBs between the same quency band signals can be noted from this figure. The DCBs between the B1I a signals are confined to a range of ±5 ns, while the DCB values between B1I an signals of the BDS3 satellites, with the exception of the PRN C33 satellite, are con a range of −31 to 55 ns. For evaluation of the stability of the BDS3 satellite-based DCBs, Figure 9 sh monthly STDs of the BDS3 satellite-based DCB estimates of the signals among th ent frequency bands. As shown in the figure, other than the C2I-C1A, C2I-C5Q, and C2I-C7A DCB estimates, which exhibit large STDs, the STDs of the DCB estim the signal combinations are less than 0.2 ns. Similar to Figures 8 and 9, the mean and STDs of the DCB estimates of each BDS2 satellite are shown in Figure 10. Th satellite-based DCB estimates are confined to a range of ±20 ns. The monthly STD BDS2 DCB estimates are mainly within 0.2 ns, showing better stability than thos BDS3 satellites. In addition, because the signal quality of B3I is better than that of the stability of the C2I-C6I DCB estimates outperforms that of the C2I-C7I DCB es For evaluation of the stability of the BDS3 satellite-based DCBs, Figure 9 shows the monthly STDs of the BDS3 satellite-based DCB estimates of the signals among the different frequency bands. As shown in the figure, other than the C2I-C1A, C2I-C5Q, C2I-C5I and C2I-C7A DCB estimates, which exhibit large STDs, the STDs of the DCB estimates of the signal combinations are less than 0.2 ns. Similar to Figures 8 and 9, the mean values and STDs of the DCB estimates of each BDS2 satellite are shown in Figure 10. The BDS2 satellite-based DCB estimates are confined to a range of ±20 ns. The monthly STDs of the BDS2 DCB estimates are mainly within 0.2 ns, showing better stability than those of the BDS3 satellites. In addition, because the signal quality of B3I is better than that of B2I [27], the stability of the C2I-C6I DCB estimates outperforms that of the C2I-C7I DCB estimates.  As seen from the statistics on the STDs of the DCB estimates of BDS2 and BDS3 shown in Figures 9 and 10, the IGSO satellites show better stability than the GEO and MEO satellites. This phenomenon has been attributed to the longer observation arcs and wider geographic coverage of IGSO satellites than those of other satellites, as outlined in previous studies [32]. To better understand the impact of the contributed observables on the stability of BDS-based DCB estimates, taking the C2I-C6I DCB as an example, the numbers of tracking stations and contributed observations for each BDS satellite are shown in Figure 11. We can see that the number of stations available for the computation of daily GEO satellites is less than that of IGSO satellites. In addition, the number of contributing  As seen from the statistics on the STDs of the DCB estimates of BDS2 and BDS3 shown in Figures 9 and 10, the IGSO satellites show better stability than the GEO and MEO satellites. This phenomenon has been attributed to the longer observation arcs and wider geographic coverage of IGSO satellites than those of other satellites, as outlined in previous studies [32]. To better understand the impact of the contributed observables on the stability of BDS-based DCB estimates, taking the C2I-C6I DCB as an example, the numbers of tracking stations and contributed observations for each BDS satellite are shown in Figure 11. We can see that the number of stations available for the computation of daily GEO satellites is less than that of IGSO satellites. In addition, the number of contributing stations for the DCB estimation of satellites launched after 2019 is usually fewer than 50, As seen from the statistics on the STDs of the DCB estimates of BDS2 and BDS3 shown in Figures 9 and 10, the IGSO satellites show better stability than the GEO and MEO satellites. This phenomenon has been attributed to the longer observation arcs and wider geographic coverage of IGSO satellites than those of other satellites, as outlined in previous studies [32]. To better understand the impact of the contributed observables on the stability of BDS-based DCB estimates, taking the C2I-C6I DCB as an example, the numbers of tracking stations and contributed observations for each BDS satellite are shown in Figure 11. We can see that the number of stations available for the computation of daily GEO satellites is less than that of IGSO satellites. In addition, the number of contributing stations for the DCB estimation of satellites launched after 2019 is usually fewer than 50, and this can be regarded as one reason for the degraded stability of these newly launched satellites, compared to the remaining satellites. Although the B2ab signal has the best pseudorange measurement accuracy and the strongest anti-multipath ability among the different BDS3 signals [3], as shown in Figure 12, the STDs of the C2I-C8X DCBs are still larger than those of the C2I-C6I and C2I-C7Z DCBs. It may be concluded that the contributed number of observables has a larger impact on the stability of DCB estimates than the signal quality.
emote Sens. 2021, 13, x FOR PEER REVIEW 1 and this can be regarded as one reason for the degraded stability of these newly laun satellites, compared to the remaining satellites. Although the B2ab signal has the pseudorange measurement accuracy and the strongest anti-multipath ability amon different BDS3 signals [3], as shown in Figure 12, the STDs of the C2I-C8X DCBs ar larger than those of the C2I-C6I and C2I-C7Z DCBs. It may be concluded that the con uted number of observables has a larger impact on the stability of DCB estimates tha signal quality.

BDS Receiver-Based DCBs among Different Frequency Bands
To check whether systematic bias existed in the receiver DCBs between BDS3 emote Sens. 2021, 13, x FOR PEER REVIEW 14 and this can be regarded as one reason for the degraded stability of these newly laun satellites, compared to the remaining satellites. Although the B2ab signal has the pseudorange measurement accuracy and the strongest anti-multipath ability among different BDS3 signals [3], as shown in Figure 12, the STDs of the C2I-C8X DCBs are larger than those of the C2I-C6I and C2I-C7Z DCBs. It may be concluded that the con uted number of observables has a larger impact on the stability of DCB estimates tha signal quality.

BDS Receiver-Based DCBs among Different Frequency Bands
To check whether systematic bias existed in the receiver DCBs between BDS3

BDS Receiver-Based DCBs among Different Frequency Bands
To check whether systematic bias existed in the receiver DCBs between BDS3 and BDS2, we designed two sets of experiments with all accessible stations from the MGEX and IGMAS networks. In the first set of experiments, both BDS2 and BDS3 satellites are used to solve for the receiver DCBs. In the second set of experiments, only BDS2 satellites are used. Since differences exist in the numbers of contributing satellites, the zero-mean condition is also different between the two solutions. Considering that the zero-mean condition imposed on all available satellites is used in the two different solutions, to unify the BDS2+BDS3-based and BDS2-based DCB solutions, the DCB estimates are only normalized by a zero-mean condition for the commonly used BDS2 satellites. Figure 13 shows the mean differences and STDs of the receiver C2I-C6I DCB estimates between BDS2-only and BDS2+BDS3 solutions during the period of DOY 214-244, 2020. There is no change in the receiver hardware device during the test period for receivers shown in this figure. The statistics shown in Table 8 reveal that the differences in receiver DCBs between the two different solutions have a strong dependence on the receiver type. This conclusion is in accordance with the study of Wang et al. [24]. Except for the SEPT and TRIMBLE receivers, the differences in receiver DCBs between the two different solutions are almost all positive. The BDS2 and BDS3 receiver DCB differences for SEPT receivers are all negative, while those for TRIMBLE receivers fluctuate around an average value at the level of 0.07 ns. The mean bias and STD of the estimated DCB differences between the BDS2-only and BDS2+BDS3 solutions for all receiver types are 0.16 ns and 0.12 ns, respectively.
Remote Sens. 2021, 13, x FOR PEER REVIEW condition is also different between the two solutions. Considering that the zero-me dition imposed on all available satellites is used in the two different solutions, to u BDS2+BDS3-based and BDS2-based DCB solutions, the DCB estimates are only ized by a zero-mean condition for the commonly used BDS2 satellites. Figure 13 shows the mean differences and STDs of the receiver C2I-C6I D mates between BDS2-only and BDS2+BDS3 solutions during the period of DOY 2 2020. There is no change in the receiver hardware device during the test period for ers shown in this figure. The statistics shown in Table 8 reveal that the differenc ceiver DCBs between the two different solutions have a strong dependence on the type. This conclusion is in accordance with the study of Wang et al. [24]. Except SEPT and TRIMBLE receivers, the differences in receiver DCBs between the two d solutions are almost all positive. The BDS2 and BDS3 receiver DCB differences fo receivers are all negative, while those for TRIMBLE receivers fluctuate around an value at the level of 0.07 ns. The mean bias and STD of the estimated DCB diff between the BDS2-only and BDS2+BDS3 solutions for all receiver types are 0.16 0.12 ns, respectively.  To further analyze the impact of the antenna type and firmware version on between the BDS2 and BDS3 receiver DCB estimates, taking the receiver type of TR NETR9 as an example, Figure 14 shows the mean differences and STDs of the  To further analyze the impact of the antenna type and firmware version on the bias between the BDS2 and BDS3 receiver DCB estimates, taking the receiver type of TRIMBLE NETR9 as an example, Figure 14 shows the mean differences and STDs of the receiver DCB estimates between BDS2-only and BDS2+BDS3 solutions for the selected stations in three groups equipped with the same types of receivers and antennas and firmware version configuration. Table 9 shows the information of the three groups of stations and the corresponding statistical results. There are obvious significant differences in the bias results between receivers from different manufacturers. In addition, the mean differences of the receiver DCB estimates between BDS2-only and BDS2+BDS3 solutions for receivers in the second and third group are all positive values, and those of the first group are almost all negative values. It may be concluded that the differences in receiver DCB estimates between BDS2 and BDS3 are related to firmware version and types of receivers and antennas, and the firmware version has a larger impact on the bias between the BDS2 and BDS3 receiver DCBs than that of the antenna type.
mote Sens. 2021, 13, x FOR PEER REVIEW the second and third group are all positive values, and those of the first grou all negative values. It may be concluded that the differences in receiver D between BDS2 and BDS3 are related to firmware version and types of rece tennas, and the firmware version has a larger impact on the bias between t BDS3 receiver DCBs than that of the antenna type.  To further research the impact of added BDS3 observations on the stabili DCB estimates, we plotted the STD of the C2I-C6I receiver DCB estimates fo lutions in Figure 15. As shown in the figure, the mean STDs of DCB estima BDS2+BDS3 and BDS2-only solutions are 0.31 and 0.37 ns, respectively, proved stability induced by the added BDS3 observations.  To further research the impact of added BDS3 observations on the stability of receiver DCB estimates, we plotted the STD of the C2I-C6I receiver DCB estimates for the two solutions in Figure 15. As shown in the figure, the mean STDs of DCB estimates based on BDS2+BDS3 and BDS2-only solutions are 0.31 and 0.37 ns, respectively, indicting improved stability induced by the added BDS3 observations.
To further research the impact of added BDS3 observations on the stability of DCB estimates, we plotted the STD of the C2I-C6I receiver DCB estimates for the lutions in Figure 15. As shown in the figure, the mean STDs of DCB estimates b BDS2+BDS3 and BDS2-only solutions are 0.31 and 0.37 ns, respectively, indict proved stability induced by the added BDS3 observations.  To investigate the stability of the BDS3 receiver DCB estimates among the different frequency bands, the monthly STD of each individual receiver DCB estimates is shown in Figure 16. This figure shows that compared to the satellite-based STD values, the receiver DCBs show higher fluctuations. The STD values of receiver DCB estimates among the different frequency bands for all receivers located at low latitudes are larger than those of receivers located at other latitudes, indicating the latitudinal dependence of BDS3 receiver DCBs. The STDs of receiver DCBs between the B2I and B1C signals and other frequency bands are smaller than 0.4 ns and 0.6 ns, respectively, for most stations. To investigate the stability of the BDS3 receiver DCB estimates among the different frequency bands, the monthly STD of each individual receiver DCB estimates is shown in Figure 16. This figure shows that compared to the satellite-based STD values, the receiver DCBs show higher fluctuations. The STD values of receiver DCB estimates among the different frequency bands for all receivers located at low latitudes are larger than those of receivers located at other latitudes, indicating the latitudinal dependence of BDS3 receiver DCBs. The STDs of receiver DCBs between the B2I and B1C signals and other frequency bands are smaller than 0.4 ns and 0.6 ns, respectively, for most stations. To study the impact of changes in types of receivers and antennas on receiver BDS DCB estimates, taking METG, SPT0 and YAR2 stations as examples, Figure 17 shows the time series of the receiver DCB estimates for the three stations during the test period. There are several gaps, due to discontinuous data or changes in the receiver hardware device. According to the information on the three receivers shown in Table 4, the receiver To study the impact of changes in types of receivers and antennas on receiver BDS DCB estimates, taking METG, SPT0 and YAR2 stations as examples, Figure 17 shows the time series of the receiver DCB estimates for the three stations during the test period. There are several gaps, due to discontinuous data or changes in the receiver hardware device. According to the information on the three receivers shown in Table 4, the receiver type for stations METG and YAR2 and the antenna type for station SPT0 changed during the period of DOY 218-227. The abnormal variations in receiver DCB estimates for the three stations are consistent with the changes in receiver hardware devices, from which we can conclude that the jumps are largely driven by changes in receiver type or antenna type. This justifies the impacts of changes in types of receivers and antennas on the variations in receiver DCB estimates.

Performance of GIMs Estimated Using Only BDS Observations
To study the performance of the BDS in modeling the global ionospheric TEC, the GIMs that were produced based on BDS-only-obtained observational data were compared with the combined final GIM products of the IGS. In addition, GIMs based only on GPS data were also calculated, and used as a reference in the performance assessment of the BDS-based GIMs.
The bias and RMS series of the APM-based GIMs with regard to the combined final GIM products of the IGS during the period from DOY 214 to 244 in 2020 are shown in Figure 18. The figure depicts that the GPS-based GIM products provide better consistency with the IGS GIMs than do the BDS-based products. Taking the IGS-based GIM products as a reference, the accuracies of the BDS-based and GPS-based GIMs are 1.86 and 1.39 total electron content units (TECUs), respectively. Considering that the IGS-provided GIMs are calculated based on GPS and GLONASS data, the GPS-based GIMs are expected to be more consistent with the IGS GIM products than the BDS-based GIMs.

Performance of GIMs Estimated Using Only BDS Observations
To study the performance of the BDS in modeling the global ionospheric TEC, the GIMs that were produced based on BDS-only-obtained observational data were compared with the combined final GIM products of the IGS. In addition, GIMs based only on GPS data were also calculated, and used as a reference in the performance assessment of the BDS-based GIMs.
The bias and RMS series of the APM-based GIMs with regard to the combined final GIM products of the IGS during the period from DOY 214 to 244 in 2020 are shown in Figure 18. The figure depicts that the GPS-based GIM products provide better consistency with the IGS GIMs than do the BDS-based products. Taking the IGS-based GIM products as a reference, the accuracies of the BDS-based and GPS-based GIMs are 1.86 and 1.39 total electron content units (TECUs), respectively. Considering that the IGS-provided GIMs are calculated based on GPS and GLONASS data, the GPS-based GIMs are expected to be more consistent with the IGS GIM products than the BDS-based GIMs.
with the IGS GIMs than do the BDS-based products. Taking the IGS-ba as a reference, the accuracies of the BDS-based and GPS-based GIMs are electron content units (TECUs), respectively. Considering that the IGScalculated based on GPS and GLONASS data, the GPS-based GIMs more consistent with the IGS GIM products than the BDS-based GIMs  The even distribution of ionospheric observables is recognized to play an important role in high-precision GIM product calculations. To study the impact of IPP distributions on the performance differences between the BDS-based and GPS-based GIMs, the IPP distributions of the BDS-and GPS-based ionospheric observables for the first 15 min of DOY 214, 2020 are shown in Figure 19. The density of the IPPs of the BDS is significantly increased when BDS3 satellites are included. Although the distributions of the IPPs tracked by both GPS and BDS satellites can cover most of the continuous regions, the IPPs for the GPS are denser than those of the BDS in the oceans. Thus, the poorer IPP distribution of BDS can be regarded as one of the reasons for the degraded precision of the BDS-based GIMs.
To investigate the performance of the BDS-based GIM in the position domain, the user positioning performance regarding single-frequency SPP has been analyzed using BDS observations from all accessible stations listed in Figure 2 for DOY 214-244, 2020. In the SPP process, we used only pseudorange observations and broadcasted navigation messages. The results of the positioning solutions without ionospheric error corrections were also analyzed. Moreover, the IGS-provided GIMs and the GIMs calculated using the APM method with only GPS data were also applied for a more comprehensive comparison. The RMS position errors of the SPP in the horizontal and vertical coordinate components, which were obtained for each 0.5-h bin over the test period with respect to local time, are illustrated in Figure 20. According to the mean values of the RMS errors summarized in the legend of this figure, the maximum differences in the horizontal and vertical position errors between the BDS-based GIMs and other GIMs products are less than 0.03 and 0.06 m, respectively. Compared to the results obtained without ionospheric error correction, the ionospheric correction obtained from the GIM products mainly reduces the SPP positioning bias in the Up component, showing a significantly better improved positioning accuracy in the vertical direction than in the horizontal direction.
tributions of the BDS-and GPS-based ionospheric observables for the first 15 min of DOY 214, 2020 are shown in Figure 19. The density of the IPPs of the BDS is significantly increased when BDS3 satellites are included. Although the distributions of the IPPs tracked by both GPS and BDS satellites can cover most of the continuous regions, the IPPs for the GPS are denser than those of the BDS in the oceans. Thus, the poorer IPP distribution of BDS can be regarded as one of the reasons for the degraded precision of the BDS-based GIMs. To investigate the performance of the BDS-based GIM in the position domain, the user positioning performance regarding single-frequency SPP has been analyzed using BDS observations from all accessible stations listed in Figure 2 for DOY 214-244, 2020. In the SPP process, we used only pseudorange observations and broadcasted navigation messages. The results of the positioning solutions without ionospheric error corrections were also analyzed. Moreover, the IGS-provided GIMs and the GIMs calculated using the APM method with only GPS data were also applied for a more comprehensive comparison. The RMS position errors of the SPP in the horizontal and vertical coordinate components, which were obtained for each 0.5-h bin over the test period with respect to local time, are illustrated in Figure 20. According to the mean values of the RMS errors summarized in the legend of this figure, the maximum differences in the horizontal and vertical position errors between the BDS-based GIMs and other GIMs products are less than 0.03 and 0.06 m, respectively. Compared to the results obtained without ionospheric error correction, the ionospheric correction obtained from the GIM products mainly reduces the

Discussion
In this contribution, the characteristics of 13 types of BDS DCBs and the accurac BDS-based GIM are analyzed for the first time. The assessment proves that the contribu number of stations and observables has a large impact on both the stability of sate DCB estimates and the performance of the GIMs. Thus, more receivers available to tr BDS data should be developed to achieve more stable DCB estimates and better model variation in the global ionospheric TEC.

Discussion
In this contribution, the characteristics of 13 types of BDS DCBs and the accuracy of BDS-based GIM are analyzed for the first time. The assessment proves that the contributed number of stations and observables has a large impact on both the stability of satellite DCB estimates and the performance of the GIMs. Thus, more receivers available to track BDS data should be developed to achieve more stable DCB estimates and better model the variation in the global ionospheric TEC.
The issue of whether there are code intersystem biases between BDS3 and BDS2 receiver DCBs is preliminarily analyzed in this study. Consistent with the study of Wang et al. [24], the difference between BDS2 and BDS2+BDS3 receiver DCBs is found to be related to the receiver type. Additionally, this study further found that the systematic bias is also related to the antenna type and firmware version. Although the average bias between the BDS2 and BDS3 receiver DCBs is close to 0, we still cannot draw the conclusion that there is no evident systematic bias between the BDS2 and BDS2+BDS3 receiver DCBs since the mean difference for JAVAD receivers reaches 1.03 ns. Therefore, the receiver DCBs of BDS3 and BDS2 should be separately estimated when BDS data collected by some specific types of receivers such as JAVAD receivers are processed.
The test in the SPP experiment shows that there are no significant differences between the GPS-and BDS-based GIMs in mitigating the ionospheric delay effects on position errors. Notably, when estimating the accuracy of ionospheric correction products, observation data during both high-and low-activity periods should be tested. However, BDS3 had not completed global constellation deployment until June 2020, which corresponded to low solar activity. Therefore, the BDS-based GIM validation in this study is based only on data collected during a low-solar activity period, in which the impact of the ionospheric delay was less than that during high-activity periods. More experiments based on data collected under high solar activity should be performed in the future.

Conclusions
This article provides an analysis of satellite and receiver DCBs among different frequency bands for BDS2 and BDS3, and an overview of the performance assessment of GIMs produced only by BDS observables for the first time. The GIMs and DCBs are estimated by the proposed APM method based on data obtained from the IGMAS and MGEX networks. In the proposed APM method, to improve the computational efficiency, GIMs based on BDS observations at the B1I, B2I and B3I signals are first estimated, and then the satellite and receiver DCBs among the other frequency bands are estimated by removing the ionospheric delay using the previously generated GIMs. The conclusions are summarized as follows: 1.
Good agreement of the APM-based BDS satellite DCB estimates with those of CAS and DLR is found, with RMS values of 0.26 and 0.18 ns, respectively. The mean STDs of the APM, DLR and CAS satellite DCB estimates are less than 0.1 ns. Larger differences between different DCB products are found for BDS3 than BDS2 satellites. 2.
For the BDS3 DCBs of the new signals, a close match of satellite DCB estimates between the same two frequency band signals is found.

3.
The BDS satellite DCB estimates among different frequency bands for IGSO satellites show better stability than those from MEO and GEO satellites for both BDS2 and BDS3. Additionally, the stability of BDS satellite DCB estimates is found to be related to the number of contributed stations and available observables applied in the process of DCB estimation.

4.
Variations in the types of receivers and antennas have impact on receiver DCB estimates. The systematic bias between BDS2 and BDS3 receiver DCBs is found to be related to the receiver type, antenna type and firmware version. 5.
The accuracies of GIMs based on BDS-only data and GPS-only data are 1.86 and 1.39 TECU, when taking the IGS-provided GIMs as references, and there is no significant difference in the positioning accuracy when the three GIM products are applied in SPP experiments during low solar activity periods.
Notably, the conclusion on the characteristics of satellite and receiver DCBs for the new BDS3 signals and the performance of GIMs based on only BDS data is drawn from the analysis based on data spanning only a month in 2020. More research should be performed with longer-term observations of BDS3 network satellites, and more evenly distributed BDS global stations in the future.