BDS-3/GNSS Undifferenced Pseudorange and Phase Time-Variant Mixed OSB Considering the Receiver Time-Variant Biases and Its Benefit on Multi-Frequency PPP

Abstract

The legacy Global Navigation Satellite System (GNSS) satellite clock offsets obtained by the dual-frequency undifferenced (UD) ionospheric-free (IF) model absorb the code and phase time-variant hardware delays, which leads to the inconsistency of the precise satellite clock estimated by different frequencies. The dissimilarity of the satellite clock offsets generated by different frequencies is called the inter-frequency clock bias (IFCB). Estimates of the IFCB typically employ epoch-differenced (ED) geometry-free ionosphere-free (GFIF) observations from global networks. However, this method has certain theoretical flaws by ignoring the receiver time-variant biases. We proposed a new undifferenced model coupled with satellite clock offsets, and further converted the IFCB into the code and phase time-variant mixed observable-specific signal bias (OSB) to overcome the defects of the traditional model and simplify the bias correction process of multi-frequency precise point positioning (PPP). The new model not only improves the mixed OSB performance, but also avoids the negative impact of the receiver time-variant biases on the satellite mixed OSB estimation. The STD and RMS of the original OSB can be improved by 7.5–60.9% and 9.4–66.1%, and that of ED OSB (it can reflect noise levels) can be improved by 50.0–87.5% and 60.0–88.9%, respectively. Similarly, the corresponding PPP performance for using new mixed OSB is better than that of using the traditional IFCB products. Thus, the proposed pseudorange and phase time-variant mixed OSB concept and the new undifferenced model coupled with satellite clock offsets are reliable, applicable, and effective in multi-frequency PPP.

1. Introduction

With GNSS modernization and development, the traditional GPS-only dual-frequency model is gradually transitioning to multi-GNSS (GPS, BDS, GLONASS, Galileo, and regional system) and multi-frequency tendencies. The GPS Block IIF satellites can transmit L1, L2, and L5 signals [1]. The global BeiDou navigation satellite system (BDS-3) satellites can broadcast the B1C, B2a, B2b, and B2 (B2a + B2b) signals as well as the legacy B1I and B3I signals [2]. Similarly, the new GLONASS satellites have begun to transmit code division multiple access (CDMA) G1a, G2a, and G3 signal and frequency division multiple access (FDMA) G1 and G2 signals [3]. Galileo can transmit E1, E5a, E5b, E5, and E6 signals [4]. The multi-frequency observations can be beneficial in precise point positioning (PPP), cycle slip detection, precise clock offset estimation (PCE), and so on [5,6,7].

Normally, the legacy precise satellite clock offsets released by GNSS analysis centers (ACs) are generated by a dual-frequency ionosphere-free (IF) combination model [8]. For example, the recommended basic frequency pair for GPS, Galileo, and BDS are L1/L2, E1/E5a, and B1I/B3I signals [1,9,10]. Due to the characteristics of multi-frequency multi-channel GNSS signals, there is an inter-frequency bias between different frequencies, and there is an intra-frequency deviation between different channels of the same frequency. The code and phase time-variant hardware delay will be absorbed into satellite clock offsets, which leads to the inconsistency of the precise satellite clock offsets estimated by different frequencies [11]. This inconsistency is known as the inter-frequency clock bias (IFCB) [1]. The IFCB can be divided into code IFCB (CIFCB) and phase IFCB (PIFCB) [12]. The CIFCB is a linear combination of code hardware delays, and it is also called the differential code bias (DCB) or pseudorange observable-specific signal bias (OSB) [13]. The DCB is usually determined by the carrier-to-code leveling (CCL) method [14,15]. Most International GNSS Service (IGS) and International GNSS Monitoring and Assessment System (iGMAS) ACs such as Centre for Orbit Determination in Europe (CODE), Deutsches Zentrum für Luft- und Raumfahrt (DLR), Chinses Academy of Sciences (CAS), WhuHan University (WHU), and Shanghai Astronomical Observatory (SHA) use the CCL method to obtain DCB [16,17,18,19]. The PIFCB is a linear combination of phase time-variant hardware delays, which is usually generated by epoch-differenced (ED) geometry-free (GF) ionosphere-free (IF) linear combination (GFIF) [12]. At present, IGS and iGMAS ACs have not yet officially released IFCB, and the OSB files obtained from GNSS AC mainly release pseudorange and phase OSB [3,20]. Nevertheless, many references have conducted research on it. Pan, et al. [21] used the ED GFIF model to obtain GPS L5 PIFCB, and researched characteristics to model and predict it [22]. Similarly, Gong, et al. [23] used the ED GFIF model to generate the PIFCB and analyzed its periodic characteristics to develop the corresponding prediction model. Su, et al. [24] converted the PIFCB obtained by the above method into carrier-phase OSB, which greatly simplified the complexity of the client algorithm. Although this method has high computational efficiency, it has certain theoretical flaws. However, the correct method assumes that the PIFCB of the same satellite at different receivers has good consistency and the receiver PIFCB can be combined into satellite PIFCB. Actually, receiver hardware delays are greatly affected by the ground environment variation [25]. If the receiver hardware delays change, it will have a negative impact on the satellite PIFCB estimation. The IFCB can also be calculated by estimating the satellite clock offset for all frequencies. Guo and Geng [26] provided the GPS L1/L2 and L5-only types of satellite clock products for the triple-frequency PPP. Jiao, et al. [27] established a multi-frequency precise satellite clock offset estimation model to obtain clock offsets and the corresponding IFCB. Fan, et al. [28] first modeled the IFCB and then brought the model coefficients into the undifferenced model to obtain the corresponding IFCB model coefficients, which can improve IFCB estimation accuracy. Geng, et al. [29] merged IFCB into phase OSB at the third frequency, which simplifies the PPPAR process. But, it is not very convenient for multi-frequency float PPP. However, these methods require the estimation of the ambiguity parameters for the basic frequency as well as the third frequency. The ambiguity parameters will increase with the increase in observations, and its corresponding computational burden is enlarged [30]. Therefore, it is very meaningful to develop an IFCB estimation method with high computational efficiency and accuracy.

With this background, we developed the new undifferenced IFCB estimation model considering the receiver time-variant biases, and further converted it into code and phase time-variant mixed OSB to improve multi-frequency PPP (MFPPP) performance. The new model avoids the negative impact of receiver time-variant biases on the satellite mixed OSB, and it has good accuracy and low noise. First, the mathematical methods of the traditional ED GFIF-phase time-variant OSB model and new undifferenced mixed OSB model coupled with satellite clock offsets and an analysis of their characteristics are presented in detail. Then, the processing datasets and the corresponding processing strategies are presented. Naturally, the experimental analysis about the advantages of the new undifferenced mixed OSB model and its benefits on MFPPP are introduced in the following section. Finally, the conclusions are summarized.

2. Methods

This section begins with the GNSS pseudorange and phase observation models. The traditional epoch-differenced phase-only GFIF OSB estimation model and the new undifferenced models coupled with satellite clock offsets are developed in detail. The IFCB is converted to OSB type for simplifying the client calculation.

2.1. General Observation Models and OSB Concept

The GNSS original observation equations with the time-variant parts of pseudorange and phase hardware delays can be written as follows [31]:

$$\left\{\begin{array}{l}{p}_{r,j}^s(t)={\rho }_r^s(t)+d{t}_r(t)-d{t}^s(t)+b_{r,j}-b_j^s+b_{r,j}(t)-b_j^s(t)+T_r^s(t)+I_{r,j}^s(t)+{\epsilon }_{p_j}(t)\\ {\varphi }_{r,j}^s(t)={\rho }_r^s(t)+d{t}_r(t)-d{t}^s(t)+B_{r,j}-B_j^s+B_{r,j}(t)-B_j^s(t)+T_r^s(t)-I_{r,j}^s(t)+N_{r,j}^s+{\epsilon }_{{\varphi }_j}(t)\end{array}\right.$$

(1)

where $p_{r,j}^s(t)$ and ${\varphi }_{r,j}^s(t)$ are the pseudorange and phase observations on frequency j in units of meter, respectively; ${\rho }_r^s(t)$ is the geometrical range from satellite s to receiver r at epoch t; $d{t}_r(t)$ and $d{t}^s(t)$ denote the original receiver and satellite clock offsets; and $b_{r,j}$ and $b_j^s$ represent the receiver and satellite time-invariant uncalibrated code delays (UCDs), while $b_{r,j}(t)$ and $b_j^s(t)$ represent the corresponding time-variant parts; similarly, $B_{r,j}$ and $B_j^s$ depict the time-invariant uncalibrated phase delays (UPDs) and $B_{r,j}(t)$ and $B_j^s(t)$ are the time-dependent parts; $T_r^s(t)$ depicts the slant tropospheric delay; $I_{r,j}^s(t)$ is the slant ionospheric delay; $N_{r,j}^s$ denotes the integer ambiguity; and ${\epsilon }_{p_j}(t)$ and ${\epsilon }_{{\varphi }_j}(t)$ are the pseudorange and phase observation noises containing multipath and unmodeled errors.

The IGS and iGMAS ACs typically employ dual-frequency undifferenced IF combination observations to generate satellite clock offsets [31,32]:

$$\left\{\begin{array}{l}{p}_{r,I{F}_{ij}}^s(t)={\rho }_r^s(t)+d{t}_{r,I{F}_{ij}}(t)-d{t}_{I{F}_{ij}}^s(t)+T_r^s(t)+{\overline{\epsilon }}_{p_{I{F}_{ij}}}(t)\\ {\varphi }_{r,I{F}_{ij}}^s(t)={\rho }_r^s(t)+d{t}_{r,I{F}_{ij}}(t)-d{t}_{I{F}_{ij}}^s(t)+T_r^s(t)+{\overline{N}}_{r,I{F}_{ij}}^s+{\epsilon }_{{\varphi }_{I{F}_{ij}}}(t)\end{array}\right.$$

(2)

with

$$\left\{\begin{array}{l}d{t}_{r,I{F}_{ij}}(t)=d{t}_r(t)+b_{r,I{F}_{ij}}+B_{r,I{F}_{ij}}(t)\\ d{t}_{I{F}_{ij}}^s(t)=d{t}^s(t)+b_{I{F}_{ij}}^s+B_{I{F}_{ij}}^s(t)\\ {\overline{\epsilon }}_{p_{I{F}_{ij}}}(t)={\epsilon }_{p_{I{F}_{ij}}}(t)+b_{r,I{F}_{ij}}(t)-b_{I{F}_{ij}}^s(t)-B_{r,I{F}_{ij}}(t)+B_{I{F}_{ij}}^s(t)\\ {\overline{N}}_{r,I{F}_{ij}}^s=N_{r,I{F}_{ij}}^s+B_{r,I{F}_{ij}}-B_{I{F}_{ij}}^s-b_{r,I{F}_{ij}}+b_{I{F}_{ij}}^s\end{array}\right.$$

(3)

where ${(\cdot )}_{I{F}_{ij}}={\alpha }_{ij}\cdot {(\cdot )}_i+{\beta }_{ij}\cdot {(\cdot )}_j$, ${\alpha }_{ij}={\displaystyle \frac{f_i^2}{f_i^2-f_j^2}}$, ${\beta }_{ij}={\displaystyle \frac{-f_j^2}{f_i^2-f_j^2}}$, $I{F}_{ij}$ denotes the dual-frequency IF-combined operation; i and j are the IGS/iGMAS recommended frequency, where GPS is L1/L2, BDS is B1I/B3I, and Galileo is E1/E5a; $d{t}_{r,I{F}_{ij}}(t)$ and $d{t}_{I{F}_{ij}}^s(t)$ denote receiver and satellite IF-combined clock offsets; and ${\overline{\epsilon }}_{p_{I{F}_{ij}}}(t)$ and ${\overline{N}}_{r,I{F}_{ij}}^s$ are the recombined pseudorange observation noises and IF ambiguity, respectively. The satellite clock offsets published by IGS and iGMAS ACs will absorb the UCD, time-variant UPD, and reference clock, which can be written as follows [11]:

$$d{\overline{t}}_{I{F}_{ij}}^s(t)=d{t}^s(t)+b_{I{F}_{ij}}^s+B_{I{F}_{ij}}^s(t)+d{t}_{ref}=d{t}_{I{F}_{ij}}^s(t)+d{t}_{ref}$$

(4)

where $d{\overline{t}}_{I{F}_{ij}}^s(t)$ is the actual estimated clock offsets and $d{t}_{ref}$ is the reference clock.

The satellite clock offsets absorbed with the hardware delays make the measurement inconsistent at the specific frequency. The satellite clock offsets for the specific frequency can be expressed as follows:

$$\left\{\begin{array}{l}d{\overline{t}}_i^s(t)=d{t}_{I{F}_{ij}}^s(t)+{\delta }_{p_i}+{\delta }_{{\varphi }_i}(t)+d{t}_{ref}=d{t}_{I{F}_{ij}}^s(t)+{\delta }_{M_i}^s(t)+d{t}_{ref}\\ d{\overline{t}}_j^s(t)=d{t}_{I{F}_{ij}}^s(t)+{\delta }_{p_j}+{\delta }_{{\varphi }_j}(t)+d{t}_{ref}=d{t}_{I{F}_{ij}}^s(t)+{\delta }_{M_j}^s(t)+d{t}_{ref}\\ d{\overline{t}}_k^s(t)=d{t}_{I{F}_{ij}}^s(t)+{\delta }_{p_k}+{\delta }_{{\varphi }_k}(t)+d{t}_{ref}=d{t}_{I{F}_{ij}}^s(t)+{\delta }_{M_k}^s(t)+d{t}_{ref}\end{array}\right.$$

(5)

where $d{\overline{t}}_i^s(t)$, $d{\overline{t}}_j^s(t)$, and $d{\overline{t}}_k^s(t)$ denote the satellite clock offsets for the basic ith, jth, and the third kth frequency signals. ${\delta }_{p_i}$ and ${\delta }_{p_j}$ denote the pseudorange OSB for the baisc frequency signals; ${\delta }_{p_k}$ denotes the pseudorange OSB for the third frequency signal; and, in analogy, ${\delta }_{{\varphi }_i}(t)$, ${\delta }_{{\varphi }_j}(t)$, and ${\delta }_{{\varphi }_k}(t)$ are carrier-phase time-variant OSB for the basic and the third frequency. Normally, the definition of satellite IFCB is the bias between the satellite clock offsets for the different frequency, which usually contains CIFCB (also known as DCB or pseudorange OSB) and PIFCB (also known as the phase time-variant OSB). The ${\delta }_{p_k}$ and ${\delta }_{{\varphi }_k}(t)$ can be combined into the mixed IFCB with OSB type ${\delta }_{M_k}^s(t)$.

Generally, the pseudorange OSB is usually determined by the CCL method. However, the phase time-variant OSB for the basic frequency cannot be calculated by the methods of obtaining pseudorange OSB. As shown in (2), it will be absorbed into the other parameters such as clock offsets and residuals. As for the dual-frequency PPP, the phase time-variant OSB for the basic frequency will not have any negative impact. Because this bias is recombined into other parameters, the server algorithms for satellite clock offsets and the client PPP algorithms are strictly self-consistent. Hence, the carrier-phase time-variant OSB values at the basic frequency signals can be set to zero. The carrier-phase time-variant OSB values at the basic and the third frequency signals can be written as follows.

$$\left\{\begin{array}{l}\left.\begin{array}{l}{\delta }_{{\varphi }_i}^s(t)=0\\ {\delta }_{{\varphi }_j}^s(t)=0\end{array}\right\}baseline frequency\\ {\delta }_{{\varphi }_k}^s(t)={\alpha }_{ij}\cdot \left(1-{\mu }_k/{\mu }_j\right)\cdot B_i^s(t)-{\beta }_{ij}\cdot \left({\mu }_k-1\right)\cdot B_j^s(t)-B_k^s(t)\end{array}\right.$$

(6)

Equation (6) indicates that the carrier-phase time-variant OSB at different frequencies, and ${\mu }_i$ and ${\mu }_k$ denote the frequency-dependent multiplier factor.

2.2. Traditional GFIF Phase Time-Variant OSB Model

The phase time-variant OSB is commonly generated by the ED GFIF model. The GFIF combination between the third and the basic frequency can be written as follows [33]:

$${\varphi }_{r,GFI{F}_{ijk}}^s(t)={\varphi }_{r,I{F}_{ik}}^s(t)-{\varphi }_{r,I{F}_{ij}}^s(t)={\delta }_{r,{\varphi }_{I{F}_{ik}}}(t)-{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)+{\overline{N}}_{r,GFI{F}_{ijk}}^s$$

(7)

where ${\overline{N}}_{r,GFI{F}_{ijk}}^s$ denotes the differential IF-phase ambiguity. If there is no cycle slip, and the epoch difference strategy can be used to eliminate ambiguity, which can be written as follows:

$$\mathsf{\Delta }{\delta }_{r,{\varphi }_{I{F}_{ik}}}(t,t-1)-\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)={\varphi }_{r,GFI{F}_{ijk}}^s(t)-{\varphi }_{r,GFI{F}_{ijk}}^s(t-1)$$

(8)

where $\mathsf{\Delta }{\delta }_{r,{\varphi }_{I{F}_{ik}}}(t,t-1)$ and $\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)$ denote the differential IF PIFCB.

Because the PIFCB of the same satellite at different stations has good consistency, the receiver PIFCB can be combined into satellite PIFCB [33], which reads as follows:

$$\mathsf{\Delta }{\delta }_{r,{\varphi }_{I{F}_{ik}}}^s(t,t-1)={\varphi }_{r,GFI{F}_{ijk}}^s(t)-{\varphi }_{r,GFI{F}_{ijk}}^s(t-1)$$

(9)

where $\mathsf{\Delta }{\delta }_{r,{\varphi }_{I{F}_{ik}}}^s(t,t-1)$ denotes the corresponding IFCB for the paired receiver and satellite.

Assuming that satellite s can be tracked by m common viewing stations in the observational network at epoch t and t − 1, the epoch-differenced PIFCB is expressed as follows:

$$\left\{\begin{array}{l}\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)={\displaystyle \frac{{\displaystyle \sum _{r=1}^m\left[\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)\cdot {\omega }_r^s(t,t-1)\right]}}{{\displaystyle \sum _{r=1}^m\left[{\omega }_r^s(t,t-1)\right]}}}\\ {\omega }_r^s(t,t-1)=\left\{\begin{array}{l}\sin (El{e}_r^s(t,t-1)),El{e}_r^s(t,t-1)<40°\\ 1,El{e}_r^s(t,t-1)\ge 40°\end{array}\right.\end{array}\right.$$

(10)

where $El{e}_r^s(t,t-1)$ represents the average elevation at epoch t and t − 1. The PIFCB can be accumulated through epoch-differenced values, as shown below:

$${\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)={\delta }_{{\varphi }_{I{F}_{ik}}}^s(t_0)+{\displaystyle \sum _{t=t_0+1}^t\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)}$$

(11)

with

$${\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)=B_{I{F}_{i\mathrm{k}}}^s(t)-B_{I{F}_{ij}}^s(t)$$

(12)

where t₀ is the initial epoch time and the initial undifferenced PIFCB, ${\delta }_{{\varphi }_{I{F}_{ik}}}^s(t_0)$, is usually set to zero. $\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)$ is the ED IFCB between epoch $t$ and $t-1$. ${\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)$ can be accumulated by the $\mathsf{\Delta }{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t,t-1)$ from the initial epoch time. The ${\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)$ is the IF-combined PIFCB, which is difficult to be corrected in the complex PPP models such as non-basic frequency PPP, undifferenced PPP, MFPPP, and so on. According to the SINEX_BIAS Version 1.0 [34], the bias in the form of OSB type can better adapt to multi-frequency observations and the complex PPP algorithms. Combining (6) and (12), the IF-combined PIFCB can be transformed into the phase time-variant OSB, which can be expressed as follows.

$${\delta }_{{\varphi }_k}^s(t)={\beta }_{ik}^{-1}\cdot {\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)$$

(13)

The phase time-variant OSB for basic and third frequency can generated by using Equations (6) and (13).

2.3. Undifferenced Mixed OSB Model Coupled with Satellite Clock Offsets

Through the comprehensive analysis of the above calculation process, it can be found that the ED GFIF model has three major flaws. The first is that it neglects the receiver PIFCB. It assumes that the PIFCB of the same satellite at different receivers has good consistency and the receiver PIFCB can be combined into satellite PIFCB. The second is that it neglects the relationship between two adjacent epochs because it obtains PIFCB by averaging and accumulating. The third and most important is the MFPPP need to correct both pseudorange OSB and PIFCB, which makes the client algorithms more complex. For traditional methods, two bias corrections are required to achieve MFPPP, which makes the client algorithms more complex. To avoid the above defects, the undifferenced pseudorange and phase time-variant mixed OSB model coupled with satellite clock offsets can be expressed as follows:

$$\left\{\begin{array}{l}{p}_{r,I{F}_{ik}}^s(t)={\rho }_r^s(t)+d{t}_{r,I{F}_{ik}}(t)-d{t}_{I{F}_{ij}}^s(t)-{\delta }_{I{F}_{ik}}^s(t)+T_r^s(t)+{\overline{\epsilon }}_{p_{I{F}_{ik}}}(t)\\ {\varphi }_{r,I{F}_{ik}}^s(t)={\rho }_r^s(t)+d{t}_{r,I{F}_{ik}}(t)-d{t}_{I{F}_{ij}}^s(t)-{\delta }_{I{F}_{ik}}^s(t)+T_r^s(t)+{\overline{N}}_{r,I{F}_{ik}}^s+{\epsilon }_{{\varphi }_{I{F}_{ik}}}(t)\end{array}\right.$$

(14)

where ${\delta }_{I{F}_{ik}}^s(t)$ is the satellite pseudorange and phase time-variant mixed IFCB.

The IGS legacy clocks are fixed throughout the undifferenced mixed OSB model; hence, we further rewrite (14) by introducing (4) as follows:

$$\left\{\begin{array}{l}{p}_{r,I{F}_{ik}}^s(t)+d{\overline{t}}_{I{F}_{ij}}^s(t)={\rho }_r^s(t)+d{\overline{t}}_{r,I{F}_{ik}}(t)-{\delta }_{I{F}_{ik}}^s(t)+T_r^s(t)+{\overline{\epsilon }}_{p_{I{F}_{ik}}}(t)\\ {\varphi }_{r,I{F}_{ik}}^s(t)+d{\overline{t}}_{I{F}_{ij}}^s(t)={\rho }_r^s(t)+d{\overline{t}}_{r,I{F}_{ik}}(t)-{\delta }_{I{F}_{ik}}^s(t)+T_r^s(t)+{\overline{N}}_{r,I{F}_{ik}}^s+{\epsilon }_{{\varphi }_{I{F}_{ik}}}(t)\end{array}\right.$$

(15)

where $d{\overline{t}}_{I{F}_{ij}}^s(t)$ is IGS legacy clock offsets and $d{\overline{t}}_{r,I{F}_{ik}}(t)$ is the recombined receiver clock offsets, which absorb the reference clock $d{t}_{ref}$. Obviously, there is a rank-deficient problem between $d{\overline{t}}_{r,I{F}_{ik}}(t)$ and ${\delta }_{I{F}_{ik}}^s(t)$. To avoid the negative impact of receiver hardware delay on the satellite IFCB estimation, the zero-mean condition (ZMC) for a constellation is used to separate the receiver and satellite IFCB.

Naturally, there are two technical routes to obtain satellite IFCB according to whether it is with or without pseudorange OSB correction for basic frequency. If the basic frequency observations are corrected by pseudorange OSB, the third satellite IFCB will not contain satellite pseudorange OSB, which can be expressed as follows.

$${\delta }_{M_{I{F}_{ik}}}^s(t)={\beta }_{ik}\cdot \left(b_k^s-b_i^s\right)+B_{I{F}_{ik}}^s(t)-B_{I{F}_{ij}}^s(t)={\delta }_{p_{I{F}_{ik}}}^s+{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)$$

(16)

The IF-combined mixed IFCB ${\delta }_{M_{I{F}_{ik}}}^s(t)$ can be further converted into OSB type, which is expressed as follows.

$$\left\{\begin{array}{l}\left.\begin{array}{l}{\delta }_{M_i}^s(t)={\delta }_{p_i}={\beta }_{ij}\cdot \left(b_i^s-b_j^s\right)\\ {\delta }_{M_j}^s(t)={\delta }_{p_j}=-{\alpha }_{ij}\cdot \left(b_i^s-b_j^s\right)\end{array}\right\}\mathrm{basic} \mathrm{frequency} \mathrm{pair}\\ {\delta }_{M_k}^s(t)=b_k^s-b_i^s\\ +{\alpha }_{ij}\cdot \left(1-{\mu }_k/{\mu }_i\right)\cdot B_i^s(t)-{\beta }_{ij}\cdot \left({\mu }_k-1\right)\cdot B_j^s(t)-B_k^s(t)\\ ={\delta }_{p_k}+{\delta }_{{\varphi }_k}(t)={\beta }_{ik}^{-1}\cdot {\delta }_{I{F}_{ik}}^s(t)\end{array}\right.$$

(17)

If the basic frequency observations are not corrected by pseudorange OSB, the estimated mixed IFCB will absorb satellite pseudorange OSB. Hence, the mixed IFCB for the third frequency can be given as follows.

$${\delta }_{M_{I{F}_{ik}}}^s(t)=b_{I{F}_{ik}}^s-b_{I{F}_{ij}}^s+B_{I{F}_{ik}}^s(t)-B_{I{F}_{ij}}^s(t)={\delta }_{p_{I{F}_{ik}}}^s+{\delta }_{{\varphi }_{I{F}_{ik}}}^s(t)$$

(18)

The mixed IFCB ${\delta }_{M_{I{F}_{ik}}}^s(t)$ can be further converted into OSB type, which is expressed as follows.

$$\left\{\begin{array}{l}\left.\begin{array}{l}{\delta }_{M_i}^s(t)=0\\ {\delta }_{M_j}^s(t)=0\end{array}\right\}\mathrm{basic} \mathrm{frequency} \mathrm{pair}\\ {\delta }_{M_k}^s(t)={\alpha }_{ij}\cdot \left(1-{\mu }_k/{\mu }_i\right)\cdot b_i^s-{\beta }_{ij}\cdot \left({\mu }_k-1\right)\cdot b_j^s-b_k^s\\ +{\alpha }_{ij}\cdot \left(1-{\mu }_k/{\mu }_i\right)\cdot B_i^s(t)-{\beta }_{ij}\cdot \left({\mu }_k-1\right)\cdot B_j^s(t)-B_k^s(t)\\ ={\delta }_{p_k}+{\delta }_{{\varphi }_k}(t)={\beta }_{ik}^{-1}\cdot {\delta }_{I{F}_{ik}}^s(t)\end{array}\right.$$

(19)

Equation (19) indicates that the OSB for the third frequency is composed of the time-invariant UCD and time-variant UPD, which differs from Equation (17) in including and excluding pseudorange OSB for basic frequency.

In summary, the traditional IFCB model is based on three assumptions: (1) the PIFCB of the same satellite at different stations has good consistency and the receiver PIFCB can be combined into satellite PIFCB; (2) pseudorange and phase IFCB are determined by CCL and ED GFIF models, respectively; (3) the PIFCB is determined by averaging multiple GNSS tracking stations. These assumptions also reflect some defects. If the time-varying phase delays of different GNSS receivers are inconsistent, it will have a negative impact on the satellite PIFCB estimation. It is undeniable that traditional methods are more efficient than undifferenced models in terms of computational efficiency because they do not require the calculation of ambiguity. However, the new method only has one set of ambiguity parameters, so the computational efficiency is still very high. In addition, the pseudorange OSB is required for basic frequency observations, and both pseudorange OSB and phase time-variant OSB corrections are required for the third frequency observation in MFPPP. Obviously, this way is too complicated. In addition, the ACs such as CODE, DLR, CAS, WHU, and SHA basically use CCL to generate satellite pseudorange OSB products. This will lead to the poor self-consistency of those biases. The new method can avoid complex bias correction. We do not need to make pseudorange OSB corrections for the basic frequency observations. We only need to make the pseudorange and phase time-variant mixed OSB corrections for the third frequency observation, which significantly optimizes the terminal positioning algorithm. The method of (19) is used to obtain the mixed OSB.

3. Experiment Setup

Approximately, the 130 stations collected from IGS MGEX for the period day of year (DOY) from 001 to 030, 2021, are used to obtain GPS, BDS-3, and Galileo satellite mixed OSB. The distribution of the selected GNSS stations for OSB estimation can be found in Figure 1.

The traditional phase time-variant ED OSB model is relatively simple. It can be calculated according to Equations (8)–(14), and its description will not be repeated again [33]. The undifferenced mixed OSB estimation processing strategies are illustrated in

Table 1. Data processing strategy for undifferenced mixed OSB estimation.

Items	Strategies
Basic frequency pair	GPS L1/L2, BDS-3 B1I/B3I, Galileo E1/E5a
Estimated frequency	GPS L5, BDS-3 B1C/B2a, Galileo E5b/E5/E6
Elevation cutoff	7 degrees
Weighting	Elevation weight [sin(elevation)]
Filter type	Bidirectional Kalman filter (forward + backward)
Satellite orbit and clock offsets	Fixed to WHU MGEX final precise products
Receiver coordinate	Fixed to IGS SINEX-file coordinates
Satellite IFCB	Estimated as random walk (104 m2/s) [5]
Reference IFCB	Zero-mean condition
Receiver clock offsets	Estimated as white noise
Tropospheric delay	Dry part: modified Hopfield model Wet part: estimated as random walk (10−9 m2/s)
Ionospheric delay	Eliminated first order by IF observations
Satellite and receiver antenna	IGS MGEX values (igs14.atx)
Phase windup effect	Corrected [39]
Relativistic effect	Corrected [40]
Earth rotation	Corrected [40]
Tide effect	Solid earth, pole, and ocean tide [40]
Ambiguity	Estimated as float
DCB	Parameter recombination

. The GPS L1/L2, BDS-3 B1I/B3I, and Galileo E1/E5a are treated as the basic frequency, while GPS L5, BDS-3 B1C and B2a, and Galileo E5, E5b, and E6 are used to estimate satellite IFCB. The satellite orbits and clock offsets are fixed to the WHU final products, and the station coordinates are fixed to IGS weekly SINEX-file coordinates [3,35]. The ZMC method for a constellation is used to eliminate rank deficiency between receiver and satellite IFCB. As for tropospheric delay, the dry parts are corrected by the modified Hopfield model based on the Global Pressure and Temperature 3 (GPT3) and Vienna Mapping Functions 3 (VMF3) models, while the wet parts are estimated by setting the parameter based on the wet part of VMF3 [36,37]. The Phase Center Offset (PCO) and Phase Center Variations (PCVs) for GPS, BDS-3, and Galileo multi-frequency observations are corrected using the IGS antenna file (igs14.atx). In the bidirectional Kalman filter (forward and backward) of IFCB processing, the ambiguities are estimated as float constants, the ZWD is estimated as a random walk process, and the receiver clock is estimated as white noise. To avoid gross errors, and retain the IFCB’s original characteristics, the satellite IFCB is estimated as a random walk process and the receiver IFCB is estimated as white noise. In the forward Kalman Filter of PPP processing, the processing strategy is basically the same as that of mixed OSB estimation. For more details for PPP processing strategies, please refer to [5,38].

4. Results

This section presents the result of the OSB obtained from the ED GFIF model and the new undifferenced model coupled with satellite clock offsets. First, the mixed OSB service system framework is introduced based on the above theoretical analysis. Then, two OSB solving methods are compared in terms of their characteristics and MFPPP performances.

4.1. Pseudorange and Phase Time-Variant Mixed OSB Service System

We established the multi-GNSS pseudorange and phase time-variant mixed OSB service system based on the above theoretical discussions. Figure 2 illustrates the flowchart of the mixed OSB service system. First, the data download module automatically downloads GNSS observations, precise orbits, and clock offsets from IGS MGEX and iGMAS. The IGS MGEX observation values and the precise orbit clock products of the WHU IGS and iGMAS analysis center are downloaded to generate the mixed OSB in this manuscript. If the precise orbit and clock offsets are missing, the ED GFIF methods are used as a backup algorithm for obtaining the phase-only time-variant IFCB. The phase-only time-variant IFCB can be transformed into OSB type. Typically, MFPPP requires both phase-only time-variant IFCB and DCB corrections [24,41]. We synthesize the additional DCB products and phase-only time-variant OSB into the mixed OSB on the server side, which can simultaneously correct both pseudorange and phase observations. This way, it can reduce the complexity of client algorithms. If the precise orbit and clock offsets are not missing, the new undifferenced IFCB model coupled with satellite clock offsets methods is used to generate the mixed IFCB. Furthermore, this IFCB can be transformed into OSB. The client only needs to correct this bias to achieve MFPPP. Note that, unless otherwise specified, the mixed OSB in the following text refers to the pseudorange and phase time-variant mixed OSB.

4.2. Characteristics of the Mixed OSB

The mixed OSB obtained from Figure 2 contains pseudorange and phase time-variant OSB. The zero-mean constraint of satellite mixed OSB obtained from the two schemes may not be consistent, which introduces systematic errors in comparing the results. In order to better illustrate the satellite mixed OSB situation for two schemes, it is necessary to convert all the mixed OSB to the same reference datum. To simplify the description, Scheme 1 and Scheme 2 represent the mixed OSB obtained from the ED GFIF model and undifferenced model, respectively.

To analyze the characteristics at each frequency, Figure 3, Figure 4 and Figure 5 show the time series of the mixed OSB for GPS L5, BDS-3 B1C, B2a, and Galileo E5b, E5, and E6 signals on DOY 011 in 2021. The time resolution is 30 s. Figure 6 shows the corresponding mixed OSB amplitudes, in which x-axes depict the PRN of each satellite (for clarity and simplicity, specific numbers are not labeled). It is not difficult to find that the OSB for GPS BLOCK-IIF and BDS-3 signals have obvious amplitudes, while the GPS BLOCK-IIA, BLOCK-IIR, and Galileo OSB amplitudes are relatively small. The amplitude of GPS BLOCK-IIF L5 OSB is in the decimeter range, while the amplitude of BDS-2 B2I and BDS-3 B1C and B2a OSB is in the centimeter range. The OSB for Galileo E5b, Galileo E5, and Galileo E6 is around 1~2 cm. The OSB for BDS-3 B1C and BDS-3 B2a is around 1~3 cm. We all know that the accuracy of phase observation values is about 0.3 cm, so the IFCB correction should be considered in PPP. By comparing the OSB of Scheme 1 and Scheme 2, it can be seen that the trend of their changes is basically the same. The results further proved that the new undifferenced method is reasonable. But, undeniably, there are some differences in some satellites between the two schemes, especially G18, G23, C44, C45, and C46. There are four reasons for this phenomenon: (1) Scheme 1 assumes that the phase time-variant OSB of the same satellite at different stations has good consistency and the receiver phase time-variant OSB can be combined into satellite phase time-variant OSB. As for Scheme 2, the zero-mean condition for all satellites is introduced to separate the receiver and satellite OSB, which avoids the negative of the receiver on satellite OSB estimation. (2) Due to the use of ZMC in Scheme 2, all OSBs obtained from Scheme 2 will absorb the datum for a satellite constellation, while Scheme 1 will not. In other words, the OSBs obtained from Scheme 2 are similar to estimating satellite clock offsets, which will absorb the OSB for the reference satellite. Therefore, there are differences between the two schemes. (3) The number of BDS-3 C44/C45/C46 observations is relatively small, which affects the estimation accuracy [5]. (4) From a theoretical analysis of the formula, some parameters of the undifferenced model have certain correlations. The undifferenced model requires satellite orbit, satellite clock offsets, PCO, and PCV to achieve the estimation of OSB. Inevitably, the accuracy of these products will affect the estimation of undifferenced OSB. Therefore, the results with the new method suffer from irrational systematic errors. However, the undifferenced OSB model has better coupling with orbit and clock errors. The corresponding self-consistency between the satellite products is stronger, which will be beneficial for PPP.

Those systematic errors introduced by precise products can be eliminated through the ED method. Therefore, the ED method can reflect the noise level. To further compare the advantages and disadvantages of the two schemes, Figure 7, Figure 8 and Figure 9 depict the corresponding mixed ED OSB for GPS L5, BDS-3 B1C, B2a, and Galileo E5b, E5, and E6 signals. It is not hard to observe that the ED OSB of Scheme 2 is better than that of Scheme 1 in terms of stability and accuracy. The noise of Scheme 1 is significantly higher than that of Scheme 2. To more objectively illustrate the performance of Schemes 1 and 2,

Table 2. Statistical results and improvement of epoch-differenced and undifferenced OSB for GPS L5, BDS-3 B1C, B2a, and Galileo E5b, E5, and E6 signals in January 2021.

Type (Unit: mm)		GPS	BDS-3		Galileo
		L5	B1C	B2a	E5b	E5	E6
STD UD OSB	Scheme 1	18.6	5.7	8.8	1.5	2.3	2.7
	Scheme 2	17.2	1.8	5.4	1.3	0.9	2.1
	Improvement	7.5%	68.4%	38.6%	13.3%	60.9%	22.2%
RMS UD OSB	Scheme 1	20.3	5.9	9.2	1.7	2.3	2.9
	Scheme 2	18.4	2.0	6.0	1.5	1.0	2.3
	Improvement	9.4%	66.1%	34.8%	11.8%	56.5%	20.7%
STD ED OSB	Scheme 1	0.5	0.5	0.8	0.4	0.2	0.4
	Scheme 2	0.2	0.1	0.1	0.1	0.1	0.1
	Improvement	60.0%	80.0%	87.5%	75.0%	50.0%	75.0%
RMS ED OSB	Scheme 1	0.5	0.6	0.9	0.4	0.3	0.5
	Scheme 2	0.2	0.1	0.1	0.1	0.1	0.1
	Improvement	60.0%	83.3%	88.9%	75.0%	66.7%	80.0%

illustrates the statistical results and improvement of Scheme 2 in January 2021. The STD improvement of the undifferenced OSB is 7.5, 68.4, 38.6, 13.3, 60.9, and 22.2%, and the RMS can be improved by 9.4, 66.1, 34.8, 11.8, 56.5, and 20.7% with respect to GPS L5, BDS-3 B1C, B2a, and Galileo E5b, E5, and E6 signals. Moreover, the STD improvement of the ED OSB is 60.0, 80.0, 87.5, 75.0, 50.0, and 75.0%, and the RMS can be improved by 60.0, 83.3, 88.9, 75.0, 66.7, and 80.0%C. There are significant differences in improvement rates among different satellite systems or frequencies. There are two main reasons: (1) the undifferenced mixed OSB for different satellite systems are frequencies that have different scales, and (2) the number of observations for different signals is different, so there are differences in the improvement.

Overall, these discoveries illustrate that the stability and accuracy of the new undifferenced OSB estimation algorithm compared to the ED IFCB algorithm.

Notably, it can be found that there is a systematic periodic bias in the OSB obtained from Scheme 1 at BDS B2a and Galileo E6 signals from Figure 4 and Figure 5. Some of the literature suggests that this systematic periodic bias is related to the number of available observations at the GNSS tracking station. Due to the small number of MGEX stations capable of tracking BDS-3 B2a and Galileo E6 signals, the accuracy and reliability of the estimated OSB are poor [24]. The OSB estimation is affected by the receiver antenna environment and receiver and antenna quality, which results in periodic errors [1]. Unfortunately, these explanations cannot explain this systematic error at the BDS B2a and Galileo E6 signals. However, it can be found that the mixed OSB value obtained from Scheme 2 does not have significant periodic bias. By comparing the two algorithms, it can be seen that the epoch-differenced algorithm is based on the assumption that the IFCB of the same satellite at different stations has good consistency and the receiver phase time-variant OSB can be combined into satellite phase time-variant OSB. However, it is necessary to separate the satellite and receiver phase time-variant OSB in the new algorithm. To confirm that this small periodic bias is caused by some receiver phase time-variant OSB, Figure 10 shows the receiver OSB for all receiver types in the observation network. It can be seen that the JAVAD TRE-3 receiver OSB for BDS-3 B2a and Galileo E6 exhibits significant periodic variations; the OSB of other receivers is relatively stable and does not have significant deviations. Figure 11 further depicts the receiver OSB for multiple stations with the JAVAD TRE-3 receiver. This is just more evidence for the theory that JAVAD TRE-3 receiver hardware delay exhibits periodic variations at the BDS-3 B2a and Galileo E6 signals. Scheme 1 ignores the receiver phase time-variant OSB, which negatively impacts the estimation of the satellite OSB. On the contrary, the new undifferenced method uses zero-mean constraints to separate the satellite and receiver OSB, which effectively solves this problem.

4.3. The Advantages of the Mixed OSB in MFPPP

To evaluate the benefits of the mixed OSB on MFPPP, we carried out GPS L1 + L2 + L5, BDS-3 B1I + B3I + B2a, BDS-3 B1I + B3I + B1C + B2a, and Galileo E1 + E5a + E5b + E5 + E6 multi-frequency PPP. First, the PPP results and phase residuals were given. Then, the statistical PPP results of all GNSS stations were statistically analyzed in positioning accuracy and initial convergence time.

The observation data for 100 selected stations in a month are used to evaluate GPS, BDS-3, and Galileo MFPPP performance. It is worth noting that the PPP solutions without OSB are the result of only correcting the DCB. Figure 12 depicts the positioning error of the GPS, BDS-3, and Galileo MFPPP. It can be clearly found that GPS, BDS-3, and Galileo MFPPP can be improved in positioning accuracy and convergence time by correcting the satellite pseudorange and phase time-variant mixed OSB. Furthermore, the OSB estimated by Scheme 2 is superior to that of Scheme 1 in stability, accuracy, and noise level, and its PPP performance is relatively more stable and better than that of Scheme 1 in positioning accuracy and initial convergence.

Figure 13 illustrates the GPS, BDS-3, and Galileo MFPPP phase residuals, where different colors identify different satellites. The carrier-phase residuals are analyzed for the 24 h result. Significant systematic errors exist in the GPS, BDS-3, and Galileo MFPPP results without OSB correction. The RMS values of the carrier-phase residuals are obviously reduced with OSB corrections.

Figure 14 depicts the boxplot of the convergence time for GPS L1 + L2 + L5, BDS-3 B1I + B3I + B2a, BDS-3 B1I + B3I + B1C + B2a, and Galileo E1 + E5a + E5b + E5 + E6 multi-frequency PPP. Moreover, Figure 15 indicates the corresponding boxplot of the positioning accuracy in the north, east, and up components. The corresponding median and mean values are also depicted in the figures. The convergence epoch is defined as the 3D positioning errors kept within 10 cm from the current epoch to the next 20 epochs. The multi-frequency PPP positioning performance in the convergence time and positioning accuracy is significantly improved with the mixed OSB correction. Specifically, the PPP performance of Scheme 2 is superior to Scheme 1 due to the low noise and high accuracy of OSB obtained from Scheme 2. For instance, the GPS L1 + L2 + L5 PPP mean convergence time is reduced from 29.40 min to 24.20 min and then to 20.30 min, the BDS-3 B1I + B3I + B2a PPP mean convergence time is reduced from 34.7 min to 30.50 min and then to 27.30 min, the BDS-3 B1I + B3I + B1C + B2a PPP mean convergence time is reduced from 29.5 min to 25.50 min and then to 21.40 min, and the Galileo E1 + E5a + E5b + E5 + E6 PPP mean convergence time is reduced from 20.90 min to 17.50 min and then to 12.30 min, respectively. The GPS L1 + L2 + L5 PPP positioning accuracy is improved from (0.78, 0.95, 1.53) cm to (0.61, 0.75, 1.30) cm and then to (0.54, 0.60, 1.11) cm in the north, east, and up components. The BDS-3 B1I + B3I + B2a triple-frequency PPP positioning accuracy is improved from (0.83, 1.03, 1.93) cm to (0.75, 0.88, 1.55) cm and then to (0.61, 0.70, 1.31) cm in the north, east, and up components. The BDS-3 B1I + B3I + B1C + B2a PPP positioning accuracy is improved from (0.68 0.94, 1.33) cm to (0.51, 0.70, 1.20) cm and then to (0.44, 0.60, 1.10) cm in the north, east, and up components. The Galileo E1 + E5a + E5b + E5 + E6 PPP positioning accuracy is improved from (0.53 0.93, 1.21) cm to (0.43, 0.68, 1.11) cm and then to (0.38, 0.60, 1.01) cm in the north, east, and up components.

In summary, the new undifferenced mixed OSB model coupled with satellite clock offsets avoids the negative impact of the receiver time-variant biases on the satellite OSB estimation. At the same time, the noise of the mixed OSB obtained from the new method is lower, and the new method combines the pseudorange and phase time-varying OSB, which greatly simplifies the client bias correction and can be better applied to multi-frequency PPP applications.

5. Conclusions

The inconsistency of the precise satellite clock estimated by different frequencies is defined as the IFCB. Unfortunately, the traditional phase time-variant OSB model has high noise and low accuracy and is affected by the receiver time-variant biases. The multi-frequency PPP needs to correct both pseudorange OSB and phase time-variant mixed OSB, which makes the client algorithms more complex. For traditional methods, two bias corrections are required to achieve multi-frequency PPP, which makes the client algorithms more complex. To overcome those flaws and simplify the bias correction process of multi-frequency PPP, the definition of pseudorange and carrier-phase time-variant mixed OSB and a new undifferenced model coupled with satellite clock offsets are presented.

The new undifferenced mixed OSB model coupled with satellite clock offsets avoids the negative impact of the receiver time-variant biases on the satellite mixed OSB estimation. For instance, the JAVAD TRE-3 receiver OSB for BDS-3 B2a and Galileo E6 exhibits significant periodic variations, which lead to periodic bias in satellite mixed OSB for the ED GFIF model. The STD improvement of the undifferenced OSB is 7.5–60.9%, and the RMS can be improved by 9.4–66.1%. Similarly, the STD improvement of the ED OSB can be improved by 50.0–87.5%, and the RMS can be improved by 60.0–88.9%. Hence, the new undifferenced mixed OSB with characteristics of low noise level and high accuracy is more suitable for MFPPP.

With the mixed OSB correction, the GPS and BDS-3 MFPPP performance is obviously improved. With the undifferenced mixed OSB correction, the mean convergence time of GPS, BDS-3, and Galileo PPP can be reduced by several minutes, and the positioning accuracy can be improved by 11.5–20.0% for GPS L1 + L2 + L5 from (0.61, 0.75, 1.30) cm to (0.54, 0.60, 1.11) cm, 15.5–20.5% for BDS-3 B1I + B3I + B2a from (0.75, 0.88, 1.55) cm to (0.61, 0.70, 1.31) cm, 8.3–14.2% for BDS-3 B1I + B3I + B1C + B2a from (0.51, 0.70, 1.20) cm to (0.44, 0.60, 1.10) cm, and 9.0–11.7% for Galileo E1 + E5a + E5b + E5 + E6 from (0.43, 0.68, 1.11) cm to (0.38, 0.60, 1.01) cm compared with the traditional IFCB, respectively.

The GNSS pseudorange and phase time-variant mixed OSB concept and the new undifferenced model coupled with satellite clock offsets demonstrate that it is reasonable and beneficial for the GNSS field.