Evaluation of Inter-System Bias between BDS-2 and BDS-3 Satellites and Its Impact on Precise Point Positioning

Abstract

The BeiDou navigation satellite system (BDS) currently has 41 satellites in orbits and will reach its full constellation following the launch of the last BDS satellite in June 2020 to provide navigation, positioning, and timing (PNT) services for global users. In this contribution, we investigate the characteristics of inter-system bias (ISB) between BDS-2 and BDS-3 and verify whether an additional ISB parameter should be introduced for the BDS-2 and BDS-3 precise point positioning (PPP). The results reveal that because of different clock references applied for BDS-2 and BDS-3 in the International GNSS Service (IGS) precise satellites clock products and the inconsistent code hardware delays of BDS-2 and BDS-3 for some receiver types, an ISB parameter needs to be introduced for BDS-2 and BDS-3 PPP. Further, the results show that the ISB can be regarded as a constant within a day, the value of which is closely related to the receiver type. The ISB values of the stations with the same receiver type are similar to each other, but a great difference may be presented for different receiver types, up to several meters. In addition, the impact of ISB on PPP has also been studied, which demonstrates that the performance of kinematic PPP could be improved when ISB is introduced.

1. Introduction

The Beidou navigation satellite system (BDS) is a global navigation satellite system (GNSS) developed in China. According to the plan, the construction of BDS is divided into three stages [1]. The first stage is called BDS-1, to build a demonstration system; the second is to provide positioning, navigation, and timing services (PNT) for the Asian-Pacific region, called BDS-2; the construction of BDS-3 at the final stage that begun from 2015 has provided the initial global service since 27, December, 2018, which consists of three geostationary earth orbit (GEO) satellites, three inclined geosynchronous orbit (IGSO) satellites, and 24 medium earth orbit (MEO) satellites [2]. After the launch of the last BDS-3 satellite planned in June 2020, the BDS will reach its full constellation to provide a global service.

Many studies on the performance of BDS-2 have been conducted, including signal characteristics of the code observations [3,4,5], satellite visibility and geometry [6,7], code pseudorange variations [8,9,10,11], accuracy assessment of the broadcast ephemerides [12], precise orbit determination (POD) and precise clock estimation (PCE) performance [13,14,15,16,17], performance of the precise point positioning (PPP) [7,18,19], and PPP ambiguity resolution (AR) [20,21,22].

In recent years, many studies have been done to assess the performance of BDS-3 with the demonstration satellites. In terms of satellite visibility and geometry, Pan et al. (2019) presented the preliminary visibility and position dilution of precision (PDOP) performance of BDS with the few BDS-3 experimental satellites [23]. Wang et al. (2019) displayed the PDOP for the full BDS-3 constellation based on simulated almanacs [2]. In the aspect of signal quality, Zhang et al. (2017) found that the signal quality of the BDS-3 satellites was significantly greater than BDS-2 and comparable to the global positioning system (GPS) and Galileo satellite navigation system (GALILEO) [24]. The satellite-induced code variations could be ignored for BDS-3 satellites [25]. In terms of the performance of POD, Tan et al. (2016) evaluated the initial performance of POD and PCE for BDS-3 experimental satellites by using the B1I, B2I, and B3I observations [26]. Xie et al. (2019) presented that the RMS improvement of orbit overlap in three-dimension was about 42% by using 16 ground stations in global, reducing from 15.9 to 9.2 cm, when combined with the Ka-band inter-satellite link (ISL) observations [27]. Li et al. (2019) found that the overlap root mean square (RMS) of BDS-3 POD for along-track, cross-track, and radial components was 24.3, 16.1, and 8.4 cm when combined with BDS-2 and GPS [28]. Further, Yang et al. (2018) analyzed the positioning performance of BDS-3 using five demonstration satellites [25]. Qu et al. (2019) illustrated that the performance of the BDS PPP ambiguity resolution was significantly improved by using the observations of BDS-2 and BDS-3 satellites [29]. Zhang et al. (2019) evaluated that both the PPP and real-time kinematic (RTK) performance could be improved, when incorporating the observations of BDS-3 satellites [30]. In addition, the performance of a BDS-3 differential code bias (DCB) was evaluated by [31].

In all studies to date, the BDS-2 and BDS-3 are generally treated as the same system. In other words, it assumes that there is no inter-system bias between measurements to BDS-2 and BDS-3 satellites. However, this assumption must be verified since it will affect the PPP performance if the significant inter-system bias is not considered in the positioning model. Li et al. (2019) tried to evaluate the inter-system bias difference between BDS-2 and BDS-3 using three international GNSS monitoring and assessment system (iGMAS) stations but concluded that no obvious systematic bias between BDS-2 and BDS-3 was found [28]. Although this seems to agree with the previous assumption, the study has evaluated only a few receiver types. Therefore, more comprehensive evaluations are required with more receiver types to be included. For high precision applications using BDS after it becomes fully operational to provide a global service, it is therefore important to investigate and verify whether there exists an inter-system bias between the BDS-2 and BDS-3 satellites. How to estimate the inter-system bias and what is the impact of it on precise positioning should also be investigated in case the inter-system bias does exist. How to tightly combine the observations of BDS-2 and BDS-3 for precise data processing becomes an urgent problem to make the best use of BDS performance.

In this contribution, we have investigated the characteristics of ISB between BDS-2 and BDS-3 and found that there exists an inter-system bias between measurements of BDS-2 and BDS-3 satellites at a receiver due to the use of different clock references for BDS-2 and BDS-3 in the IGS precise satellites clock products and further due to inconsistent code hardware delays of BDS-2 and BDS-3 for some receiver types. The ISB values of the stations with the same receiver type are similar to each other but a significant difference up to several meters exists for different receiver types. This conclusion indicates that an additional ISB parameter should be introduced in the positioning model using BDS-2 and BDS-3, which however could be considered as a constant within a day. We have also evaluated the impact of ISB on the convergence time and accuracy for PPP using BDS-2 and BDS-3. The rest of the paper is organized as follows. Firstly, Section 2 describes the PPP mathematical model considering ISB between BDS-2 and BDS-3. Secondly, Section 3 displays the information of signal frequency and precise products for BDS. Then, the analysis of ISB between BDS-2 and BDS-3 is detailed in Section 4. Afterwards, the impact of ISB on the BDS PPP performance is presented in Section 5. In addition, a discussion is made in Section 6. Finally, the conclusions are summarized in Section 7.

2. Methods

The code and carrier phase observations of the ionosphere-free (IF) combination can be linearized as [32]:

$$P_{r,IF}^{s,Q}={\rho }_r^{s,Q}+c\left(d{t}_r-d{t}^{s,Q}\right)+T_r^{s,Q}+d_{r,IF}^Q-d_{IF}^{s,Q}+{\epsilon }_{P_{r,IF}^{s,Q}}$$

(1)

$$L_{r,IF}^{s,Q}={\rho }_r^{s,Q}+c\left(d{t}_r-d{t}^{s,Q}\right)+T_r^{s,Q}+{\lambda }_{IF}\left(N_{r,IF}^{s,Q}+b_{r,IF}^Q-b_{IF}^{s,Q}\right)+{\epsilon }_{L_{r,IF}^{s,Q}}$$

(2)

where superscript s and subscript r denote the satellite and receiver, respectively; superscript Q denotes the satellite system, i.e., C2 for BDS-2 satellites and C3 for BDS-3 satellites; $P_{r,IF}^{s,Q}$ and $L_{r,IF}^{s,Q}$ denote the observations of the IF code and carrier-phase in meters, respectively; ${\rho }_r^{s,Q}$ is the geometric distance between satellite s and receiver r in meters; c is the speed of light in meters; $d{t}_r$ and $d{t}^{s,Q}$ represent the receiver and satellite clock biases in seconds; $T_r^{s,Q}$ is the slant troposphere delay between satellite s and receiver r in meters; $d_{r,IF}^Q$ and $d_{IF}^{s,Q}$ are the hardware delays of the IF code observation for receiver r and satellite s in meters, respectively. Similarly, $b_{r,IF}^Q$ and $b_{IF}^{s,Q}$ are the hardware delays of the IF carrier-phase observation for receiver r and satellite s in cycles, respectively; ${\lambda }_{IF}$ denotes the IF wavelength in meters; $N_{r,IF}^{s,Q}$ is the ambiguity of the IF carrier-phase observation in cycles; ${\epsilon }_{P_{r,IF}^{s,Q}}$ and ${\epsilon }_{L_{r,IF}^{s,Q}}$ represent the IF code and carrier-phase measurement noises in meters. For simplicity, the other biases are supposed to be corrected well in Equations (1) and (2), such as the effects of phase winding in the satellite and receiver antenna, the phase center variation (PCV), relativistic effect, and tide load.

Further, $d_{r,IF}^Q$ and $d_{IF}^{s,Q}$ can be expressed as:

$$\left\{\begin{array}{c}{d}_{r,IF}^Q={\alpha }_{1,2}{d}_{r,1}^Q+{\beta }_{1,2}{d}_{r,2}^Q\\ d_{IF}^{s,Q}={\alpha }_{1,2}{d}_1^{s,Q}+{\beta }_{1,2}{d}_2^{s,Q}\end{array}\right.$$

(3)

$${\alpha }_{1,2}=\frac{f_1^2}{f_1^2-f_2^2} , {\beta }_{1,2}=-\frac{f_2^2}{f_1^2-f_2^2}$$

(4)

where f_i is the frequency of the ith signal.

As is known, the expression of the clock bias estimated in Equations (1) and (2) is determined by Equation (1) [33], that is, the hardware delays of the IF code observation will be absorbed by the clock bias. In addition, in order to study the characteristics of ISB between BDS-2 and BDS-3, the receiver clock bias of BDS-3 is chosen as a reference and an ISB parameter is introduced to BDS-2 observations in this paper, which is mainly because most BDS-2 satellites can only be observed in the Asian-Pacific region. Then, the PPP mathematical model of BDS-2 and BDS-3 can be described as follows:

$$P_{r,IF}^{s,C2}={\rho }_r^{s,C2}+c\left(dt{\prime }_r^{C3}-dt{\prime }^{s,C2}\right)+T_r^{s,C2}+IS{B}_r+{\epsilon }_{P_{r,IF}^{s,C2}}$$

(5)

$$P_{r,IF}^{s,C3}={\rho }_r^{s,C3}+c\left(dt{\prime }_r^{C3}-dt{\prime }^{s,C3}\right)+T_r^{s,C3}+{\epsilon }_{P_{r,IF}^{s,C3}}$$

(6)

$$L_{r,IF}^{s,C2}={\rho }_r^{s,C2}+c\left(d{t^{\prime }}_r^{C3}-dt{\prime }^{s,C2}\right)+T_r^{s,C2}+{\lambda }_{IF}{B}_{r,IF}^{s,C2}+IS{B}_r+{\epsilon }_{L_{r,IF}^{s,C2}}$$

(7)

$$L_{r,IF}^{s,C3}={\rho }_r^{s,C3}+c\left(dt{\prime }_r^{C3}-dt{\prime }^{s,C3}\right)+T_r^{s,C3}+{\lambda }_{IF}{B}_{r,IF}^{s,C3}+{\epsilon }_{L_{r,IF}^{s,C3}}$$

(8)

where

$$cd{t^{\prime }}_r^{C3}=cd{t}_r+d_{r,IF}^{C3}$$

(9)

$$cdt{\prime }^{s,C2}=cd{t}^{s,C2}+d_{IF}^{s,C2}$$

(10)

$$cdt{\prime }^{s,C3}=cd{t}^{s,C3}+d_{IF}^{s,C3}$$

(11)

$${\lambda }_{IF}{B}_{r,IF}^{s,C2}={\lambda }_{IF}\left[N_{r,IF}^{s,C2}+\left(b_{r,IF}^{C2}-b_{IF}^{s,C2}\right)\right]-\left(d_{r,IF}^{C2}-d_{IF}^{s,C2}\right)$$

(12)

$${\lambda }_{IF}{B}_{r,IF}^{s,C3}={\lambda }_{IF}\left[N_{r,IF}^{s,C3}+\left(b_{r,IF}^{C3}-b_{IF}^{s,C3}\right)\right]-\left(d_{r,IF}^{C3}-d_{IF}^{s,C3}\right)$$

(13)

$$IS{B}_r={\alpha }_{1,2}\left(d_{r,1}^{C2}-d_{r,1}^{C3}\right)+{\beta }_{1,2}\left(d_{r,2}^{C2}-d_{r,2}^{C3}\right)$$

(14)

Since the receiver and satellite clock parameters cannot be separated from each other in the precise satellite clock estimation [34], then a clock datum must be selected for each system [35]. Fortunately, the precise clock products of all satellites in a GNSS, e.g., BDS-2 and BDS-3, contain the same clock datum, which can be compensated by the receiver clock parameters in users. In this paper, the clock datum of IGS precise satellite clock products for BDS-2 and BDS-3 is set to $D_{IGS}^{C2}$ and $D_{IGS}^{C3}$, respectively. Equations (10) and (11) can be re-written as:

$$dt{\prime }_{IGS}^{s,C2}=d{t}^{s,C2}+d_{IF}^{s,C2}+D_{IGS}^{C2}$$

(15)

$$dt{\prime }_{IGS}^{s,C3}=d{t}^{s,C3}+d_{IF}^{s,C3}+D_{IGS}^{C3}$$

(16)

When the IGS precise products are applied, Equations (9) and (14) can be re-expressed as:

$$d{t^{\prime }}_r^{C3}=d{t}_r+d_{r,IF}^{C3}+D_{IGS}^{C3}$$

(17)

$$ISB{\prime }_{r,IGS}={\alpha }_{1,2}\left(d_{r,1}^{C2}-d_{r,1}^{C3}\right)+{\beta }_{1,2}\left(d_{r,2}^{C2}-d_{r,2}^{C3}\right)+\left(D_{IGS}^{C2}-D_{IGS}^{C3}\right)$$

(18)

At current, the IGS analysis center (AC) providing global users with BDS-3 satellites precise orbit and clock products includes Wuhan University (WUM) and German Research Centre for Geosciences (GBM). Similarly, the expressions of ISB when using GBM and WUM precise products can be presented as:

$$ISB{\prime }_{r,GBM}={\alpha }_{1,2}\left(d_{r,1}^{C2}-d_{r,1}^{C3}\right)+{\beta }_{1,2}\left(d_{r,2}^{C2}-d_{r,2}^{C3}\right)+\left(D_{GBM}^{C2}-D_{GBM}^{C3}\right)$$

(19)

$$ISB{\prime }_{r,WUM}={\alpha }_{1,2}\left(d_{r,1}^{C2}-d_{r,1}^{C3}\right)+{\beta }_{1,2}\left(d_{r,2}^{C2}-d_{r,2}^{C3}\right)+\left(D_{WUM}^{C2}-D_{WUM}^{C3}\right)$$

(20)

For station r, the signal delay of receiver is the same in Equations (19) and (20). Then, a new parameter $∆ISB{\prime }_{r,GBM-WUM}$ can be derived as:

$$∆ISB{\prime }_{r,GBM-WUM}=ISB{\prime }_{r,GBM}-IS{B^{\prime }}_{r,WUM}=\left(D_{GBM}^{C2}-D_{GBM}^{C3}\right)-\left(D_{WUM}^{C2}-D_{WUM}^{C3}\right)$$

(21)

If we would assume that the same clock datums are applied for BDS-2 and BDS-3 in IGS, that is:

$$D_{GBM}^{C3}=D_{GBM}^{C2}$$

(22)

$$D_{WUM}^{C3}=D_{WUM}^{C2}$$

(23)

then

$$∆ISB{\prime }_{r,GBM-WUM}=0$$

(24)

It means that the values of $∆ISB{\prime }_{r,GBM-WUM}$ for all stations should be zero in case of BDS-2 and BDS-3 with the same clock references in IGS.

According to Equation (18), by constructing the single difference between two stations, the part of clock datums for BDS-2 and BDS-3 in the ISB can be completely eliminated, that is, only the code hardware delays at the receiver are retained.

$$∆ISB{\prime }_{r,r_0,IGS}={\alpha }_{1,2}\left[\left(d_{r,1}^{C2}-d_{r,1}^{C3}\right)-\left(d_{r_0,1}^{C2}-d_{r_0,1}^{C3}\right)\right]+{\beta }_{1,2}\left[\left(d_{r,2}^{C2}-d_{r,2}^{C3}\right)-\left(d_{r_0,2}^{C2}-d_{r_0,2}^{C3}\right)\right]$$

(25)

where r₀ represents the reference station.

If we would make the assumption that the code hardware delays of BDS-2 and BDS-3 at the receiver are consistent with each other:

$$d_{r,1}^{C2}=d_{r,1}^{C3}$$

(26)

$$d_{r,2}^{C2}=d_{r,2}^{C3}$$

(27)

then

$$∆ISB{\prime }_{r,r_0,IGS}=0$$

(28)

It indicates that the values of $∆ISB{\prime }_{r,r_0,IGS}$ for all stations should be zero if the above assumptions are true. Otherwise, the code hardware delays of BDS-2 and BDS-3 at the receiver may be different.

3. BDS-3 Precise Products Description

At present, BDS-3 is close to a full constellation. In the near future, BDS will become a true global positioning system such as the GPS. As of May 2020, up to 41 BDS satellites are available to global users, including 26 BDS-3 satellites and 15 BDS-2 satellites. Figure 1 presents the ground tracks of BDS-2 and BDS-3 satellites at day of year (DOY) 110, 2020. It can be seen that the ground trajectories of GEO, IGSO, and MEO satellites are different. The trajectories density of BDS has been increased with both BDS-2 and BDS-3 satellites. The ground trajectories of the MEO satellites confine from 56.1° S to 56.1° N in latitude around the world. The ground tracks of the IGSO satellites present several figure-of-eight loops centered at the Asian-Pacific region, approximately confined from 57.5° S to 57.5° N in latitude. The GEO satellites are evenly distributed from 50° to 160° E along the equator, which are almost stationary.

In addition to the three signals shared with the BDS-2 satellites, i.e., B1I (1561.098 MHz), B2I (1207.140 MHz), and B3I (1268.520 MHz), the BDS-3 satellites have also broadcasted three new signals, i.e., B1C (1575.420 MHz), B2a (1176.450 MHz), and B2a+b (1191.795 MHz) [36]. The information of the precise orbit and satellites clock products for BDS-2 and BDS-3 currently released by GBM and WUM is listed in

Table 1. The information of IGS precise products for BDS-3.

Institution	ID	Observation Types	PCO/PCV	Interval	Availability
WHU	WUM0MGXFIN	B1I, B3I for BDS-2, and BDS-3	igs14_2101.atx	SP3: 15 min CLK: 30 s	Since DOY 265, 2019
GFZ	GBM0MGXRAP			SP3: 5 min CLK: 30 s/5 min	Since DOY 332, 2019

.

4. The Analysis of ISB Between BDS-2 and BDS-3

This section analyzes the characteristics of ISB between BDS-2 and BDS-3 based on the observations from multi-GNSS experiment (MGEX) tracking stations, and confirms whether it is necessary to introduce an ISB parameter when combining the observations of BDS-2 and BDS-3.

4.1. Processing Strategy for ISB

Usually, the ISB parameter can be regarded as white noise or constant parameter. For the white noise, the ISB parameter is independent with each other among different epochs and uncorrelated over time, in the case of which an ISB parameter needs to be introduced for each epoch. For the case of constant parameter, the ISB parameter is estimated as a constant parameter over time, that is, only one ISB parameter needs to be introduced among different epochs during a day. Due to the uncertain characteristics of ISB, both white noise (WN) and constant parameter (CON) schemes will be applied to decide the suitable scheme for ISB in this paper. This is because the white noise scheme can more realistically reflect the variation of ISB parameter over time. Obviously, if the ISB results under the WN and CON schemes are consistent within one day, it indicates that the ISB are stable within a day, which can be seen as a constant parameter.

For the satellite-induced code biases in BDS-2 satellites, the elevation-dependent correction model according to Wanninger and Beer (2015) is applied for BDS-2 IGSO and MEO satellites, but it is not suitable for BDS-2 GEO satellites [10]. It is worth noting that the correction is a somewhat arbitrary absolute level, which may cause inconsistent code hardware delays for GEO and IGSO/MEO at the receiver. For a station in the Asian-Pacific region, up to seven BDS-2 IGSO and five BDS-2 GEO satellites could be visible, however the number of visible BDS-2 MEO is about one in most time because of only three MEO satellites in the BDS-2 constellation. Thus, the consistency of receiver code hardware delays for BDS-2 IGSO and GEO is to be studied in this paper by comparing the ISB values of two different combinations, e.g., BDS-2 IGSO+BDS-3 and BDS-2 GEO+BDS-3. Finally, three stations in Australia, such as CEDU (SEPT POLARX5 5.3.2), MRO1 (TRIMBLE NETR9 5.44), and NNOR (SEPT POLARX5TR 5.3.2) are selected, which is mainly because the satellite visibility of them are similar to each other. The distribution of three stations are shown in Figure 2.

With the precise products of GBM, Figure 3 shows the ISB results with both WN and CON schemes for two different combinations, e.g., BDS-2 IGSO+BDS-3 and BDS-2 GEO+BDS-3, respectively at DOY 110, 2020. As the figure shows, for the BDS-2 IGSO+BDS-3 and BDS-2 GEO+BDS-3 combinations, the ISB results of each station tends to be stable under the WN and CON schemes with a range less than 0.3 m, which are very close to each other. It indicates that the ISB parameters for these two combinations are stable within a day. Further, the ISB result of the BDS-2 IGSO+BDS-3 combination is close to that of BDS-2 GEO+BDS-3 during a day, with a range of less than 0.3 m after convergence, especially for station NNOR. That is, there is no obvious systematic bias between BDS-2 IGSO and GEO, which indicates that the BDS-2 IGSO, MEO, and GEO can be still treated as one system after the corrections of the satellite-induced code biases for IGSO and MEO.

In order to verify the suitable stochastic model for ISB between BDS-2 and BDS-3, four MGEX stations with different receiver types are selected, including JFNG (TRIMBLE NETR9 5.43), CEDU (SEPT POLARX5 5.3.2), NNOR (SEPT POLARX5TR 5.3.2), WTZZ (JAVAD TRE_3 DELTA 3.7.9), and the distribution of these stations is presented in Figure 4. As a comparison, both GBM and WUM precise products are used to estimate ISB.

Figure 5 shows the series of ISB within a day for JFNG, CEDU, NNOR, and WTZZ with both WN and CON schemes at DOY 110, 2020. It can be seen that the ISB of each station estimated by the two strategies can remain stable after a period of time within a range of less than 0.1 m, but the ISB values of four stations are different. The ISB result of GBM is greater than that of WUM for each station, which may be caused by the different clock datums applied in the two ACs. The stability of ISB for JFNG, NNOR, and CEDU in the Asian-Pacific region is better than that of WTZZ. This is mainly because the station WTZZ is located in the European region, where the number of observable BDS-2 satellites is limited, so more time is needed by WTZZ to be stable compared to others. It is worth noting that the ISB series estimated with the CON scheme is very close to that of the WN scheme for both GBM and WUM precise products in each station, which indicates that the ISB parameter is stable within a day, and can be considered as a constant.

Figure 6 shows the series of $∆ISB{\prime }_{r,GBM-WUM}$ calculated by Equation (21) within a day for the above four MGEX stations at DOY 110, 2020. As can be seen from the figure, the series of $∆ISB{\prime }_{r,GBM-WUM}$ for each station is very stable within one day, and the standard deviations (STD) are better than 0.14 m. The mean values of JFNG, CEDU, NNOR, and WTZZ during a day are 1.07, 1.03, 1.03, and 0.86 m, respectively, which are very close to a constant. According to Equation (21), only clock datums of precise clock products for BDS-2 and BDS-3 in GBM and WUM are contained in parameter $∆ISB{\prime }_{r,GBM-WUM}$, which are the same for all stations within one day. In theory, the value of $∆ISB{\prime }_{r,GBM-WUM}$ for each station should be also the same, which is consistent with the results shown in Figure 6. However, the fact that different clock datums are applied to BDS-2 and BDS-3 in IGS precise products is presented, because if the BDS-2 and BDS-3 with the same clock references in IGS precise products, then Equations (22) and (23) will be true. That is, the value of $∆ISB{\prime }_{r,GBM-WUM}$ for each station should be close to zero, but it is obviously inconsistent with the results shown in Figure 6. According to Equation (18), different clock datums applied to BDS-2 and BDS-3 in IGS precise products may cause the non-zero ISB value in each station.

Based on the above analysis, the value of ISB between BDS-2 and BDS-3 is non-zero, which is caused by the fact that different clock datums for BDS-2 and BDS-3 are used in IGS, so it is necessary to introduce an additional ISB parameter when using both BDS-2 and BDS-3 observations. Further, the result shows that the ISB is very stable within a day and can be regarded as a constant. Therefore, the CON scheme will be applied to estimate the ISB parameter in this paper.

4.2. The Consistency of BDS-2 and BDS-3 Code Hardware Delays

According to Equation (18), the estimated ISB includes two parts: The difference of clock datums between BDS-2 and BDS-3, and the difference of hardware delays of BDS-2 and BDS-3 code observations at the receiver. Although the clock references of BDS-2 and BDS-3 in IGS precise products are the same for all stations in a day, it may be different between different days. Therefore, in order to only analyze the part of receiver hardware delays of code observations for BDS-2 and BDS-3, the clock references included in ISB should be eliminated. In this paper, the single-difference (SD) ISB value between two stations is constructed according to Equation (25) to eliminate the clock references of BDS-2 and BDS-3.

Then, up to 85 MGEX stations (Figure 7) supporting BDS-2 and BDS-3 observations in global are selected to calculate the values of $∆ISB{\prime }_{r,r_0,IGS}$, covering several receiver manufacturers, such as JAVAD, SEPT, and TRIMBLE.

Table 2. Receiver information of the selected MGEX stations.

Manufacturer	Receiver Types	Receiver Version	Supported BDS-3 Satellites
JAVAD	JAVAD TRE_3	3.7.6	C19-C30, C32-C37
	JAVAD TRE_3 DELTA	3.7.9	C19-C30, C32-C46
SEPTENTRIO	SEPT ASTERX4	4.4.3	C19-C30, C32-C37
	SEPT POLARX5	5.3.0	C19-C20, C32-C37
		5.3.2
	SEPT PLOARX5TR	5.3.2
TRIMBLE	TRIMBLE ALLOY	5.37	C19-C30
		5.43	C19-C30, C32-C46
		5.44
	TRIMBLE NETR9	5.37	C19-C30
		5.42
		5.43
		5.44

details the information of the selected receiver types. Station NNOR located in Australia is set as the reference. The observations collected from DOY 091 to 110, 2020 for 85 MGEX stations are used. The sampling interval is 30 s, and the ISB parameter is estimated as a constant during a day for each station.

The $∆ISB{\prime }_{r,r_0,IGS}$ results for all selected MGEX stations by using GBM and WUM precise products are shown in Figure 8 and Figure 9, from DOY 091 to 110, 2020. The $∆ISB{\prime }_{r,r_0,IGS}$ of each station is very stable over time, but one of the stations has a unique result compared to the others for GBM and WUM. It is confirmed that this station is MCHL with the receiver type TRIMBLE ALLOY 5.37 in Australia. Except for the station MCHL, in the 20-day time span, the values of $∆ISB{\prime }_{r,r_0,IGS}$ for the other stations approximately range from −3.3 to 1.3 m with GBM precise products, and −3.3 to 1.5 m with WUM precise products, which are in a range of about 4.8 m for both GBM and WUM. Obviously, not all the stations are with the same ISB values, which are inconsistent with the expected value zero by Equation (28). It means that Equations (26) and (27) may be not true for some receiver types listed in Table 2, that is, the code hardware delays of BDS-2 and BDS-3 are different for these receiver types, especially the TRIMBLE ALLOY 5.37.

Further, from the above 85 MGEX stations, four receiver types with more stations are used to evaluate the relationship between the value of $∆ISB{\prime }_{r,r_0,IGS}$ and the receiver type, including SEPT POLARX5 5.3.2 with 32 stations; SEPT POLARX5TR 5.3.2 with 9 stations; TRIMBLE NETR9 5.43 with 9 stations; and TRIMBLE NETR9 5.44 with 16 stations. The distribution of the four selected receiver types is shown in Figure 10.

Figure 11, Figure 12, Figure 13 and Figure 14 present the values of $∆ISB{\prime }_{r,r_0,IGS}$ for SEPT POLARX5 5.3.2, SEPT POLARX5TR 5.3.2, TRIMBLE NETR9 5.43, and TRIMBLE NETR9 5.44, respectively, from DOY 091 to 110, 2020. Compared to the results shown in Figure 8 and Figure 9, it can be seen that the variation range of $∆ISB{\prime }_{r,r_0,IGS}$ values for the stations with the same receiver type is significantly small, especially the receiver types, such as SEPT POLOARX5TR 5.3.2 and TRIMBLE NETR9 5.44. The STDs of the $∆ISB{\prime }_{r,r_0,IGS}$ estimated by using WUM precise products are slightly better than that of GBM. With the WUM precise products, the mean values of SEPT POLARX5 5.3.2, SEPT POLARX5TR 5.3.2, TRIMBLE NETR9 5.43, and TRIMBLE NETR9 5.44 are −0.88, −0.08, −0.55, and 0.13 m, among which SEPT POLARX5 5.3.2 is with the largest absolute value. The stations with receiver type TRIMBLE NETR9 5.44 display the best consistency, followed by SEPT PLOARX5TR 5.3.2 and SEPT POLARX5 5.3.2, while the TRIMBLE NETR9 5.43 is the worst. The STD of SEPT POLARX5TR 5.3.2 and TRIMBLE NETR9 5.44 with the WUM precise products are 0.27 and 0.2 m, which indicates that the values of $∆ISB{\prime }_{r,r_0,IGS}$ may be related to the receiver type and the stations with the same receiver type may have similar values. However, the $∆ISB{\prime }_{r,r_0,IGS}$ values of different receiver types may be quite different, such as TRIMBLE NETR9 5.44 and SEPT POLARX5 5.3.2, which are also inconsistent with the expected result with a value of zero in Equation (28) for all receiver types. It shows again that Equations (26) and (27) may not be true for some receiver types.

The above research shows that different clock references are applied for BDS-2 and BDS-3 in the IGS precise satellites clock products, and for some receiver types, the code hardware delays of BDS-2 and BDS-3 are inconsistent, which both cause the non-zero value of ISB between the BDS-2 and BDS-3 at the receiver. In other words, when combining BDS-2 and BDS-3 observations for PPP, an additional ISB parameter needs to be introduced. Further, the results show that the values of ISB are closely related to the receiver type. The stations with the same receiver type may be with the similar value, but a great difference may be presented for different receiver types, up to several meters.

5. The Impact of ISB on BDS PPP

According to Section 4, it can be known that an additional ISB parameter is necessary for BDS-2 and BDS-3, which can be regarded as a constant within a day. As is well known, the position accuracy is actually determined by the carrier-phase observations when combining the code and phase observations. Since the weight values of the carrier-phase observations are much larger than the code observations, if no ISB parameters were introduced to Equations (5)–(8), it can actually be completely absorbed by the ambiguity parameter and will not affect the accuracy of the carrier-phase observation, which is mainly due to the constant characteristics of ISB. Unfortunately, the bias caused by the ISB parameter in the pseudorange observations may be absorbed by the residuals. Therefore, the impact of ISB on BDS-2 and BDS-3 PPP will be evaluated in the section.

5.1. Data Description and Processing Strategy

Datasets collected from six MGEX stations with different receiver types, i.e., JFNG (TRIMBLE NETR9 5.43), MRO1 (TRIMBLE NETR9 5.44), MCHL (TRIMBLE ALLOY 5.37), CEDU (SEPT POLARX5 5.3.2), NNOR (SEPT POLARX5TR 5.3.2), and WTZZ (JAVAD TRE_3 DELTA 3.7.9) are employed to assess the performance of BDS-2 and BDS-3 PPP. The observations of the six MGEX stations are from DOY 091 to 110, 2020, with a sampling interval of 30 s. Figure 15 displays the distribution of the six MGEX stations.

The observations of each station will be processed under two schemes: With an additional ISB parameter considered as a constant and without the ISB parameter, respectively. For the sake of brevity, ISB_CON denotes the solutions with the ISB parameter, and ISB_NO denotes the solutions without the ISB parameter. The precise products of GBM are applied for PPP processing. The observations of B1I and B3I are used for BDS-2 and BDS-3 satellites. The cut-off elevation is set to 7°. The weight value of observations is set by the weighted strategy based on elevation, i.e., ${\sin }^2\left(E\right)$ [37]. The zenith tropospheric wet delay is considered as a random walk process with the empirical power density $15mm\sqrt{h}$ [38] and initialized every 2 h. A float carrier-phase ambiguity parameter is estimated for each arc. Since the PCO and PCV information of a receiver antenna is not available to BDS at current, the corresponding GPS values will be used instead. In addition, the station coordinates from the IGS weekly SNX solution are considered as the truth.

5.2. Results and Analysis

In this paper, the performance of PPP is assessed in terms of convergence time and position accuracy. In this study, the convergence time is defined as the time needed to obtain a position accuracy less than the threshold for the consecutive twenty epochs [37]. The threshold of the east, north, and up directions is set to 5, 5, 10 cm for static PPP [39], and 10, 10, 20 cm for kinematic PPP [30].

Figure 16 shows the ISB values of the six selected MGEX stations from DOY 91 to 110, 2020. Meanwhile, the mean value and STD of the ISB for each station are listed in

Table 3. The mean ISB of six stations with GBM precise products from DOY 091 to 110, 2020, unit—m.

	CEDU	JFNG	MCHL	MRO1	NNOR	WTZZ
mean	2.89	2.09	−4.84	−0.29	1.87	−1.23
std	0.06	0.08	0.08	0.06	0.06	0.07

. It can be seen that the ISB series of each station during 20 consecutive days is very stable with a range of less than 0.3 m, and the STD is better than 0.08 m. However, the ISB values of all stations are not same. The ISB value of station MCHL is with the largest absolute value, -4.84 m, and the ISB absolute value of MRO1 is the smallest, only −0.29 m. It is worth noting that the ISB values of the stations JFNG and MRO1 with receiver type TRIMBLE NETR9 are 2.09 and −0.29 m respectively, and the difference is up to 2.38 m, which is mainly because different receiver versions are used. The station MRO1 with the receiver version 5.44 is significantly smaller than the station JFNG with the receiver version 5.43, which indicates that the latest receiver version may be better compatible with BDS-2 and BDS-3 signals.

Since station MCHL has the largest absolute value of ISB, Figure 17 presents the static PPP solution of station MCHL at DOY 97, 2020. For the static PPP, under the ISB_CON and ISB_NO schemes, station MCHL can both achieve a three-dimensional convergence in about 30 min. However, the convergence of ISB_CON is different than that of ISB_NO in the east, north, and up directions at the beginning of the convergence, which lasted about 30 min. The position accuracy of the east and up directions under the ISB_CON strategy is better than the ISB_NO strategy. It is worth noting that after a short time, these differences in three directions between ISB_CON and ISB_NO could be ignored, which have little impact on the final convergence time of static PPP. After about 30 min of convergence, the position accuracy under the ISB_CON and ISB_NO schemes all tend to coincide with each other in the east, north, and up directions, respectively. It indicates that the ISB parameter may only affect the initial convergence of static PPP, but the effect on the position accuracy after convergence is negligible.

In addition, the kinematic PPP solution of station MCHL is displayed in Figure 18 at DOY 97, 2020. In the initial convergence of about 60 min, the coordinate accuracy under the ISB_CON strategy in the east, north, and up directions is better than that of the ISB_NO strategy, especially in the east and up directions. After about 60 min of convergence, the position under the ISB_CON and ISB_NO strategies in the east, north, and up directions also tend to coincide, which is similar to the static PPP. It means that the ISB parameter will not affect the position accuracy after PPP convergence. Compared with the static PPP, the influence of ISB on the initial convergence of kinematic PPP is more obvious, and the introduction of ISB can speed up the convergence.

Figure 19 and Figure 20 show the statistics of the mean convergence time of the static and kinematic PPP for the six MGEX stations under the ISB_CON and ISB_NO schemes from DOY 91 to 110, 2020. For the static PPP, under the ISB_CON and ISB_NO schemes, the six MGEX stations can achieve three-dimensional convergence in about 30 min, and there is no significant difference in the east, north, and up directions. It indicates that the ISB parameter only affects the position accuracy during the initial convergence of static PPP, and the impact on the convergence time is negligible. For the kinematic PPP, it is worth noting that the three-dimensional convergence time under ISB_CON is shorter than that of ISB_NO for the stations CEDU, JFNG, MCHL, NNOR, and WTZZ, which is mainly because the ISB absolute value of these stations is more than 1.2 m. The improvement is about 5 min. However, for station MRO1, the three-dimensional convergence time under the ISB_CON strategy is longer than ISB_NO. The reason may be that the ISB value of MRO1 is only −0.29 m, and an additional ISB parameter into the kinematic PPP model may slightly weaken the model strength.

Table 4. RMS (unit—cm) of static PPP solutions under the ISB_CON and ISB_NO schemes at different convergence epochs of 15 min, 30 min, 1 h, 2 h and 6 h, from DOY 091 to 110, 2020.

		15 min		30 min		1 h		2 h		6 h
		CON	NO	CON	NO	CON	NO	CON	NO	CON	NO
CEDU	east	15.4	17.0	6.7	6.8	4.1	3.9	2.8	2.7	0.5	0.5
	north	3.1	3.1	1.9	1.9	1.2	1.3	0.8	0.9	0.6	0.6
	up	10.6	11.9	3.3	3.8	2.4	2.3	2.2	2.2	1.5	1.5
JFNG	east	26.5	30.6	8.7	10.0	3.3	3.3	1.4	1.4	0.9	0.9
	north	9.0	9.0	2.0	2.0	1.1	1.1	0.6	0.6	0.6	0.6
	up	47.9	44.5	7.3	7.1	4.1	3.9	3.2	3.1	2.4	2.3
MCHL	east	23.2	25.3	10.0	10.4	3.7	3.7	3.6	3.6	0.6	0.7
	north	6.4	7.3	2.3	2.5	1.5	1.6	1.1	1.1	0.4	0.4
	up	10.3	12.0	5.3	6.0	3.3	3.2	2.5	2.4	1.8	1.8
MRO1	east	19.6	19.2	5.5	5.4	2.6	2.6	1.6	1.6	0.5	0.5
	north	7.1	6.9	1.8	1.7	1.0	1.0	0.6	0.6	0.4	0.4
	up	36.7	36.4	12.3	12.3	4.0	4.0	1.6	1.6	1.4	1.4
NNOR	east	14.2	14.5	6.6	6.5	3.0	2.9	2.1	2.1	0.7	0.7
	north	6.5	6.3	2.8	2.7	1.2	1.2	0.5	0.5	0.4	0.4
	up	16.0	17.4	5.9	6.1	3.8	3.7	1.7	1.7	2.1	2.1
WTZZ	east	9.7	9.3	4.5	4.4	2.1	2.1	1.1	1.1	0.4	0.4
	north	7.3	7.0	3.3	3.3	1.7	1.7	1.0	1.0	0.3	0.3
	up	16.8	17.0	6.7	6.7	2.5	2.5	1.6	1.6	1.5	1.5

details the RMS of the static PPP under the ISB_CON and ISB_NO schemes at different convergence epochs for six MGEX stations, from DOY 091 to 110, 2020. For station MRO1, the RMS of the ISB_CON and ISB_NO schemes are basically the same for different session lengths. However, for the other five stations, after 15 min of convergence, the RMS of the east and up directions with the ISB_CON scheme is better than that of ISB_NO, which means that the ISB could improve the position accuracy during the initial convergence period. When the convergence time exceeds 30 min, there is no obvious difference in the position accuracy between ISB_CON and ISB_NO, which further indicates that the impact of ISB on the position accuracy is basically negligible after convergence.

6. Discussion

The results in Section 5 indicate that if the ISB parameter is ignored in Equations (5)–(8), the bias will not affect the accuracy of the carrier-phase observations, which may be absorbed by the residuals of the pseudorange observations. To verify this, the pseudorange residuals will be discussed in detail.

In general, the residuals of the pseudorange observations is considered as white noise, and the mathematical expectation is zero. However, due to factors such as the multipath, the value of the mathematical expectation may be not zero. For station r, the pseudorange residuals expectations under ISB_CON for BDS-2 and BDS-3 are supposed as:

$$\left\{\begin{array}{c}E\left({\epsilon }_{P_{r,IF}^{s,C2}}\right)=e_{r,IF}^{C2}\\ E\left({\epsilon }_{P_{r,IF}^{s,C3}}\right)=e_{r,IF}^{C3}\end{array}\right.$$

(29)

At the same time, we assume that the constant part of the ISB absorbed by the pseudorange residuals of the BDS-2 and BDS-3 is ${\mathsf{\Delta }}_r^{C2}$ and ${\mathsf{\Delta }}_r^{C3}$, respectively when no ISB parameters are introduced for BDS-2 an BDS-3. In this case, Equations (5)–(8) could be re-expressed as:

$$P_{r,IF}^{s,C2}={\rho }_r^{s,C2}+c\left(\overline{d}{t}_r^{C3}-dt{\prime }^{s,C2}\right)+T_r^{s,C2}+\epsilon {\prime }_{P_{r,IF}^{s,C2}}$$

(30)

$$P_{r,IF}^{s,C3}={\rho }_r^{s,C3}+c\left(\overline{d}{t}_r^{C3}-dt{\prime }^{s,C3}\right)+T_r^{s,C3}+\epsilon {\prime }_{P_{r,IF}^{s,C3}}$$

(31)

$$L_{r,IF}^{s,C2}={\rho }_r^{s,C2}+c\left(\overline{d}{t}_r^{C3}-d{t}^{s,C2}\right)+T_r^{s,C2}+{\lambda }_{IF}\left(B_{r,IF}^{s,C2}+{\mathsf{\Delta }}_r^{C2}/{\lambda }_{IF}\right)+{\epsilon }_{L_{r,IF}^{s,C2}}$$

(32)

$$L_{r,IF}^{s,C3}={\rho }_r^{s,C3}+c\left(\overline{d}{t}_r^{C3}-d{t}^{s,C3}\right)+T_r^{s,C3}+{\lambda }_{IF}\left(B_{r,IF}^{s,C3}+{\mathsf{\Delta }}_r^{C3}/{\lambda }_{IF}\right)+{\epsilon }_{L_{r,IF}^{s,C3}}$$

(33)

where

$$\epsilon {\prime }_{P_{r,IF}^{s,C2}}={\epsilon }_{P_{r,IF}^{s,C2}}+{\mathsf{\Delta }}_r^{C2}$$

(34)

$$\epsilon {\prime }_{P_{r,IF}^{s,C3}}={\epsilon }_{P_{r,IF}^{s,C3}}+{\mathsf{\Delta }}_r^{C3}$$

(35)

$$c\overline{d}{t}_r^{C3}=cdt{\prime }_r^{C3}-{\mathsf{\Delta }}_r^{C3}$$

(36)

$$ISB{\prime }_r={\mathsf{\Delta }}_r^{C2}-{\mathsf{\Delta }}_r^{C3}$$

(37)

Obviously, the pseudorange residuals expectations under ISB_NO for BDS-2 and BDS-3 could be derived as:

$$\left\{\begin{array}{c}E\left(\epsilon {\prime }_{P_{r,IF}^{s,C2}}\right)=e{\prime }_{r,IF}^{C2}=e_{r,IF}^{C2}+{\mathsf{\Delta }}_r^{C2}\\ E\left(\epsilon {\prime }_{P_{r,IF}^{s,C3}}\right)=e{\prime }_{r,IF}^{C3}=e_{r,IF}^{C3}+{\mathsf{\Delta }}_r^{C3}\end{array}\right.$$

(38)

Further, the value of ISB calculated based on Equations (29), (37), and (38) is presented as:

$$ISB{\prime }_{r,CAL}=\left[e{\prime }_{r,IF}^{C2}-e{\prime }_{r,IF}^{C3}\right]-\left[e_{r,IF}^{C2}-e_{r,IF}^{C3}\right]$$

(39)

Figure 21 and Figure 22 show the residuals of station MCHL under the ISB_CON and ISB_NO schemes, respectively. It can be seen that under the ISB_CON and ISB_NO schemes, the carrier-phase residuals are both randomly distributed along the zero Y-axis, and the mean value of the residuals is approximately 0 m, which indicates that the bias caused by ISB has no effect on the carrier-phase residuals. However, it is worth noting that the pseudorange residuals of the ISB_CON is obviously different from that of ISB_NO. The pseudorange residuals of ISB_CON is randomly distributed along the Y-axis with a near zero value, especially for BDS-3, while the pseudorange residuals of ISB_NO appear to deviate from the zero value in the Y-axis. What is more, for BDS-2 and BDS-3, the mean value and RMS under ISB_NO are greater than those of ISB_CON. Obviously, the pseudorange residuals of ISB_NO have been contaminated by the ISB.

As shown in Figure 21 and Figure 22, in the case of ISB_CON and ISB_NO, the pseudorange residuals expectations (unit—m) of station MCHL for BDS-2 and BDS-3 are:

$$e_{r,IF}^{C2}=-0.65, e_{r,IF}^{C3}=0.05$$

(40)

$$e{\text{'}}_{r,IF}^{C2}=-2.32, e{\text{'}}_{r,IF}^{C3}=3.34$$

(41)

According to Equations (39)–(41), the ISB value calculated by the pseudorange residuals is derived as (unit—m):

$$ISB{\prime }_{r,CAL}=-4.96$$

(42)

Based on Figure 16, the estimated value of ISB under ISB_CON is about -4.96 m for MCHL at DOY 97, 2020, which is consistent with the result in Equation (42). It verifies that the bias will be completely absorbed by the pseudorange residuals of BDS-2 and BDS-3 when no ISB parameter is introduced for BDS-2 and BDS-3, that is, the carrier-phase observation will not be contaminated by the bias. Additionally, it further proves that the ISB parameter between BDS-2 and BDS-3 is constant within a day, because according to the principle of least squares, it can be completely absorbed by the ambiguity parameter in the carrier-phase observation.

7. Conclusions

At present, up to 41 BDS satellites have been in orbits to provide positioning, navigation, and timing services to global users. In order to verify whether an additional ISB parameter needs to be introduced for BDS-2 and BDS-3 PPP, the characteristics of an inter-system bias parameter between BDS-2 and BDS-3 has been investigated by using observations of 85 MGEX stations in global with a range from DOY 91 to 110, 2020. The results show that the value of ISB between BDS-2 and BDS-3 at the receiver is non-zero, which is mainly because: Different clock references are applied for BDS-2 and BDS-3 in the IGS precise satellites clock products, and for some receiver types, the code hardware delays of BDS-2 and BDS-3 are inconsistent. Thus, an additional ISB parameter needs to be introduced for BDS-2 and BDS-3 PPP. Moreover, the ISB could be considered as a constant within a day, and is closely related to the receiver type. The stations with the same receiver type have a similar ISB value to each other, but a great difference may be presented for different receiver types, up to several meters.

In addition, the impact of the ISB parameter on the BDS-2 and BDS-3 static and kinematic PPP performance have also been studied in terms of convergence time and position accuracy. The results show that the impact of ISB on PPP is related to the absolute value of ISB, especially for the kinematic PPP. For the static PPP, the position accuracy with ISB is better than the case without the ISB parameter in the period of the early 30 min convergence, especially in east and up directions, but the effect of ISB on the position after convergence is negligible. For the kinematic PPP, in the initial convergence of about 60 min, the position accuracy of the stations with a large ISB absolute value in the east, north, and up directions is obviously better than the case without ISB, especially in the east and up directions, and the convergence performance could be improved at the same time. Similar to the static PPP, the effect of ISB on the improvement of position accuracy is not obvious after convergence. The reason for the above results is that the position accuracy is actually determined by the carrier-phase observations when combining the code and phase observations. Since the weight values of the carrier-phase observations are much larger than the code observations, if no ISB parameters were introduced for BDS-2 and BDS-3, the bias will be completely absorbed by the pseudorange residuals, which has no effect on the accuracy of carrier-phase observations.