BDS-3 differential code bias estimation with undifferenced uncombined model based on triple-frequency observation

Since 2015, the new generation global BDS system, i.e., BDS-3, has started its development with five experimental satellites demonstration system and has announced its initial global service officially on December 27, 2018. Among the various characteristics to be analyzed for the new generation BDS satellites, the differential code bias (DCB) is of special attention since that it has a direct dependence on the new signals, i.e., B1C and B2a, and it is one of the most intricacy problems in the ionosphere sensing and positioning with multi-GNSS and multi-frequency observations. To take the full capability of the triple-frequency BDS signals, this paper proposed a new method for the DCB estimation in which the undifferenced uncombined observations are processed in PPP mode. In addition, with the intention to estimate all the unknowns, including the DCB, in a single filter, the DESIGN (deterministic plus stochastic ionosphere model for GNSS) method is applied for the ionospheric delay constrains in this method. In the formula derivation, special attention is given to the DCB and clock parameters due to different frequencies for B1I/B1C, etc. Then, the efficiency of the new method is assessed with observations of 23 iGMAS stations capable for BDS triple-frequency tracking and 21 IGS stations capable for GPS triple-frequency tracking during DOY 001 to 090, 2019. Moreover, the traditional DCB estimation method by employing the geometry-free (GF) combination with the ionospheric delay removed by global ionosphere map product is also performed for comparison purpose. The experimental results suggest that by using the undifferenced uncombined solution, rather than the GF combination, the BDS-2 DCB on B1IB2I and B1CB3I can be improved, especially for the MEO satellites. Regarding to the DLR products, the undifferenced uncombined DCB solution presents a RMS of 0.32 ns and 0.27 ns for B1IB2I and B1CB3I, respectively. Concerning the BDS-3 satellites DCB, it is GF combination that performs better by a factor of 12.7%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$12.7\%$$\end{document} and 15.2%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$15.2\%$$\end{document} for B1CB2a and B1CB3I, respectively. This is mainly due to the fact that the undifferenced uncombined DCB solution is sensitive to the limited precision of the BDS-3 orbit and clock. This conclusion is further confirmed by the improvement in the GPS DCB solution with the new method. Compared with the GF combination solution, the STD for daily repeatability improves from 0.088 to 0.061 ns and 0.119 to 0.090 ns for satellite on C1WC2W and C1WC5X, respectively, by using the undifferenced uncombined model.


Introduction
To establish a global navigation satellite system independently, China has been working on the development of BDS since 1990 and enables the operational service of centimeter level relative and standalone positioning throughout the Asia-Pacific region since 2012, which is referred to as BDS-2 (Shi et al. 2012a;Yang et al. 2014). Currently, BDS is under the construction of the global service system as scheduled by the three-step development strategy, and this new generation BDS is simply called BDS-3 .
The BDS-3 starts its development with a five experimental satellites demonstration system since 2015, including two IGSO satellites, I1-S (C31) and I2-S (C32), and three MEO satellites, M1-S (C33), M2-S (C34) and M3-S (C35) (Xie et al. 2017;Yang et al. 2018). The goal of this demonstration system is to comprehensively test and analyze its new features, e.g., the new types of the navigation signals, the inter-satellite links (ISL), and the uploading injection and downloading reception signal, etc. After the validation of the designed performance, the first two BDS-3 satellites (C19 and C20) were launched on November 5, 2017. With the intention to achieve fully operation by 2020, the deployment of BDS-3 is speed up since then and launched 11 new satellites up to August 25, 2018(China Satellite Navigation Office 2017. Now, the BDS-3 constellation is consisted of 1 GEO, 2 IGSO and 20 MEO satellites, transmitting new navigation signals, referred to as B1C (1575.42 MHz) and B2a(1176.45 MHz). As a global system is to operate compatibly with other GNSS, these new signals are designed to share the same frequencies and similar modulation method to GPS and Galileo (Xiao et al. 2016).
Along with the deployment of BDS-3, a wide range of valuable literatures have been published. Attributed to the project of International GNSS Monitoring and Assessment Service (iGMAS), there are over 10 stations that are capable for BDS-3 satellite tracking. Though the number of stations is rather limited, Tan et al. (2016) presented the first attempt of precise orbit determination (POD) for the experimental satellites with the old B1, B2 and B3 signals. Xie et al. (2017) presents the POD results with the new navigation signals, and the orbit overlap and the satellite laser ranging (SLR) evaluation both suggested that the POD precision is about 2-6 dm in 3D and 1-2 dm in the radial component. As a great advantage of the new generation BDS, the contribution of the ISL in terms of POD was analyzed and confirmed by Ren et al. (2017) with different data processing strategies. Taking the limited measurements of BDS-3 into consideration, Hu et al. (2018) refined the cycle slip detection and repair strategy in their study to improve the usability of the observations as well as the precision of POD. Besides the frequency compatibility of the new signals, a better performance in terms of measurement noise and anti-multi-path property is also expected. Specifically, the satellite-induced pseudo-range multi-path error that elevation angle related in BDS-2 has been investigated in great details and has been modeled by Wanninger and Beer (2015) for IGSO and MEO and Lou et al. (2016) for GEO, respectively. As indicated by Xie et al. (2017), Zhang et al. (2017b) and Yang et al. (2018), this kind of systematic errors of the BDS-3 signal is obviously reduced to a negligible level. Concerning the new signals, it is stated that the RMS of B1C signal is higher than those of other signals at the same elevation angle. The RMS of the pseudo-range noises of B2a+b is the smallest ). In addition, equipped with the rubidium atomic frequency standard and the passive hydrogen maser frequency standard, Xie et al. (2017) argued that the clock of BDS-3 performs better than that of BDS-2. It should be emphasized that all the above-mentioned studies were mainly focused on the experimental satellites since that only two BDS-3 satellites were available by the end of 2017. For more details on the performance of BDS-3, including the qualities of signals in space, the time synchronization and timing precision, etc., we refer to Yang et al. (2018).
Though the elevation angle-dependent variation of satellite code bias is removed for BDS-3, the constant part still exists in both the receiver and satellite hardware, which is usually referred to as differential code bias (DCB) since its absolute value is actually a nonestimable parameter (Warnant 1997;Schaer and Dach 2010;Montenbruck et al. 2014). Typically, the DCB can be classified into two categories: The first type is the inter-frequency bias, which represents the bias between code observations at two different frequencies, e.g., P1P2 DCB. The second type is the intra-frequency bias, which represents the bias between two code observations at the same frequency, e.g., P1C1 DCB (Leandro et al. 2007;Wang et al. 2016). In this paper, we focus our attention on the inter-frequency DCB. For high-precision GNSS applications, e.g., single-frequency precise point positioning (PPP), ionospheric delay modeling, the effect of DCB should be carefully calibrated. Thus, a number of methods have been proposed for the estimation of GNSS DCB (McCaffrey et al. 2017). Notably the work of IGS Ionosphere Working Group (Iono-WG), in which the GNSS DCB is generated routinely as a by-product of global ionosphere map (GIM), and typically the sphere harmonic function (SHF) is employed for the ionospheric delay modeling to separate the DCB (Schaer 1999). In this case, the precision of DCB estimation mainly depends on the density and distribution of the global GNSS tracking network. Thus, for the BDS satellite DCB solution with a sparse network in the early stage, Li et al. (2012) proposed a new approach, i.e., IGGDCB, in which the ionospheric vertical TEC is modeled by generalized triangular series (GTS) function with adaptive parameter set for each individual station, and the demonstration suggested that the precision of IGGDCB is around 0.1 and 0.2 ns for GPS and BDS, respectively (Li et al. 2012;Shi et al. 2016). Alternatively, Montenbruck et al. (2014) suggested to estimate the multi-GNSS DCB by removing the ionospheric delay with the GIM products, and depending on the number of stations involved, a daily stability of 0.05-0.3 ns has been reported. As the method is rather easy to implement, it has been widely used for code bias analysis since then, given that the high precise ionosphere correction, for instance, GIM is available. In addition, it should be noted that the geometry-free (GF) combination observation is typically enabled in these studies. To access the full capabilities of the multi-GNSS multi-frequency signals, the undifferenced uncombined observation model has drawn increasing interest in the GNSS community (Schönemann et al. 2011;Gu et al. 2015b). In addition, it has also been demonstrated as an efficient way for simultaneous estimation of ionospheric delay and DCB even with low-cost receivers (Gu et al. 2013;Zhang et al. 2017a).
Obviously, the performance of the undifferenced uncombined observation model mainly depends on the capability of the particular ionospheric parameterization method involved. Thus, the DESIGN (deterministic plus stochastic ionosphere model for GNSS) model is developed by Shi et al. (2012a) and further refined by Zhao et al. (2018) for high-precision undifferenced uncombined GNSS data processing, and the PPP results constrained with DESIGN are promising for both single-and dual-frequency observation (Lou et al. 2015;Zhao et al. 2018).
Concerning the characteristics of BDS-3 DCB, especially the DCB for new navigation signals on both side of satellite and receiver, as well as the performance of DCB estimation based on the undifferenced uncombined model with a limited stations, this contribution performs a comprehensive analysis based on 3-month BDS-2 and BDS-3 triple-frequency observations collected with 23 iGMAS. In addition, the GPS triple-frequency DCB solution is also carried out as a reference for the 3 months with 21 IGS stations. This paper is organized as follows: First, some conventions are adopted in Sect. 2. Afterward, the DCB estimation methods are introduced in Sect. 3, including the undifferenced uncombined observation model and GF combination method based. In Sect. 4, the DCB estimation experiments are carried out and analyzed. Finally, the conclusion is drawn in Sect. 5.

Notation
In this paper, we adopt the following conventions: Matrices and vectors are denoted in bold form, while scalars are denoted in regular form. The term in its bold form stands for the corresponding vector, e.g., P r = P 1 r P 2 r . . . P n s r T and t = t 1 t 2 . . . t n s T are the pseudo-range observed minus computed (OMC) vector for receiver r and satellite clock vector for all the n s satellites. For more details concerning these vectors, we refer to Sect. 3.1. In addition, a few notations are defined for future reference: as defined, z s is a s × 1 vector with elements of zeros and u s is a s × 1 vector with elements of ones, while Z s is a s × s zero matrix and U s is a s × s identity matrix. The dimensions and lengths of such vectors will generally be obvious from context. J i j is the IF transformation matrix for observation on frequency i, j(i = j), and i, j ∈ ( 1 2 3 ) as for BDS observation in this study.

Methods
The undifferenced uncombined observation of the GNSS pseudo-range and carrier phase is generally expressed as in which P s r , f , Φ s r , f are pseudo-range and carrier phase from receiver r to satellite s(s = 1, . . . , n s ) on frequency f in length units, respectively. ρ is the geometric distance, while antenna phase center corrections as well as the phase windup should be applied to P, Φ before ρ becomes unassociated with the frequency; t( s r ) is the receiver and satellite clock error in length units, respectively; T z is the zenith tropospheric delay that can be converted to slant with the mapping function α; I denotes the zenith TEC with the frequency and mapping function M s -dependent factor β f = 40.3 f 2 ·M s ; N is the float ambiguity in length units; b s r , f is the frequencydependent code bias, which can be expressed as where b r , f and b s f are the constant code bias that receiveronly dependent and satellite-only dependent, respectively; b s f (E) is the satellite-induced elevation-dependent variation; b s type, f depends on specific receiver hardware and tracking implementation of satellite systems as argued by Gong et al. (2018); thus, it is the code bias both receiver and satellite related. For BDS-2, b s f (E) and b s type, f of different satellites can be corrected according to Lou et al. (2016) and Gong et al. (2018), respectively. While, for BDS-3, b s f (E) is negligible, and b s type, f is ignored in this study since that the number of BDS-3 stations with new signal is too limited for a reliable analysis which requires sufficient number of tracking stations equipped with the same type of receivers. Thus, b s r , f is further simplified as

Undifferenced uncombined observation model based on DESIGN
In the code bias estimation based on undifferenced uncombined observation model, the ρ s r is assumed to be known with precise satellite orbit and receiver coordinate, while the IGS final satellite clock product is also applied. Concerning the precise orbit and clock of BDS satellites, especially for the BDS-3, the products of Wuhan University are utilized (ftp:// igs.gnsswhu.cn/pub/gnss/products/mgex). And the BDS-2 and BDS-3 are treated as one system by using B1I/B3I measurements for both BDS-2 and BDS-3 satellites. As the precise clock is estimated with B1/B3 IF combination in which the satellite clock and code bias lumped together (Xie et al. 2017), i.e., the symbol := here means "is replaced by." Similarly, denote i.e., estimate receiver clock for BDS-2 and BDS-3 separately. Then, by substituting precise satellite clock (9) into Eq. (6) and recall the IF transformation matrix J i j , we have the observation equations in which P s r and Φ s r are the OMC for pseudo-range and carrier phase, respectively, i.e., P s r = P s r − ρ s r and Φ s r = Φ s r − ρ s r with precise coordinate of satellite and receiver; as presented, since the satellite products used in this study are generated simultaneously with the same signals, i.e., B1I/B3I for both BDS-2 and BDS-3. When applying these products in the BDS-3 data processing with the new signals, special attention should be paid to the clock and DCB parameters, as well as the value of frequencies for different signals.
For the ionosphere delay I s , we adopt the DESIGN method as suggested in the authors' previous studies (Shi et al. 2012a;Lou et al. 2015;Zhao et al. 2018) while a 0 is the average value of ionosphere delay over the station; a 1 , a 2 and a 3 , a 4 are the coefficients of the two second-order polynomials along east-west and south-north direction, respectively; r s r is the residual ionospheric effect for each satellite; dL( s r ), dB( s r ) are the longitude and latitude differences between the ionospheric pierce point (IPP) and the approximate location of station, respectively.
Furthermore, the daily variation of deterministic part is expressed as Fourier series with fixed frequency components f 1 = 1/24, f 2 = 1/12, f 3 = 1/8. Then, we estimate the amplitude, i.e., x 0 , y j , z j , as daily constant that updated every 24 h while the stochastic part is estimated as random walk with the variogram summarized as Eqs. (15) in which a s is the maximum correlation distance, which is about 9000 s (Zhao et al. 2018), and c s varies with geomagnetic latitude B as Gaussian function as while the variation of ionosphere activity is modeled by c s,min and c s,max with Epstein function as with x dependent on the sunspot number s n x(s n ) = (s n − 100)/20 In addition, as suggested by Shi et al. (2012a), the a priori ionosphere delay correction I s corr retrieved from, e.g., GIM (Global Ionosphere map), can usually be applied as constrains to further improve the solution I s corr = a 0 + a 1 dL + a 2 dL 2 + a 3 dB + a 4 dB 2 + r s r + ε I When constrained with DESIGN model, the ionosphere delay can be separated from the code bias parameters as demonstrated by Zhao et al. (2018). However, the receiver clock error and the code bias are still linear dependent, and the following conditions should be introduced (Gu et

Geometry-free combination model based on GIM
Alternatively, another efficient way for the code bias estimation is based on the GF combination as suggested by Montenbruck et al. (2014), in which the GIM is introduced to remove the ionospheric delay. From the between-frequency difference measurements, we have the GF combination of frequency 1 and frequency 2, as well as the GF combination of frequency 1 and frequency 3 as follows: in which the superscript i j represents the frequency. Note that the satellite ambiguity and code bias remain constant in continuous observation arc, the sum of the satellite ambiguity and code bias can be obtained by moving average method: in which nArc represents the epoch number of the tracking arc for satellite s. For simplicity, (23) and (22), the slant ionospheric delay together with code bias can be obtained: the zenith ionospheric delay can be eliminated by GIM product. As a result, satellite code bias, together with receiver code bias, can be extracted as: however, it is noted that b r ,i j and b s i j are linearly dependent; thus, a zero-mean condition is introduced as: as we mentioned before, the receiver code bias of BDS-2 and BDS-3 differs from each other. Therefore, two zero-mean conditions should be set for BDS-2 satellites and BDS-3 satellites, respectively.

Data and strategy
To cope with the increasing demand of high-precision GNSS applications, such as the multi-frequency data processing, real-time applications, the GNSS/INS integrated navigation, as well as the processing of new signals transmitted by BDS-3, we have launched the FUSING (FUSing IN Gnss) project. At present, the FUSING software is capable for the precise orbit determination (POD) of GNSS, high-rate satellite clock estimation, ionosphere and troposphere modeling, satellite DCB estimation, as well as multi-frequency multi-GNSS precise positioning Zhao et al. 2018;Yang et al. 2019). To further confirm the above analysis with numerical experiment, we realized both the undifferenced uncombined method with DESIGN model and the geometry-free combination method for DCB estimation, based on FUSING platform. DCB on both BDS-2, BDS-3 and GPS are generated with different methods. Observations of 23 iGMAS stations during DOY 001 to 090, 2019 with an interval of 30 s are collected for BDS-2 and BDS-3 DCB estimation. From our experience, the data quality of these 23 iGMAS stations is not as stable as that of IGS stations. Thus, observations from 21 IGS stations of the same period are collected to assess   Fig. 1 Signals Table 2 Receiver types Table 3 Receiver Daily repeatability of estimated DCB values a Abbreviated as UDUC for simplicity in the following discussion b Actually the GF combination is insensitive to the precision of the geometry terms, e.g., receiver coordinate, etc   (19), the GIM products are utilized as pseudo-range for our method, while, for GF combination, the GIM products of CODE are utilized to correct the ionosphere delay directly. In addition, the signal and receiver details are presented in Tables 2  and 3, respectively. As we can see, all the receivers of the iGMAS network used in this study are produced by Chinese manufacturers.

BDS DCB estimation
Since the variation of BDS-2 satellite DCB is eliminated with the model suggested by Lou et al. (2016), and the variation of BDS-3 satellite DCB is negligible as argued by Yang et al. (2018), the satellite DCB is estimated as daily constant in this study. The DCB estimation performance is evaluated by STD for each satellite over the experimental period. Then, the averaged STD of BDS-2, for instance, is derived by taking the mean value of STD for all the BDS-2 satellites. The statistics of the BDS satellite DCB estimation for the two strategies are shown in Fig. 2. Note that the STD was shown as the error bar with its value magnified ten times, i.e., in unit of 0.1 ns. Over all, the averaged STD of BDS B1B2 DCB is 0.27 ns and 0.29 ns for undifferenced uncombined and GF combination solution, respectively, as presented in Fig. 2, while for BDS B1B3 DCB, the averaged STD is 0.17 ns and 0.16 ns for undifferenced uncombined model and GF combination model, respectively. A detailed analysis suggested that the undifferenced uncombined model outperforms the GF combination model in the BDS-2 DCB solution by an improvement of about 28.1%. Concerning the BDS-3 DCB solution, however, the GF combination model has a better precision by a factor of about 13.5%. This is reasonable since In the upper part shows the B1B2 DCB results, and the B1B3 DCB results are shown in the bottom part. Besides, BDS-2 satellites are arranged in the left part, while BDS-3 satellites are arranged in the right part with background painted gray. The STD was shown as the error bar with its value magnified ten times, i.e., in unit of 0.1 ns. The MEO satellites of BDS-2 are highlighted in red in the x-axis label that the DCB solution based on the undifferenced uncombined model relies on the precision of satellite orbit and clock products. Unfortunately, the BDS-3 products have a limited precision due to its imperfect models and insufficient number of tracking stations in its current state. Similarly, the most significant improvement in the BDS-2 DCB estimation can be found for the MEO satellites, i.e., C11, C12, C14, the orbit and clock of which are typically better determined than the GEO and IGSO satellites. By the comparison between the results on different frequencies, it is illustrated that the B1B3 DCB generally outperforms that of B1B2 for both the undifferenced uncombined solution and the GF combination solution. This is not surprising since that the signal quality on B3 is better than that of B2.   Figure 3 further presents the series of BDS DCB over DOY 001 to DOY 090, 2019 with the undifferenced uncombined solution in black and the GF combination solution in red. It is noted that for some satellites there is no DCB generated since that the satellite is failed to be tracked by all the 23 stations.
Apart from satellite DCB, the results of receiver DCB are also presented in Tables 4, 5 and 6 for BDS-2 B1B2, BDS-2 B1B3, BDS-3 B1B2 and BDS-3 B1B3, respectively. It is found that the UNICORE stations, i.e., BYNS, CNYR, KUN1 and SHA1, cannot track the new signals of BDS-3, while the other receivers cannot track B1I/B2I signals of BDS-2, and all these iGMAS stations are capable for BDS-2 B1I/B3I signal tracking. Unlike the satellite BDS-2 DCB solutions, the receiver DCB solutions with different methods have the consistent precision with an averaged STD of about 0.8ns and 0.47ns for B1B2 and B1B3, respectively. Actually, the result is not surprising since that though the MEO satellite performs better with the undifferenced uncombined solution, the result of GEO and IGSO is roughly the same for both solutions as presented in Fig. 2. Concerning the receiver BDS-3 DCB solution in Fig. 6, we find that the precision can be improved from 0.56ns to 0.46ns and from 0.42 to 0.37 ns for B1B2 and B1B3, respectively, by the GF combination solution. The comparison between B1B2 and B1B3 receiver DCB of BDS-3 again confirms that the B3I signal performs better than B2a. In addition, the satellite DCB performs much more stable than that of receiver no matter which method or which signal was involved.

GPS DCB estimation
Taking the imperfect BDS-3 satellite orbit and clock products, as well as the poor quality of BDS-3 measurements into consideration, we carried out the GPS satellite DCB estimation in this section to further confirm that with the high-precision satellite products, the undifferenced uncombined model outperforms the traditional GF combination in DCB estimation.
The results of DCB estimates over DOY 001 to DOY 090, 2019 for each GPS satellites are plotted in Fig. 4 for different methods. We can find that compared with the GF combination DCB solution, the averaged STD over the experimental  Tables 7 and 8 for P1P2 and P1P5, respectively. Note that the stations equipped with the LEICA receivers, i.e., KRGG, MYVA, THTG, and the stations equipped with the SEPT receivers, i.e., AREG, HARB, KOUG, NRMG, cannot track the C5X signals that analyzed in this study. Over all, the averaged STD for undifferenced uncombined DCB combination solution performs slightly better than that of GF combination DCB solution. It should be emphasized that in the calculation of the averaged STD over different stations, the stations REUN, THTG and ULAB are not included since the receiver code bias of these stations exhibits remarkable variability during the experimental period as illustrated in Fig. 5. For example, there are jumps around DOY 063, 2019 for the station REUN on P1P2 and P1P5, and the station THTG on P1P2, while a significant fluctuation is also observed for the station ULAB on P1P2 and P1P5. The study of Zhang et al. (2016) also pointed out that code bias of some receiver may exhibit instability due to the hardware configuration and the environment and thus should be given special attention.

Comparison with the MGEX DCB product
In the above discussion, the performance of DCB for BDS and GPS is mainly assessed in terms of daily stability, i.e.,  Before the evaluation, it should be noted that there are no DCB products available for the new BDS-3 signals up to now. Thus, only the DCB for BDS-2 and GPS is analyzed. In addition, though DLR does not provide the C1WC2W and C1WC5X DCB products directly, it provides the daily C1CC1W, C1CC2W, C1CC5X products. Then the reference value of C1WC2W and C1WC5X DCB can be derived from b C1W C2W = b C1CC2W − b C1CC1W and b C1W C5X = b C1CC5X − b C1CC1W , respectively. However, the DCBs b C1W C2W and b C1W C5X derived in this way no longer satisfy the constrains of Eq. (20) due to different number of satellites involved. Thus, there is a constant bias between DLR DCB and ours for all the satellites, which is further removed in the following statistics.
By the comparison between the DLR DCB and ours day by day for each satellite each product, the corresponding residual distribution on different frequencies is plotted in Figs. 6 and 7 for BDS-2 and GPS, respectively. As illustrated in Fig. 6, the undifferenced uncombined solution has a better accuracy by an improvement of 23.8% and 41.3% for BDS-2 B1B2 and BDS-2 B1B3, respectively. The improvement is mainly attributed to the advantage of undifferenced uncombined model in DCB solution with high-accuracy MEO satellite orbit and clock products as discussed in Sect. 4.2. Concerning the GPS results in Fig. 7, the RMS decreased  from 0.23 to 0.17 ns for GPS P1P2 by using the undifferenced uncombined model, which is consistent with result of Fig. 4. As presented in right panel of Fig. 7, however, the GF combination has better consistency with DLR products and the residual for C1WC5X performs abnormal since that the distribution is far from the normal distribution. The reason may lie in the fact that the C1WC5X stations are rather limited compared with that of C1WC2W for GPS. Among these stations, two stations suffer from exhibiting remarkable variability as shown in Fig. 5.

Conclusion
Concerning the estimation of DCB for the new BDS-3 signals with a limited number of stations, this contribution proposed a new method by utilizing the undifferenced uncombined observation with the ionospheric delay constrained by DESIGN. The experiments are based on the observations collected from 23 iGMAS stations for DOY 001-090, 2019. As a comparison, the GF combination method in which the ionospheric delay is removed by the GIM products is also implemented for DCB estimation. Moreover, the GPS triplefrequency observations from 21 IGS stations of this period are also collected for DCB estimation as a reference.
The DCB of BDS-2 and BDS-3 on both satellites and receivers is estimated. The daily repeatability is utilized as the statistics for the evaluation. Concerning the performance of DCB estimates based on our new method, the STD is 0.205 ns and 0.132 ns for BDS-2 satellite DCB on B1IB2I and B1IB3I, respectively, while the STD is 0.331 ns and 0.204 ns for BDS-3 satellite DCB on B1CB2a and B1CB3I, respectively. Compared with the DCB solution with the GF combination method, the new method performs better in BDS-2 DCB solution by a factor of about 28.1%. A detailed analysis suggested that the improvement is mainly attributed to the high accuracy of BDS-2 MEO satellite orbit and clock products. For the same reason, the BDS-2 DCB outperforms that of BDS-3 significantly with the undifferenced uncombined solution since the precision of BDS-3 satellite ephemeris is rather limited due to imperfect models and insufficient number of tracking stations. In addition, regarding to the DLR DCB products, it is confirmed that our new method gives a better accuracy by an improvement of about 23.8% and 41.3% for BDS-2 B1IB2I and BDS-2 B1IB3I, respectively. Concerning the receiver BDS-2 DCB, the solutions with different methods have the consistent precision with an averaged STD of about 0.8 ns and 0.47 ns for B1IB2I and B1IB3I, respectively, while for receiver BDS-3 DCB, by using the GF combination performs better, and the precision can be improved from 0.56 to 0.46 ns and from 0.42 to 0.37 ns for B1CB2a and B1CB3I, respectively.
To avoid the effects of limited accuracy of ephemeris for BDS, the efficiency of DCB solution based on undifferenced uncombined solution is further assessed with 21 IGS GPS stations. By using the undifferenced uncombined solution, rather than the GF combination, the STD for daily repeatability improves from 0.088 to 0.061 ns and 0.119 to 0.090 ns for satellite on C1WC2W and C1WC5X, respectively. The statistics further confirms that the undifferenced uncombined model performs better in terms of precision when high-accuracy satellite orbit and clock products are applied.
Finally, as mentioned above, due to the limited BDS-3 observations and limited stations capable for BDS-3 tracking, this contribution is still an initial study of BDS-3 DCB estimation. Meanwhile, efforts should be made to improve the BDS-3-related models as well as BDS-3 precise products. With the development of the BDS-3 constellation as well as the tracking networks, it is expected that a more comprehensive and stable result will be achieved in BDS-3 DCB estimation, which will definitely benefit high-precision GNSS applications based on BDS-3.