Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.020106

ARTICLE

Modeling and Prediction of Inter-System Bias for GPS/BDS-2/BDS-3 Combined Precision Point Positioning

Zejie Wang1, Qianxin Wang1,* and Sanxi Li2

1School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou, 221116, China
2Beijing Railway Electrification School, Beijing, 102202, China
*Corresponding Author: Qianxin Wang. Email: wqx@cumt.edu.cn
Received: 04 November 2021; Accepted: 11 January 2022

1  Introduction

The Precise Point Positioning (PPP) technology was first proposed and realized by Zumberge et al. in 1997 to achieve efficient and robust analysis of GPS data from large networks [1]. In the past two decades, because of its high precision, efficiency, and flexibility, the PPP technology has been widely applied to precise orbit determination of low earth orbiters, earthquake monitoring, tsunami warning, atmospheric retrieval, and other fields [2–7]. However, the traditional PPP technique has limitations of long initialization time and limited positioning precision [8]. The advent of multi-GNSS provides a great opportunity to reduce the convergence time and improve the positioning precision because multi-GNSS has advantages in the “richness” of constellations, signals, and frequencies [9–11]. The PPP with multi-GNSS is called MGPPP. However, the Inter-System Bias (ISB) measurement of Multi-GNSS, including the time system offset, the coordinate system difference, and the inter-system hardware delay bias, must be considered for Multi-GNSS data fusion processing. Generally, ISB is regarded as a new variable in GNSS data processing, and it is an important issue in the fusion of multi-GNSS data [12–16].

To control the effect of ISB, three types of methods are widely used in MGPPP: (1) The Parameter Removing Method (PRM) eliminates the inter-system hardware delay differences by Between-Satellite Single-Difference (BSSD). In addition, the method eliminates the time and coordinate system differences by adopting unified satellite and clock products or performing time and coordinate system transformations. Currently, the PRM method is widely used in the loosely combined PPP model [17]. (2) The Parameter Estimation Method (PEM) takes the ISB as an additional parameter and makes estimations or calibrations through observation equations. The PEM method is widely used in the un-differenced PPP model [18,19] and the tightly combined PPP model (a single reference satellite for all observations from all GNSS systems [20]). (3) The Parameter Simplification Method (PSM) assimilates the ISB into the receiver clock, ambiguity, and pseudo-range residuals parameters to simplify MGPPP data processing and preserve the same parameter solution model as that of the traditional single GNSS PPP [21,22].

These methods of ISB processing are proposed and used for different application scenarios. Therefore, they have inherent disadvantages and advantages. In the PRM method, the effect of ISB can be eliminated, but the magnitude of ISB cannot be obtained and there is a requirement on the satellite numbers in each GNSS system (e.g., at least two satellites). In the PEM method, the nominal ISB can be estimated as an additional parameter, but the absolute ISB cannot be obtained because it is strongly related to the clock and ambiguity parameters. Also, additional ISB need to be estimated, which is impractical for the positioning under limited satellite availability. In the PSM method, there is no need to estimate the ISB, and the model of MGPPP is consistent with that of a single GNSS PPP, but the characteristics of some parameters are changed. In this case, the users need to choose a suitable method to process the ISB.

Considering that PEM is the most popular method for processing ISB, this study focuses on the investigation of the PEM-based ISB processing strategy for GPS/BDS-2/BDS-3 combined PPP. It is acknowledged that there are obvious ISB between BDS-2 and GPS observations [23,24], and the “inter-satellite-type biases” (ISTBs) have been found in different types of BDS including Geostationary Earth Orbit (GEO), Inclined GeoSynchronous Orbit (IGSO), and Medium Earth Orbit (MEO) [25]. Besides the newly designed satellite signals, different satellite buses, updated rubidium atomic frequency standards, and new passive hydrogen masers have been used in BDS-3 satellites [26]. Compared with BDS-2 satellites, great improvements have been made for BDS-3 satellites. However, some problems need to be further investigated: (1) Are there obvious ISB between BDS-2 and BDS-3 observations (ISBBDS-2, BDS-3)? (2) What are the characteristics of the ISB time series if they exist? (3) How to model and predict these ISB parameters accurately? (4) How to effectively utilize the predicted ISB to improve the positioning precision and convergence time of MGPPP?

To investigate the above problems, the following work has been done in this study. Firstly, the estimation method of ISBBDS-2, BDS-3 is discussed based on GPS/BDS-2/BDS-3 combined PPP. Then, the ISB are estimated by using one month of the GPS/BDS-2/BDS-3 observations from 11 Multi-GNSS Experiment (MGEX) stations and 10 international GNSS Monitoring and Assessment System (iGMAS) stations. The significances of estimating ISBBDS-2, BDS-3 are tested and analyzed in Section 2. Then, these ISB time series are modeled and predicted by a newly proposed method called QM+ARIMA+SW, and the PISB precision of each station is analyzed in Section 3. Subsequently, a processing strategy of the PISB is proposed for MGPPP, where the PISB is used by adding some virtual observation equations and adaptive factors. Meanwhile, several experimental schemes are designed and implemented to validate the effectiveness of the proposed method. The positioning precision and convergence of these schemes are compared under different satellite visibility conditions in Section 4. Finally, this study is concluded in Section 5.

2  ISB Estimation

2.1 ISB Estimation Method

Generally, the ionosphere-free linear combinations of GPS, BDS-2, and BDS-3 observations can be expressed as the following Eqs. (1)–(6):

PGIF=ρG+c[dtr−dtGs]+brGP+bGsP+TG+εGIFP(1)

PB2IF=ρB2+c[dtr−GB−dtB2s]+brB2P+bB2sP+TB2+εB2IFP(2)

PB3IF=ρB3+c[dtr−GB−dtB3s]+brB3P+bB3sP+TB3+εB3IFP(3)

ΦGIF=ρG+c[dtr−dtGs]+brGΦ+bGsΦ+TG+NGIF+εGIFΦ(4)

ΦB2IF=ρB2+c[dtr−GB−dtB2s]+brB2Φ+bB2sΦ+TB2+NB2IF+εB2IFΦ(5)

ΦB3IF=ρB3+c[dtr−GB−dtB3s]+brB3Φ+bB3sΦ+TB3+NB3IF+εB3IFΦ(6)

In the above equations, the subscripts G, B2, B3 respectively refer to GPS, BDS-2, BDS-3 satellite systems; P and Φ respectively indicate the pseudo-range and carrier phase measurements (unit: meter); IF indicates the ionosphere-free combination; ρ indicates the actual geometric distance from the receiver to the satellite, which includes the receiver coordinates (xr, yr, zr) and the satellite coordinates (xs, ys, zs); c indicates the speed of light in vacuum (unit: meter per second); dtr and dts respectively indicate the clock errors of the receiver and the satellite (unit: second); GB indicates the time offset of GPS to BDS; brPand bsP respectively indicate the frequency-dependent code hardware delays of the receiver and the satellite (unit: second); brΦ and bsΦ respectively indicate the frequency-dependent carrier phase hardware delays of the receiver and the satellite (unit: second); T indicates the tropospheric delay (unit: meter); N indicates the ambiguity parameters of the ionosphere-free linear combinations (unit: meter); εp and εΦ respectively indicate the noise and the un-modeled errors of the ionosphere-free linear combinations of the code and carrier phase measurements (unit: meter).

To obtain precise and consistent products, the satellite orbit and clock of GPS, BDS-2, and BDS-3 are produced in advance, and the precision of these products is validated. Specifically, the GNSS observations for one month (from June 01 to 30, 2018) of 93 IGS-MGEX stations and 22 iGMAS stations are used to produce the satellite orbit and clock of GPS, BDS-2, and BDS-3. This study adopts the data processing strategy proposed in [27]. The IGS-MGEX products are taken as references to validate the precision of GPS and BDS-2 satellite orbit and clock determination results. Meanwhile, the precision of GPS, BDS-2, and BDS-3 satellite orbit and clock determination results is validated by the discrepancy of overlapping arcs. The signal tracking capability of each station is illustrated in Fig. 1. The statistics of precision results are presented in Fig. 2 and Table 1.

Figure 1: The signal tracking capability of each station used to determine the satellite orbit and clock of GPS, BDS-2, and BDS-3 [28]

Figure 2: The precision of satellite orbit (on the left side, the mean RMS of all satellites in one day) and clock (on the right side, the mean standard deviation of all satellites) determination results in one month (from June 01 to 30, 2018)

It can be seen from Table 1 that the precision of GPS/BDS-2/BDS-3 satellite orbit (1D RMS) and clock (STD) determination results is better than 2 cm/10 cm/15 cm and 0.1 ns/0.3 ns/0.5 ns, respectively. Also, the precision of the BDS-2 satellite orbit and clock determination results is slightly better than that of BDS-3. This is mainly because the number of BDS-2 tracking stations (80 stations) is larger than that of BDS-3 tracking stations (19 stations). It is expected that the precision of the BDS-3 satellite orbit and clock determination results can be further improved when more stations are updated to track the BDS-3 satellite signal.

After the precise satellite orbit and clock, productions are substituted into Eqs. (1)–(6), the dts and (xs, ys, zs) will become known parameters. In our experiments, the time and coordinate systems of the GPS, BDS-2, and BDS-3 satellite orbit and clock products are referred to as GPST and WGS-84, respectively. Thus, the time offset GB can be eliminated if the above orbit and satellite clock products are used. Besides, the satellite clock corrections include the effect of the ionosphere-free linear combination of the satellite code hardware delays (bsP). Therefore, Eqs. (1)–(6) can be rewritten as Eqs. (7)–(12):

PGIF=ρG+c[dtr−dtGprecs]+brGP+TG+εGIFP(7)

PB2IF=ρB2+c[dtr−dtB2precs]+brB2P+TB2+εB2IFP(8)

PB3IF=ρB3+c[dtr−dtB3precs]+brB3P+TB3+εB3IFP(9)

ΦGIF=ρG+c[dtr−dtGprecs]−bGsP+brGΦ+bGsΦ+TG+NGIF+εGIFΦ(10)

ΦB2IF=ρB2+c[dtr−dtB2precs]−bB2sP+brB2Φ+bB2sΦ+TB2+NB2IF+εB2IFΦ(11)

ΦB3IF=ρB3+c[dtr−dtB3precs]−bB3sP+brB3Φ+bB3sΦ+TB3+NB3IF+εB3IFΦ(12)

where dtsprec= dts+GB-bsp/c. Generally, the code hardware delays of the receiver (brP) can be assimilated into the receiver clock correction (c×dtr). Meanwhile, the differences between the code and the carrier phase hardware delays of the receiver (brΦ-brp) and satellite (bsΦ-bsp) can be assimilated into the ambiguities (N). If the new symbols dt r’ and N’ are adopted to denote the new receiver clock correction (c×dtr+brGP) and the new ambiguities (N+brΦ-brp+bsΦ-bsp), an additional parameter (i.e., the difference of receiver code hardware delays between different satellite systems, such as brGP-brB2P or brGP-brB3P) must be added, and the added parameter is ISB. Based on this, Eqs. (7)–(12) can be rewritten as Eqs. (13)–(18):

PGIF=ρG+cdtr′−cdtGprecs+TG+εGIFP(13)

PB2IF=ρB2+cdtr′−cdtB2precs+ISBG,B2+TB2+εB2IFP(14)

PB3IF=ρB3+cdtr′−cdtB3precs+ISBG,B3+TB3+εB3IFP(15)

ΦGIF=ρG+cdtr′−cdtGprecs+TG+NGIF′+εGIFΦ(16)

ΦB2IF=ρB2+cdtr′−cdtB2precs+ISBG,B2+TB2+NB2IF′+εB2IFΦ(17)

ΦB3IF=ρB3+cdtr′−cdtB3precs+ISBG,B3+TB3+NB3IF′+εB3IFΦ(18)

According to Eqs. (13)–(18), the ISB parameters (ISBG, B2 and ISBG, B3) and the other parameters (xr, yr, zr, dt’, T, N’) can be estimated by the same estimator, e.g., the least square (LSQ) adjustment method, or other adjustment methods. Then, the value of ISBB2, B3 can be obtained by subtracting ISBG, B2 with ISBG, B3.

2.2 ISB Estimation Results

To estimate ISBG, B2, ISBG, B3, and ISBB2, B3, the observations are collected from 21 stations (include 11 MGEX stations marked in black and 10 iGMAS stations marked in red) for one month (from June 1 to 30, 2018). All the 21 stations can receive the satellites signals of GPS (L1, L2), BDS-2 (B1I, B3I), and BDS-3 (B1I, B3I) simultaneously; the 10 iGMAS stations can receive the satellites signals of new BDS-3 (B1C, B2a). Because the ISB is highly dependent on the receiver type [29], the receiver types of the above 21 stations are distinguished with different colors, as shown in Fig. 3. Table 2 lists the ISB estimation strategy used in our experiments. Fig. 4 illustrates the MGPPP precision of each station in the north (N), east (E), and up (U) components when using the GPS (L1, L2)/BDS-2 (B1I, B3I)/BDS-3 (B1I, B3I) observation data. Fig. 5 shows the MGPPP error distribution of the 10 iGMAS stations when using the GPS (L1, L2)/BDS-2 (B1I, B3I)/BDS-3 (B1I, B3I) and GPS (L1, L2)/BDS-2 (B1I, B3I)/BDS-3 (B1C, B2a) observation data. The positioning precision is evaluated by comparing the coordinates with external values obtained from more than 200 stations at CUMAC (GNSS Analysis Center at China University of Mining and Technology) based on the double-different (DD) solution. For the same station, the coordinate difference between the external values of the DD solution and the IGS final products [30] is just about 3 mm in the horizontal direction and 6 mm in the vertical direction. This indicates the high precision of the DD solution. Thus, the external values obtained by the DD solution can be regarded as the actual values for the MGPPP of the 11 MGEX stations and 10 iGMAS stations. Fig. 6 shows the ISB estimation results of the 21 stations on 01 June 2018. Table 3 lists the statistical estimation results of ISB BDS-2 (B1I, B3I), BDS-3(B1I, B3I) and ISB BDS-2 (B1I, B3I), BDS-3(B1C, B2a) from all the stations that are classified by the receiver type. Fig. 7 illustrates the ISB estimation results of the CLGY, WUH1, DWIN, and ZHON stations for one month.

Figure 3: The geographical distribution and receiver types of the 11 IGS MGEX stations and the 10 iGMAS stations used for ISB estimation

Figure 4: The PPP precision (the mean RMS of 30 days) of 21 stations (10 iGMAS stations and 11 MGEX stations) in the N, E, and U components when using GPS (L1, L2)/BDS-2 (B1I, B3I)/BDS-3 (B1I, B3I) observation data

Figure 5: The distribution of the MGPPP positioning errors for the 10 iGMAS stations. Top panels: the error distribution when using GPS (L1, L2)/BDS-2(B1I, B3I)/BDS-3(B1I, B3I) observation data from June 1 to 30, 2018. Bottom panels: the error distribution when using 10 iGMAS stations using GPS (L1, L2)/BDS-2(B1I, B3I)/BDS-3(B1C, B2a) observation data from June 01 to 30, 2018

Figure 6: The ISB estimation results of the 21 stations on 01 June 2018 Top panels: ISB GPS (L1, L2), BDS-2 (B1I, B3I); Bottom panels: ISB GPS (L1, L2), BDS-3 (B1C, B2a)

Figure 7: The ISB estimation results of the CLGY, WUH1, DWIN, and ZHON stations for one month. Left panels: the ISB GPS (L1, L2), BDS-2 (B1I, B3I) of CLGY, WUH1, DWIN, and ZHON stations from June 01 to 30, 2018; Right panels: the ISB GPS (L1, L2), BDS-3 (B1C, B2a) of CLGY, WUH1, DWIN, and ZHON stations from June 01 to 30, 2018

It can be seen from Figs. 4 and 5 that the daily solution precision of GPS/BDS-2/BDS-3 combined PPP is better than 2 cm for most of the stations, and the positioning precision of IGS stations is better than that of iGMAS stations. This is because IGS stations have a better receiver quality than iGMAS stations. The following observations can be made from the results in Figs. 6 and 7 and Table 3: (1) The ISB are highly dependent on the receiver type, and the ISB obtained from the IGS stations are better than that obtained from the iGMAS stations; (2) there is no apparent difference in the ISB of BDS-2 and BDS-3 observations if they all use B1I and B3I signals, but there is an obvious difference in the ISB of BDS-2 and BDS-3 observations if they use B1I/B3I and B1C/B2a signals; (3) the ISB have high stability (about 0.1–0.2 ns) for the same receiver in one day, but the ISB vary greatly (about 1–2 ns) day by day, and systematic variations are presented. Therefore, the ISB time series can be modeled and predicted, and the prediction results will be conducive to reducing the convergence time and improving the positioning precision of MGPPP, especially in the condition of limited satellites.

3  ISB Modeling and Prediction

The investigation of ISB modeling and prediction has been performed in [31,32]. However, the precision of predicting ISB next day is about 0.5–1.5 ns, which is not high enough. To improve the precision of short-term ISB prediction, a combination method based on quadratic model (QM), autoregressive integrated moving average model (ARIMA), and sliding window strategy (SW), i.e., QM+ARIMA+SW, is proposed. The principle of the proposed method is described as follows.

First, a quadratic function is employed to fit the original ISBoi series by the least square method:

ISBtf=At2+Bt+C(19)

where ISBft indicates the fitting value of the ISB at the tth epoch; A, B, and C are the coefficients of the quadratic, first-order, and constant terms of the quadratic function, respectively. Fig. 8 illustrates the fitting results of the stations based on the quadratic function. Due to the limited paper space, only the fitting results of WUH1 and ZHON stations are shown here.

Figure 8: The fitting results of ISB GPS (L1, L2), BDS-2 (B1I, B3I) and ISB GPS (L1, L2), BDS-3 (B1C, B2a) of WUH1, and ZHON stations from June 01 to 30, 2018 based on QM

Then, the fitting residuals of Eq. (19) are further modeled by ARIMA [33], which combines the difference method and the ARMA model. The modeling operation can be expressed as Eq. (20):

ISBtr=∑i=1paiISBt−ir+∑i=1qbiνt−i(20)

where ISBrt indicates the difference between ISBot and ISBft; (a, b) and (p, q) are the coefficients and the orders of the ARMA model, respectively; v is the error term, {vt} ∼ WN (0, σ2), and WN indicates white noise. Noted that the time series modeled by the ARMA model must meet the requirement for stability. Thus, the stability of {ISBrt} should be tested first. If {ISBrt} is unstable, the difference operation will be performed to make it stable.

The combination of different methods and the ARMA model can be expressed as {ISBrt}∼ARIMA (p, d, q), where d indicates the times of different operations. Fig. 9 shows the fitting results of ISBrt based on the ARIMA model. Table 4 lists the fitting precision of ISBft+ISBrt based on QM+ARIMA. The Mean Absolute Error (MAE) is adopted as the precision index of ISB modeling and prediction [34], and it is calculated as Eq. (21):

MAEj=1n∑t=1n|(ISBtf+ISBtr)−ISBto|(21)

where ISBot is the original value of the ISB at the tth epoch, and (ISBft+ISBrt) are replaced with the prediction value ISBpt to evaluate the prediction precision of ISB; n is the total number of epochs for ISB modeling or prediction; j indicates the jth station.

Finally, the ISB at the (t+1)th epoch is predicted according to Eq. (22):

ISBt+1p=A(t+1)2+B(t+1)+C+∑i=1paiISB(t+1)−ir+∑i=1qbiv(t+1)−i(22)

Figure 9: The fitting results of ISBrt of WUH1 and ZHON stations from June 01 to 30, 2018 based on ARIMA

In Eq. (22), A, B, C, (a, b), and (p, q) are known parameters, and they are obtained by fitting the models according to Eqs. (19) and (20). To reduce the influence of historical data and strengthen the effect of new data, a slide window strategy is adopted in this study. Specifically, the data at the 1st epoch are deleted and the data at the (t+1)th epoch are added. The window size is always a constant. To validate the effectiveness of the proposed method, the original ISB of each station are taken as basic data series to model and predict the ISB based on the proposed QM+ARIMA+SW method, and the window size is always set to 23. The ISB prediction results and the statistics of the prediction precision are respectively shown in Fig. 10 and Table 5.

Figure 10: The ISB prediction results of WUH1 and ZHON stations from June 24 to 30, 2018 based on QM+ARIMA+SW

It can be seen from Table 5 that based on the proposed QM+ARIMA+SW method, the prediction precision of the ISB in the next day is about 0.1−0.6 ns, which is better than that of the existing research results (0.5−1.5 ns). However, there is a significant difference in the predicted ISB (PISB) precision of each station. For example, the PISB precision of the ZHON station is high (about 0.1 ns), and that of the DWIN station is poor (about 0.6 ns). The results indicate that it is difficult to find an effective method to accurately predict the ISB of all stations. This is because the time-varying characteristic of ISB is closely related to the capability (stability) of the receiver hardware.

4  Selective Utilization of the Predicted ISB

4.1 Selective ISB Prediction Method

The predicted ISB of each station can be obtained by the proposed QM+ARIMA+SW method, but the prediction precision varies. To address this issue, a better approach needs to be developed to handle the PISB of each station separately. However, the PISB precision of each station cannot be known in advance. Thus, the reliability of the PISB needs to be determined before they are used. To effectively utilize PISB, a new PISB processing strategy for GPS/BDS-2/BDS-3 combined PPP is proposed in this study. The fundamentals of the proposed strategy are described as follows.

Here, Eqs. (13)–(18) are still used to estimate the ISB parameters (ISBG, B2 and ISBG, B3) and the other parameters (xr, yr, zr, dt’, T, N’). Meanwhile, two virtual observation equations are added to estimate the ISB parameters, which are shown in Eqs. (23) and (24):

ISBG,B2P=ISBG,B2+εG,B2P(23)

ISBG,B3P=ISBG,B3+εG,B3P(24)

where ISBP and εP denote the predicted PISB and its prediction error, respectively.

Besides, an adaptive factor α is introduced, and the setting of its value refers to [35]:

α={1|ΔX|≤C0C0|ΔX|(C1−|ΔX|C1−C0)2C0<|ΔX|≤C10|ΔX|>C1(25)

where C0 and C1 are constants, C0 = 1.0–1.5, C1 = 3.0–4.5. The calculation of |ΔX| is shown as follows:

|ΔX|=||Xe−Xp||tr(∑Xp)(26)

where Xe is the state vector of estimated ISB, and it is calculated according to Eqs. (13)–(18); and Xp is the state vector of predicted ISB, and it is calculated according to Eqs. (23) and (24); Σ and tr denote the covariance matrix and the trace of a matrix, respectively. Finally, the unknown parameters (xr, yr, zr, dt’, T, N’, ISBG, B2, ISBG, B3) can be solved by combining Eqs. (13)–(18) and (23)–(26). The calculation formulas are shown as follows:

X=(ATPA+αPp)−1(ATPL+αPPXP)(27)

∑X=(ATPA+αPp)−1σ02(28)

where A, P, and L are respectively the design matrix, weight matrix, and observation vector, and they are calculated according to the original Eqs. (13)–(18); Pp and Xp are respectively the weight matrix and the predicted state vector of ISB, and they are calculated according to the added virtual observation Eqs. (23) and (24); X and ΣX are the final estimated state vector and its covariance matrix, respectively; σ20 is a scale factor, and α is an adaptive factor. Since the matrices Pp and Xp are obtained through matrix extension of P and X, the dimensions of Pp and Xp should be the same as those of P and X.

Note that the adaptive factor α ranges between 0 and 1, and it balances the contribution of the OISB and the PISB to the final estimation of ISB. If α = 0, the solutions to Eqs. (27) and (28) are the same as those of Eqs. (13)–(18); in this case, the PISB is not trusted and ignored. If α = 1, the solutions to Eqs. (27) and (28) are the same as those of directly combining Eqs. (13)–(18) and Eqs. (23) and (24); in this case, the PISB is fully trusted and treated the same as the OISB. If 0 < α < 1, the weight of the virtual observations obtained according to Eqs. (23) and (24) are adjusted by the difference between the PISB and the OISB. This indicates that the method proposed in this study has a general form of PISB processing, and the methods of ignoring and directly utilizing PISB are the special cases of the proposed method.

4.2 Experiments and Analysis

To validate the proposed method, four experimental schemes are designed in this study:

(1)   Scheme 1: α is set to 0. That is, the information of the PISB is ignored, and only GPS/BDS-2 observations are used for MGPPP;

(2)   Scheme 2: α is set to 0, but the observations of GPS/BDS-2/BDS-3 are used for MGPPP;

(3)   Scheme 3: α is set to 1. That is, the information of the PISB is treated the same as the OISB, and the observation is the same as that of scheme 2;

(4)   Scheme 4: α is determined by Eq. (25). That is, the information of the PISB is utilized selectively, and the observation is the same as that of scheme 2.

Then, the four experimental schemes are tested under different satellite availability conditions. Based on the test results, the impacts of different PISB processing strategies on MGPPP in different satellite visibility environments are analyzed.

In the open sky environment, the observation data of the CLGY, WUH1, DWIN, and ZHON stations for one week are used. Meanwhile, the satellite cut-off elevation angle is set to 5∘, and the satellite signals of GPS (L1+L2), BDS-2 (B1I+B3I), and BDS-3 (B1C+B2a) are adopted. The parameter processing strategy shown in Table 2 is adopted, but the ISB are processed by schemes 1, 2, 3, and 4, respectively. Then, the convergence time and positioning precision of the four schemes are compared. Figs. 11 and 12 respectively show the averaged positioning precision and convergence time of each station under the four schemes from June 24 to 30, 2018 in the open sky environment. Table 6 lists the statistical results of the positioning precision and convergence time of each day under the four schemes.

Figure 11: The averaged positioning precision of each station under the four schemes in the open sky environment

Figure 12: The averaged convergence time of each station under the four schemes in the open sky environment

From Figs. 11, 12, and Table 6, the following observations can be made: (a) Compared with the results of scheme 1, in the open sky environment, the additional BDS-3 observations in scheme 2 improve the positioning precision by 11%, 12%, and 5% and reduce the convergence time by 8%, 10%, and 8% in the N, E, and U components, respectively; (b) compared with the results of scheme 2, the additional PISB in scheme 3 lead to a slight improvement (0%, 6%, and 5% in the N, E, and U components) of the positioning precision and an obvious improvement (5%, 14%, and 3%, in the N, E, and U components) of the convergence time; (c) compared with the results of scheme 3, the selective utilization of PISB in scheme 4 does not improve the positioning precision, but the convergence time can be further reduced (5%, 4%, and 3%, in the N, E, and U components); (d) it should be noted that the precision of the PISB is relatively high (0.1∼0.6 ns) in our experiments. Therefore, the difference in the positioning precision between the results of scheme 3 and scheme 4 is not significant. However, if the precision of PISB is very poor, a system error will be introduced to the positioning results, and the PISB processing strategy of scheme 3 will be aggressive.

In the blocked sky environment, the same observations and data processing strategy as those in the open sky environment are used. Meanwhile, the satellite cut-off elevation angle is set to 30∘ to simulate the blocked sky environment. Fig. 13 shows the number of visible satellites and the position dilution of precision (PDOP) at the CLGY, WUH1, DWIN, and ZHON stations on June 24, 2018. Figs. 14 and 15 respectively illustrate the averaged positioning precision and convergence time of the four schemes at each station from June 24 to 30, 2018 in the blocked sky environment. Table 7 lists the statistical results of the averaged positioning precision and convergence time of all stations in the N, E, and U components from June 24 to 30, 2018 in the blocked sky environment.

Figure 13: The number of visible satellites and PDOP values at the CLGY, WUH1, DWIN, and ZHON stations on June 24, 2018, in the blocked sky environment

Figure 14: The averaged positioning precision of each station under the four schemes in the blocked sky environment

Figure 15: The averaged convergence time of each station under the four schemes in the blocked sky environment

From Figs. 14, 15, and Table 7, the following observations can be made: (a) the use of additional BDS-3 observations contribute more to the positioning precision and the convergence time in the blocked sky environment than that in the open sky environment. Compared with scheme 1, scheme 2 improves the positioning precision and the convergence time by 22%, 18%, 9% and 16%, 9% and 9% in the N, E, and U components, respectively; (b) the effect of the PISB on the positioning precision is still not obvious, but the improvements of the positioning precision by scheme 3 (0%, 6%, and 6% in the N, E, and U components) and scheme 4 (0%, 6%, and 3% in the N, E, and U components) in the blocked sky environment are larger than those in the open sky environment. This is because the number of observed satellites is significantly decreased in the blocked sky environment; (c) the effect of the PISB on the convergence time in the blocked sky environment is more significant than that in the open sky environment. Compared with scheme 2 and scheme 3, the convergence time of scheme 3 and scheme 4 is improved by 8%, 17%, 10%, and 8%, 9%, 4% in the N, E, and U components, respectively; (d) generally, the use of additional BDS-3 observations, the PISB, and the selective utilization of PISB in the blocked sky environment contributes to more significant improvement of the positioning precision and convergence time of MGPPP than that in the open sky environment.

5  Conclusions

The observations of Multi-GNSS can significantly improve the availability, continuity, reliability, and precision of navigation and positioning services. However, the influence of ISB in multi-GNSS data fusion processing must be considered and investigated, and it has become a hotspot in the field of GNSS research. This study focuses on the short-term ISB modeling and prediction methods, as well as the processing strategy of predicted ISB for GPS, BDS-2, and BDS-3 combined PPP. Based on the theoretical analysis and actual computation and comparison results, the following conclusions can be made:

(1)   There is no obvious difference of ISB between BDS-2 and BDS-3 observations when B1I and B3I signals are used. However, an obvious difference of ISB can be found between BDS-2 and BDS-3 observations when B1I/B3I and B1C/B2a signals are used.

(2)   The ISB highly depend on the receiver type, and they are relatively stable (about 0.1–0.2 ns) for the same receiver for one day. However, they vary greatly day by day (about 1–2 ns) and present some systematic variations. Therefore, the short-term ISB of each station can be modeled and predicted.

(3)   However, since the time-varying characteristic of ISB is closely related to the capability and stability of the receiver hardware, it is difficult to find an effective method to accurately predict the ISB of all stations. In the experiments, based on the proposed QM+ARIMA+SW method, the prediction precision of ISB for the next day is about 0.1−0.6 ns for different stations.

(4)   To effectively utilize the predicted ISB, this study proposes a new data processing strategy for GPS/BDS-2/BDS-3 combined PPP. The predicted ISB is used by adding two virtual observation equations, and an adaptive factor is introduced to balance the contribution of the OISB and the PISB to the final estimations of ISB. In the open sky environment, the QM+ARIMA+SW method achieves almost the same positioning precision of MGPPP as the direct utilization of the PISB and the shortest convergence time in all schemes above.

(5)   In the blocked sky environment, the proposed method achieves better performance than that in the open sky environment. Also, the proposed method improves the positioning precision and convergence time significantly. Therefore, the selective utilization of the PISB contributes to better performance than the method that does not use PISB or use PISB directly. However, the effectiveness of the proposed method should be further validated by more real GNSS kinematic data (e.g., vehicle-borne, and airborne measurements) in some challenging environments.

Funding Statement: This work was supported by “The National Key Research and Development Program of China (No. 2020YFA0713502)”, “The National Natural Science Foundation of China (No. 41874039)”, “Jiangsu National Science Foundation (No. BK20191342)”, “Fundamental Research Funds for the Central Universities (No. 2019ZDPY-RH03)”. We acknowledge TopEdit LLC for the linguistic editing and proofreading during the preparation of this manuscript.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M., Webb, F. H. (1997). Precise point positioning for the efficient and robust analysis of GPS data from large networks. Journal of Geophysical Research, 102(B3), 5005–5017. DOI 10.1029/96JB03860. [Google Scholar] [CrossRef]

 2.  Kouba, J., Héroux, P. (2001). Precise point positioning using IGS orbit and clock products. GPS Solutions, 5(2), 12–28. DOI 10.1007/PL00012883. [Google Scholar] [CrossRef]

 3.  Li, X., Ge, M., Zhang, Y., Wang, R., Xu, P. et al. (2013). New approach for earthquake/tsunami monitoring using dense GPS networks. Scientific Report, 3, 1–10. DOI 10.1038/srep02682. [Google Scholar] [CrossRef]

 4.  Yuan, Y., Zhang, K., Rohm, W., Choy, S., Norman, R. et al. (2014). Real-time retrieval of precipitable water vapor from GPS precise point positioning. Journal of Geophysical Research: Atmospheres, 119(16), 10044–10057. DOI 10.1002/2014JD021486. [Google Scholar] [CrossRef]

 5.  Sun, Z., Xu, T. (2019). The quality assessment of non-integer-hour data in gps broadcast ephemerides and its impact on the accuracy of real-time kinematic positioning over the South China Sea. Computer Modeling in Engineering & Sciences, 119(2), 263–280. DOI 10.32604/cmes.2019.04425. [Google Scholar] [CrossRef]

 6.  Zhao, Q., Pan, S., Gao, C. (2020). BDS/GPS/LEO triple-frequency uncombined precise point positioning and its performance in harsh environments. Measurement, 151, 107216. DOI 10.1016/j.measurement.2019.107216. [Google Scholar] [CrossRef]

 7.  Xiang, W., Jing, X. (2021). Remote sensing monitoring method based on BDS-based maritime joint positioning model. Computer Modeling in Engineering & Sciences, 127(2), 801–818. DOI 10.32604/cmes.2021.013568. [Google Scholar] [CrossRef]

 8.  Bisnath, S., Gao Y, B. (2007). Current state of precise point positioning and future prospects and limitations. Observing our Changing Planet, 133, 615–624. DOI 10.1007/978-3-540-85426-5. [Google Scholar] [CrossRef]

 9.  Li, X., Zhang, X., Ren, X., Fritsche, M., Wickert, J. et al. (2015). Precise positioning with current multi-constellation global navigation satellite systems: GPS, GLONASS, Galileo and BeiDou. Scientific Report, 5, 8328. DOI 10.1038/srep08328. [Google Scholar] [CrossRef]

10. Montenbruck, O., Steigenberger, P., Prange, L., Deng, Z., Zhao, Q. et al. (2017). The multi-GNSS experiment (MGEX) of the international GNSS service (IGS)–achievements, prospects, and challenges. Advances in Space Research, 59(7), 1671–1697. DOI 10.1016/j.asr.2017.01.011. [Google Scholar] [CrossRef]

11. Zhao, X., Ge, Y., Ke, F. (2020). Investigation of real-time kinematic multi-GNSS precise point positioning with the CNES products. Measurement, 166, 108231. DOI 10.1016/j.measurement.2020.108231. [Google Scholar] [CrossRef]

12. Yang, Y., Li, J., Xu, J., Tang, J. (2011). Generalized DOPs with consideration of the influence function of signal-in-space errors. Journal of Navigation, 64(S1), 3–18. DOI 10.1017/S0373463311000415. [Google Scholar] [CrossRef]

13. Odijk, D., Teunissen, P. J. G. (2013). Characterization of between-receiver GPS-Galileo inter-system biases and their effect on mixed ambiguity resolution. GPS Solutions, 17(4), 521–533. DOI 10.1007/s10291-012-0298-0. [Google Scholar] [CrossRef]

14. Paziewski, J., Wielgosz, P. (2015). Accounting for Galileo-GPS inter-system biases in precise satellite positioning. Journal of Geodesy, 89(1), 81–93. DOI 10.1007/s00190-014-0763-3. [Google Scholar] [CrossRef]

15. Lou, Y., Zheng, F., Gu, S., Wang, C., Guo, H. et al. (2016). Multi-GNSS precise point positioning with raw single-frequency and dual-frequency measurement models. GPS Solutions, 20(4), 849–862. DOI 10.1007/s10291-015-0495-8. [Google Scholar] [CrossRef]

16. Zhang, F., Liu, C., Xiao, G., Zhang, X., Feng, X. (2020). Estimating and analyzing long-term multi-GNSS inter-system bias based on uncombined PPP. Sensors, 20(5), 1499. DOI 10.3390/s20051499. [Google Scholar] [CrossRef]

17. Rabbou, M., EI-Rabbany, A. (2017). Performance analysis of precise point positioning using multi-constellation GNSS: GPS, GLONASS, Galileo and BeiDou. Survey Review, 49(352), 39–50. DOI 10.1080/00396265.2015.1108068. [Google Scholar] [CrossRef]

18. Liu, T., Yuan, Y., Zhang, B., Wang, N., Tan, B. et al. (2017). Multi-GNSS precise point positioning (MGPPP) using raw observations. Journal of Geodesy, 91, 253–268. DOI 10.1007/s00190-016-0960-3. [Google Scholar] [CrossRef]

19. Liu, G., Guo, F., Wang, J., Du, M., Qu, L. (2020). Triple-frequency GPS un-differenced and uncombined PPP ambiguity resolution using observable-specific satellite signal biases. Remote Sensing, 12(14), 2310. DOI 10.3390/rs12142310. [Google Scholar] [CrossRef]

20. Akram, A., Ahmed, E. (2016). Precise point positioning using triple GNSS constellations in various modes. Sensors, 16(6), 2–14. DOI 10.3390/s16060779. [Google Scholar] [CrossRef]

21. Chen, J., Zhang, Y., Wang, J., Yang, S., Dong, D. (2015). A simplified and unified model of multi-GNSS precise point positioning. Advances in Space Research, 55(1), 125–134. DOI 10.1016/j.asr.2014.10.002. [Google Scholar] [CrossRef]

22. Chen, J., Wang, J., Zhang, Y., Yang, S., Chen, Q. (2016). Modeling and assessment of GPS/BDS combined precise point positioning. Sensors, 16(7), 1151–1164. DOI 10.3390/s16071151. [Google Scholar] [CrossRef]

23. Choi, B., Cho, C., Cho, J., Lee, S. (2015). Multi-GNSS standard point positioning using GPS, GLONASS, BeiDou and QZSS measurements recorded at MKPO reference station in South Korea. Journal of Positioning, Navigation and Timing, 4(4), 205–211. DOI 10.11003/JPNT.2015.4.4.205. [Google Scholar] [CrossRef]

24. Jiang, N., Xu, Y., Xu, T., Xu, G., Sun, Z. (2017). GPS/BDS short-term ISB modeling and prediction. GPS Solutions, 21(1), 163–175. DOI 10.1007/s10291-015-0513-x. [Google Scholar] [CrossRef]

25. Nadarajah, N., Teunissen, P. J. G., Raziq, N. (2013). BeiDou Inter-satellite-type bias evaluation and calibration for mixed receiver attitude determination. Sensors, 13(7), 9435–9463. DOI 10.3390/s130709435. [Google Scholar] [CrossRef]

26. Zhao, Q., Wang, C., Guo, J., Wang, B., Liu, J. (2018). Precise orbit and clock determination for BeiDou-3 experimental satellites with yaw attitude analysis. GPS Solutions, 22(4), 1–13. DOI 10.1007/s10291-017-0673-y. [Google Scholar] [CrossRef]

27. Hu, C., Wang, Q., Wang, Z., Moraleda, A. (2018). New-generation BeiDou (BDS-3) experimental satellite precise orbit determination with an improved cycle-slip detection and repair algorithm. Sensors, 18(5), 1402. DOI 10.3390/s18051402. [Google Scholar] [CrossRef]

28. Jiao, G., Song, S., Jiao, W. (2020). Improving BDS-2 and BDS-3 joint precise point positioning with time delay bias estimation. Measurement Science and Technology, 31(2), 025001. DOI 10.1088/1361-6501/ab41cf. [Google Scholar] [CrossRef]

29. Mi, X., Sheng, C., El, M., Zhang, B. (2021). Characteristics of receiver-related biases between BDS-3 and BDS-2 for five frequencies including inter-system biases, differential code biases, and differential phase biases. GPS Solutions, 25(3), 113. DOI 10.1007/s10291-021-01151-w. [Google Scholar] [CrossRef]

30. The IGS final products (2020). http://www.igs.org/products. [Google Scholar]

31. Zhang, H., Hao, J., Liu, W., Yu, H., Tian, Y. (2017). A kalman filter method for BDS/GPS short-term ISB modeling and Prediction. Acta Geodaetica et Cartographics Sinica, 45(S2), 31–38. DOI 10.11947/j.AGCS.2016.F023. [Google Scholar] [CrossRef]

32. Shang, R., Gao, C., Gao, W., Zhang, R., Peng, Z. (2021). Multi-GNSS inter-system model for complex environments based on optimal state estimation. Measurement Science and Technology, 32(5), 054006. DOI 10.1088/1361-6501/abdae5. [Google Scholar] [CrossRef]

33. George, E. P. B., Gwilym, M. J. (1970). Time series analysis forecasting and control. Journal of Time Series Analysis, 3(3228), 88–94. DOI 10.1111/jtsa.12194. [Google Scholar] [CrossRef]

34. Wang, Q., Hu, C., Zhang, K. (2019). A BDS-2/BDS-3 integrated method for ultra-rapid orbit determination with the aid of precise satellite clock offsets. Remote Sensing, 11(15), 1758. DOI 10.3390/rs11151758. [Google Scholar] [CrossRef]

35. Yang, Y., He, H., Xu, G. (2001). Adaptively robust filtering for kinematic geodetic positioning. Journal of Geodesy, 75(2–3), 109–116. DOI 10.1007/s001900000157. [Google Scholar] [CrossRef]

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.