The Preliminary Results for Five-System Ultra-Rapid Precise Orbit Determination of the One-Step Method Based on the Double-Difference Observation Model

Abstract

The predicted parts of ultra-rapid orbits are important for (near) real-time Global Navigation Satellite System (GNSS) precise applications; and there is little research on GPS/GLONASS/BDS/Galileo/QZSS five-system ultra-rapid precise orbit determination; based on the one-step method and double-difference observation model. However; the successful development of a software platform for solving five-system ultra-rapid orbits is the basis of determining and analyzing these orbits. Besides this; the different observation models and processing strategies facilitate to validate the reliability of the various ultra-rapid orbits. In this contribution; this paper derives the double-difference observation model of five-system ultra-rapid precise orbit determination; based on a one-step method; and embeds this method and model into Bernese v5.2; thereby forming a new prototype software platform. For validation purposes; 31 days of real tracking data; collected from 130 globally-distributed International GNSS Service (IGS) multi-GNSS Experiment (MGEX) stations; are used to determine a five-system ultra-rapid precise orbit. The performance of the software platform is evaluated by analysis of the orbit discontinuities at day boundaries and by comparing the consistency with the MGEX orbits from the Deutsches GeoForschungsZentrum (GFZ); between the results of this new prototype software platform and the ultra-rapid orbit provided by the International GNSS Monitoring and Assessment System (iGMAS) analysis center (AC) at the Institute of Geodesy and Geophysics (IGG). The test results show that the average standard deviations of orbit discontinuities in the three-dimension direction are 0.022; 0.031; 0.139; 0.064; 0.028; and 0.465 m for GPS; GLONASS; BDS Inclined Geosynchronous Orbit (IGSO); BDS Mid-Earth Orbit (MEO); Galileo; and QZSS satellites; respectively. In addition; the preliminary results of the new prototype software platform show that the consistency of this platform has been significantly improved compared to the software package of the IGGAC.

1. Introduction

The launches of the first Galileo In-Orbit Validation Element (GIOVE) satellite in 2005 and the first experimental satellite of the Chinese regional BeiDou navigation satellite system (BDS-2) in 2007 marked the beginning of a new era in satellite navigation [1]. Currently, the navigation satellite systems that already provide position, navigation, and timing services include American GPS, Russian GLONASS, European Galileo, Chinese BDS, Japanese Quasi-Zenith Satellite System (QZSS), and so on. Montenbruck et al. have also analyzed the signal-in-space ranging errors of the broadcast ephemeris messages of these five navigation satellite systems [2]. By the end of January 2018, the GPS/GLONASS/BDS/Galileo/QZSS constellation was composed of 31/24/14/15/3 satellites, respectively, with operational capability. Each navigation satellite system also has its own future constellation construction plan. For users, there will be more satellites available in the future, and these satellites will also have better performance (e.g., in 2020, the Chinese global BeiDou navigation satellite system (BDS-3) will be composed of five geostationary earth orbit (GEO), three inclined geosynchronous orbit (IGSO), and twenty-seven medium earth orbit (MEO) satellites). Recently, BDS-3 has entered the intensive launch phase, and two BDS-3 satellites were launched on October 15, 2018. At the end of October 2018, there are 18 BDS-3 satellites being tested (http://www.csno-tarc.cn/system/constellation). The construction and in-orbit operation of multiple navigation satellite systems can provide more reliable and accurate results for global navigation satellite system (GNSS) applications (e.g., atmospheric remote sensing [3,4,5], precise point positioning [6,7,8], precise relative positioning [9,10], and the calculation of Earth rotation parameters [11]). At the same time, the multiple constellations also bring great challenges to the data processing of multi-GNSS. In order to promote the scientific and engineering applications of multi-GNSS in high-precision fields, the International GNSS Service (IGS) initiated the multi-GNSS experiment (MGEX) in 2011 (ftp://igs.org/pub/resource/pubs/IGS M-GEX VF.pdf). Generally, seven MGEX analysis centers (ACs) provide special multi-GNSS precise orbit products to MGEX. These ACs are the Centre National d’Etudes Spatiales (CNES) and Collecte Localisation Satellites (CLS) in France; the Center for Orbit Determination in Europe (CODE) in Switzerland; the Deutsches GeoForschungsZentrum (GFZ), the Technische Universität München (TUM), and the European Space Agency (ESA) in Germany; Wuhan University (WU) in China, and the Japan Aerospace Exploration Agency (JAXA) in Japan. Among them, CODE, GFZ, and WU can provide five-system final precise orbit products, and ESA have only provided two-week final precise orbit products in the past.

Because the accuracy of satellite orbits directly affects the accuracy of GNSS user position, navigation, and timing applications, GNSS precise orbit determination (POD) technology has become one of the key technologies in some GNSS research institutes (e.g., MGEX ACs). The techniques for the precise orbit determination of GNSS satellites are mainly based on microwave GNSS observations and satellite laser ranging (SLR) data [12]. For these reasons, these institutes have independently developed or redeveloped corresponding GNSS POD software platforms. In general, the GNSS POD software platforms used by CNES/CLS, CODE, GFZ, TUM, ESA, WU, and JAXA are POD GINS [13], Bernese v5.3 [14], EPOS.P8 [15], Bernese v5.0 [16], NAPEOS 3.8 [17], PANDA [18], and MADOCA [19], respectively. In this paper, we mainly consider the processing method of microwave GNSS observations, and the process of GNSS POD; the error equations are established by using non-difference observations in all MGEX ACs, except for the double-difference observations used by Bernese v5.3. In addition, the solution strategies used by MGEX ACs in calculating multi-GNSS orbit parameters are not exactly the same: CNES/CLS, CODE [20], and GFZ use the one-step method strategy; TUM use the two-step method strategy [21]; and WU use the three-step method strategy. Studies have shown that more available satellites benefit the estimation of common parameters (such as station coordinates and tropospheric parameters) when using the one-step method. Theoretically, the one-step method strategy is a more rigorous solution, which is more conductive to multi-GNSS POD [22]. IGS MGEX also encourages MGEX ACs to use the one-step method. At present, IGS provides GPS/GLONASS two-system ultra-rapid orbits, and the International GNSS Monitoring and Assessment System (iGMAS) provides four-system (GPS/GLONASS/BDS/Galileo) ultra-rapid orbits (http://www.igmas.org/). Besides, MGEX provides final precise orbits, and all MGEX ACs are testing their own ultra-rapid orbits. However, these ultra-rapid orbits are not fully released for GNSS users, and the published research and application of ultra-rapid orbits are mainly focused on a single or dual system [23,24]. Two other examples also show this problem: Kevin K. Choi et al. evaluated the prediction accuracy of IGS ultra-rapid orbit products [25]; and IGS ACs, such as CODE, introduced their solution strategies for IGS ultra-rapid orbit products [26,27]. Recently, based on the one-step method and non-difference observations, Li et al. proposed a method to solve the five-system ultra-rapid orbit problem [28]. There is relatively little data available at present regarding the introduction and research of a multi-GNSS (especially for five systems) ultra-rapid orbit based on the one-step method and the error equation of double-difference observations. On the other hand, a diversity of observation models and processing strategies will facilitate to validate the ultra-rapid orbits offered by various POD software platforms. Because the fact that the predicted part of the ultra-rapid orbits is important for a (near) real-time precise service, it is necessary to analyze these five-system ultra-rapid precise orbits provided by the one-step method and the error equation of double-difference observations. Developing a multi-GNSS POD software platform based on this method and error equation will be the premise of the corresponding analysis. However, Bernese v5.3 has not been released publicly, and Bernese v5.2 can only handle GPS/GLONASS/Galileo three-system data.

Therefore, this paper will derive the double-difference observation model of a multi-GNSS ultra-rapid POD based on the one-step method. A new prototype software platform for a five-system ultra-rapid POD will be developed on the basis of the Bernese v5.2 framework, and this model and method. The structure of this paper is as follows: First, the multi-GNSS observation model and the processing strategies of the software platform are introduced in the second section. Then, the preliminary results of this software platform are shown in the third section. Finally, the discussion and conclusions are given.

2. Materials and Methods

2.1. The Multi-GNSS Observation Model of the New Prototype Software Platform

In this paper, during the processing of the ultra-rapid POD solution, dual-frequency ionosphere-free (IF) combined observations were used to eliminate the first order of ionospheric delays. Hadas et al. showed that the impact of higher order ionosphere delays can reach several millimeters for the radial, along-track, and cross-track orbit components [29]. Considering that this paper only plans to get a preliminary result for five-system ultra-rapid orbits, the impact of higher order ionosphere delays was not considered. The observation equation of the IF combined observations can be linearized and expressed as follows:

$$\begin{array}{l}{l}_{r,IF}^s={\mu }_r^s·\left(\Phi {\left(t,t_0\right)}^s·o_o^s-R_{r,erp}·r_r\right)-\Delta t_{IF}^s+\Delta t_{r,IF}+{\partial }_r^s·{\delta }_{erp}\\ +{\lambda }_{IF}\left(b_{r,IF}-b_{IF}^s+N_{r,IF}^s\right)+m_{r,trop}^s·{\delta }_{r,trop}+{\epsilon }_{r,IF}^s\end{array}$$

(1)

$$\begin{array}{l}{p}_{r,IF}^s={\mu }_r^s·\left(\Phi {\left(t,t_0\right)}^s·o_o^s-R_{r,erp}·r_r\right)-\Delta t_{IF}^s+\Delta t_{r,IF}+{\partial }_r^s·{\delta }_{erp}\\ +c\left(d_{r,IF}-d_{IF}^s\right)+m_{r,trop}^s·{\delta }_{r,trop}+{\omega }_{r,IF}^s\end{array}$$

(2)

where

$$o_o^s={\left(a,e,i,\Omega ,\omega ,u_0,{\mathrm{p}}_{1~9}^s\right)}^T$$

(3)

$${\delta }_{erp}={\left(x_0,y_0,{\theta }_0\right)}^T$$

(4)

in addition, $l_{r,IF}^s$ and $p_{r,IF}^s$ denote the “observed minus computed” IF phase and pseudo-range observables from satellite $s$ to receiver $r$, respectively; ${\mu }_r^s$ represents the unit vector in the direction from receiver $r$ to satellite $s$; and $\Phi {\left(t,t_0\right)}^s$ is the state transition matrix from the initial epoch $t_0$ to current epoch $t$. $o_o^s$ denotes the initial orbit state vector of satellite $s$, including the semimajor axis of the orbit $a$, the eccentricity of the orbit $e$, the inclination of the orbital plane with respect to the equatorial plane $i$, the right ascension of the ascending node $\Omega$, the argument of perigee $\omega$, the argument of latitude $u_0$, and the solar radiation pressure (SRP) parameters of the empirical CODE orbit model (ECOM) developed by the Astronomical Institute, University of Bern $\left({\mathrm{p}}_{1~9}^s\right)$ [30,31]. $R_{r,erp}$ denotes the satellite coordinates rotation matrix from the Earth-fixed frame to the inertial frame, according to the Earth rotation parameters (ERP); $r_r$ is the receiver position increments relative to the a priori position; $\Delta t_{IF}^s$ is the IF combined clock bias of satellite $s$; $\Delta t_{r,IF}$ is the IF combined clock bias of receiver $r$; ${\partial }_r^s$ is the partial derivative of the IF combined observation to ERP; ${\delta }_{erp}$ is the ERP, including polar motion $\left(x_0,y_0\right)$ and length of day change ${\theta }_0$; $c$ is the speed of light in vacuo; ${\lambda }_{IF}$ is the wave-length of the IF combination; $b_{r,IF}$ and $b_{IF}^s$ are the receiver $r$ and satellite $s$ phase delays of IF combination, respectively; $N_{r,IF}^s$ is the ambiguity of the IF combination; $d_{r,IF}$ and $d_{IF}^s$ are the receiver $r$ and satellite $s$ code biases of the IF combination; $m_{r,trop}^s$ is a related mapping function with respect to the zenith tropospheric wet delays ${\delta }_{r,trop}$; and ${\epsilon }_{r,IF}^s$ and ${\omega }_{r,IF}^s$ are the sum of measurement noise and multipath effects for the IF combined phase and pseudo-range observations, respectively.

In the same navigation satellite system, the observation equation of the linearized double-difference IF combination observation is as follows:

$$\begin{array}{l}\Delta \nabla l_{r1r2,IF}^{s1s2}=\left(\Delta {\mu }_{r1r2}^{s1}·\Phi {\left(t,t_0\right)}^{s1}·o_o^{s1}-\Delta {\mu }_{r1r2}^{s2}·\Phi {\left(t,t_0\right)}^{s2}·o_o^{s2}\right)\\ +{\lambda }_{IF}·\Delta \nabla N_{r1r2,IF}^{s1s2}+\Delta \nabla {\partial }_{r1r2}^{s1s2}·{\delta }_{erp}\\ -\left(\Delta {\mu }_{r1}^{s1s2}·R_{r1,erp}·r_{r1}-\Delta {\mu }_{r2}^{s1s2}·R_{r2,erp}·r_{r2}\right)\\ +\left(\Delta m_{r1,trop}^{s1s2}·{\delta }_{r1,trop}-\Delta m_{r2,trop}^{s1s2}·{\delta }_{r2,trop}\right)+\Delta \nabla {\epsilon }_{r1r2,IF}^{s1s2}\end{array}$$

(5)

$$\begin{array}{l}\Delta \nabla p_{r1r2,IF}^{s1s2}=\left(\Delta {\mu }_{r1r2}^{s1}·\Phi {\left(t,t_0\right)}^{s1}·o_o^{s1}-\Delta {\mu }_{r1r2}^{s2}·\Phi {\left(t,t_0\right)}^{s2}·o_o^{s2}\right)\\ +\Delta \nabla {\partial }_{r1r2}^{s1s2}·{\delta }_{erp}\\ -\left(\Delta {\mu }_{r1}^{s1s2}·R_{r1,erp}·r_{r1}-\Delta {\mu }_{r2}^{s1s2}·R_{r2,erp}·r_{r2}\right)\\ +\left(\Delta m_{r1,trop}^{s1s2}·{\delta }_{r1,trop}-\Delta m_{r2,trop}^{s1s2}·{\delta }_{r2,trop}\right)+\Delta \nabla {\omega }_{r1r2,IF}^{s1s2}\end{array}$$

(6)

where $\Delta \nabla l_{r1r2,IF}^{s1s2}$ is the double-difference IF combined phase of $l_{r,IF}^s$ between receivers ($r1,r2$) and satellites ($s1,s2$), and $\Delta \nabla p_{r1r2,IF}^{s1s2}$ has an analogous meaning; $\Delta {\mu }_{r1}^{s1s2}$ is the computed vector ${\mu }_{r1}^{s1}$ minus ${\mu }_{r1}^{s2}$, and $\Delta {\mu }_{r2}^{s1s2}$ has an analogous meaning; $\Delta \nabla N_{r1r2,IF}^{s1s2}$ is the double-difference IF combined ambiguity; $\Delta \nabla {\partial }_{r1r2}^{s1s2}$ is the double-difference partial derivative of ${\partial }_r^s$; $\Delta m_{r1,trop}^{s1s2}$ is the computed mapping function given by $\Delta m_{r1,trop}^{s1}$ minus $\Delta m_{r1,trop}^{s2}$, and $\Delta m_{r2,trop}^{s1s2}$ has an analogous meaning; $\Delta \nabla {\epsilon }_{r1r2,IF}^{s1s2}$ is the double-difference value of ${\epsilon }_{r,IF}^s$ between receivers ($r1,r2$) and satellites ($s1,s2$), and $\Delta \nabla {\omega }_{r1r2,IF}^{s1s2}$ has an analogous meaning.

Based on the linearized double-difference IF combined observation equation, the five system (GPS/GLONASS/BDS/Galileo/QZSS) ultra-rapid POD observation model using the one-step method is obtained as follows:

$$\begin{array}{l}\Delta \nabla l_{r1r2,IF}^{S1S2}=\left(\Delta {\mu }_{r1r2}^{S1}·\Phi {\left(t,t_0\right)}^{S1}·o_o^{S1}-\Delta {\mu }_{r1r2}^{S2}·\Phi {\left(t,t_0\right)}^{S2}·o_o^{S2}\right)\\ +{\lambda }_{IF}·\Delta \nabla N_{r1r2,IF}^{S1S2}+\Delta \nabla {\partial }_{r1r2}^{S1S2}·{\delta }_{erp}\\ -\left(\Delta {\mu }_{r1}^{S1S2}·R_{r1,erp}·r_{r1}-\Delta {\mu }_{r2}^{S1S2}·R_{r2,erp}·r_{r2}\right)\\ +\left(\Delta m_{r1,trop}^{S1S2}·{\delta }_{r1,trop}-\Delta m_{r2,trop}^{S1S2}·{\delta }_{r2,trop}\right)+\Delta \nabla {\epsilon }_{r1r2,IF}^{S1S2}\end{array}$$

(7)

$$\begin{array}{l}\Delta \nabla p_{r1r2,IF}^{S1S2}=\left(\Delta {\mu }_{r1r2}^{S1}·\Phi {\left(t,t_0\right)}^{S1}·o_o^{S1}-\Delta {\mu }_{r1r2}^{S2}·\Phi {\left(t,t_0\right)}^{S2}·o_o^{S2}\right)\\ +\Delta \nabla {\partial }_{r1r2}^{S1S2}·{\delta }_{erp}\\ -\left(\Delta {\mu }_{r1}^{S1S2}·R_{r1,erp}·r_{r1}-\Delta {\mu }_{r2}^{S1S2}·R_{r2,erp}·r_{r2}\right)\\ +\left(\Delta m_{r1,trop}^{S1S2}·{\delta }_{r1,trop}-\Delta m_{r2,trop}^{S1S2}·{\delta }_{r2,trop}\right)+\Delta \nabla {\omega }_{r1r2,IF}^{S1S2}\end{array}$$

(8)

where

$$S=\left(G or R_k or C or E or J\right),$$

(9)

besides, $G$, $R_k$, $C$, $E$, and $J$ refer to GPS, GLONASS, BDS, Galileo, and QZSS, respectively; and $k$ is the frequency factor in the frequency division multiple access (FDMA) mode.

Thus, in the five-system ultra-rapid POD solution, the unknown parameters, $X_{pod}$, which need to be estimated include the initial orbit state vector $o_o^s$, the receiver position increments $r_r$, the double-difference IF combined ambiguity $\Delta \nabla N_{r1r2,IF}^{s1s2}$, the ERP ${\delta }_{erp}$, and the zenith tropospheric wet delays ${\delta }_{r,trop}$, which are summarized in the following vector:

$$X_{pod}=\left(o_o^{S1}\dots o_o^{Sm},r_r,\Delta \nabla N_{r1r2,IF}^{S1S2}\dots \Delta \nabla N_{r1r2,IF}^{SmSn},{\delta }_{erp},{\delta }_{r,trop}\right)$$

(10)

In addition, in the data processing of Multi-GNSS observations in this paper, elevation-dependent weighting was adopted, and the weight function $w\left(z\right)$ is defined as follows [32]:

$$w\left(z\right)=\{\begin{array}{cc}\hfill 10000\times {\cos }^2\left(z\right)for& phase\hfill \\ \hfill {\cos }^2\left(z\right)for& code\hfill \end{array}$$

(11)

where $z$ is the zenith angle of the satellite.

The above models are the important basis of the new prototype software platform for the five-system ultra-rapid POD. The following shows the processing strategies of the new prototype software platform.

2.2. The Processing Strategies

The new prototype software platform was developed on the basis of the Bernese v5.2 framework, and the above model and method, so the platform should be developed first to be able to simultaneously handle both five satellite systems and GNSS measurements stored in version 3 of the receiver independent exchange data format (RINEX) [33].

2.2.1. Coordinate and Time Systems of the New Prototype Software Platform

A unified coordinate system and time system are needed before five-system ultra-rapid POD. For the coordinate system, GPS, GLONASS, BDS, Galileo, and QZSS use the WGS84, PZ-90, CGCS2000, GTRF, and GRS80 ellipsoid systems [34,35,36,37,38], respectively. These coordinate systems are defined to approach the International Terrestrial Reference System (ITRS). According to Section 2.1, the double-difference observation model was used in our software platform, and so the differences between WGS84 and other coordinates can be ignored, except PZ-90 [39]. The transformation between WGS84 and PZ-90 can be realized by the Helmert seven-parameter model, and the transformation model is shown below:

$$\left(\begin{array}{c}{X}_{GPS}\\ Y_{GPS}\\ Z_{GPS}\end{array}\right)=\left(\begin{array}{c}\Delta X\\ \Delta Y\\ \Delta Z\end{array}\right)+\left(1+\alpha \right)\left(\begin{array}{ccc}1& \Delta \omega & -\Delta \phi \\ -\Delta \omega & 1& \Delta \epsilon \\ \Delta \phi & -\Delta \epsilon & 1\end{array}\right)\left(\begin{array}{c}{X}_{GLO}\\ Y_{GLO}\\ Z_{GLO}\end{array}\right)$$

(12)

where ${\left(\begin{array}{ccc}{X}_{GPS}& Y_{GPS}& Z_{GPS}\end{array}\right)}^T$ and ${\left(\begin{array}{ccc}{X}_{GLO}& Y_{GLO}& Z_{GLO}\end{array}\right)}^T$ are the coordinates in the WGS84 and PZ-90 systems, respectively. Furthermore, $\left(\begin{array}{ccccccc}\Delta X& \Delta Y& \Delta Z& \alpha & \Delta \omega & \Delta \phi & \Delta \epsilon \end{array}\right)$ are the seven transformation parameters (ftp://ftp.aiub.unibe.ch/BSWUSER52/GEN/DATUM).

For the time system, GPS, GLONASS, BDS, Galileo, and QZSS use the GPST, GLNT, BDT, GST, and QZSST systems, respectively. These time systems conform to UTC (USNO), UTC (SU), UTC (BSNC), UTC (PTB), and UTC (NICT), respectively. Their transformation model is shown below [34,35,36,37,38,39]:

$$\{\begin{array}{l}GPST=GLNT-3h+\Delta t_{LS}+{\delta }_{GPST-GLNT}\\ GPST=BDT+1356weeks+14s+{\delta }_{GPST-BDT}\\ GPST=GST+{\delta }_{GPST-GST}\\ GPST=QZSST+{\delta }_{GPST-QZSST}\end{array},$$

(13)

where $\Delta t_{LS}$ is the leap second of $GLNT$. ${\delta }_{GPST-GLNT}$, ${\delta }_{GPST-BDT}$, ${\delta }_{GPST-GST}$, and ${\delta }_{GPST-QZSST}$ are the micro-deviations between GPST and the GLNT, BDT, GST, and QZSST systems. In our software platform, these micro-deviations are ignored.

2.2.2. Partial Available Strategies and the Flowchart of the New Prototype Software Platform

There is research that has showed that the change of satellite orbit model can remarkably change the results of Galileo satellite POD [40], and Montenbruck et al. summarized the five-system satellite attitude models [41]. In this paper, the yaw-steering is assumed for all satellites, which is not yet correct for BDS and QZSS satellites. Besides, the ECOM 5-parameter model will be used for all satellites in the next section, but the model was not designed for Galileo, BDS, and QZSS satellites. The quality of the ultra-rapid orbits mostly depends on the satellite attitude models and the SRP model; therefore, the state of the new prototype software platform can be considered as preliminary.

Furthermore, during the processing of the five-system ultra-rapid POD in our new prototype software platform, the length of the arc and the setting of the stochastic pulses referred to the method proposed by Simon Lutz et al. [27]. The 06 (UTC), 12 (UTC), 18 (UTC), and 24 (UTC) ultra-rapid orbits used 54-h, 60-h, 66-h, and 72-h observations, respectively. Besides this, 3–5 sets of stochastic pulses were used in each arc. The cut-off angle was set to 3 degrees, and the data sampling rate was 180 s. The available partial processing strategies are shown in

Table 1. The partial strategies of the new prototype software platform. IF: Ionosphere-free.

Items	Description
Processing scheme	IF double-difference network processing
Solar radiation	ECOM 5-parameter model
Sampling interval	180 s
Cut-off angle	3 degrees
Tropospheric delay	GPT model, GMF mapping function, two-hourly ZTDs [42,43]
A priori ERP	IERS Bulletin A and continuous piecewise linear function [44]
Attitude	Yaw-steering always assumed for all sat
Ambiguity	Real constant value for each ambiguity arc
PCO/PCV(Sat)	GPS, GLONASS: PCO and PCV from IGS; Galileo: PCO and PCV from ESA; BeiDou, QZSS: PCO from IGS MGEX
PCO/PCV(Rec.)	GPS, GLONASS: IGS PCO and PCV; new GNSS: adopted from GPS L1 and L2
Update frequency	6 h
A priori orbits	Broadcast ephemeris
A priori coordinates	Previous 3-day solution, strong constraint (0.01, 0.01, 0.01) m
Earth gravity	EGM2008 12 × 12 [45]
Solid earth tide and pole tide	IERS2010 [46]
Ocean tide	FES2004 [47]
Nutation model	IAU2000 [48]
N-body gravitation	Sun, moon, and other planets (DE405)
Relativity effect	IERS2010 [46]

, and

Table 2. The observation types selected by the new prototype software platform. Sys: System; G: GPS; R: GLONASS (k: the frequency factor); C: BDS; E: Galileo; J: QZSS.

Sys	Frequency		Signal
	Band	MHz	Phase/Code
G	L1	1575.42	L1C/C1C	L1W/C1W
	L2	1227.60	L2W/C2W L2S/C2S L2C/C2C	L2X/C2X L2L/C2L L2D/C2D
$R_k$	G1	1598.0625	L1C/C1C
		−1608.75	L1P/C1P
	G2	1242.9375	L2P/C2P
		−1251.25	L2C/C2C
E	E1	1575.42	L1X/C1X	L1C/C1C
	E5a	1176.45	L5X/C5X	L5Q/C5Q
C	B1	1561.098	L1I/C1I	L2I/C2I
			L2X/C2X
	B2	1207.14	L7I/C7I	L7X/C7X
J	L1	1575.42	L1C/C1C	L1X/C1X
	L2	1227.60	L2X/C2X	L2L/C2L
			L2S/C2S

summarizes the observation types of each satellite navigation system selected from the observation files. In our software platform, the five-system observation and navigation files were processed simultaneously in response to MGEX’s call for the joint processing of multi-GNSS, which is also needed for the five-system ultra-rapid POD using the one-step method. In addition, the flowchart of the new prototype software platform is shown in Figure 1.

The above section introduces the models, methods, and processing strategies used in the new prototype software platform. Although the double-difference IF combined observation equations are constructed in the same navigation satellite systems, all satellites and these equations are solved together in a least square estimator.

3. Preliminary Results

In order to test the basic performance of the new prototype software platform for the five-system ultra-rapid POD, one month’s data collected from 130 globally-distributed IGS MGEX stations were used to determine and analyze the ultra-rapid orbit. Figure 2 shows the distribution of the selected tracking stations (the triangles show stations that can observe the four systems (GPS/GLONASS/BDS/Galileo), and the circles show stations that can observe the QZSS). Figure 2 shows that there are 130 stations which can observe four systems. Among them, 67 stations can observe all five systems. Figure 2 also shows that the tracking stations of QZSS have regional distribution. To facilitate the comparison with the MGEX orbit products, we unified the coordinate system to WGS84 and we unified the time system to GPST.

The test period was from 1 January 2018 to 31 January 2018. In this paper, the ultra-rapid orbits of MEO and IGSO satellites were tested. The methods to assess the predicted parts of multi-GNSS ultra-rapid orbits are mainly the analysis of the orbit discontinuities at day boundaries, comparison with IGS/MGEX orbit, and SLR validation, etc. [28,49,50]. In this section, we used the first two validation methods, and the assessment steps were carried out as follows:

• Five-system ultra-rapid orbits were estimated, and the time latency for estimated ultra-rapid orbits is shown in Figure 3;
• The six-hour predicted part of each arc were compared to the final MGEX orbit products of GFZ (the accuracy codes are 0 for all satellites in gbm19861.sp3), and the root mean square (RMS) of the position deviation in the radial, along-track, and cross-track directions orbit components was calculated. Figure A1, Figure A2, Figure A3, Figure A4 and Figure A5 of Appendix A represent the differences between the predicted orbits of five-system ultra-rapid orbits and the orbits of GFZ for GPS, GLONASS, BDS, Galileo, and QZSS satellites. Figure 4, Figure 5, Figure 6 and Figure 7 represent the root mean square (RMS) of these satellites;
• The analysis of the orbit discontinuities at day boundaries. Figure A6, Figure A7, Figure A8, Figure A9 and Figure A10 of Appendix A represent the residuals of the fitting procedure for GPS, GLONASS, BDS, Galileo, and QZSS satellites in the radial, along-track, and cross-track directions orbit components. Figure 8, Figure 9, Figure 10 and Figure 11 represent the standard deviations (SD) of these satellites.
• The Institute of Geodesy and Geophysics (IGG), as one of the iGMAS ACs, has provided ultra-rapid orbits estimated by the two-step method [51]; therefore, the six-hour predicted part of each arc for these orbits (named IGGAC in this paper) were compared to the final MGEX orbits of GFZ in the radial, along-track, and cross-track directions orbit components. In this test period, the available ultra-rapid orbits of IGGAC included satellites of GPS and GLONASS from 1 January 2018 to 31 January 2018, and of Galileo from 14 January 2018 to 31 January 2018. Figure A11, Figure A12 and Figure A13 of Appendix A represent the differences between the predicted orbits of IGGAC and the orbits of GFZ for GPS, GLONASS, and Galileo satellites. Figure 12, Figure 13, Figure 14 and Figure 15 represent the root mean square (RMS) of these satellites.

  Table A1. The differences between the results of our new prototype software platform and IGGAC in consistency with the final MGEX orbit products of GFZ.

  Items		The New Prototype Software Platform (m)	IGGAC (m)	Improved Percentage (%)
  GPS	Radial	0.010	0.034	70.59
  	Along	0.015	0.080	81.25
  	Cross	0.013	0.056	76.79
  	3D	0.022	0.103	78.51
  GLONASS	Radial	0.014	0.040	65.00
  	Along	0.023	0.113	79.65
  	Cross	0.016	0.067	76.12
  	3D	0.031	0.137	77.19
  Galileo	Radial	0.014	0.062	77.42
  	Along	0.020	0.170	88.24
  	Cross	0.015	0.090	83.33
  	3D	0.029	0.202	85.82

  of Appendix A shows the differences between the results of our new prototype software platform and IGGAC in consistency with the final MGEX orbit products of GFZ.

Based on the above calculations, we obtained the following test results. Figure 3 shows that the average time consumption for calculation was 2.02 h; from Figure 4, Figure 5, Figure 6 and Figure 7, compared to the final MGEX orbit products of GFZ in the radial, along-track, and cross-track directions orbit components, for the six-hour predicted orbits of GPS satellites, the average RMSs were 0.027, 0.064, and 0.052 m, respectively; for the six-hour predicted orbits of GLONASS satellites, the average RMSs were 0.033, 0.098, and 0.054 m, respectively; for the six-hour predicted orbits of BDS IGSO satellites, the average RMSs were 0.106, 0.234, and 0.184 m, respectively; for the six-hour predicted orbits of BDS MEO satellites, the average RMSs were 0.074, 0.180, and 0.114 m, respectively; for the six-hour predicted orbits of Galileo satellites, the average RMSs were 0.049, 0.093, and 0.102 m, respectively; for the six-hour predicted orbits of QZSS satellites, the average RMSs were 0.594, 1.338, and 0.804 m, respectively; from Figure 8, Figure 9, Figure 10 and Figure 11, analysis of the orbit discontinuities at day boundaries in the radial, along-track, and cross-track directions orbit components, for the six-hour predicted orbits of GPS satellites, the average SDs were 0.010, 0.015, and 0.013 m respectively; for the six-hour predicted orbits of GLONASS satellites, the average SDs were 0.014, 0.023, and 0.016 m respectively; for the six-hour predicted orbits of BDS IGSO satellites, the average SDs were 0.070, 0.108, and 0.053 m, respectively; for the six-hour predicted orbits of BDS MEO satellites, the average SDs were 0.032, 0.047, and 0.030 m, respectively; for the six-hour predicted orbits of Galileo satellites, the average SDs were 0.014, 0.020, and 0.015 m, respectively; for the six-hour predicted orbits of QZSS satellites, the average SDs were 0.241, 0.349, and 0.190 m, respectively; from Figure 12, Figure 13, Figure 14 and Figure 15, comparison between IGGAC and the final MGEX orbit products of GFZ in the radial, along-track, and cross-track directions orbit components, for the six-hour predicted orbits of GPS satellites, the average RMSs were 0.034, 0.080, and 0.056 m respectively; for the six-hour predicted orbits of GLONASS satellites, the average RMSs were 0.040, 0.113, and 0.067 m, respectively; for the six-hour predicted orbits of Galileo satellites, the average RMSs were 0.062, 0.170, and 0.090 m, respectively.

4. Discussion

According to the above results of the orbit discontinuities at day boundaries, for the three-dimensional SDs, we can know that they are about 3, 6, 14, and 47 cm for GPS/GLONASS/Galileo, BDS MEO, BDS IGSO, and QZSS satellites, respectively. This shows that the performance to process IGSO satellites of this new prototype software platform needs to be improved compared with MEO satellites. Besides, for the results of step (2) in Section 3, considering the orbit assessment results of GFZ, which was based on the one-step method and no-difference observation model [52], we can know that the double-difference observation model of five systems and the no-difference observation model of five systems have good consistency. In addition, by analyzing the results of step (2) and (4) in Section 3, compared to the orbits provided by the software package of IGGAC, we can know that the consistency of the results calculated by this new prototype software platform is significantly improved in the radial, along-track, cross-track, and 3D directions.

As mentioned above, our new prototype software platform is only a preliminary state. In view of the fact that only the orbits of non-GEO satellites were tested in this paper, the performance of our new prototype software platform in processing the orbit of GEO satellites needs to be further verified. Besides this, the new prototype software platform did not build a special attitude model and SRP model for the orbit normal mode, which is worth further improvement.

5. Conclusions

In this study, we derived the double-difference observation model of the GPS/GLONASS/BDS/Galileo/QZSS five-system ultra-rapid POD based on the one-step method. Additionally, using the framework of Bernese v5.2 and this model and method, we developed a new prototype software platform for five-system ultra-rapid precise orbit calculation. Subsequently, one month of real tracking data collected from 130 globally distributed IGS MGEX stations were used to test the software platform. The results show that this new prototype software platform has the preliminary function of processing five-system ultra-rapid orbits. In terms of the analysis of the orbit discontinuities at day boundaries, the average SDs in 3D were 0.022, 0.031, 0.139, 0.064, 0.028, and 0.465 m for GPS, GLONASS, BDS IGSO, BDS MEO, Galileo, and QZSS satellites, respectively. Besides, based on the one-step method, the double-difference observation model of five systems and the no-difference observation model of five systems have good consistency. Compared to the ultra-rapid orbits of IGGAC, in 3D direction of orbit consistency with the final MGEX orbit products of GFZ, the improved percentages of this new prototype software platform were 78.51%, 77.19%, and 85.82% for GPS, GLONASS, and Galileo satellites, respectively. The consistency in other directions has also improved significantly.