Epoch-by-epoch estimation and analysis of BeiDou Navigation Satellite System (BDS) receiver differential code biases with the additional BDS-3 observations

The differential code bias (DCB) of the Global Navigation Satellite System (GNSS) is an important error source in ionospheric modeling, which was generally estimated as constants every day. However, the receiver DCB may be changing due to the varying spatial environments and temperatures. In this paper, a method based on the global ionospheric map (GIM) of the Center for Orbit Determination in Europe (CODE) is presented to estimate the BeiDou Navigation Satellite System (BDS) receiver DCB with epoch-by-epoch estimates. The BDS receiver DCBs are analyzed from 30 d of Multi-GNSS Experiment observations. The comparison of estimated receiver DCB of BDS with the DCB provided by the German Aerospace Center (DLR) and the Chinese Academy of Sciences (CAS) shows a good agreement. The root-mean-square (rms) values of receiver DCB are 0.43 and 0.80 ns with respect to the DLR and CAS estimates, respectively. In terms of the intraday variability of receiver DCB, most of the receiver DCBs show relative stability within 1 d with the intraday standard deviation (SD) of less than 1 ns. However, larger fluctuations with more than 2 ns of intraday receiver DCB are found. Besides, the intraday stability of receiver DCB calculated by the third-generation BDS (BDS-3) and the second-generation BDS (BDS-2) observations is compared. The result shows that the intraday stability of BDS-3 receiver DCB is better than that of BDS-2 receiver DCB.

1 Introduction

The Global Positioning System (GPS) has been widely used as a useful tool in ionospheric monitoring and modeling (Hernández-Pajares et al., 1999; McCaffrey et al., 2017; Choi and Lee, 2018). With its rapid development in recent years, the BeiDou Navigation Satellite System (BDS) begins to play an important role in global navigation, positioning, timing, and related applications (Su and Jin, 2019; Wang et al., 2019). Since the third-generation BeiDou Navigation Satellite System (BDS-3) provided global services at the end of 2018, more and more stations have been constructed to track the BDS signals. In addition, BDS can also provide available observations globally for ionospheric modeling and research. The ionospheric total electron contents (TECs) can be obtained by the geometry-free linear combination of the dual-frequency Global Navigation Satellite System (GNSS) code observations, in which the differential code bias (DCB) is one of the main error terms to be estimated and removed (Sanz et al., 2017; Jin and Su, 2020).

In general, DCB is considered to be an unknown parameter together with ionospheric parameters in global or regional ionospheric modeling (Jin et al., 2012). However, only a limited number of tracking stations could be used for the earlier regional system (BDS-2) because of its regional coverage (Montenbruck et al., 2013). In a method called IGGDCB (where IGG stands for the Institute of Geodesy and Geophysics, Chinese Academy of Sciences), which was proposed to estimate the BDS DCB, the single-station ionosphere modeling method is used to overcome a limited number of global tracking stations for BDS (Li et al., 2012). Furthermore, in order to make full use of the useful information of the raw measurement, a method based on triple-frequency uncombined precise point positioning (UPPP) is used to estimate the BDS DCB (Shi et al., 2015; Fan et al., 2017; Liu et al., 2019). All the above methods need to estimate the ionospheric TEC simultaneously for DCB estimation. Therefore, an existing global ionospheric map (GIM) is used to eliminate ionospheric TEC for improving the BDS DCB estimation efficiency (Jin et al., 2016; Li et al., 2019; Wang et al., 2019). At present, two International GNSS Service (IGS) analysis centers can provide DCB products of BDS. One is the German Aerospace Center (DLR), in which the DCB products are generated by the GIM based on the Multi-GNSS Experiment (MGEX) observational data (Montenbruck et al., 2014). The other IGS analysis center is the Chinese Academy of Sciences (CAS), where the DCB products are estimated by the IGG method based on the multi-GNSS observations (Wang et al., 2016). The BDS DCB products provided by DLR and CAS can be obtained at ftp://cddis.gsfc.nasa.gov/pub/gps/products/mgex/dcb/ (last access: 20 May 2020). It should be noted that these MGEX stations are selectively used in the DCB estimation by DLR and CAS, and the stations used by DLR and CAS are different. In other words, only the receiver DCB of stations used in DLR and CAS can be provided in their DCB products.

In the DCB estimation, the satellite and receiver DCBs of BDS are generally estimated as constants every day (Li et al., 2014; Xue et al., 2016). However, the receiver DCB may vary within 1 d due to varying spatial environments and temperatures. The estimation of receiver DCB as constant every day may cause errors in ionospheric modeling if the receiver DCB has significant intraday fluctuations. It would have been better to analyze the intraday variation of receiver DCB before the estimation of receiver DCB as constant over a day. As can be seen from previous studies, the stability of BDS receiver DCB is not as good as that of the satellite (Jin et al., 2016; Xue et al., 2016). Thus, a study on estimation and analysis of short-term variations of multi-GNSS DCB was carried out by Li et al. (2018), who used the GIM by IGS and the satellite DCB by DLR to estimate the receiver DCB with an hourly resolution. As shown by the results, an obvious short-term variation was found in the receiver DCB within 1 d. However, the receiver DCB may exhibit variations within 1 h. Furthermore, more global available observations from BDS-3 can be used in DCB estimation instead of only using BDS-2 as in previous studies. Since 2013, the BDS measurements have been collected by the MGEX network (Montenbruck et al., 2017). With the rapid development of MGEX, the BDS-3 measurements can be tracked by more MGEX stations.

In order to better analyze the short-term variations of BDS receiver DCB with the additional BDS-3 observations, this paper presents a DCB estimation method based on the GIM by IGS to estimate the BDS satellite and receiver DCB. The receiver DCB that is treated as the changing parameter within 1 d is estimated epoch by epoch. The intraday stability of receiver DCBs are analyzed through BDS-2+BDS-3 observations. Besides, the intraday stability of the receiver DCB obtained by BDS-2-only and BDS-3-only observations are compared. The rest of the paper is organized as follows. In Sect. 2, the DCB estimation method and the experiment are introduced. In Sect. 3, the experimental results and corresponding analysis are presented. Finally, conclusions are given in Sect. 4.

2 Method and data

In this section, the method used to estimate receiver DCB (epoch by epoch) is firstly introduced. The corresponding equations for the DCB estimation method are described in detail. Then, the experimental data used in this study are described.

2.1 DCB estimation method

As we all know, the ionospheric observation equation can be obtained by the geometry-free linear combination of dual-frequency GNSS code observations. The code observations smoothed by carrier phase can be expressed as follows (Wellenhof et al., 1992; Hernández-Pajares et al., 2009, 2017):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{P}}_{\mathrm{4},\mathrm{r}}^{\mathrm{s}}\left(t\right)& ={\stackrel{\mathrm{̃}}{P}}_{\mathrm{2},\mathrm{r}}^{\mathrm{s}}\left(t\right)-{\stackrel{\mathrm{̃}}{P}}_{\mathrm{1},\mathrm{r}}^{\mathrm{s}}\left(t\right)=\mathrm{40.28}\left({\displaystyle \frac{\mathrm{1}}{f_{\mathrm{1}}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{f_{\mathrm{2}}^{\mathrm{2}}}}\right)\\ & \cdot {\mathrm{STEC}}_{\mathrm{r}}^{\mathrm{s}}\left(t\right)+c\cdot \left({\mathrm{DCB}}_{\mathrm{r}}+{\mathrm{DCB}}^{\mathrm{s}}\right),\end{array}\end{array}$$

where “s”, “r”, and “t” denote the satellite, receiver, and epoch, respectively; f₁ and f₂ are the first and second frequencies of the observations, respectively; ${\stackrel{\mathrm{̃}}{P}}_{\mathrm{1},\mathrm{r}}^{\mathrm{s}}$ and ${\stackrel{\mathrm{̃}}{P}}_{\mathrm{2},\mathrm{r}}^{\mathrm{s}}$ are smoothed code observations at frequencies f₁ and f₂, respectively; ${\stackrel{\mathrm{̃}}{P}}_{\mathrm{4},\mathrm{r}}^{\mathrm{s}}$ is the smoothed geometry-free linear combination of code observations; c is the speed of light; DCB_r and DCB^s are the receiver and satellite DCB, respectively; and STEC${}_{\mathrm{r}}^{\mathrm{s}}$ is the slant total electron content along the signal-propagation path between GNSS satellite “s” and ground receiver “r”. The slant total electron content can be converted to the vertical TEC (VTEC) by a mapping function, which can be defined as follows (Schaer, 1999):

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}\mathrm{STEC}=\mathrm{VTEC}\cdot M\left(z\right);\\ M\left(z\right)=\frac{\mathrm{1}}{\cos z^{\prime }},\sin z^{\prime }=\frac{R}{R+H}\sin z;\end{array}\right.\end{array}$$

where R (6371 km) is the mean radius of the Earth; H (450 km) is the height of the assumed single-layer ionosphere; and z and z^′ denote the satellite zenith angles of a satellite at the receiver and the corresponding ionospheric pierce point (IPP), respectively.

According to Eq. (1), the satellite and receiver DCBs can be obtained by estimation together with the ionospheric parameters in ionospheric modeling or by elimination of the STEC with the existing GIM. The satellite and receiver DCBs are generally estimated as constants every day. In order to consider the intraday variability of the receiver DCB, the receiver DCB that is treated as the changing parameter within 1 d is estimated epoch by epoch in this paper. The satellite DCB is still estimated as a constant over 1 d, which can be defined as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathrm{DCB}}_{\mathrm{r}}\left(t\right)& +{\mathrm{DCB}}^{\mathrm{s}}=\left({\stackrel{\mathrm{̃}}{P}}_{\mathrm{4},{\mathrm{r}}^{\mathrm{s}}}\right(t)-\mathrm{40.28}\left({\displaystyle \frac{\mathrm{1}}{f_{\mathrm{1}}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{f_{\mathrm{2}}^{\mathrm{2}}}}\right)\\ & \cdot {\mathrm{VTEC}}_{\mathrm{r}}^{\mathrm{s}}\left(t\right)\cdot M\left(z\right))/c,\end{array}\end{array}$$

where VTEC values can be eliminated with the GIM provided by IGS. Therefore, Eq. (1) can be further described as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{r}\left[\begin{array}{cccc}{\mathbf{IR}}_{\mathrm{1}}& \mathrm{\cdots }& \mathrm{\cdots }& {\mathbf{IS}}_{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{\cdots }& \mathrm{\cdots }& {\mathbf{IR}}_{\mathrm{r}}& {\mathbf{IS}}_{\mathrm{r}}\end{array}\right]\times \left[\begin{array}{c}{\mathrm{DCB}}_{\mathrm{1}}\left(t\mathrm{1}\right)\\ \mathrm{⋮}\\ {\mathrm{DCB}}_{\mathrm{r}}\left(t_{\mathrm{r}}\right)\\ {\mathrm{DCB}}^{\mathrm{1}}\\ \mathrm{⋮}\\ {\mathrm{DCB}}^{\mathrm{s}}\end{array}\right]\\ =\left[\begin{array}{c}\left({\stackrel{\mathrm{̃}}{P}}_{\mathrm{4},\mathrm{r}}^{\mathrm{s}}\right(\mathrm{1})-k\cdot {\mathrm{VTEC}}_{\mathrm{r}}^{\mathrm{s}}(\mathrm{1}\left)\right)/c\\ \mathrm{⋮}\\ \left({\stackrel{\mathrm{̃}}{P}}_{\mathrm{4},\mathrm{r}}^{\mathrm{s}}\right(m)-k\cdot {\mathrm{VTEC}}_{\mathrm{r}}^{\mathrm{s}}(m\left)\right)/c\end{array}\right]\end{array},\end{array}$$

where IR_r and IS_r are the coefficient matrices of receiver DCB and satellite DCB for the rth station, respectively, which consist of 0 and 1. The sizes of IR_r and IS_r are m_r×t_r and m_r×s, respectively; m_r is the number of available observations for the rth station; $m=m_{\mathrm{1}}+\mathrm{\cdots }{m}_{\mathrm{r}}$ is the total of available observations for all stations; s is the total number of all available satellites; t_r is the number of available epoch; DCB_r(t_r) is the vector of 1×t_r; DCB^s is the DCB of the sth satellite; and $k=\mathrm{40.28}(\frac{\mathrm{1}}{f_{\mathrm{1}}^{\mathrm{2}}}-\frac{\mathrm{1}}{f_{\mathrm{2}}^{\mathrm{2}}})\cdot M\left(z\right)$. As can be seen from Eq. (4), the numbers of parameters to be estimated are $t_{\mathrm{1}}+\mathrm{\cdots }{t}_{\mathrm{r}}+s$, m_r>t_r, and $m>t_{\mathrm{1}}+\mathrm{\cdots }{t}_{\mathrm{r}}+s$. The number of observations is much larger than the number of parameters to be estimated. Therefore, the receiver and satellite DCB values can be obtained based on the least square method. Notably, the zero-mean condition should be added to the satellite DCB in least square adjustment due to the correlation between satellite and receiver DCB, which can be defined as below:

$$\begin{array}{}\text{(5)}& \sum _{i=\mathrm{1}}^s{\mathrm{DCB}}^i=\mathrm{0}.\end{array}$$

In order to reduce the influence of multipath noise and mapping function error, we set a 20^∘ cutoff elevation angle. In addition, the weight of observations is used in DCB estimation, which can be expressed as follows:

$$\begin{array}{}\text{(6)}& P={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+{\cos }^{\mathrm{2}}E)}},\end{array}$$

where P is the weight of DCB observations and E is elevation angle of the satellite.

Based on Eqs. (4)–(6), the DCB values of the epoch-by-epoch receiver and daily satellite can be obtained. In this study, the BDS receiver and satellite DCBs with the additional BDS-3 observations are estimated by the method mentioned above. The intraday stability of receiver DCBs are analyzed by using BDS-2+BDS-3 observations. Besides, the intraday stability of receiver DCB obtained by using BDS-2-only and BDS-3-only observations are compared. To implement the DCB estimation, the MATLAB processing program is used, which is developed based on the M_DCB software (Jin et al., 2012).

2.2 Experimental data

In order to verify the performance of the DCB estimation method and analyze the intraday stability of BDS receiver with the additional BDS-3 observations, we use 30 d (January 2019) of Multi-GNSS Experiment observations, which include 109 stations. The MGEX observations are downloaded from ftp://cddis.gsfc.nasa.gov/pub/gps/data/ (last access: 20 May 2020), and the GIM is used from ftp://cddis.gsfc.nasa.gov/pub/gps/products/ionex/(last access: 20 May 2020). Besides, the precise satellite ephemeris is provided by Wuhan University, which is available at ftp://cddis.gsfc.nasa.gov/pub/gps/products/mgex/ (last access: 20 May 2020). It should be noted that the DCB estimation and analysis in this study are mainly for C2I–C6I DCB since the C2I and C6I code observations from BDS-2 and BDS-3 can be tracked simultaneously by more stations. Figure 1 presents the distribution of selected MGEX stations in January 2019, in which different types of receivers are distinguished by colors. The corresponding information of the receiver is listed in Table 1. As can be seen from Fig. 1 and Table 1, the Trimble receiver is used in most stations. In addition, the average number of available observation epochs for the selected MGEX stations in January 2019 is counted, as shown in Fig. 2. It is clear that greater quantities of available observations are shown in the stations located in the Asia Pacific region.

Figure 1 Distribution of selected MGEX stations in January 2019. Six colors represent six different receiver types, which are listed in Table 1.

Figure 2 Distribution of average number of available observation epochs for selected MGEX stations in January 2019.

Table 1 Summary of selected MGEX stations in January 2019.

3 Results and analysis

The section begins by verifying the performance of the DCB estimation method, in which the satellite and receiver DCBs are compared with the DCB products by CAS and DLR. Then, the intraday stability of the receiver DCB is analyzed. The intraday stability of receiver DCB obtained by using BDS-2-only and BDS-3-only observations are finally compared.

3.1 Validation of DCB estimation method.

To verify the performance of the DCB estimation method, we take the DCB products from CAS and DLR as references, and we use the mean difference and rms values to evaluate the estimated satellite and receiver DCBs. Figure 3 shows the mean difference and rms of the estimated satellite DCB with respect to CAS and DLR values. It can be seen that the estimated satellite DCB shows a good agreement with CAS and DLR. The mean differences are mostly within ±0.2 ns, and their mean values with respect to DLR and CAS are −0.004 and −0.006 ns, respectively. The rms values are less than 0.4 ns, and the mean rms values with respect to DLR and CAS are 0.18 and 0.23 ns, respectively. It can be concluded that the DCB estimation method can obtain a satellite DCB that is in good agreement with DLR and CAS. In other words, it has little influence on the estimation result of satellite DCB when the receiver DCB is estimated epoch by epoch.

Figure 3 The mean difference and rms of estimated satellite DCB with respect to CAS and DLR.

Figure 4 The difference between the estimated BDS receiver DCB with respect to CAS and DLR.

Figures 4 and 5 present the mean difference and rms of the estimated BDS receiver DCB with respect to CAS and DLR values. It should be noted that the receiver DCB for some of our selected MGEX stations cannot be provided by CAS and DLR since these stations were not used in their DCB estimation. In this study, 101 station reference values can be provided by CAS and 61 station reference values can be provided by DLR. As it can be seen, when DLR is taken as a reference, the mean difference mostly ranges from −0.5 to 0.5 ns, and their mean value is 0.05 ns. The corresponding rms values are mostly less than 0.5 ns, and the mean rms is 0.43 ns. However, the difference between our estimated receiver DCB and CAS shows a larger value. The difference ranges from −2 to 2 ns; the values for some stations located in low latitudes are more than 1 ns. In terms of rms, the values for most stations are less than 1 ns, and the corresponding mean rms is 0.80 ns. The reason may be related to the single-station ionosphere modeling adopted by CAS, since the modeling accuracy of stations in low latitudes is low. The evaluation results show that the DCB estimation method can obtain the receiver DCB value with certain accuracy.

Figure 5 The rms of estimated receiver DCB with respect to CAS and DLR.

3.2 Intraday stability of receiver DCB

In this study, the time series values of the receiver DCB within 1 d can be obtained through the epoch-by-epoch estimation for the receiver DCB. In order to better understand the intraday variations and stability of receiver DCB, we calculate and analyze the intraday standard deviation (iSD) and the intraday fluctuation as

$$\begin{array}{}\text{(7)}& \begin{array}{l}\mathrm{iSD}=\sqrt{\frac{\sum \left({\mathrm{RDCB}}_i-\stackrel{\mathrm{‾}}{\mathrm{RDCB}}\right)}{N-\mathrm{1}}},\\ {\mathrm{fDCB}}_i={\mathrm{RDCB}}_i-{\mathrm{RDCB}}_{\mathrm{1}},\end{array}\end{array}$$

where iSD is the intraday SD, RDCB_i is the receiver DCB at epoch i, $\stackrel{\mathrm{‾}}{\mathrm{RDCB}}$ is the corresponding mean value of 1 d, N is the number of the available epochs, fDCB_i is the corresponding value of fluctuation at epoch i, and RDCB₁ is the receiver DCB at the first available epoch.

Figure 6 The intraday SD of receiver DCB for all stations on 4 consecutive days (2, 3, 4, and 5 January 2019).

Four consecutive days (2, 3, 4, and 5 January 2019) are taken as an example. The intraday SD of receiver DCB for all stations on the 4 consecutive days is shown in Fig. 6, and the statistical results of intraday SD are listed in Table 2. As can be seen from Fig. 6, the intraday SD of receiver DCB for all stations shows similar values on 4 consecutive days. This result can also be shown in Table 2, with similar minimum, maximum, and mean values of intraday SD on the 4 consecutive days. It indicates that the intraday stability of the receiver DCB for most stations is almost the same on the 4 consecutive days.

Table 2 Statistics of intraday SD of receiver DCB for all stations on 4 consecutive days.

As can also be seen from Fig. 6, the intraday SD of the receiver DCB for most stations is less than 1 ns on the 4 consecutive days. In particular, the intraday SD of the receiver DCB for the stations located in the Asia Pacific region are significantly less than that for other regions. The reason may be related to greater quantities of available observations of these stations in the Asia Pacific region. In addition, some stations with larger intraday SD of the receiver DCB can be found in low-latitude regions, which may be related to the fact that the ionosphere is more active at low latitudes (Li et al., 2018). It can be concluded that the available observations of the station and the level of ionospheric activity of the area where the station is located may be two important factors affecting the intraday stability of the receiver DCB.

Figure 7 The intraday fluctuation (blue line) and intraday SD (red line) of the receiver DCB for the nine selected stations within 30 d.

Figure 7 shows the intraday fluctuation and intraday SD of the receiver DCB for the nine selected stations located in the Asia Pacific region within 30 d. It can be seen that the intraday fluctuation of the receiver DCB is mostly within ±1 ns, and the corresponding intraday SD is mostly less than 0.4 ns. It is clear that the receiver DCB with larger fluctuation shows lower intraday stability and vice versa. Although intraday receiver DCB of the nine selected stations are relatively stable, larger fluctuations with more than 2 ns of intraday receiver DCB can be found.

Figure 8 Statistical distribution of intraday SD of receiver DCB for all station in 30 d.

In order to better explain the intraday stability of the receiver DCB, Fig. 8 shows the statistical distribution of intraday SD of the receiver DCB for all stations on 30 d, where $\left|\mathrm{0.5}\right|$, $\left|\mathrm{1}\right|$, and $\left|\mathrm{2}\right|$ represent the percentage of the absolute value of intraday SD is less than 0.5, 1, and 2 ns, respectively. More than 94 % of the intraday SDs are less than 1 ns, indicating that most of the receiver DCBs show relative stability within 1 d with the intraday SD of less than 1 ns.

3.3 Comparison between BDS-2-only and BDS-3-only solutions

As we all know, BDS-3 is the next-generation satellite navigation system. The additional BDS-3 observations used in this study contribute to the increasing available observables in DCB estimation. To further understand the contribution of BDS-3, we compare the intraday stability of the receiver DCB obtained by using BDS-2-only and BDS-3-only observations. Figures 9 and 10 show the mean intraday SD of receiver DCB and the corresponding statistical distribution for all stations within 30 d by using BDS-2-only and BDS-3-only observations. Obviously, the result shows that the intraday stability of the BDS-3 receiver DCB is better than that of the BDS-2 receiver DCB.

Figure 9 The mean intraday SD of receiver DCB for all station in 30 d by using BDS-2-only and BDS-3-only observations.

Figure 10 Statistical distribution of intraday SD of receiver DCB by using BDS-2-only and BDS-3-only observations.

4 Conclusions

The short-term variations of epoch-by-epoch BDS receiver DCB are estimated and analyzed with additional BDS-3 observations. The BDS receiver DCBs with the additional BDS-3 observations are analyzed through 30 d of Multi-GNSS Experiment observations. To verify the performance of the DCB estimation method, we take the DCB products from CAS and DLR as references, and we use the mean difference and rms values to evaluate the estimated satellite and receiver DCBs. The estimated satellite DCB shows a good agreement with CAS and DLR, since the mean differences are mostly within ±0.2 ns, and their mean values with respect to DLR and CAS are −0.004 and −0.006 ns, respectively. The rms values are less than 0.4 ns, and the mean rms values with respect to DLR and CAS are 0.18 and 0.23 ns, respectively. In terms of the intraday variability of the receiver DCB, most of the receiver DCBs show relative stability within 1 d with the intraday SD of less than 1 ns. However, larger fluctuations with more than 2 ns of intraday receiver DCB can be found. It can be concluded that the available observations of the station and the level of ionospheric activity of the area where the station is located may be two important factors affecting the intraday stability of the receiver DCB. As shown by the results of the intraday stability of receiver DCB obtained by using BDS-2-only and BDS-3-only observations, the intraday stability of the BDS-3 receiver DCB is better than that of the BDS-2 receiver DCB.