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

Abstract

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 BDS-2 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%.

1. 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. The BeiDou navigation satellite system (BDS) satellite’s information (www.csno-tarc.cn/en/system/constellation, up to June 2020).

PRN	SVN	NORAD ID	SVN	System	Manufacture	Notes(UTC)	Civil Signal
C01	C020	44231	GEO-8	BDS-2	CAST	17 May, 2019	B1I/B2I/B3I
C02	C016	38953	GEO-6	BDS-2	CAST	25 Oct, 2012	B1I/B2I/B3I
C03	C018	41586	GEO-7	BDS-2	CAST	12 June, 2016	B1I/B2I/B3I
C04	C006	37210	GEO-4	BDS-2	CAST	1 Nov, 2010	B1I/B2I/B3I
C05	C011	38091	GEO-5	BDS-2	CAST	25, Feb, 2012	B1I/B2I/B3I
C06	C005	36828	IGSO, 1	BDS, 2	CAST	1, Aug, 2010	B1I/B2I/B3I
C07	C007	37256	IGSO, 2	BDS, 2	CAST	18, Dec, 2010	B1I/B2I/B3I
C08	C008	37384	IGSO, 3	BDS, 2	CAST	10, Apr, 2011	B1I/B2I/B3I
C09	C009	37763	IGSO, 4	BDS, 2	CAST	27, Jul, 2011	B1I/B2I/B3I
C10	C010	37948	IGSO, 5	BDS, 2	CAST	2, Dec, 2011	B1I/B2I/B3I
C11	C012	38250	MEO, 3	BDS, 2	CAST	30, Apr, 2012	B1I/B2I/B3I
C12	C013	38251	MEO, 4	BDS, 2	CAST	30, Apr, 2012	B1I/B2I/B3I
C13	C017	41434	IGSO, 6	BDS, 2	CAST	30, Mar, 2016	B1I/B2I/B3I
C14	C015	38775	MEO, 6	BDS, 2	CAST	19, Sep, 2012	B1I/B2I/B3I
C16	C019	43539	IGSO, 7	BDS, 2	CAST	10, Jul, 2018	B1I/B2I/B3I
C19	C201	43001	MEO, 1	BDS, 3	CAST	5, Nov, 2017	B1I/B3I/B1C/B2a/B2b
C20	C202	43002	MEO, 2	BDS, 3	CAST	5, Nov, 2017	B1I/B3I/B1C/B2a/B2b
C21	C206	43208	MEO, 3	BDS, 3	CAST	12, Feb, 2018	B1I/B3I/B1C/B2a/B2b
C22	C205	43207	MEO, 4	BDS, 3	CAST	12, Feb, 2018	B1I/B3I/B1C/B2a/B2b
C23	C209	43581	MEO, 5	BDS, 3	CAST	29, Jul, 2018	B1I/B3I/B1C/B2a/B2b
C24	C210	43582	MEO, 6	BDS, 3	CAST	29, Jul, 2018	B1I/B3I/B1C/B2a/B2b
C25	C212	43603	MEO, 11	BDS, 3	SECM	24, Aug, 2018	B1I/B3I/B1C/B2a/B2b
C26	C211	43602	MEO, 12	BDS, 3	SECM	24, Aug, 2018	B1I/B3I/B1C/B2a/B2b
C27	C203	43107	MEO, 7	BDS, 3	SECM	12, Jan, 2018	B1I/B3I/B1C/B2a/B2b
C28	C204	43108	MEO, 8	BDS, 3	SECM	12, Jan, 2018	B1I/B3I/B1C/B2a/B2b
C29	C207	43245	MEO, 9	BDS, 3	SECM	29, Mar, 2018	B1I/B3I/B1C/B2a/B2b
C30	C208	43246	MEO, 10	BDS, 3	SECM	29, Mar, 2018	B1I/B3I/B1C/B2a/B2b
C32	C213	43622	MEO, 13	BDS, 3	CAST	19, Sep, 2018	B1I/B3I/B1C/B2a/B2b
C33	C214	43623	MEO, 14	BDS, 3	CAST	19, Sep, 2018	B1I/B3I/B1C/B2a/B2b
C34	C216	43648	MEO, 15	BDS, 3	SECM	15, Oct, 2018	B1I/B3I/B1C/B2a/B2b
C35	C215	43647	MEO, 16	BDS, 3	SECM	15, Oct, 2018	B1I/B3I/B1C/B2a/B2b
C36	C218	43706	MEO, 17	BDS, 3	CAST	18, Nov, 2018	B1I/B3I/B1C/B2a/B2b
C37	C219	43707	MEO, 18	BDS, 3	CAST	18, Nov, 2018	B1I/B3I/B1C/B2a/B2b
C38	C220	44204	IGSO, 1	BDS, 3	CAST	20, Apr, 2019	B1I/B3I/B1C/B2a/B2b
C39	C221	44337	IGSO, 2	BDS, 3	CAST	24, Jun, 2019	B1I/B3I/B1C/B2a/B2b
C40	C224	44709	IGSO, 3	BDS, 3	CAST	4, Nov, 2019	B1I/B3I/B1C/B2a/B2b
C41	C227	44864	MEO, 19	BDS, 3	CAST	16, Dec, 2019	B1I/B3I/B1C/B2a/B2b
C42	C228	44865	MEO, 20	BDS, 3	CAST	16, Dec, 2019	B1I/B3I/B1C/B2a/B2b
C43	C226	44794	MEO, 21	BDS, 3	SECM	23, Nov, 2019	B1I/B3I/B1C/B2a/B2b
C44	C225	44793	MEO, 22	BDS, 3	SECM	23, Nov, 2019	B1I/B3I/B1C/B2a/B2b
C45	C223	44543	MEO, 23	BDS, 3	CAST	22, Sep, 2019	B1I/B3I/B1C/B2a/B2b
C46	C222	44542	MEO, 24	BDS, 3	CAST	22, Sep, 2019	B1I/B3I/B1C/B2a/B2b
C59	C217	43683	GEO, 1	BDS, 3	CAST	1, Nov, 2018	B1I/B3I
C60	C229	45344	GEO, 2	BDS, 3	CAST	9, Mar, 2020	B1I/B3I
C61	C230	45807	GEO, 3	BDS, 3	CAST	23, Jun, 2020	B1I/B3I

PRN, Pseudo, Random Noise; SVN, space vehicle number; GEO, Geostationary Earth Orbit; IGSO, Inclined Geosynchronous Satellite Orbit; MEO, Medium Earth Orbit; CAST, China Academy of Space Technology (CAST); SECM, Shanghai Engineering Center for Microsatellites.

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.

PRN, Pseudo, Random Noise; SVN, space vehicle number; GEO, Geostationary Earth Orbit; IGSO, Inclined Geosynchronous Satellite Orbit; MEO, Medium Earth Orbit; CAST, China Academy of Space Technology (CAST); SECM, Shanghai Engineering Center for Microsatellites.

2. Methodology

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

2.1. Observations

GNSS dual, frequency, ionospheric, free, combined observations considering satellite, and receiver, dependent ISCBs can be expressed as follows [26]:

$$\left\{\begin{array}{l}{P}_{r,IF}^s={\rho }_r^s+d{t}_r+b_{r,P,ISCB}^s-\left(dt{0}^s+d{t}_{}^s\right)-B_P^s+M_r^s\cdot zwd+{\epsilon }_P^s\\ {\mathsf{\Phi }}_{r,IF}^s={\rho }_r^s+d{t}_r+b_{r,\mathsf{\Phi }}^s-\left(dt{0}^s+d{t}_{}^s\right)-B_{\mathsf{\Phi }}^s+N_r^s+M_r^s\cdot zwd+{\epsilon }_{\mathsf{\Phi }}^{}\end{array}\right.$$

(1)

where the superscript s identifies the satellite and the subscript r identifies the receiver; P and $\mathsf{\Phi }$ represent the ionospheric, free pseudorange and carrier, phase functions corrected for dry tropospheric delay using a model, respectively; N is the ambiguity (in meters); $\rho$ is the traveled geometric distance between the satellite to the receiver for the instances of transmission and reception of the signal; $d{t}_r$ represents the receiver clock errors; the zenith wet delay zwd is estimated from the observations using the wet mapping function M; $\epsilon$ 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; $d{t}^s$ represents the satellite clock errors, while $dt{0}^s$ 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 $dt{0}^s$ 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 $dt{0}^s$ and receiver, dependent inter, satellite code bias $b_{r,P,ISCB}^s$ in (1) are the same, which makes it difficult to separate the two parameters if there is no extra constraint. The carrier, phase biases, $b_{r,\mathsf{\Phi }}^s$ and $B_{\mathsf{\Phi }}^s$, are generally called uncalibrated phase delays (UPDs) [12], which are usually involved in ambiguity parameters. The pseudorange biases $b_{r,P,ISCB}^s$ 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.

2.2. ISCB Estimation Algorithm

From (1), the coefficients of the satellite, dependent delay errors, such as $B_P^s$, $dt{0}^s$, and $d{t}^s$, are the same. Among them, $B_P^s$, the satellite, dependent code bias, is relatively stable [27], which is actually included in satellite clock products provided by the International GNSS Service (IGS) based on ionospheric, free, combined observations [5], i.e., $d{\tilde{t}}^s=d{t}^s+B_P^s$. In addition, broadcast ephemeris provides timing group delay (TGD) correction parameters, which are another representation of the satellite, dependent code bias. $dt{0}^s$ can be regarded as the clock offset bias relative to datum clock time. In a positioning solution, $dt{0}^s$ can be absorbed by receiver, dependent ISCB parameters. Considering known fixed station coordinates, the linearized undifferenced observation (1) of the single, station, multi, GNSS, real, time ISCB estimation model with the same datum is as follows:

$$\begin{array}{l}{v}_P^s=d{t}_r+\left(b_{r,P,ISCB}^s-dt{0}_{}^s\right)+M_r^s\cdot zwd-l_P^s\\ v_{\mathsf{\Phi }}^s=d{t}_r+\left(b_{r,P,ISCB}^s-dt{0}_{}^s\right)+{\stackrel{⌢}{N}}_r^s+M_r^s\cdot zwd-l_{\mathsf{\Phi }}^s\end{array}$$

(2)

where v represents the residual; ${\stackrel{⌢}{N}}_r^s=N_r^s+\left(b_{\mathsf{\Phi }}^s-b_{r,P,ISCB}^s\right)-\left(B_{\mathsf{\Phi }}^s-B_P^s\right)$. 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:

$$ISC{B}_{}^s=b_{r,P,ISCB}^s-dt{0}_{}^s$$

(3)

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:

$${\displaystyle {\sum }_{i=1}^{n_S}ISC{B}_{}^{s_i}}=0$$

(4)

where $n_S$ 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:

$$X=\left[\begin{array}{c}d{t}_{1\times 1}\\ zw{d}_{1\times 1}\\ ISC{B}_{n_S\times 1}\\ N_{n_S\times 1}\end{array}\right]$$

(5)

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:

$$\left[\begin{array}{l}{v}_t\\ 0\end{array}\right]=A_t\cdot X_t+\left[\begin{array}{l}{l}_t\\ 0\end{array}\right]$$

(6)

where $l_t$ is the OMC matrix:

$$A=\left[\begin{array}{cccc}\left[1_{2n_S\times 1}\right]& M_{2n_S\times 1}& \left[\begin{array}{c}{I}_{n_S\times n_S}\\ I_{n_S\times n_S}\end{array}\right]& \left[\begin{array}{c}{0}_{n_S\times n_S}\\ I_{n_S\times n_S}\end{array}\right]\\ 0& 0& \left[1_{1\times n_S}\right]& \left[0_{1\times n_S}\right]\end{array}\right]$$

(7)

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; $M_{2n_S\times 1}$ is the projection coefficient matrix of the troposphere parameter.

3. 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.

3.1. 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. Experimental site information (October to December 2019).

Manufacturer	Type	Firmware	Antenna	Site	BDS Signal
Septentrio	POLARX5TR	5.3.0	SEPCHOKE_B3E6	CEBR	B1I/B2I/B3I
	POLARX5	5.3.0	SEPCHOKE_B3E6	KIRU, KOUR, REDU
			TRM59800.00	NKLG
Trimble	ALLOY	5.37 5.37	TRM57971.00	BRST	B1I/B2I/B3I
			TRM59800.00	MCHL
		5.42	TRM57971.00	UNB3
		5.43	LEIAR25.R3	KIR8, MAR7
	NETR9	5.42	TRM59800.00	METG
			TRM115000.00	POAL, POVE, SALU, SAVO, TOPL, UFPR

. 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. Processing strategy, model, and parameters of the inter, satellite code bias (ISCB) estimation.

Parameters		Model
Observations		UD IF pseudorange/carrier phase
Prior information		Pseudorange 1.0 m; carrier phase 0.02 cycles
Cut, off elevation		7°
Observation weights		p = 1, elev. >30° p = 2sin(elev.), elev. ≤30°
Interval		30s
PCO/PCV	Satellite	GPS, Galileo PCO: IGS14.atx [28]; BDS, 2 GEO PCO: IGS M, GEX; BDS, 2 IGSO/MEO PCO: ESA Mode [29]; BDS, 3 MEO: www.beidou.gov.cn GPS PCV: IGS14.atx; BDS, Galileo PCV: Uncorrected
	Receiver	GPS PCO: IGS14.atx; BDS, Galileo PCO: Same as GPS; GPS PCV: IGS14.atx; BDS, Galileo PCV: Same as GPS
Adjustment method		Square root information filtering [30]
Troposphere delay		Saastamoinen model + GMF mapping function random, walk process for each epoch
Receiver clock offset		Estimated as white noise
Ambiguity		Float
ISCB		Constant estimation with zero, mean condition every day [25]

UD, Undifferenced; IF, Ionospheric-Free; PCO, Phase Center Offset; PCV, Phase Center Variation; ESA, Europe Space Agency; ISCB, Inter-Satellite Code Bias.

.

UD, Undifferenced; IF, Ionospheric-Free; PCO, Phase Center Offset; PCV, Phase Center Variation; ESA, Europe Space Agency; ISCB, Inter-Satellite Code Bias.

3.2. 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).

$$\begin{array}{l}AV{E}^s=\frac{{\displaystyle \sum _{i=1}^n\left(ISC{B}_i^s\right)}}{n}\\ STD=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(ISC{B}_i^s-AV{E}^s\right)}^2}}{n-1}}\end{array}$$

(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).

$$IS{B}^{C2-C3}=\frac{{\displaystyle \sum _{i=1}^{n_{C2}}\left(ISC{B}_i^{C2}\right)}}{n_{C2}}-\frac{{\displaystyle \sum _{i=1}^{n_{C3}}\left(ISC{B}_i^{C3}\right)}}{n_{C3}}$$

(9)

Table 4. The average ISCBs of BDS, 2 and BDS, 3 and the inter, system bias (ISB) and standard deviation (STD) of the Septentrio and Trimble receivers (m).

Manufacturer	Station	ISCB Average		ISB between BDS, 2 and BDS, 3	ISCB STD
		BDS, 2	BDS, 3		BDS, 2	BDS, 3
Septentrio	CEBR	0.906	−0.476	1.382	0.477	0.191
	KIRU	0.853	−0.534	1.387	0.210	0.153
	KOUR	1.267	−0.248	1.515	0.417	0.197
	REDU	0.905	−0.612	1.517	0.383	0.225
	NKLG	1.105	−0.515	1.620	0.523	0.191
Average		1.007	−0.477	1.484	0.402	0.191
Trimble	BRST	−1.664	1.484	−3.148	0.217	0.340
	MCHL	−2.979	3.276	−6.255	0.232	0.202
	UNB3	−0.766	0.142	−0.907	0.325	0.163
	KIR8	−0.168	0.091	−0.259	0.208	0.135
	MAR7	−0.106	0.069	−0.176	0.245	0.166
	METG	0.497	−0.467	0.964	0.233	0.176
	POAL	−0.885	0.223	−1.108	0.418	0.255
	POVE	−0.686	0.172	−0.858	0.416	0.289
	SALU	−0.669	0.167	−0.836	0.551	0.427
	SAVO	−0.339	0.085	−0.424	0.478	0.244
	TOPL	−0.633	0.158	−0.792	0.408	0.211
	UFPR	−0.683	0.171	−0.854	0.429	0.225
Average		−0.757	0.464	−1.221	0.347	0.236

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 Figure 2 and Figure 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].

3.3. ISCB Time Variation Characteristics

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.

4. Results and Discussion

4.1. 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.

4.2. 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:

$$\left\{\begin{array}{l}{v}_{i,P}^{s,C2}\left(t\right)=u_i^{s,C2}\cdot dx+d{t}_r\left(t\right)+\left(b_{i,P,ISCB}^{s,C2}-dt{0}_i^{s,C2}\right)+M_i^{s,C2}\cdot zwd\left(t\right)-l_{i,P}^{s,C2}\left(t\right)\\ v_{i,\mathsf{\Phi }}^{s,C2}\left(t\right)=u_i^{s,C2}\cdot dx+d{t}_r\left(t\right)-dt{0}_i^{s,C2}+b_{i,\mathsf{\Phi }}^{s,C2}-B_{j,\mathsf{\Phi }}^{s,C2}+N_i^{s,C2}+M_i^{s,C2}\cdot zwd\left(t\right)-l_{i,\mathsf{\Phi }}^{s,C2}\left(t\right)\\ v_{j,P}^{s,C3}\left(t\right)=u_j^{s,C2}\cdot dx+d{t}_r\left(t\right)+\left(b_{j,P,ISCB}^{s,C3}-dt{0}_j^{s,C3}\right)+M_j^{s,C3}\cdot zwd\left(t\right)-l_{j,P}^{s,C3}\left(t\right)\\ v_{j,\mathsf{\Phi }}^{s,C3}\left(t\right)=u_j^{s,C2}\cdot dx+d{t}_r\left(t\right)-dt{0}_j^{s,C3}+b_{j,\mathsf{\Phi }}^{s,C3}-B_{j,\mathsf{\Phi }}^{s,C3}+N_j^{s,C3}+M_j^{s,C3}\cdot zwd\left(t\right)-l_{j,\mathsf{\Phi }}^{s,C3}\left(t\right)\end{array}\right.$$

(10)

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:

$$\left\{\begin{array}{l}{v}_{i,P}^{s,C2}\left(t\right)=u_i^{s,C2}\cdot dx+d{t}_r\left(t\right)+M_i^{s,C2}\cdot zwd\left(t\right)-{\tilde{l}}_{i,P}^{s,C2}\left(t\right)\\ v_{i,\mathsf{\Phi }}^{s,C2}\left(t\right)=u_i^{s,C2}\cdot dx+d{t}_r\left(t\right)+{\tilde{N}}_i^{s,C2}+M_i^{s,C2}\cdot zwd\left(t\right)-l_{i,\mathsf{\Phi }}^{s,C2}\left(t\right)\\ v_{j,P}^{s,C3}\left(t\right)=u_j^{s,C2}\cdot dx+d{t}_r\left(t\right)+M_j^{s,C3}\cdot zwd\left(t\right)-{\tilde{l}}_{j,P}^{s,C3}\left(t\right)\\ v_{j,\mathsf{\Phi }}^{s,C3}\left(t\right)=u_j^{s,C2}\cdot dx+d{t}_r\left(t\right)+{\tilde{N}}_j^{s,C3}+M_j^{s,C3}\cdot zwd\left(t\right)-l_{j,\mathsf{\Phi }}^{s,C3}\left(t\right)\end{array}\right.$$

(11)

where ${\tilde{l}}_P^s=l_P^s+ISC{B}^s$, ${\tilde{N}}_i^s=N_i^s-dt{0}_i^s+b_{i,\mathsf{\Phi }}^s-B_{\mathsf{\Phi }}^s$. 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:

$$X=\left[\begin{array}{c}d{x}_{3\times 1}\\ d{t}_{1\times 1}\\ zw{d}_{1\times 1}\\ N_{n_S\times 1}\end{array}\right]$$

(12)

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:

$$v_t=A_t\cdot X_t+{\tilde{l}}_t$$

(13)

where:

$${\tilde{l}}_t=l_t+ISCB$$

(14)

where $ISCB$ is the vector of the BDS, 2 and BDS, 3 satellites ISCBs:

$$A=\left[\begin{array}{cccc}\left[B_{2n_S\times 3}\right]& \left[1_{2n_S\times 1}\right]& M_{2n_S\times 1}& \left[\begin{array}{c}{0}_{n_S\times n_S}\\ I_{n_S\times n_S}\end{array}\right]\end{array}\right]$$

(15)

where $B_{2n_S\times 3}$ is the unit vector of each satellite in the position direction $d{x}_{3\times 1}$.

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.

4.3. 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:

$$\begin{array}{l}AV{E}^s=\frac{{\displaystyle \sum _{i=1}^{n_{sta}}\left(RE{S}_i^s\right)}}{n_{sta}}\\ ST{D}^s=\sqrt{\frac{{\displaystyle \sum _{i=1}^{n_{sta}}{\left(RE{S}_i^s-AV{E}^s\right)}^2}}{n_{sta}-1}}\\ RM{S}^s=\sqrt{\frac{{\displaystyle \sum _{i=1}^{n_{sta}}{\left(RE{S}_i^s\right)}^2}}{n_{sta}}}\end{array}$$

(16)

where $n_{sta}$ 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.

4.4. Improvements in BDS, 2+3 Joint PPP

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.

Based on the PPP results of DOY 359–365 of 2019 of the stations in 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. Convergence comparisons and improvements in the different convergence conditions of BDS, 2+3 joint PPP with and without ISCB self, calibration.

Model		Convergence Condition				Average
		40 cm	30 cm	20 cm	10 cm
Convergence time in the horizontal direction (min)	BDS	52	68	91.5	133	/
	BDS+ISCB	18.5	38	58.5	99.5	/
	Improvement rate	64%	44%	36%	25%	42%
Convergence time in the vertical direction (min)	BDS	31.5	40	64	/	/
	BDS+ISCB	18.5	32.5	48.5	/	/
	Improvement rate	41%	19%	24%	/	28%

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.

5. Conclusions

For reasons mostly related to chip shape distortions, GNSS observations are corrupted by receiver, dependent code biases. Receiver, dependent ISCB and clock offset biases are brought into observation residuals, which degrades the pseudorange correction accuracy. Aiming to determine the receiver, dependent code biases, we presented a navigation signal ISCB estimation algorithm and analyzed the clustering code bias characteristics of BDS, 2 and BDS, 3. Then, we proposed a BDS, 2 and BDS, 3 ISCB self, calibration method for single, station and effects analysis of the BDS, 2+3 joint solution. The results were as follows:

• The ISCB real, time estimation algorithm for single stations presented in this paper accurately stripped off the receiver, dependent ISCB and clock offset biases from the pseudorange observations simultaneously.
• We analyzed the ISCB characteristics of BDS, 2 and BDS, 3 based on the B1I + B3I signal and found there to be an obvious clustering code bias phenomenon between the BDS, 2 and BDS, 3 satellites, leading to systematic biases existing in the BDS, 2+3 joint solution.
• 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.
• 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%.