Assessment of Different Stochastic Models for Inter-System Bias between GPS and BDS

Abstract

Inter-system bias (ISB) will affect accuracy and processing time in integrated precise point positioning (PPP), and ISB stochastic models will largely determine the quality of ISB estimation. Thus, the impacts of four different stochastic models of ISB processing will be assessed and studied in detail to further reveal the influence of ISB in positioning. They are ISB-PW considering ISB as a piece-wise constant, ISB-RW considering ISB as random walk, ISB-AD considering ISB as an arc-dependent constant, and ISB-WN considering ISB as white noise. Together with the model without introducing ISB called ISB-OFF, i.e., five different schemes, ISB-OFF, ISB-PW, ISB-RW, ISB-AD, and ISB-WN, will be designed and tested in this study. From the results of pseudorange residuals, it can be noticed that when considering ISB, the Root-Mean-Square (RMS) of ionosphere-free combined pseudorange residuals are much smaller than without ISB (ISB-OFF). The results of convergence time and positioning accuracy analysis show that PPP performance with ISB-AD is even worse than ISB-OFF, when using the precise products from the German Research Centre for Geosciences (GFZ) named as GBM products here; while the strategies of ISB-RW, and ISB-WN achieve the best results. For the products from Wuhan University called WUM products, a completely different result is achieved. PPP with the stochastic models of ISB-PW and ISB-AD perform best. The most likely reason is the ISB stochastic models applied by the analysis centers are consistent with those used in the PPP on the user side. So, ISB-RW, or ISB-WN is recommended when GBM products are used, and for the WUM products, ISB-PW, or ISB-AD is chosen. From the statistics of PPP precision during the convergence period, it can be concluded that considering ISB also has a great improvement on combined PPP accuracy during the initialization phase.

1. Introduction

Along with the modernization of the Global Positioning System (GPS) [1], the revitalization of the Russian Global Navigation Satellite System (GLONASS) [2], the launch of more European Global Navigation Satellite System (Galileo) satellites [3], and the fast development of Chinese BeiDou Navigation Satellite System (BDS) [4,5], more and more satellites can be observed at one place. Together with the initiation of the Multi-GNSS Experiment (MGEX) in 2012 [6], precise point positioning (PPP) has been stepped up into the stage of multi-system combined mode. However, because of different time systems, reference frames, and hardware delays within receivers used in different systems [7,8], biases exist between different systems when we are making multi-GNSS combined PPP. If we blindly apply observations from all systems without any calibration of these biases, a poor accuracy result will be achieved, and sometimes even worse than that of the single-system PPP. Hence, biases between different systems need to be carefully considered and calibrated, which is the basis of multi-system integrated positioning. Often, an extra parameter named inter-system bias (ISB) is introduced to cover the biases between different satellite systems [9,10,11]. Currently, the processing strategy for ISB is usually to estimate it together with other unknown items, such as coordination, receiver clock, and tropospheric delay, etc. As for the estimation stochastic models of ISB, normally a strategy of considering ISB as an arc-dependent constant is applied in Paziewski and Wielgosz [12], Li, et al. [13], Guo, et al. [14], and Lou, et al. [15]. Besides applying the code ISB estimation, the stochastic model of considering ISB as an arc-dependent constant is also used in the phase ISB estimation for PPP ambiguity resolution (PPP-AR). Because day-to-day variation of phase ISB is less than 0.1 cycles, the phase ISB is usually estimated as an arc-dependent constant every 24 h [16,17]. Jiang, et al. [9] regard ISB as a piece-wise parameter, which can increase the sample size and is good for ISB modeling and prediction. Besides these two strategies, another two stochastic models can be used, which are considering ISB as white noise or as random walk. However, up to now, no specific analysis has been made to assess the difference between these stochastic models and their impact on ISB estimation, therefore in this study, we will evaluate and compare the influence of different stochastic models on ISB, positioning accuracy, and processing time. Based on the results, the recommendations on the selections for the ISB stochastic models will be given.

The paper is organized as follows: Section 2 gives the introductions on the PPP model for ISB estimation and processing strategies. In addition, the detailed descriptions of four different ISB stochastic models and five designed schemes are made. Then, the experimental data used in this study are declared in Section 3. The results and analysis are shown in Section 4, inside pseudorange observation residuals. The estimated ISBs, convergence time, positioning accuracy with GBM and WUM precise products, and accuracy improvement during the convergence time are discussed and compared with five different strategies. Finally, conclusions of this study are made in Section 5.

2. Model for Inter-System Bias Estimation and Descriptions of the Stochastic Models

The main biases between two satellite systems are related to the satellite constellations, the signals broadcast, and time and coordinate frame conventions [18]. Thus, in our multi-GNSS combined processing model, we introduce an extra ISB parameter to cover the effect of the biases between the two systems, and the ISB is estimated with other unknowns using the ionospheric-free (IF) PPP processing.

2.1. Model for Inter-System Bias Estimation and Processing Strategies

The ionosphere-free (IF) pseudorange and phase observation functions between the receiver and a satellite can be written as:

$$\begin{array}{l}P=\rho +c·(d{t}_r-d{t}^s)+m·ZTD+(b_r-b^s)+\varsigma \\ L=\rho +c·(d{t}_r-d{t}^s)+m·ZTD+N+(B_r-B^s)+\epsilon \end{array}$$

(1)

where $P$, $L$ are pseudorange and carrier-phase IF combined observations, respectively; $\rho$ is the range between the receiver and the satellite; $c$ is the light speed; $d{t}_r$ is the receiver clock error and $d{t}^s$ is the satellite clock error; $m$ and $ZTD$ are mapping function and zenith tropospheric delay; the first-order error of ionospheric delay is eliminated by ionosphere-free combination; $b_r$, $b^s$ and $B_r$, $B^s$ are the pseudorange and carrier-phase hardware delay biases of the ionospheric-free function for receivers and satellites; $N$ is ambiguity in distance; $\varsigma$ and $\epsilon$ denote unmodeled parameters (e.g., the measurement noise and multipath errors for pseudorange and carrier phase).

The pseudorange hardware delay biases $b_r$, $b^s$ in the GNSS observations can be assimilated into the clock errors $c·(d{t}_r-d{t}^s)$. The carrier-phase hardware delay biases $B_r$, $B^s$ are not considered in most GNSS data processing, since the carrier-phase hardware delay biases is satellite dependent and stable over time; thus, it can be absorbed by the ambiguity [19,20]. Using the precise satellite orbit and clock products, Equation 1 can be rewritten as:

$$\begin{array}{l}P=\rho +c·d{\tilde{t}}_r+m·ZTD+\varsigma \\ L=\rho +c·d{\tilde{t}}_r+m·ZTD+\tilde{N}+\epsilon \end{array}$$

(2)

where $d{\tilde{t}}_r$ and $\tilde{N}$ are redefined receiver clock error and ambiguity:

$$\begin{array}{l}c·d{\tilde{t}}_r=c·d{t}_r+b_r\\ \tilde{N}=N+B_r-b_r\end{array}$$

(3)

Because the satellite hardware delays are contained in the precise satellite clock and can be eliminated on the user side when applying the precise products in Equation (2) [21]. In Equation (3), the ambiguity term is not an integer property anymore since it contains the bias term, which is the uncalibrated hardware delay biases, showing a fractional cycle property. Thus, the term $B_r-b_r$ is called the uncalibrated phase delay [22].

The tight combined PPP model requires the estimation of an extra ISB parameter. Hence, the combined observation model can be expressed as:

$$\begin{array}{l}{P}^G={\rho }^G+c·d{\tilde{t}}_r+m^G·ZTD+{\varsigma }^G\\ L^G={\rho }^G+c·d{\tilde{t}}_r+m^G·ZTD+{\tilde{N}}^G+{\epsilon }^G\\ P^O={\rho }^O+c·(d{\tilde{t}}_r+ISB)+m^O·ZTD+{\varsigma }^O\\ L^O={\rho }^O+c·(d{\tilde{t}}_r+ISB)+m^O·ZTD+{\tilde{N}}^O+{\epsilon }^O\end{array}$$

(4)

where in the notation of the receiver clock error, mapping functions, the ambiguity, and some other parameters, we use the superscript G as GPS to distinguish them from the parameters of other systems, which are denoted with the superscript O. $d{\tilde{t}}_r$ is a common item between different systems, being the receiver clock term relating to the GPS. The $ISB$ is the inter-system bias parameter between GPS and other systems in the unit of time. It will be modeled as an unknown parameter and estimated together with the other parameters such as coordinates and the receiver clock error. So, the corresponding estimated parameters vector $X$ is,

$$X={[\Delta r,d{\tilde{t}}_r,ISB,ZTD,{\tilde{N}}^G,{\tilde{N}}^O]}^T$$

(5)

where $\Delta r$ denotes the vector of the receiver position increments, ${\tilde{N}}^G$ and ${\tilde{N}}^O$ are the modified ambiguity parameters for different systems. Here, we take GPS and BDS systems as an example for data processing, and a Kalman filter is employed to estimate the unknown parameters in the processing [23,24,25].

Suppose the initial state value is $X_0={\widehat{X}}_0$, the variance is ${\sum }_{X_0}$, so the state estimation ${\widehat{X}}_i$ at epoch $i$ can be obtained recursively by Kalman filter:

(a)

One-step state prediction:

$$X_{i,i-1}={\mathsf{\Phi }}_{i,i-1}{\widehat{X}}_{i-1}$$

(6)

(b)

One-step prediction of covariance:

$${\sum }_{X_{i,i-1}}={\mathsf{\Phi }}_{i,i-1}{\sum }_{X_{i-1}}{\mathsf{\Phi }}_{i,i-1}^T+Q_i$$

(7)

(c)

Filtering gain:

$$K_i={\sum }_{X_{i,i-1}}{A}_i^T{(A_i{\sum }_{X_{i,i-1}}{A}_i^T+R_i)}^{-1}$$

(8)

(d)

Status update:

$${\widehat{X}}_i=X_{i,i-1}+K_i(L_i-A_i{X}_{i,i-1})$$

(9)

(e)

State covariance update:

$${\sum }_{X_i}=(I-K_i{A}_i){\sum }_{X_{i,i-1}}$$

(10)

where ${\mathsf{\Phi }}_{i,i-1}$ is transition matrix; $A_i$ is the coefficient matrix of observation equation; $Q_i$ and $R_i$ is the variance matrix of system noise sequence and variance matrix of observation noise, respectively; $X_{i,i-1}$ and ${\sum }_{X_{i,i-1}}$ is one-step prediction value and its variance-covariance matrix, respectively; $K_i$ is gain matrix; ${\widehat{X}}_i$ and ${\sum }_{X_i}$ is filtered estimation and its variance-covariance matrix, respectively; $I$ denotes the unit matrix.

In this experiment, during the preprocessing we detect and repair clock jumps to avoid misidentifying an observation jump caused by the receiver clock jumps as the cycle slip. Afterwards, the Geometry-Free (GF) and Hatch-Melbourne-Wubbena (HMW) [26,27,28] combinations are used to detect cycle slips. The previously described Kalman filter is employed to estimate the unknown parameters with a sliding window of 3 h station by station, which means the 24-h observation data for each station will be divided into 8 processing sessions. GPS and BDS combined PPP in static and simulated kinematic modes are applied in our study. In static PPP mode, the solution strategy for position coordinates is to consider as contents, while for the kinematic PPP mode, they are estimated as white noise at each epoch. To analyze the impacts of different precise products from different analysis centers (AC) on the ISB estimation, the MGEX precise products from two ACs, GFZ (GBM) and Wuhan University (WUM) are used to fix the satellite orbit and clock in the combined PPP model.

The dry component of the tropospheric delay is corrected with the Saastamoinen model, while the residual zenith wet delay is estimated as a random-walk process. We apply the Global Mapping Function (GMF) [29] to convert zenith delay to slant directions. The cutoff elevation angle is set as 7° and an elevation-dependent weighting method is used. Moreover, the phase-wind up effects [30], the solid earth tide, the ocean loading tide [31] and relativistic effects are corrected by the model. Up to now, BDS-specific receiver antenna Phase Center Offsets (PCOs) and Phase Center Variations (PCVs) are not provided in the igs14.atx. In this case, GPS-specific receiver PCOs and PCVs are applied instead. The initial standard deviation for the carrier-phase observations in GPS and BDS systems are set to 0.003 m, while that for the pseudorange observations is set to 0.3 m. The details of models and strategies related to the data processing are listed in

Table 1. Models and strategies for data processing.

Item	Models/Strategies
Data	GPS + BDS, 8 stations
Processing-arc-window	3 h, 8 stations in one-week period have totally 448 tests
Signal selection	GPS:L1 and L2; BDS:B1 and B2
Estimator	Kalman filter
Elevation cut off	7°
Interval rate	30 s
Satellite orbit and clock	Fixed to MGEX (GBM or WUM) products
Tropospheric delay	Saastamoinen model corrected for the dry component; estimate the residual wet component with random-walk processing
Mapping function	Global Mapping Function (GMF)
Ionospheric delay	First-order effect eliminated by ionospheric-free linear combination
Receiver phase center	PCO and PCV for GPS from igs14.atx are used; Corrections of BDS is applied the same with GPS
Satellite phase center	PCO and PCV are used with igs14.atx
ISB	Schemes ISB-OFF, ISB-PW, ISB-RW, ISB-AD, ISB-WN
Phase ambiguities	Estimated as a constant for each arc

.

The estimation of the ISB parameter can be performed with different stochastic models, such as white noise, random walk, piece-wise constant, or processing-arc-dependent constant, which will be discussed and studied in detail in the following.

2.2. The Different Stochastic Models

To not only verify the importance of ISB, and its influence on the accuracy of tight combined PPP model, but also compare the impacts of different ISB estimation stochastic models on the positioning, five different schemes are proposed in this study. The differences between different schemes mainly focus on the configurations on the ISB. Here, we will take one 3-h processing-arc-window containing 360 observation epochs as an example.

• Neglecting inter-system bias between GPS and BDS is named ISB-OFF, which means the positioning model of BDS and GPS is separated, and the two systems are processed without any relation. The BDS observations are processed the same as GPS ones. However, the ISB between GPS and BDS will remain in the BDS processing model and will be reflected in the pseudorange residuals. The configuration equation can be expressed as:

  $$IS{B}_{OFF}={\left[\begin{array}{ccccc}NUL{L}_1& NUL{L}_2& \dots & NUL{L}_{359}& NUL{L}_{360}\end{array}\right]}^T$$

  (11)

  where $IS{B}_{OFF}$ is the ISB in the scheme ISB-OFF. $NUL{L}_1$ to $NUL{L}_{360}$ is an empty matrix which means without considering ISB.
• Estimating the ISB as a piece-wise constant every 30 min (ISB-PW). ISB parameter is initialized at the first epoch, then re-initialized, and updated every 30 min [9]. The configuration equation can be written as:

  $$\begin{array}{l}IS{B}_{PW}={\left[\begin{array}{cccccc}ISB{P}_1& ISB{P}_2& ISB{P}_3& ISB{P}_4& ISB{P}_5& ISB{P}_6\end{array}\right]}^T\\ ISB{P}_i=\left[\stackrel{60}{\overbrace{\begin{array}{ccccc}IS{B}_i& IS{B}_i& \dots & IS{B}_i& IS{B}_i\end{array}}}\right],i=1,2,3,4,5,6\end{array}$$

  (12)

  where $IS{B}_{PW}$ means the ISB in the scheme ISB-PW. The ISB in one 3-h processing-arc-window can be divided into six parts according to size of piece-wise (every 30 min), which are denoted as $ISB{P}_1$ to $ISB{P}_6$. $ISB{P}_i$ is one piece-wise ISB, and $i$ is from 1 to 6.
• Estimating the ISB as a random-walk processing (ISB-RW). With this scheme, ISB is initialized at the first epoch, and then with a time-related spectral density of 0.001 $m/\sqrt{s}$. The configuration equation can be formed as:

  $$\begin{array}{l}IS{B}_{RW}={\left[\begin{array}{ccccc}ISB{R}_1& ISB{R}_2& \dots & ISB{R}_{359}& ISB{R}_{360}\end{array}\right]}^T\\ ISB{R}_i=ISB{R}_1+sd\times 30\times (i-1),i=1,2,\dots ,359,360\end{array}$$

  (13)

  where $IS{B}_{RW}$ is the ISB in the scheme ISB-RW. $ISB{R}_i$ is the ISB for one epoch, with $i$ from 1 to 360. $sd$ means the time-related spectral density.
• Estimating the ISB as a processing-arc-dependent constant (ISB-AD). During the whole processing-arc period, ISB is only initialized at the first epoch. If the processing-arc window is one day, this means the ISB is estimated as a constant each day. In this paper, we make the period of our processing tests as 3 h, so the ISB is considered to be a 3-h arc-dependent constant. The configuration equation can be indicated as:

  $$IS{B}_{AD}={\left[\stackrel{360}{\overbrace{\begin{array}{ccccc}ISB& ISB& \dots & ISB& ISB\end{array}}}\right]}^T$$

  (14)

  where $IS{B}_{AD}$ means the ISB in the scheme ISB-AD, which is estimated as a constant during the 3-h processing-arc window.
• Estimating the ISB as white noise (ISB-WN). In this way, the ISB is initialized for each epoch, which is the same as the normal processing mode as the receiver clock error. The configuration equation can be shown as:

  $$IS{B}_{WN}={\left[\begin{array}{ccccc}IS{B}_1& IS{B}_2& \dots & IS{B}_{359}& IS{B}_{360}\end{array}\right]}^T$$

  (15)

  where $IS{B}_{WN}$ means the ISB in the scheme ISB-WN. The ISB in each epoch will be initialized and estimated.

For convenience, the five schemes mentioned above marked as ISB-OFF, ISB-PW, ISB-RW, ISB-AD, and ISB-WN, respectively, are summarized in

Table 2. The description of GPS/BDS ISB estimating schemes.

Schemes	Descriptions
ISB-OFF	Neglecting inter-system bias
ISB-PW	Estimating the ISB as piece-wise constant
ISB-RW	Estimating the ISB as a random-walk processing
ISB-AD	Estimating the ISB as a processing-arc-dependent constant
ISB-WN	Estimating the ISB as white noise

.

3. Experimental Data

To validate the different proposed ISB estimation stochastic models and the influence of the extra parameter ISB on PPP performance, observations from 8 stations in the MGEX network, named as CUT0, KARR, MRO1, PERT, MOBS, NNOR, STR1, and YAR2, were chosen, covering a one-week period which is GPS week 1965 with the DOY (day of year) from 246-252, 2017. These 8 stations are all located in Australia, where we can usually observe more than 7 BDS satellites [5]. The distribution and the receiver types equipped in the stations are shown in Figure 1 and

Table 3. Details of 8 selected stations location and receiver types.

Station ID	Location			Receiver Type
	Latitude (°)	Longitude (°)	Height (m)
CUT0	−32.0039	115.8948	24.000	TRIMBLE NETR9
KARR	−20.9814	117.0972	109.247	TRIMBLE NETR9
MRO1	−26.6966	116.6375	354.069	TRIMBLE NETR9
PERT	−31.8019	115.8852	12.920	TRIMBLE NETR9
MOBS	−37.8294	144.9753	40.578	SEPT POLARX4TR
NNOR	−31.0487	116.1927	234.984	SEPT POLARX4
STR1	−35.3155	149.0109	800.032	SEPT POLARX5
YAR2	−29.0466	115.3470	241.291	SEPT POLARX4TR

.

The data of these 8 stations will be processed with the tight combined PPP model, and with the different stochastic models for ISB estimation mentioned above. The results in static and kinematic mode using GBM and WUM precise products will be analyzed and explained in the following sections.

4. Results and Discussion

This section presents the results and analysis performed to assess the difference between diverse ISB stochastic models and their impacts on positioning accuracy and convergence time. As mentioned before, processing is re-initialized eight times per day to evaluate the influence on convergence time. One initialization is considered to be a cold start. A time slide window of 3 h is chosen to ensure enough time for converging into a certain accuracy in nearly all the cases and offer enough test samples. Considering the entire experiment over all stations, it gives 448 cold starts (8 stations × 7 days × 8 initializations).

Here, we use the station coordinates from the IGS weekly solutions as the references, and the positioning errors are defined as the increments between PPP results and the references. The statistical analysis of positioning errors are made with the methods of 68% and 95% quantiles. These two statistical parameters are chosen instead of mean and standard deviation due to the possible remaining biases that might cause results do not follow a normal distribution [32]. The convergence time is defined as when the positioning error is lower than 0.2 m (95%) and 0.1 m (68%) in the North, East, and UP components [15]. This statistical approach is also used in Bree and Tiberius [33], Jr, et al. [32], and Zhou, et al. [34].

Since different MGEX analysis center (GBM, or WUM) has its own processing strategy to obtain precise multi-GNSS orbit and clock products, during the product generation, different ISB stochastic models will be employed with different ACs. In this study, we will apply the precise products from GBM and WUM to validate the difference between these two kinds of products and analyze the potential influence of the different ISB estimation strategies used in different ACs on the positioning and convergence time.

4.1. Pseudorange and Carrier-Phase Observation Residuals

As is already known, the remaining measurement noises and other unmodeled errors will be assimilated by the observation residuals, which is a very important indicator to assess the accuracy of the positioning model. Thus, here the pseudorange and carrier-phase observation residuals will be assessed and compared in the five different schemes.

The pseudorange residuals at station STR1 (with the GBM products as an example) of each satellite with different colors in GPS and BDS system, respectively, are shown in Figure 2 and Figure 3. Figure 2 is the result of static mode, while Figure 3 is the kinematic mode. From the figures, we can see that both in the static and kinematic modes, compared to considering ISB, neglecting ISB (ISB-OFF) will lead to a higher pseudorange residual. After introducing the ISB parameter (ISB-PW, ISB-RW, ISB-AD, and ISB-WN), the ionospheric-free pseudorange (PC) residuals in GPS system will on average improve 4.8 cm (i.e., 5.4%) in static mode, and 5.1 cm (i.e., 5.7%) in kinematic mode. The improvement is not obvious, which is in accordance with the ISB processing strategy that the GPS receiver clock error is considered to be the reference clock error, and ISB covers the biases between BDS and GPS. Thus, the ISB has no effect on the pseudorange residuals in GPS system. As for the BDS system, the PC residual of ISB-OFF is largest; because the remaining errors and difference are not absorbed by ISB, they go into the residuals. The statistical results of the other four schemes are similar, and compared to ISB-OFF, the RMS of pseudorange residuals for ISB-PW, ISB-RW, ISB-AD, and ISB-WN can reduce by 21.3 m (95.9%) both in static mode and kinematic mode. Hence, we can conclude that it is necessary to consider the ISB parameter in multi-GNSS combined positioning. In addition, the carrier-phase observation residuals of GPS and BDS satellites in station STR1 with the GBM products are also taken as an example, the results of static and kinematic modes are displayed in Figure 4 and Figure 5, respectively.

It can be noted from Figure 4 and Figure 5 that the RMS of carrier-phase observation residuals is two orders of magnitude lower than that of the pseudorange residuals, which is caused by the different measurement accuracy levels of pseudorange and carrier-phase observations. Compared with neglecting ISB, the RMS of carrier-phase residuals considering ISB has no significant improvement regarding the pseudorange residuals. That is because in the carrier-phase PPP performance except ISB, the ambiguity is also an unknown parameter. If the ISB is neglected, the impact of unmodeled ISB will be partly or even totally assimilated by the ambiguity, which means the effect of considering ISB or not will not obviously reflect in the carrier-phase residuals. Besides, we can see that when using GBM products, the RMS of carrier-phase residuals in ISB-RW and ISB-WN are the best. Compared with the residuals in ISB-OFF and ISB-AD, the RMS can on average improve by 25.5% and 33.8% in GPS and BDS satellites in the static mode, and for the kinematic mode the improvement rate is 20.2% and 36.3% in GPS and BDS satellites. The carrier-phase observation residuals mainly reflected the accuracy of PPP performance, so we can conclude that the ISB stochastic models impact the carrier-phase observation residuals, and will have a further influence on the PPP performance. In the following parts, the impact of ISB stochastic models on the ISB estimation and PPP performance will be studied in detail.

4.2. Comparison of Estimated Inter-System Bias

The estimated ISB is the most direct factor to reflect which the stochastic model is the optimal one for the PPP performance and how to select an ideal stochastic model for different precise products. Thus, the comparisons of the estimated ISBs from the schemes ISB-PW, ISB-RW, ISB-AD, and ISB-WN using GBM and WUM precise products are made. Figure 6 shows the estimated ISBs using the GBM products in static and kinematic modes. The first 3-h ISB results on DOY 251, 2017 in the station STR1 are taken as an example.

From Figure 6, we can see that the four estimated ISBs from PPP performance with different stochastic models using GBM products have completely different results. The ISBs in ISB-AD have the largest values. The ISBs in ISB-PW are coincident with those in ISB-AD in the first 30 min, but they show different values in the following periods because ISB in ISB-PW was initialized every 30 min. The ISBs in ISB-RW and ISB-WN demonstrate the consistent trend but with a systematic difference, which will be absorbed by the ambiguity in the PPP solution. To deeply show the feature ISBs with different precise products, PPP performance using the WUM products is made. The estimated ISBs with the WUM product is shown in Figure 7.

It can be noticed from Figure 7 that entirely different from the results in Figure 6, the ISBs in ISB-PW, ISB-RW, and ISB-AD have the same trend and very close values. A consistent trend is displayed between the ISB in ISB-WN and the other three ISBs, but the ISB in ISB-WN is rougher than others. The impacts of the different features of ISBs on PPP performance will be analyzed and assessed in the next step.

4.3. Convergence Time and Positioning Accuracy

Positioning accuracy is evaluated by comparing the coordinates with external reference values, which are from the final IGS weekly solution. As for the statistics of convergence time, we use the statistical method of 68% and 95% quantiles mentioned above. In the 68% quantile case, we make the convergence time as the first epoch time whose error is less than 0.1 m and the errors in the following epochs are all less than 0.1 m. As for the 95% quantile, the criterion is less than 0.2 m.

4.3.1. The Case with GBM Precise Products

Figure 8 presents the errors in North (N), East (E), and UP (U) components of GPS/BDS combined static PPP in different schemes. Meanwhile,

Table 4. Convergence times of GPS and BDS combined static and kinematic PPP in five different schemes (unit: min).

	Static (95%)			Static (68%)			Kinematic (95%)			Kinematic (68%)
	N	E	U	N	E	U	N	E	U	N	E	U
OFF	8.5	-	20.0	7.5	74.5	16.0	14.5	-	-	11.5	-	25.0
PW	5.0	63.0	31.0	4.5	45.5	31.0	7.5	72.5	34.5	6.0	63.5	25.5
RW	5.0	13.0	13.5	4.0	11.5	12.0	7.0	20.5	18.0	5.5	16.0	14.0
AD	5.0	-	-	4.5	-	-	7.5	-	-	6.0	-	25.5
WN	5.0	13.5	14.0	4.5	11.5	12.0	7.5	19.5	17.5	5.5	17.5	13.5

shows the convergence time derived in the N, E, and U components, based on the statistics over all the tests.

From Figure 8 and Table 4, it can be noticed that PPP performance at some components in ISB-OFF, and ISB-AD schemes cannot converge in one 3-h processing-arc window. It also can be seen that using the GBM products, the convergence time of ISB-OFF is 8.5 min and 20.0 min in N, and U components at 95% quantile. It takes 7.5 min, 74.5 min and 16.0 min to converge into the 68% quantile in three components. As for the N component, the convergence time of four schemes with ISB (ISB-PW, ISB-RW, ISB-AD, and ISB-WN) is on average improved by 41.2% from 8.5 min to 5.0 min, and 41.7% from 7.5 min to 4.4 min at 95% and 68% quantile, respectively. That means after considering the ISB, the convergence rate speeds up, compared to ISB-OFF. The convergence time of ISB-PW at the U component is longer than ISB-OFF, this may be caused by the stochastic model of ISB-PW here being inconsistent with the model of ISB estimation applied by GFZ for precise product generation, which is not mentioned in any reference, but we have confirmed with the colleagues responsible for the matter in GFZ. From the results, this inconsistency will slow down the ISB convergence rate, and PPP performance time will be increased. This finding corresponds with the results in Figure 6. At E component, ISB-PW performs faster than ISB-OFF. It also can be seen that ISB-RW and ISB-WN have similar results and perform best among the proposed schemes. ISB-PW performs better than ISB-AD, and at the end of the 3-h period, its accuracy can reach closer to ISB-RW, and ISB-WN. It is quite incredible that compared with ISB-OFF, ISB-AD has a poorer result. In other words, GPS/BDS combined PPP considering the ISB with GBM products performs worse than neglecting the ISB, which is largely due to the inconsistency of ISB stochastic model applied in PPP processing on the user side and the precise products generation in GFZ. Therefore, we recommend PPP users adopt the same stochastic model of ISB estimation as the analysis center.

Figure 9 indicates the errors in N, E, and U components of kinematic PPP with different schemes, and Table 4 also shows the convergence time of kinematic PPP at the 95% and 68% quantile in three different components. It can be seen that the convergence times of all schemes in kinematic mode are longer than the static mode, because of different coordinate estimation strategies. However, similar to the static PPP, both at 95% and 68% quantile, ISB-AD performs worst in the convergence performance, while ISB-RW and ISB-WN perform best. PPP performance cannot be converged in E and U components for ISB-OFF, and ISB-AD at 95% quantile. For the 68% quantile, there are no convergence time values in E component. Compared with ISB-PW, the convergence time of ISB-RW at 95% quantile can be improved by 71.7% from 72.5 min to 20.5 min in E component, and 47.8% from 34.5 min to 18.0 min in U component, respectively. As for the 68% quantile, the improvement can reach 74.8% from 63.5 min to 16.0 min in E component, and 45.1% from 25.5 min to 14.0 min. ISB-WN has a similar convergence time as ISB-RW. It can also be noted that the convergence time in U component is less than that in E component in ISB-PW, ISB-RW, and ISB-WN. The first reason may be the phase ambiguities are more correlative with the E component than U component. Thus, in some cases, the convergence time of E component will be enlarged because of their correlation with phase ambiguities [35]. The second reason may be caused by the GPS and BDS satellites with north-south ground tracking in the earth-fixed reference frame [27].

GPS/BDS combined positioning accuracy in static and kinematic modes for five different schemes is listed in

Table 5. Positioning accuracy of GPS/BDS combined static and kinematic PPP with five different schemes (unit: cm).

	Static (95%)			Static (68%)			Kinematic (95%)			Kinematic (68%)
	N	E	U	N	E	U	N	E	U	N	E	U
OFF	9.5	31.9	23.6	1.3	9.3	4.3	14.6	34.6	31.8	4.6	12.7	10.7
PW	1.4	3.6	3.9	0.7	1.6	1.8	2.9	6.7	8.6	1.4	3.0	3.8
RW	1.2	1.6	2.7	0.7	0.7	1.3	2.0	3.0	5.8	0.8	1.3	2.6
AD	12.0	36.3	34.3	3.6	18.9	12.8	19.1	50.6	43.8	7.9	24.8	19.8
WN	1.2	1.6	2.8	0.7	0.7	1.3	2.0	3.0	5.8	0.8	1.3	2.6

. In static PPP, we make the statistics using the results in the last epoch. As for the kinematic mode, since PPP performance cannot converge in some schemes, we make statistics with the same processing period, here using two hours later for each test. We can see, from the table, ISB-RW and ISB-WN have a similar processing accuracy, and perform best, while ISB-AD performs worst. Compared with ISB-AD, in static PPP at 95% quantile, the accuracy of ISB-RW and ISB-WN are improved by 90% from 12.0 cm to 1.2 cm, 95.6% from 36.3 cm to 1.6 cm, and 91.8% from 34.3 cm to 2.8 cm in N, E, and U components, respectively. At 68% quantile, the remarkable improvement can also be achieved with 80.6%, 96.3%, and 89.8% in three components. For kinematic mode at 95% quantile, the accuracy improvements are 89.5% from 19.1 cm to 2.0 cm, 94.1% from 50.6 cm to 3.0 cm, and 86.8% from 43.8 cm to 5.8 cm in N, E, and U components, respectively, while at 68% quantile, 89.9%, 94.8%, and 86.9% accuracy increment in three components can be acquired. Compared to ISB-OFF, schemes with the ISB parameter (i.e., ISB-PW, ISB-RW, and ISB-WN) can get an obvious accuracy improvement, which indicates the importance of considering ISB. Therefore, we conclude that from the validation of convergence time and performance accuracy using the GBM precise products, it is necessary to introduce the ISB parameter, and estimate it as a random-walk processing, or white noise, if you want to get a higher processing precision in a short time period (e.g., 3 h in this study). The scheme of ISB estimated as an arc-dependent constant (ISB-AD) is not recommended, because of its longest convergence time and worst positioning accuracy.

4.3.2. The Case with WUM Precise Products

The analysis of the convergence time and positioning accuracy in GPS/BDS combined kinematic PPP using WUM precise products is shown in Figure 10. Meanwhile,

Table 6. Convergence time and positioning accuracy of GPS/BDS combined kinematic PPP with the WUM products in five different schemes.

	Convergence Time (95%)/min			Convergence Time (68%)/min			Accuracy (95%)/cm			Accuracy (68%)/cm
	N	E	U	N	E	U	N	E	U	N	E	U
OFF	7.5	17.5	17.0	6.5	15.5	13.0	4.9	6.1	7.6	2.0	2.7	3.2
PW	7.0	14.0	14.0	5.0	12.5	12.0	4.8	6.6	7.6	2.0	2.9	3.2
RW	7.0	21.0	14.5	5.5	17.5	12.5	5.0	7.0	7.8	2.1	3.3	3.3
AD	7.0	14.0	14.0	5.0	12.5	12.0	4.8	6.2	7.4	2.0	2.7	3.1
WN	7.0	23.5	14.5	5.5	19.0	12.5	5.2	6.8	8.0	2.2	3.3	3.3

presents the statistics of convergence time and positioning accuracy. The ISB stochastic model applied in AC Wuhan University is considered it as an arc-dependent constant [14]. From the figure and the table, we can see that different from the tests with the GBM products, the performance using the WUM products have similar accuracy level whether considering the ISB parameter or not. The possible reason is that ISB estimation strategy used in WUM products is the same as the ambiguity processing strategy, and both are estimated as a constant each arc window. Thus, the ISB is highly correlative with the ambiguity; in this case, if the ISB is not introduced, the impact of its error will be absorbed by ambiguity. This finding corresponds with the statements from Figure 7. For the convergence time, in N component the improvements are small (6.7% from 7.5 min to 7.0 min) at the 95% quantile, and at 68% quantile, it is on average 19.2%, from 6.5 min to 5.25 min. For the E component, the schemes ISB-PW and ISB-AD perform better than ISB-OFF, and ISB-RW is worse than ISB-OFF, while ISB-WN performs worst. Compared with ISB-WN, the convergence time of ISB-PW and ISB-AD can be improved by 40.4% from 23.5 min to 14.0 min, and 34.2% from 19.0 min to 12.5 min at 95% quantile and 68% quantile, respectively. Similar to N component, the improvements of convergence time in U component considering ISB are low, which can on average reach 16.2%, and 5.8% at 95% quantile, and 68% quantile, respectively.

From these two cases with GBM and WUM precise products, we can conclude that GPS/BDS combined PPP processing with orbit and clock products from different ACs can get a completely different result, because during the generations of multi-GNSS orbit and satellite clock products, different ISB resolution strategies (stochastic model) are applied in different ACs, and if the users use the same ISB strategies as the ACs (GBM and WUM) during the PPP performance, a better result can be obtained. So, based on the tests in this study, we recommend users consider the ISB parameter as random-walk processing, or white noise if the GBM products are used; however, when the WUM products are applied, it is better to consider the ISB as a piece-wise constant, or an arc-dependent constant.

4.4. Accuracy Improvement during the Convergence Period

During the convergence period, the accuracy of GPS/GLONASS combined PPP can be improved after considering the inter-frequency biases (IFB) [36]. So, similar analysis of accuracy improvement in GPS/BDS fusion PPP when considering ISB during the initialization phase is made in this study.

Table 7. Positioning accuracy of GPS/BDS combined kinematic PPP during the convergence period (8 min and 16 min) in five different schemes with GBM precise products (unit: cm).

	8 min (95%)			16 min (95%)			8 min (68%)			16 min (68%)
	N	E	U	N	E	U	N	E	U	N	E	U
OFF	71.4	78.4	191.8	76.1	94.7	178.7	49.3	57.9	106.1	35.6	45.5	76.4
PW	58.2	69.3	135.6	43.3	58.8	100.0	24.6	31.4	58.1	18.1	27.1	42.7
RW	58.7	71.2	135.9	42.9	53.1	99.0	24.3	30.6	57.4	17.6	23.7	41.5
AD	58.2	69.3	135.6	43.4	58.8	100.0	24.6	31.4	58.1	18.1	27.1	42.7
WN	60.5	69.9	137.1	43.8	54.3	99.2	24.6	31.1	57.5	17.8	24.3	41.5

and

Table 8. Positioning accuracy of GPS/BDS combined kinematic PPP during the convergence period (8 min and 16 min) in five different schemes with WUM precise products (unit: cm).

	8 min (95%)			16 min (95%)			8 min (68%)			16 min (68%)
	N	E	U	N	E	U	N	E	U	N	E	U
OFF	104.8	122.8	247.2	51.3	59.7	137.7	32.1	38.7	82.0	23.0	29.1	58.6
PW	59.2	68.7	136.9	42.0	52.0	97.6	25.3	30.5	57.0	18.2	22.9	41.0
RW	59.3	68.2	136.8	42.5	55.5	98.0	25.4	31.4	57.0	18.4	24.3	41.1
AD	59.2	68.7	136.9	42.0	52.0	97.6	25.3	30.5	57.0	18.2	22.9	41.0
WN	59.6	73.4	137.0	43.2	57.6	98.8	25.7	32.5	57.6	18.6	25.6	41.5

show the statistics of the positioning accuracy during the convergence period. It can be seen from the tables that no matter which product is applied, if ISB is neglected, the accuracy in N, E, and U components is lower than the schemes with considering the ISB during 8 min and 16 min convergence period both at 95% and 68% quantiles. Therefore, we can conclude that considering ISB also has a great improvement on the PPP positioning accuracy during the convergence period. With GBM precise products, the average improved rate can be 40.1%, 35.4%, and 41.0% in N, E, and U components, respectively. In addition, for the WUM products, this rate can be 25.4%, 22.3%, and 33.4% in three components, respectively.

5. Conclusions

To study the effect of stochastic models for the ISB estimation and the importance of considering ISB, five different ISB processing strategies are proposed. From the comparison of pseudorange observation residuals between five strategies, it can be found that the ISB-OFF strategy has the largest residuals, and the other four schemes have much lower results. Compared to ISB-OFF, the RMS of pseudorange residuals in ISB-PW, ISB-RW, ISB-AD, and ISB-WN can reduce by 21.3 m (improvement rate is 95.9%) both in static mode and kinematic mode. That is to say, it is quite necessary to consider the effect of ISB parameter in multi-GNSS combined positioning.

As for the analyses on the convergence time and positioning accuracy, the cases with GBM and WUM precise products are carried out to reveal the potential impact on ISB caused by the different solution strategies in the analysis center. The result of the GBM case shows that ISB-RW and ISB-WN strategies have a similar processing accuracy, and perform best, while ISB-AD performs worst, even worse than the ISB-OFF, which is against the conclusion of pseudorange observation residuals results, the most probable reason being that the ISB stochastic model during GBM precise product generation is not set as an arc-dependent constant. The result of the WUM case is opposite; the schemes ISB-PW and ISB-AD perform best, and ISB-RW is worse than ISB-OFF, while ISB-WN performs worst. Based on the above findings, it can be concluded that the ISB processing strategies of GBM and WUM products are quite different. Thus, the recommendation of the ISB stochastic model selection can summarized as follows: it is better to consider the ISB parameter as random walk, or white noise if the GBM products are used; however, when the WUM products are applied, it is better to consider the ISB as a piece-wise constant, or an arc-dependent constant.

The test of accuracy improvement during the convergence period is performed as well. From the results, the conclusion can be drawn that just as after a converged period, considering ISB or not also has a strong impact on PPP positioning accuracy during the convergence period.