Characteristics of receiver-related biases between BDS-3 and BDS-2 for five frequencies including inter-system biases, differential code biases, and differential phase biases

It is foreseeable that the BeiDou navigation satellite system with global coverage (BDS-3) and the BeiDou navigation satellite (regional) system (BDS-2) will coexist in the next decade. Care should be taken to minimize the adverse impact of the receiver-related biases, including inter-system biases (ISBs), differential code biases (DCB), and differential phase biases (DPB) on the positioning, navigation, and timing (PNT) provided by global navigation satellite systems (GNSS). Therefore, it is important to ascertain the intrinsic characteristics of receiver-related biases, especially in the context of the combination of BDS-3 and BDS-2, which have some differences in their signal level. We present a method that enables time-wise retrieval of between-receiver ISBs, DCB, and DPB from multi-frequency multi-GNSS observations. With this method, the time-wise estimates of the receiver-related biases between BDS-3 and BDS-2 are determined using all five frequencies available in different receiver pairs. Three major findings are suggested based on our test results. First, code ISBs are significant on the two overlapping frequencies B1II and B2b/B2I between BDS-3 and BDS-2 for a baseline with non-identical receiver pairs, which disrupts the compatibility of the two constellations. Second, epoch-wise DCB estimates of the same type in BDS-3 and BDS-2 can show noticeable differences. Thus, it is unreasonable to treat them as one constellation in PNT applications. Third, the DPB of BDS-3 and BDS-2 may have significant short-term variations, which can be attributed to, on the one hand, receivers composing baselines, and on the other hand, frequencies.


Introduction
The BeiDou navigation satellite system with global coverage (BDS-3) has been fully operational since July 2020 and has the potential to enable a wide range of applications for positioning, navigation, and timing (PNT) all over the world (Yang et al. 2020). Its constellation comprises 30 satellites, including three satellites in GEostationary Orbit (GEO), three in Inclined GeoSynchronous Orbit (IGSO), and 24 in Medium-altitude Earth Orbit (MEO) . To achieve compatibility and interoperability with other global navigation satellite systems (GNSS), and backward compatibility with the BeiDou navigation satellite (regional) system (BDS-2), BDS-3 transmits five navigational signals in space, namely B1I at 1561.098 MHz, B2b at 1207.140 MHz, B3I at 1268.520 MHz, B1C at 1575 MHz and B2a at 1176.450 MHz (Yang et al. 2018;Yuan et al. 2020). Meanwhile, attention should be paid to the fact that, although BDS-3 has already been operational, BDS-2 will still be in service for at least another decade (Yang et al. 2019;Mi et al. 2020b) and forms an important part of BDS globalization. As a regional system serving the Asia-Pacific, BDS-2 constellation comprises five GEO satellites, seven IGSO satellites, and three MEO satellites (CSNO 2019; Montenbruck et al. 2013). Signals at three frequencies are used in BDS-2, namely B1I at 1561.098 MHz, B2I at 1207.140 MHz, and B3I at 1268 MHz, ensuring the PNT service of BDS-2 (Odolinski et al. 2014b;Yang et al. 2014).
It has long been recognized that multiple constellations and multiple frequencies are becoming available that benefit PNT services in accuracy, integrity, and availability (Odijk and Teunissen 2012;Tian et al. 2017). In such cases, in addition to the unification of coordinate and time reference frames, differences in receiver hardware delays related to using signals from different systems should also be considered (Gioia and Borio 2016). These biases are called intersystem biases (ISBs) and are caused by the correlation process within the GNSS receiver (Gao et al. 2017a;Paziewski and Wielgosz 2014). In general, the effect of the receiver ISBs is considered a major source of error in the combined processing of data from different GNSS (Odijk et al. 2016). For example, ISBs have to be considered not only in realtime kinematic (RTK) positioning and precise point positioning (PPP) Geng et al. 2019), but also in integer ambiguity resolution enabled PPP (PPP-RTK) (Khodabandeh and Teunissen 2016). In addition, applications based on multi-GNSS observations, such as time and frequency transfer (Tu et al. 2019;Verhasselt and Defraigne 2019), and atmospheric retrieval (Lu et al. 2020;Pan and Guo 2018), also need to pay attention to the impact of ISBs.
In addition to ISBs, the hardware delay differences experienced by different frequencies in a single GNSS constellation, namely differential code biases (DCB) and differential phase biases (DPB), are also important sources of errors limiting GNSS-based PNT applications (Odolinski and Teunissen 2017b;Sanz et al. 2017). For example, the lumped effect of satellite and receiver DCB and DPB is generally considered a major source of error in ionospheric retrieval from GNSS observables (Brunini and Azpilicueta 2010). Fortunately, the satellite DCB and DPB are fairly stable over a considerable time for each GNSS constellation (Sardon et al. 1994). The variability of receiver DCB and DPB may be evident over a relatively short period, e.g., two hours to one day, due to temperature perturbations around the receivers (Zha et al. 2019;Zhang and Teunissen, 2015;Zhang et al. 2016). Thus, the handling of receiver DCB and DPB is important to ensure GNSS-derived ionospheric retrieval accuracy and reliability. Furthermore, precision and reliability of positioning and time and frequency transfer are also constrained by receiver DCB and DPB (Dach et al. 2002;Odolinski et al. 2015). With high-precision and high-reliability DCB and DPB, PPP and RTK, as well as PPP-RTK can reach even higher levels (Gao et al. 2017b;Odolinski and Teunissen 2017a;Psychas et al. 2019). Moreover, accurate calibration of DCB and DPB is an important prerequisite for time and frequency transfer based on GNSS (Defraigne and Baire 2011;Huang and Defraigne 2016).
Therefore, it is important to ascertain the intrinsic characteristics of receiver-related biases between BDS-3 and BDS-2, since they are treated, quite naturally, as one system. For this purpose, we first present a method that allows DCB, DPB, and ISBs to be estimated simultaneously and continuously. With this method, the characteristics of code and phase ISBs between BDS-3 and BDS-2, and the DCB and DPB of BDS-3 and BDS-2 are determined.

Calibration of receiver-related biases using multi-GNSS observables
The system of code and phase observation equations based on single-differenced (SD), serving as a starting point of developing the algorithm, reads (Odolinski et al. 2015) where p s * ab,j (i) and s * ab,j (i) are the SD code and phase observations associated with two receivers a and b . GNSS constellation is * , which is distinguished by different letters ( * = A, B, … ). The satellite identifier is s * , the frequency is j , and i denotes the epoch. x ab (i) is the column vector of geometric unknowns and the corresponding coefficient g s * ab (i) is the receiver-to-satellite unit vector. dt ab (i) is the SD receiver clock between a and b . The symbols d * ab,j (i) and * ab,j (i) denote, respectively, the SD receiver code biases and phase biases.l s * ab is the SD slant ionospheric delay, and j = f 2 1 f 2 j is the frequency-dependent factor.T s * ab denotes the SD tropospheric delay. z s * ab,j denotes the SD ambiguity and j denotes the wavelength, respectively. Note that all of the variables involved in (1) are in meters, except z s * ab,j is in cycles. s * ab,j and e s * ab,j denote the SD random observation noise and unmodeled effects such as multipath for the code and phase observations, respectively. Consider a zero baseline or a short one of less than a few kilometers, we can assume no differential ionospheric and tropospheric effects. In this case, the model can also be referred to as an ionospheric-fixed and tropospheric-fixed model. Then, the SD code and phase observation is rewritten as (Odolinski et al. 2014a) (1) However, even though (2) does not depend on ionospheric and tropospheric delays, the model is still not of full rank. A rank deficit occurs in two ways. One is the linear dependency between the columns of the receiver clock and the code/phase biases, and the other is the column dependency between phase biases and the SD ambiguities (Mi et al. 2019a, b). By applying the S-system transformation, full rank can be achieved by constraining a minimum set of parameters, or S-basis (Zhang et al. 2018). It should be noted that the choice of S-basis is not unique, which dictates the estimability and the interpretation of parameters. For example, the rank deficiency of the first type can be solved by fixing the receiver code biases of one-or multi-constellation, which results in two different models: classical differencing and inter-system differencing. A detailed explanation and comparison of classical differencing and inter-system differencing can be found in Mi et al. (2020a), which will not be repeated here. In our study, inter-system differencing is adopted. Thus, the rank deficiency between the columns of the receiver clock and the code/ phase biases are solved by fixing the code biases on the first frequency of only one of the constellations. Also, we have the rank defects between phase biases and ambiguities, which are solved by fixing the SD ambiguities of one reference satellite per constellation. Once the rank defects have been solved, we have the following full-rank model where d A ab,j (i) and ̃A ab,j (i) are the DCB and DPB, and their counterpart d * ab,j (i) and ̃ * ab,j (i) are the code and phase ISBs, and their interpretations are given in Table 1.
Recall that in the estimable form of DPB ab,1 correspond to the reference satellites s A = 1 A and s * = 1 * . As is well known that one reference cannot be visible for a long time (usually less than a few hours). In general, when the observation period exceeds a few tens of hours, it is inevitable to change the reference satellite more than once. However, in this case, abrupt jumps will be introduced in the epoch-wise estimates of DPB and phase ISBs, which is not helpful in restoring their characteristics.
The datum ambiguities can be considered time-invariant as long as the reference satellites are tracked continuously without cycle slip. Thus, for multi-epochs, the number of rank defects between phase biases and the ambiguities remains unchanged and is independent of the epoch number. Hence, when the SD ambiguities of one satellite ( j z 1 * ab,j ) are selected as a datum at epoch i , the SD ambiguities of the remaining satellites ( j z s * ab,j , s = 2, 3, ⋯ , m ) will absorb j z 1 * ab,j , and thus have the DD form ( j z 1 * s * ab,j ). Then, j z 1 * s * ab,j can be transferred to epoch i + 1 , even though the reference satellite is no longer visible, j z 1 * ab,j in j z 1 * s * ab,j can still serve as a datum. In this case, the SD code and phase observation equations at epoch i + 1 read: Receiver phase bias of the first frequencỹ . In this way, the DPB and phase ISBs are estimated continuously without changing the reference satellite.
To accurately calibrate these biases, the baseline and DD ambiguities are precisely estimated in advance using a strategy that does not change the reference satellite. Then, the baseline and DD ambiguities are subtracted from (3). In this case, the SD code and phase observation equations with fixed baseline and ambiguities are expressed as, Since the datum in DD ambiguities has not changed, no discontinuity will be present in the estimates of DPB and phase ISBs.

Experimental setup
In our analysis, we selected two sets of GNSS data, measured by three and four collocated receivers, respectively, with five observation types (B1I, B1C, B2a, B2b, B3I) of BDS-3 and three types (B1I, B2I, B3I) of BDS-2. See Table 2 for detailed characteristics. Two points deserve noted from the table. First, the receivers with IDs APM1, APM2, APM3, and APM4, which comprise two Trimble and two Septentrio receivers, are connected to a common antenna, implying that they can create a total of six zero baselines. Second, the receivers IGG1, IGG2, and IGG3 are each equipped with a single antenna, creating three baselines, with lengths of 1.8 m, 5.6 m, and 6.7 m. The sampling interval of the first set of receivers (AMP1 to AMP4) was 30 s and that of the second set (IGG1 to IGG3) was 10 s.
The cutoff elevation angle was set to 10° to discard particularly noisy code and phase observations. The elevationdependent weighting function used can be expressed as s * = u ∕ sin(E s * ) and s * p = u p∕ sin(E s * ) (Euler and Goad 1991), where s * and s * p donate the standard deviations of the phase and code observations of satellite s * , E s * r is the elevation of satellite, u and u p are the undifferenced zenith-referenced a priori phase and code standard deviations, assumed here as 3 mm and 0.3 m, respectively. The satellite positions that are required for elevation angle determination are computed using the broadcast ephemeris. The Detection, Identification and Adaptation (DIA) procedure is used to detect and eliminate the effect of outliers (Teunissen 2018), and the LAMBDA method was used for integer ambiguity resolution (Teunissen 1995;Chang et al. 2005). For the sake of brevity, only partial baselines selected during some of the experimental days are reported here in the analysis of ISBs, DCB and DPB estimates. These results are representative of all the experimental results that we obtained. In addition, in order to reduce the impact of the multipath effect caused by the tall buildings around the receivers, sidereal filtering is implemented (Wang et al. 2018). In sidereal filtering, the multipath model needs to be shifted by a certain period, usually close to a sidereal day, thus, this method is only applicable to APM1 to APM4 with multiple days of observations available.

Characterization of BDS-3 and BDS-2 ISBs estimates
This section first describes the B1I and B3I code and phase ISBs estimates between BDS-3 and BDS-2, as those frequencies are identical in the two constellations. Then, the characterization of ISBs estimates between BDS-3 B2b and BDS-2 B2I, which are overlapping frequencies but with different signal modulations, is analyzed separately. Figure 1 shows the estimates of B1I and B3I code ISBs between BDS-3 and BDS-2. When comparing the righthand panels showing the B3I code ISBs, with the left-hand panels of B1I ISBs, we can see that the B3I code ISBs estimates fluctuate randomly around zero. In other words, for B3I signals, the code ISBs are not shown to be present for  Figure 2 depicts the B1I code ISBs estimated for each of the three baselines composed of different receivers connected to separate antennas. As might have been expected, the B1I code ISBs between BDS-3 and BDS-2 is not close to zero for any of the three baselines, indicating a significant bias. The trend in the ISBs estimates, particularly for IGG1-IGG3 and IGG2-IGG3, can be attributed to multipath effects. However, due to the short observation time, sidereal filtering cannot be used to weaken the influence of multipath. Even so, these results confirm the previous observation that the B1I code ISBs based on a mixed-receiver combination are nonzero, and this suggests that we should consider the difference between B1I code of BDS-3 and BDS-2 in practice. Figure 3 is analogous to Fig. 1 but illustrating the results of B1I and B3I phase ISBs between BDS-3 and BDS-2. Unexpectedly, unlike the code ISBs, we find that, for each baseline considered, using both identical and mixed-receiver pairs, the B1I and B3I phase ISBs are always randomly distributed around a mean value almost zero. In other words, there is no reason to expect the presence of phase ISBs for B1I and B3I between BDS-3 and BDS-2. Figure 4 shows the estimates of code and phase ISBs between BDS-3 B2b and BDS-2 B2I, from which two conclusions can be drawn. First, similar to the B1I and B3I  phase ISBs between BDS-3 and BDS-2, interoperability can be achieved for the phase observations between BDS-3 B2b and BDS-2 B2I. Second, the mean of the code ISBs between BDS-3 B2b and BDS-2 B2I based on APM1-APM3, a mixed-receiver combination, are estimated as nonzero, thereby suggesting the code ISBs should be considered when mixing BDS-3 B2b and BDS-2 B2I for a baseline with nonidentical receiver pairs. Similar to Fig. 2, but for different signals, Fig. 5 shows the code ISBs between BDS-3 B2b and BDS-2 B2I for the short three baselines. Likewise, there is a certain trend in each panel, which we believe is also due to the multipath effect. Overall, it can be confirmed that there are nonzero mean code ISBs between BDS-3 B2b and BDS-2 B2I, which breaks the interoperability between BDS-3 and BDS-2.
To summarize, two conclusions can be drawn. First, there are no ISBs on the phase observations of the three overlapping frequencies of BDS-3 and BDS-2, so they can be treated as one constellation. Second, and more importantly, the code observations of BDS-3 and BDS-2 can be compatible on B3I but not on B1I and B2I/B2b, for a baseline with non-identical receiver pairs. Thus, the differences in code observations of B1I and B2I/B2b between the two constellations must be considered when mixing BDS-3 and BDS-2 observations.

Characterization of BDS-3 and BDS-2 DCB estimates
In this section, the DCB of the common frequencies (B1I-B3I, B1I-B2b/B2I) of BDS-3 and BDS-2 are first analyzed. Then, using B1I as a reference, DCB of the new frequencies of BDS-3 (B1I-B1C and B1I-B2a) will be reported. Figure 6 shows the B1I-B3I DCB of BDS-3 (left) and BDS-2 (right), with each color representing a different baseline. Normally, one would expect DCB of the same frequency combination to be completely consistent for BDS-3 and BDS-2, as they are considered one constellation. The expected results can be seen in the baselines with identical receiver pairs (blue and red lines), where the DCB of B1I-B3I barely differs between BDS-3 and BDS-2. However, the situation becomes different when referring to a mixedreceiver combination. Here, we see there is a clear distinction of B1I-B3I DCB between BDS-3 and BDS-2, with a difference of 1.128 m. In this case, the difference of B1I-B3I DCB between BDS-3 and BDS-2 must be taken into account and ignoring it may adversely affect PNT applications. Moreover, focusing on each panel, we see that these estimates fluctuate randomly around their mean value, with no apparent trend over time, indicating that B1I-B3I DCB has no significant short-term change for each of the above three baselines. Figure 7 is analogous to Fig. 6, except that it shows the DCB of BDS-3 B1I-B2b and BDS-2 B1I-B2I. These results confirm the previous observations that for a baseline   comprising the same receiver type, no difference will occur between the DCB of the overlapping frequency combination of BDS-3 and BDS-2. However, at the same time, the characteristics of DCB between BDS-3 and BDS-2 based on the mixed-receiver combination are interesting. The bottom two panels (yellow lines) show a difference between BDS-3 B1I-B2b and BDS-2 B1I-B2I of 1.535 m, which indicates the DCB difference between BDS-3 B1I-B2b and BDS-2 B1I-B2I should be carefully considered. Figure 8 shows the BDS-3 B1I-B1C and B1I-B2a DCB estimates, from which we can see that the values of DCB are significant, but the intra-day stability is also noticeable. After sidereal filtering, although the multipath effect is obviously weakened, we can still see that the DCB estimation has a certain trend, e.g., see bottom right. This phenomenon can be attributed to the combination of residual multipath effects and other unmodeled errors.
An important conclusion can be drawn that inconsistencies in the characterization of B1I-B3I and B1I-B2b/B2I DCB in BDS-2 and BDS-3 must be fully considered in PNT applications.

Characterization of BDS-3 and BDS-2 DPB estimates
Concerning the DPB, Fig. 9 shows the estimates of B1I-B3I DPB of BDS-3 (left) and BDS-2 (right) for three pairs of , where the datum ambiguities j z 1 A ab,j − 1 z 1 A ab,1 are included. Thus, we can only analyze the fractional part of DPB without obtaining its absolute value. However, from the compatibility of BDS-3 and BDS-2 in phase ISBs, DPB should also be compatible for the two constellations since DPB and phase ISBs have a linear relationship that can be inferred from each other.
We pay attention here to the three cases depicted in Fig. 9, from which two conclusions can be drawn. First, it is found that the B1I-B3I DPB of BDS-3 and BDS-2 are slightly different in magnitude due to the datum ambiguities, but remarkably similar in the trend (compare left and right). This also illustrates the compatibility of DPB between BDS-3 and BDS-2. Second, the short-term temporal variations of B1I-B3I DPB are significant in BDS-3 and BDS-2. See the second case (red lines), except for the obvious intraday variation, there is an obvious jump around 12:00 (see green ellipse). We think this is probably due to a sudden change in the ambient temperature (Zhang et al. 2016;Mi et al. 2020b). The same phenomenon occurs in the third case (yellow lines). From the comparison of the three baselines, it can be concluded that the response of different receivers to ambient temperature is inconsistent.
Similar to Fig. 9, Fig. 10 shows the DPB but for the BDS-3 B1I-B2b and BDS-2 B1I-B2I. These results confirm the previous conclusions that the short-term temporal  variations of DPB must be considered in BDS-3, BDS-2, and their combination. In addition, interestingly, the shortterm variations are present in all three cases, especially the jump at around 12:00. Combined with the DPB characteristics of B1I-B3I, this can be attributed to the fact that the B2I and B2b signals of the receivers involved in the experiment are more sensitive to the environment. Thus, it can be concluded that different frequencies have different response mechanisms to ambient temperature. Figure 11 is analogous to Fig. 8, except it shows the DPB estimates for three baselines. Pay attention to the differences between the receivers used; we can see that the DPB estimates of B1I-B1C and that of B1I-B2a are slightly different in their magnitude and mean values, but noticeably similar in their trend. This also indicates that the receiver-related biases of B1C and B2a of the used receivers have similar mechanisms in response to the environment. Back to Fig. 9, a similar phenomenon also exists in the DPB estimates of B1I-B2b, which means B2b has an environment response similar to that of B1C and B2a.
One important conclusion can be drawn from this analysis. The receiver-related DPB of BDS-3 and BDS-2 may have short-term variations, which can be estimated by analyzing the performance of different receivers and frequencies.

Conclusions
We have presented a method for simultaneously estimating the receiver-related biases, including inter-system biases (ISBs), differential code biases (DCB), and differential phase biases (DPB). This method has the following characteristics, which make it well suited for use in retrieving receiver-related biases. First, an advantage has been taken of not changing reference satellites, thereby enabling the continuity of DPB and phase ISBs estimates. Second, use has been made of a single-differenced (SD) full-rank model. This ensures compatibility of all kinds of cases from singlefrequency single-constellation to multi-frequency multiconstellation data, and more importantly, the reasonable simultaneously estimation of DCB, DPB, and ISBs.
Special care should be taken when using BeiDou navigation satellite system with global coverage (BDS-3) and the BeiDou navigation satellite (regional) system (BDS-2), as they are considered compatible and thus treated as one constellation. With this in mind, we applied the method detailed above to several sets of GNSS data with all five frequencies of BDS-3 and three frequencies of BDS-2, covering a range of receiver types and observation periods. The time-wise estimates of the DCB, DPB, and ISBs between BDS-3 and BDS-2, using all five frequencies available, were presented.
It was experimentally shown that the phase observations of the three overlapping frequencies between BDS-3 and BDS-2 are indeed compatible. In other words, the phase observations of the three overlapping frequencies can be processed as one constellation when mixing BDS-3 and BDS-2. However, when we referred to code observations, the situation got complicated. The code ISBs of B1I and B2b/B2I between BDS-3 and BDS-2 are estimated as nonzero but of B3I are not. That is to say, for code observations, only B3I can achieve full compatibility between BDS-3 and BDS-2, and the differences of B1I and B2b/B2I between BDS-3 and BDS-2 must be taken into account. In addition, care should be taken to the difference of DCB involving B1I and B2I/ B2b between BDS-3 and BDS-2, as they were found to be significant in our experiment. Moreover, interestingly, we found that DPB of BDS-3 and BDS-2 may have significant short-term variations, which are closely related to receiver and frequency.

Data availability
The raw data were provided by the Innovation Academy for Precision Measurement Science and Technology, Chinese Academy of Sciences, Wuhan, China. The raw data used during the Fig. 11 Time series of BDS-3 B1I-B1C (left) and B2I-B2a (right) DPB estimates on DOY 340 of 2020. Three zero baselines are used, including APM1-APM2 with two ALLOY receivers (top), APM3-APM4 with two POLARX5 receivers (middle), APM1-APM3 with ALLOY and POLARX5 receivers (bottom). The numbers in each panel represent the mean and STD in cycles study are available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.