Clustering Code Biases between BDS-2 and BDS-3 Satellites and Effects on Joint Solution

: China's BeiDou navigation satellite system (BDS) has finished global constellation construction and can achieve joint solution, simultaneously relying on the B1I + B3I signals of the BDS2 and BDS-3 satellites. For reasons mostly related to chip shape distortions, navigation satellite observations are corrupted by receiver-dependent code biases. Those biases are brought into observation residuals and degrade the pseudorange correction accuracy. Herein, we present a code bias estimation algorithm, using what we found to be an obvious clustering code bias phenomenon between the BDS-2 and BDS-3 satellites, leading to the systematic biases existing in the BDS-2+3 joint solution. Therefore, we propose a BDS-2+3 joint solution with code bias self-calibration, which can accurately strip off clustering code biases between the BDS-2 and BDS-3 satellites, and can greatly improve precise point positioning (PPP) convergence speed and accuracy. The statistics showed that the residual biases and root mean square (RMS) improved by 36% and 15% and the convergence time improved by approximately 35%. In the convergence stage, the positioning accuracy improved by approximately 38% and 21% in the horizontal and vertical directions, respectively. Meanwhile, in the post-convergence stage, the accuracy improved by approximately 10%.


Introduction
Early ionosphere studies based on Global Navigation Satellite System (GNSS) identified that the measurements of satellite signals are affected by differential code biases (DCBs) [1,2]. When aiming to estimate the total electron content from GNSS data, these code biases must be calibrated so as to not degrade the estimation results [3]. Generally, the satellite-dependent code bias is included in the satellite clock products provided by International GNSS Service (IGS) based on ionospheric-free (IF) combined observations [4], and is generally considered ignorable in ionospheric-free combined precise point positioning (PPP) using IGS's precise products [5]. In addition, broadcast ephemeris provides timing group delay (TGD) correction parameters, which are another representation of satellite-dependent code biases. However, in 2019, Hauschild et al. found that clock offset biases exist in precise clock solutions based on different independent global networks and can reach several nanoseconds, which also introduces an unequal bias into observations [6].
On the receiver side, similar to the effect of multipath signal reflections, the individual chip shape distortions in the ranging signals generated by receivers with different correlators and/or front-end designs lead to different shifts in the correlator's tracking point for each pseudorange number (PRN), which then introduces different inter-satellite code biases (ISCBs) into the observations [6,7]. In GNSS positioning, unequal pseudorange biases cannot be completely absorbed by receiver clock offset and are therefore brought into observation residuals, which greatly reduces pseudorange correction accuracy, and then affects positioning performance. More recent studies have shown that residual pseudorange biases exist between manufacturers with different correlator types and front-end designs [3]. Aiming to enable PPP with ambiguity resolution (AR), fractional carrier-phase biases (FCBs) comprise a more recent branch of research [8][9][10][11]. Code biases are typically lumped into these carrier-phase bias parameters and are thus still relevant [3], and also affect the efficiency of the AR and convergence rate of PPP [12,13].
Hauschild and Montenbruck verified the existence of inter-receiver code biases based on different receiver types for zero baselines and found that the bias inconsistencies are largest for the legacy signals of Russia's GLObal NAvigation Satellite System (GLONASS) due to the inter-channel biases of the frequency-division multiple access (FDMA) signals, but Global Positioning System (GPS) and BeiDou navigation satellite system (BDS) signals are also affected by significant biases [14]. The difference in GLONASS can reach 4-8 m [15]. Notably, the receiver-dependent ISCB of each manufacturer remains stable over a long period of time; the standard deviations (STDs) for GPS and Galileo are at the centimeter level, while that for GLONASS is at the decimeter level [15,16].
In the initial construction phase, Galileo launched four In-Orbit Validation (IOV) satellites, designed as prototypes of Full Operational Configuration (FOC) satellites, to carry out some experiments from 2011 to 2012 [17]. Currently, there are three IOV satellites and 21 FOC satellites in-orbit service (https://www.gsc-europa.eu/system-service-status/constellation-information). Some analyses have found that some manufacturers show larger code bias inconsistencies for Galileo IOV satellites compared the rest of the FOC satellites [6]. For the BeiDou system, China's government finished the construction of the BDS-2 regional navigation system in December 2012 and provided B1I/B2I/B3I civil signals [18,19]. All BDS-2 satellites are based on the DongFangHong-3A satellite platform manufactured by the China Academy of Space Technology (CAST) [20]. From November 2017 to June 2020, the BDS-3 global system launched 24 Medium Earth Orbit (MEO), three Geostationary Earth Orbit (GEO), and three Inclined Geosynchronous Satellite Orbit (IGSO) satellites successfully, which provide B1I/B3I/B1C/B2a/B2b civil signals [21,22]. BDS-3 satellites are based on two different satellite platforms developed by CAST and the Shanghai Engineering Center for Microsatellites (SECM), China Academy of Science. Table 1 shows the BDS satellites' information. In the table, we can be observed that the global users could achieve a joint solution by simultaneously relying on the B1I + B3I signals of the BDS-2 and BDS-3 satellites. Compared with BDS-2, the BDS-3 satellite's signal quality, satellite orbit, and signal-in-space accuracy have been greatly improved [22,23]. Even though B1I and B3I are backward-compatible signals and use the same modulation mode, such as the Galileo IOV and FOC satellites, there are also significant code bias inconsistencies between the BDS-2 and BDS-3 satellites [24].
In an undifferenced (UD) model of a multi-GNSS fusion solution, the GPS observations are usually selected as the reference datum, while the inter-system bias (ISB) parameters actually mean that all computed code biases of the other systems are obtained relative to the biases for the GPS observations [25]. Thus, the code bias inconsistencies of Galileo (IOV and FOC) and BDS (BDS-2 and BDS-3) destroy the understanding of belonging to same system.
Currently, most research focuses on the code bias of the GPS, GLONASS, and Galileo satellites, but the code biases of the BDS satellites have not been analyzed in depth, especially for BDS-3. Aiming to elucidate the receiver-dependent code biases, herein, we present a navigation signal ISCB estimation algorithm and analyze the clustering code bias characteristics of BDS-2 and BDS-3. Then, we propose a BDS-2 and BDS-3 ISCB self-calibration method for single-station and effect analysis of the BDS-2+3 joint solution.

Methodology
For convenience of methodological description, we first introduce the observations considering ISCBs and then design the ISCB estimation algorithm.

Observations
GNSS dual, frequency, ionospheric, free, combined observations considering satellite, and receiver, dependent ISCBs can be expressed as follows [26]: where the superscript s identifies the satellite and the subscript r identifies the receiver; P and Φ represent the ionospheric, free pseudorange and carrier, phase functions corrected for dry tropospheric delay using a model, respectively; N is the ambiguity (in meters); ρ is the traveled geometric distance between the satellite to the receiver for the instances of transmission and reception of the signal; r dt represents the receiver clock errors; the zenith wet delay zwd is estimated from the observations using the wet mapping function M; ε is the observation noise of the respective ionospheric, free function; b and B represent the ionospheric, free, receiver, and satellite, dependent signal delay biases, respectively; s dt represents the satellite clock errors, while 0 s dt represents the clock offset biases. Theoretically, after correcting by the precision clock product, the satellite clock offset residuals should just include the random error terms, and not the deviation terms. Hauschild et al. found that there were different clock offset biases in precise clock solutions based on different independent global networks, which can reach several nanoseconds [6]. Thus, the clock offset bias 0 s dt can be regarded as the clock offset bias relative to datum clock time. In a positioning solution based on precise clock products, the coefficients of the satellite, dependent 0 s dt and receiver, dependent inter, satellite code bias  [12], which are usually involved in ambiguity parameters. The pseudorange biases , , s r P ISCB b are called receiver, dependent ISCBs. As analyzed above, the individual chip shape distortions in ranging signals lead to different shifts in the correlator's tracking point for each PRN and then introduce different ISCBs into the observations [6]. Thus, we added the superscript s to express that ISCB b is related to satellites.

ISCB Estimation Algorithm
From (1) where v represents the residual; . l is the observed minus calculated (OMC) difference from the satellite to the receiver, where the observed part is the observations of the pseudorange and carrier phase, and the calculated part is the geometric distance between the satellite and the receiver calculated by the precise orbit position at signal transmission and station coordinates at signal reception. The ISCB parameters are added to the linearized equations for every epoch as follows: The ISCB parameters are estimated as constant, and a zero constraint is applied to the sum of all ISCBs at one station to separate the ISCBs from every epoch. Then, the ISCBs can be obtained as follows: where S n is the number of satellites involved in each epoch of each system. From Equation (3), it can be found that the estimated ISCB includes not only the receiver, dependent ISCB, but also the clock offset biases in the clock products. Then, ISCB real, time estimation can be realized in single, station mode.
In the ISCB single, station, real, time estimation model, the parameters to be estimated mainly include clock offset, troposphere, the ISCB of each satellite, and ambiguity parameters, as shown in the following equation: According to Equations (2) and (4), the pseudorange and carrier, phase observations are added into the observation equations, and then the t, th epoch linearized observation equation can be expressed as follows: where t l is the OMC matrix: where 1 is the vector with all elements 1; 0 is the matrix with all elements 0; I is the identity matrix; A is the observation coefficient matrix of one station;

Characteristics of BDS, 2 and BDS, 3 ISCB
The BeiDou system has finished global constellation construction and can achieve a joint solution by simultaneously relying on the B1I+B3I signals of the BDS, 2 and BDS, 3 satellites. In this part, we analyze the ISCB characteristics of BDS, 2 and BDS, 3.

Data Preparation
Currently, some of the manufacturers of IGS can track BDS, 3 satellites, such as Septentrio and Trimble. In order to analyze and compare the ISCB characteristics of BDS, 2 and BDS, 3 signals and the relationship with the receiver brand, type, firmware, and antenna, 17 IGS stations of Septentrio and Trimble were selected. The distribution of the stations is shown in Figure 1 and the information of these stations is shown in Table 2. The multi, GNSS precise orbit and clock products provided by the German Research Center for Geoscience (GFZ) were selected [4]. It should be noted that GFZ began to provide BDS, 2 and BDS, 3 precise products based on B1I+B3I signals from December 2019 (day of year-DOY 335); therefore, the data from 1 to 31 December, 2019 (DOY 335-365) were selected as the test arc.  The processing strategy of the code, division multiple access (CDMA) signal ISCB estimation is shown in Table 3.

Clustering Biases between BDS, 2 and BDS, 3
The ISCBs of those stations listed in Table 2 were estimated and obtained every day following the processing strategy in Table 3. In order to analyze the characteristics of the BDS, 2 and BDS, 3 ISCBs, we statistically calculated the average and stability of all of the BDS satellites' ISCBs from DOY 335 to 365 of 2019 following (8).
where n is the number of ISCBs of the satellite s, AVE is the average, and STD is the standard deviation. Figure 2 shows the average and stability of all of the BDS, 2 and BDS, 3 satellites' ISCBs of Septentrio and Trimble in Figure 1 from DOY 335 to 365 of 2019, where the horizontal axis refers to the BDS, 2 and BDS, 3 satellites, the different marking lines represent the different receivers, and the error bars represent the ISCB STD values for those receivers. Based on Figure 2, it can be concluded that the ISCBs have a good consistency for the Septentrio devices, while the differences of the Trimble brand are relatively large and can be up to 6-8m, i.e., BRST and MCHL. It can be also found in Figure 2 that there is an obvious clustering phenomenon for the BDS, 2 and BDS, 3 satellites for all Septentrio receivers and some of the Trimble ones. In an undifferenced model of a multi, GNSS fusion solution, the GPS observations are usually selected as the reference datum, while the ISB parameters actually mean that all of the computed code biases of the other systems are obtained relative to the biases for the GPS observations [25]. Thus, one receiver's systematic bias of two clusters of BDS, 2 and BDS, 3 ISCBs are actually the ISB and can be calculated by the average of the BDS, 2 ISCBs and the BDS, 3 ISCBs, as per Equation (9). Table 4 provides the statistics of the receivers of Septentrio and Trimble in Figure 1 and Table 2, including the average ISCB, the ISB between BDS, 2 and BDS, 3, and the STD. From the statistics provided by Table 4, the clustering code bias phenomenon become clearer. The ISCB cluster difference of the Septentrio devices were close, and the ISB of BDS2 and BDS, 3 was approximately 1.5 m. Meanwhile, there were significant differences among the Trimble devices, especially BRST and MCHL. The MCHL ISB of BDS2 and BDS, 3 was up to 6.3 m, which would seriously affect the positioning performance of the BDS, 2+3 joint solution based on B1I+B3I. From Table 4, we can also conclude that the STD of the BDS ISCB remained at decimeter level, worse than the centimeter level of GPS and Galileo [16], and the ISCBs of the BDS, 3 satellite were more stable than those of the BDS, 2 satellite. Aiming to identify the significant clustering code biases between BDS, 2 and BDS, 3, Figure 3 shows all of the BDS satellites' ISCB time series of the part receivers of Septentrio (NKLG and REDU) and Trimble (BRST and MCHL) from DOY 335 to 365 of 2019. Among them, the blue, tone lines are the satellites of BDS, 2, while the red, tone lines are the satellites of BDS, 3. The two different color clusters show the clustering bias characteristics of BDS, 2 and BDS, 3 more clearly. Combining Figures 2 and 3 and Table 4, it can be observed that there are some differences in the different manufacturers, receiver types, firmware, and antennas in terms of the BDS code bias characteristics. As Septentrio showed, some manufacturers have good consistency among devices, while some are irregular, such as Trimble. In order to verify the universality of the clustering code biases between BDS, 2 and BDS, 3, approximately 90 global BDS, 2 and BDS, 3 tracking stations were also selected, as shown in Figure 4. Figure 5 shows the distribution of these stations' ISCBs for BDS, 2 and BDS, 3, as estimated with the additional constraint (4). Similar to Figure 3, the blue, tone marks in Figure 5 are the satellites of BDS, 2, while the red, tone marks are the BDS, 3 satellites. The results show the ISCB clustering phenomenon between BDS, 2 and BDS, 3 is universal in the individual stations. As per the analysis above, the differences in the average ISCB of the BDS, 2 and BDS, 3 satellites are also the ISB parameters estimated in the undifferenced model; therefore, systematic clustering biases should be considered in the BDS, 2+3 joint solution.
In combination with Figure 3, it is interesting to note that the code biases between BDS, 2 and BDS, 3 are opposite for the Septentrio and Trimble stations. The inconsistency in the phenomenon of inter, manufacturer code biases also exists in Figure 4, which was caused by manufacturers with different correlator types and front, end designs [3].   Figure 6 shows the time series of the BDS, 2 and BDS, 3 ISCBs for the MCHL station in the first three hours. As with the ambiguities for the carrier, phase observations, there was a significant convergence period (about 0.5 h) with the accumulation of observations, after which all of the ISCBs stabilized. Thus, these ISCBs can be applied to the self, calibrating code biases of the subsequent epoch.

Systematic Biases in the BDS, 2+3 Joint Solution
Those results above show that there is an obvious clustering code bias phenomenon between the BDS, 2 and BDS, 3 satellites, which leads to systematic biases existing in the BDS, 2+3 joint solution. Thus, based on the B1I+B3I signals, BDS, 2 and BDS, 3 should be regarded as two individual navigation satellite systems. According to multi, GNSS positioning theory, the ISB parameters should be added into observations [25,31].
By introducing ISB parameters into observations, they can absorb the shared systematic biases between BDS, 2 and BDS, 3, but cannot describe the dispersed clustering code biases between the BDS, 2 and BDS, 3 satellites accurately, which leads to unequal code biases in residuals, while individual ISCB parameters can avoid this. On account of relying on its own observations and ephemeris products with no other external products, the ISCB estimation algorithm above can realize ISCB self, calibration at a single station independently. Based on the stable ISCB, we propose the BDS, 2+3 joint solution with code bias self, calibration, using which we analyzed the residuals and positioning improvements.

BDS, 2+3 Joint Model with Code Bias Self-Calibration
Considering the existence of unequal receiver, dependent clustering code biases between BDS, 2 and BDS, 3, the linearized undifferenced observation (1) of the BDS, 2+3 joint solution is as follows: where C2 and C3 represent the BDS, 2 and BDS, 3 satellites, respectively; i and j are the i, th BDS, 2 satellite and j, th BDS, 3 at the t, th epoch; u is the unit vector from the receiver to the satellite direction; dx represents the receiver position incremental vector relative to the previous epoch.
Similarly to the ISCB estimation process, in the BDS, 2+3 joint model, considering clustering code biases, the parameters to be estimated mainly include position, clock offset, troposphere, the ISCB of each satellite, and the ambiguity parameters. The analyses above show that the receiver, dependent ISCB remains stable over a long time. Thus, the ISCBs estimated in the previous epoch can be applied to self, calibrating the clustering code biases of the subsequent epoch to achieve the synchronization of BDS, 2 and BDS, 3. Thus, (10) can be further simplified as follows: , . ISCB can be obtained from (3).
From (11), the BDS, 2+3 joint model with code bias self, calibration just estimates the general parameters, such as the position, clock offset, troposphere, and ambiguity parameters, as follows: The pseudorange and carrier, phase observations are added into the observation equations, and then the t, th epoch linearized observation equation can be expressed as follows: where: where ISCB is the vector of the BDS, 2 and BDS, 3 satellites ISCBs: is the unit vector of each satellite in the position direction From (11) and (14), the ISCBs obtained in previous epoch are used to correct the OMC of the BDS, 2 and BDS, 3 satellites, which not only resolves the systematic biases existing in BDS, 2+3, but also improves the pseudorange correction accuracy. Theoretically, with the improvement of the pseudorange correction accuracy, positioning performance will also be elevated.

Residual Analysis
In order to analyze the effects of ISCB self, calibration on the BDS, 2+3 joint solution, we analyzed the pseudorange residuals based on the PPP results of the stations in Figure  4 from DOY 359 to 365 of 2019. Figure 7 shows a comparison of the PPP pseudorange residuals before and after ISCB self, calibration, where the horizontal axis refers to the BDS, 2 and BDS, 3 satellites, and the different marks in upper two sub, figures represent the residual distribution of all of the receivers in Figure 4 for one satellite. The blue, tone marks are the residuals with no ISCB self, calibration, while the red, tone marks are the residuals after ISCB self, calibration. The histograms of the third and fourth subfigures are the root mean square (RMS) and STD statistics of the pseudorange residuals with and without ISCB correction. The RMS and STD can be obtained as follows: where sta n refers to the stations number of the residuals of the satellite s, AVE is the residuals average, STD is the standard deviation, and RMS is the root mean square.
From the results, compared to the uncalibrated BDS, 2+3 joint PPP, the ISCB self, calibration can further correct the pseudorange residuals to close to 0. The statistics show that the biases and RMS improved by 36% and 15%, respectively. However, the STD did not change significantly. It was found from Equation (3) that the multi, GNSS ISCB real, time estimation algorithm for single stations proposed in this paper can accurately strip off the clustering code biases of BDS, 2 and BDS, 3 and the clock offset biases from pseudorange observations. We interrupted the BDS, 2+3 joint PPP experiment every 2 h to simulate the reconvergence process and used the ISCB estimated in the previous two hours to self, calibrate the pseudoranges of the subsequent epoch. Then, we analyzed the influence of the ISCB self, calibration on the BDS, 2+3 joint PPP. Figure 8 shows one station's horizontal and vertical convergence comparison of the BDS, 2+3 joint PPP with and without ISCB self, calibration. The black line indicates the PPP results with no ISCB self, calibration, while the red line shows the PPP results after ISCB self, calibration. It can be seen from the figure that the convergence speed of the BDS, 2+3 joint PPP self, calibrated by ISCB significantly improved.  Figure 6, the 68% accuracy of each epoch in the convergence stage was analyzed, and a comparison of the accuracy is shown in Figure 9 (convergence condition was 10 cm in the horizontal direction and 20 cm in the vertical direction). Table 5 shows convergence comparisons and improvements in the different convergence conditions of BDS, 2+3 joint PPP with and without ISCB self, calibration. It can be seen from Figure 9 and Table 5 that ISCB self, calibration improved the convergence speed and accuracy of BDS, 2+3 joint PPP greatly. The convergence time to 10 cm in the horizontal direction reduced from 133 min to approximately 100 min, while the convergence time to 20 cm in the vertical reduced from 64 min to approximately 48.5 min. As can be seen in Table 5, the convergence time improved by approximately 42% and 28% in the horizontal and vertical directions, respectively. In addition, ISCB self, calibration significantly improved the BDS, 2+3 joint PPP accuracy in the convergence and post, convergence stages. In the convergence stage, the accuracy improved by approximately 38% and 21% in the horizontal and vertical directions, respectively. Meanwhile, in the post, convergence stage, the accuracy improved by approximately 10% in both the horizontal and vertical directions.

Model
Convergence Condition Average 40 cm 30 cm 20 cm 10 cm 3. We proposed the BDS, 2+3 joint solution with code bias self, calibration, which can accurately strip off clustering code biases between the BDS, 2 and BDS, 3 satellites and can greatly improve the PPP convergence speed and accuracy. 4. The statistics showed that the residual biases and RMS of BDS, 2+3 joint PPP improved by 36% and 15%, respectively, and the convergence time improved by approximately 35%. In the convergence stage, the positioning accuracy improved by approximately 38% and 21% in the horizontal and vertical directions, respectively. Meanwhile, in the post, convergence stage, the accuracy improved by approximately 10%.