Precise Orbit and Clock Products of Galileo, BDS and QZSS from MGEX Since 2018: Comparison and PPP Validation

Abstract

In recent years, the development of new constellations including Galileo, BeiDou Navigation Satellite System (BDS) and Quasi-Zenith Satellite System (QZSS) have undergone dramatic changes. Since January 2018, about 30 satellites of the new constellations have been launched and most of the new satellites have been included in the precise orbit and clock products provided by the Multi Global Navigation Satellite System (Multi-GNSS) Experiment (MGEX). Meanwhile, critical issues including antenna parameters, yaw-attitude models and solar radiation pressure models have been continuously refined for these new constellations and updated into precise MGEX orbit determination and precise clock estimation solutions. In this context, MGEX products since 2018 are herein assessed by orbit and clock comparisons among individual analysis centers (ACs), satellite laser ranging (SLR) validation and precise point positioning (PPP) solutions. Orbit comparisons showed 3D agreements of 3–5 cm for Galileo, 8–9 cm for BDS-2 inclined geosynchronous orbit (IGSO), 12–18 cm for BDS-2 medium earth orbit (MEO) satellites, 24 cm for BDS-3 MEO and 11–16 cm for QZSS IGSO satellites. SLR validations demonstrated an orbit accuracy of about 3–4 cm for Galileo and BDS-2 MEO, 5–6 cm for BDS-2 IGSO, 4–6 cm for BDS-3 MEO and 5–10 cm for QZSS IGSO satellites. Clock products from different ACs generally had a consistency of 0.1–0.3 ns for Galileo, 0.2–0.5 ns for BDS IGSO/MEO and 0.2–0.4 ns for QZSS satellites. The positioning errors of kinematic PPP in Galileo-only mode were about 17–19 mm in the north, 13–16 mm in the east and 74–81 mm in the up direction, respectively. As for BDS-only PPP, positioning accuracies of about 14, 14 and 49 mm could be achieved in kinematic mode with products from Wuhan University applied.

1. Introduction

The Multi Global Navigation Satellite System (Multi-GNSS) Experiment (MGEX) [1] project was initiated by the International GNSS Services (IGS) in 2011 to prepare the IGS for multi-GNSS processing in view of the emergence of BeiDou Navigation Satellite System (BDS), Galileo and Quasi-Zenith Satellite System (QZSS). With several years’ development, MGEX has now built up a global multi-GNSS network of more than 300 stations to collect all available GNSS signals. Based on the observations of the MGEX network, precise orbit and clock products are now routinely generated by seven MGEX analysis centers (ACs).

Among all the MGEX ACs, the Center for Orbit Determination in Europe (CODE) [2], GeoForschungsZentrum Potsdam (GFZ) [3], Technische Universität München (TUM) [4] and Wuhan University (WHU) [5] provide MGEX products for all the emerging constellations including BDS, Galileo and QZSS. The Centre National d’Etudes Spatiales (CNES) and Collecte Localisation Satellites (CLS) only provide Galileo products [6], and the Japan Aerospace Exploration Agency (JAXA) [7] only provides QZSS products. In addition, the Shanghai Observatory (SHAO) recently joined the group of MGEX ACs by providing rapid products covering Galileo and BDS [8].

In view of the variety of MGEX orbit and clock products and the absence of combined MGEX products [9], it could be confusing for multi-GNSS users to choose the most appropriate products and to understand their positioning performances. In this context, assessments of the MGEX products—particularly the new BDS, Galileo, and QZSS—are of great importance and necessity in order to properly utilize the current MGEX products and adequately prepare for the upcoming generation of combined MGEX products. Steigenberger et al. previously performed an assessment of Galileo products from CNES/CLS, CODE, GFZ and TUM for a 20-week period in 2013. They reported a general orbit consistency of 5–30 cm among these four ACs [10]. Guo et al. further assessed the quality of MGEX orbit and clock products during 2013–2015 for the emerging Galileo, BDS and QZSS constellations. From this study, satellite laser ranging (SLR) validations demonstrated an orbit accuracy of about 10 cm for the Galileo and BDS-2 inclined geosynchronous orbit (IGSO) and medium earth orbit (MEO) satellites, and about 20 cm for QZSS IGSO J01 [11]. Thereafter, in 2017, Montenbruck et al. assessed the orbit and clock products of the new constellations for the first half of 2016. Signal-in-space range error values of 3–7 cm were reported for Galileo and BDS-2 IGSO/MEO satellites [9].

Since January 2018, BDS has launched 28 satellites including 26 BDS-3 satellites and two BDS-2 satellites. At the end of 2018, BDS-3 began a global service. In parallel, following the launch of three QZSS satellites in June 2017, QZSS completed a four-satellite constellation and announced an official start to QZSS services in November 2018. Galileo also launched four new satellites into orbit in 2018, and all of them have since started signal transmission. Most of these new satellites have been included in the MGEX products, but their orbit and clock performances have not been fully investigated in current literature.

In addition, several critical issues have been continuously studied for BDS, Galileo and QZSS precise orbit determination (POD) and precise clock estimation (PCE), including antenna parameters [12], yaw-attitude models [13,14], solar radiation pressure (SRP) models [15,16,17], inter system bias (ISB) [14] and ambiguity fixing [18,19]. Moreover, the European GNSS Service Centre (EGSC) and Cabinet Office (CAO), Government of Japan, published the metadata of the Galileo and QZSS satellites, including the phase center offsets and variations (PCOs and PCVs), attitude laws and physical characteristics [20,21,22,23,24]. Satellite metadata are helpful for advanced POD and PCE processing [25,26]. With these refined models and new information, several updates have been implemented by MGEX ACs. As a result, further detailed assessments of the MGEX products are in need to identify potential improvements and deficiencies.

In this study, we present an assessment of the MGEX Galileo, BDS and QZSS products since 2018. After the introduction, an overview of the product availability as well as processing strategies is given. Then, quality of the precise orbit and clock products is analyzed in terms of orbit/clock comparisons and SLR validation. Thereafter, precise point positioning (PPP) performances with MGEX orbit and clock products are assessed. Finally, conclusions and discussions are given.

2. Precise MGEX Orbit and Clock Products

Currently, multi-GNSS precise orbit and clock products are generated by seven MGEX ACs, including CODE, GFZ, CNES/CLS, JAXA, SHAO, TUM and WHU. In addition, the European Space Operations Centre of the European Space Agency (ESA) has been providing multi-GNSS orbit and clock products since the end of 2017 (http://navigation-office.esa.int/products/gnss-products/) [27]. In this section, an overview on these products is given in terms of product availability and processing strategies.

Table 1. Overview of the precise MGEX orbit and clock products as of November 2019.

Institutions	ID	Constellations 9	Orbit	Clock	Remarks
CODE 1	COD0MGXFIN	GRECJ	5 min	30 s	GEOs excluded
ESA 2	ESM	GRECJ	5 min	30 s	BDS-3 included; GEOs excluded
GFZ 3	GFZ0MGXRAP	GRECJ	5 min	30 s
CNES/CLS 4	GRG0MGXFIN	GRE	15 min	30 s
JAXA 5	JAX0MGXFIN	GRJ	5 min	30 s	J02, J03 and J07 excluded
SHAO 6	SHA0MGXRAP	GREC	15 min	5 min
TUM 7	TUM0MGXRAP	ECJ	5 min	/
WHU 8	WUM0MGXFIN	GRECJ	15 min	30 s	BDS-3 included

¹ CODE stands for Center for Orbit Determination in Europe. ² ESA stands for European Space Agency. ³ GFZ stands for GeoForschungsZentrum. ⁴ CNES/CLS stands for Centre National d’Etudes Spatiales and Collecte Localisation Satellites. ⁵ JAXA stands for Japan Aerospace Exploration Agency. ⁶ SHAO stands for Shanghai Observatory. ⁷ TUM stands for Technische Universität München. ⁸ WHU stands for Wuhan University. ⁹ The letters G, R, E, C and J denote Global Positioning System (GPS), Global Navigation Satellite System (GLONASS), Galileo, BDS and QZSS, respectively.

summarizes the precise orbit and clock products provided by the seven MEGX ACs and the ESA. As of November 2019, all the seven MGEX ACs utilize the long filename convention for their products (http://mgex.igs.org/IGS_MGEX_Products.php). The initial three-character field of each filename specifies the analysis center and thus is used to denote the orbit/clock products from different ACs throughout the rest of this article (see “ID” in Table 1). Specifically, “COD”, “ESM”, “GFZ”, “GRG”, “JAX”, “SHA”, “TUM” and “WUM” are used to denote orbit/clock products from CODE, ESA, GFZ, CNES/CLS, JAXA, SHAO, TUM and WHU, respectively. For now, COD, ESM, GFZ, JAX and TUM orbit products are provided with 5-min sampling, while the GRG, SHA and WUM orbit products are provided at 15 min intervals. The complementary clock products are made available by seven ACs (excepting TUM). In 2017, only GRG and GFZ clock products were generated with 30 s sampling [9], but now 30 s clock products are generated by six ACs (excepting SHAO).

The availability of the MGEX products is depicted in Figure 1. Since day of year (DOY) 001 of 2018, COD and ESM orbit and clock products were provided without interruption. Nevertheless, TUM products were interrupted for more than five months from DOY300 2018 to DOY088 2019, and frequent interruptions were observed for SHA products as well. As for the new constellations, Galileo was included in COD, ESM, GFZ, GRG, SHA, TUM and WUM. BDS-2 was included in COD, ESM, GFZ, SHA and WUM, as well as TUM, since DOY217 2018. Please note that BDS-2 geostationary orbit (GEO) satellites were excluded in COD and ESM products. BDS-3 was only included in ESM and WUM since DOY340 2018 and DOY001 2019, respectively. SHA included BDS-3 products for four weeks from DOY196 2019 to DOY223 2019, but BDS-3 satellites are excluded in SHA products for now. QZSS satellites were included in COD, ESM, GFZ, JAX, TUM and WUM. The JAX products only included J01, while the other five products included all the QZSS IGSO satellites (i.e., J01, J02 and J03). Note that the JAX products once included J02 and J03 for a short time from DOY070 2019 to DOY083 2019. In addition, GFZ, TUM and WUM products also included the GEO satellite J07.

The main strategies of the MGEX POD and PCE processes are summarized in

Table 2. Comparison of the processing strategies of MGEX ACs (as of November 2019).

Item		CODE	ESA	GFZ	CNES/CLS	JAXA	SHAO	TUM	WHU
Differencing	Orbit	DD 1	UD	UD	UD	UD	UD	UD	UD
	Clock	UD 2	UD	UD	UD	UD	UD	UD	UD
Frequencies	Galileo	E1/E5a	E1/E5a	E1/E5a	E1/E5a	/	E1/E5a	E1/E5a	E1/E5a
	BDS	B1/B2	B1/B3	B1/B2	/	/	B1/B3	/	B1/B3
	QZSS	L1/L2	L1/L2	L1/L2	/	L1/L2	/	L1/L5	L1/L2
Elevation cutoff	Orbit	3°	10°	7°	8°	10°	7°	5°	7°
	Clock	5°	10°	7°	8°	10°	7°	/	7°
Arc length	/	3 d	1 d	3 d	30 h	7 d	3 d	3 d	1 d

¹ DD stands for double differenced. ² UD stands for un-differenced.

. Different software packages were employed by different ACs: Bernese by CODE and TUM; NAPEOS by ESA; EPOS by GFZ; POD GINS by CNES/CLS; MADOCA by JAXA; and PANDA by SHAO and WHU. CODE used double differenced (DD) observations for the POD process, while un-differenced (UD) observations are commonly used by other ACs. The POD arc durations ranged from 24 h to several days. CNES/CLS used an additional 3 h of data from the neighboring days (3 + 24 + 3 h) to get smaller discontinuities at the day boundaries [10]. The satellite PCO and PCV values of Galileo and QZSS published by the EGSC and CAO were updated into the Antenna Exchange (ANTEX) format file [28,29,30,31] and were commonly used by most ACs. As for BDS-2 satellites, MGEX conventional PCO values [32] were used by most ACs, while ESA and GFZ adopted PCO and PCV values estimated by Dilssner et al. [33].

The five-parameter Extended CODE Orbit Model (ECOM) [34] was adopted by GFZ for the new BDS, Galileo and QZSS satellites [3]. According to the IGS Technical Report 2018 [8], GFZ switched the SRP model from ECOM to ECOM2 [35] for IGS final, rapid and ultra-rapid processing in November 2018. However, it is uncertain whether ECOM2 was also applied by GFZ to generate MGEX products. CODE and CNES/CLS used the ECOM2 model. Please note that CODE adopted a so-called terminator reference frame since DOY198 2018 to describe the SRP accelerations for QZSS and BDS-2 satellites during orbit normal (ON) mode [36]. ESA adopted an a priori box-wing model to enhance the ECOM model, and estimated box-wing parameters were used for BDS and QZSS. In addition, one constant and two first-order harmonic terms were also estimated in along-track (http://navigation-office.esa.int/products/gnss-products/esm.acn). TUM also improved SRP models for Galileo and QZSS by using the a priori box-wing models, while ECOM is still used for BDS POD [4]. WHU has been applied for different SRP models for different satellites. The ECOM model has been used for the BDS-3, QZSS J02, J03 and J07 satellites. The cuboid model proposed by Zhao et al. [37] has been used for QZSS J01 POD. For BDS-2 IGSO and MEO satellites, a tightly constrained acceleration bias in the along-track direction has been added to the ECOM model [5]. For BDS-2 GEO satellites, an improved a priori model has been adopted to reduce the perturbation caused by the communication antenna [17]. As for JAXA, it has utilized a more complicated ECOM model with 13 parameters for QZSS POD [7].

3. Assessment of Precise MGEX Orbit and Clock Products

In this section, multi-GNSS products from 1 January 2018 to 1 November 2019 are collected. The consistency is assessed by orbit and clock comparisons among different ACs, and the orbit accuracy is further evaluated by the external SLR validations.

3.1. Orbit Comparisons

Satellite orbits from different ACs were compared at 15 min intervals for the along-track, cross-track and radial components. The comparison results will be discussed according to the individual Galileo, BDS and QZSS constellation.

3.1.1. Galileo

For Galileo, the ESM products were taken as the reference due to its smallest standard deviation (STD) values of SLR residuals among all ACs, which will be discussed in the next section. The daily RMS values of Galileo orbit differences with respect to ESM were averaged over all the available satellites for each AC, and the average RMS values in along-track, cross-track and radial directions are depicted in Figure 2. Outliers that were larger than the range of the vertical axes are represented as points on the horizontal axes. For most of the ACs, the RMS values of Galileo orbit differences were between 2 and 5 cm in the three components. COD showed the best agreement with regard to ESM, with the mean RMS values of 1.8, 1.4 and 1.6 cm in the along-track, cross-track and radial components, respectively. TUM showed a small mean RMS value of 1.8 cm in the radial component, while the along-track and cross-track components had the largest mean RMS values of 4.6 and 3.4 cm, respectively. It is worth mentioning that the orbit consistency between GRG and ESM was obviously improved since DOY280 2018, and the mean RMS values descended from 4.0, 3.5 and 3.0 cm to 2.4, 2.0 and 2.5 cm in the along-track, cross-track and radial components, respectively. This is because the GRG Galileo products have been computed via an undifferenced ambiguity fixing strategy since Global Positioning System (GPS) week 2022 (October 2018) [6].

Table 3. 3D RMS of Galileo orbit comparisons (units: cm).

	COD	ESM	GFZ	GRG	SHA	TUM	WUM
COD	-	2.8	4.2	4.0	6.6	6.2	4.9
ESM		-	3.5	3.9	6.1	6.1	4.7
GFZ			-	4.5	6.4	6.7	5.3
GRG				-	7.2	5.8	5.1
SHA					-	7.3	6.6
TUM						-	7.1
WUM							-

summarizes the averaged 3D RMS values of Galileo orbit comparison results among all ACs. Differences that were three times larger than the medians of the series were excluded from the statistics of comparison results. Note that only comparison results after DOY280 2018 were counted for GRG. As shown in the table, the consistency among COD, ESM, GFZ, GRG and WUM was at the 3–5 cm level.

3.1.2. BDS

Due to the different characteristics among the GEO, IGSO and MEO satellites, BDS orbit quality was assessed according to individual satellites rather than by averaging all satellites. In addition, BDS-2 and BDS-3 are discussed separately. WUM was taken as the reference for comparison because it was the only AC providing orbit and clock products for all the available BDS satellites. Although BDS-2 GEO satellites were included in TUM products since DOY281 2018, C02 and C03 products were still absent and C01 products were only available for about two months. In addition, COD and ESM products excluded BDS-2 GEOs from their POD solutions. As a result, only BDS-2 GEO orbits from GFZ, SHA and TUM products were compared with WUM.

Figure 3 shows the RMS of BDS-2 orbit differences with respect to WUM in the along-track, cross-track and radial components, and the average 3D RMS values of IGSO and MEO comparison results are summarized in

Table 4. 3D RMS of BDS-2 IGSO/MEO orbit comparisons (units: cm).

		COD	ESM	GFZ	SHA	TUM	WUM
COD	IGSO	-	21.4	20.6	17.9	24.0	18.0
	MEO	-	10.6	12.7	11.2	13.4	11.1
ESM	IGSO		-	11.8	15.8	23.1	17.4
	MEO		-	8.1	8.8	11.2	8.4
GFZ	IGSO			-	11.8	22.4	14.5
	MEO			-	7.9	12.5	8.3
SHA	IGSO				-	24.4	13.8
	MEO				-	12.0	8.7
TUM	IGSO					-	24.2
	MEO					-	12.1

. As can be seen, GEO satellites exhibited large orbit differences among ACs, which were up to a few meters in the along-track component and about 1.0 and 0.5 m in the cross-track and radial components, respectively. The orbits of IGSO/MEO satellites agreed well with each other among different ACs. For IGSOs, most of the RMS orbit values were within 0.1 m in the three components. Note that the newly launched IGSO satellite C16 showed larger RMS values than the other IGSO satellites. For MEOs, the COD, ESM, GFZ, SHA and TUM solutions showed an agreement of 6–8 cm in the along-track, 4–8 cm in the cross-track and 2–4 cm in the radial component with respect to WUM, respectively.

Figure 4 shows the RMS values of BDS-3 orbit differences between ESM and WUM in along-track, cross-track and radial components. One may notice that the orbit differences of satellites manufactured by the Shanghai Engineering Center for Microsatellites (SECM) were obviously larger than those of satellites manufactured by the China Academy of Space Technology (CAST) in the along-track and cross-track components. The mean RMS values of CAST satellites were 12.3, 9.2 and 4.1 cm, while the SECM satellites had larger mean RMS values of 25.0, 15.2 and 5.7 cm. This indicates that the orbit consistency of CAST satellites was better than that of SECM satellites. The different consistencies of CAST and SECM satellites may have been caused by the different PCO values of BDS-3 satellites used in ESM and WUM solutions. For the ESM solution, PCO values based on inhouse estimation were used, while for the WUM solution, conventional MGEX PCO values were used. The better orbit consistencies of CAST satellites may have been attributable to a better consistency between the two sets of PCO values used for CAST satellites.

3.1.3. QZSS

The qualities of QZSS orbit products were assessed by individual satellites. The orbits of three IGSO satellites (i.e., J01, J02 and J03) were compared with respect to ESM due to its best completeness among the three IGSO satellites. The orbits of the GEO satellite J07 are included by GFZ, WUM, and TUM. However, the J07 products of TUM were only available for 36 days, due to a long interruption. As a result, J07 orbits were compared only between GFZ and WUM.

Figure 5 shows the RMS values of QZSS orbit differences in along-track, cross-track and radial components. For J01, the periods with ON mode are marked as the grey area in the top left subfigure. It is clear that the RMS values increased significantly during and after switching the attitude mode to ON mode. As a result, only the RMS values in the yaw-steering (YS) period account for J01 in

Table 5. 3D RMS of QZSS IGSO orbit comparisons (units: cm).

		COD	ESM	GFZ	JAX	TUM	WUM
COD	J01	-	14.4	75.6	14.7	36.9	23.3
	J02	-	15.2	39.3	-	44.5	36.0
	J03	-	10.8	32.6	-	24.6	34.2
ESM	J01		-	78.6	15.6	31.4	19.2
	J02		-	41.0	-	41.3	31.6
	J03		-	38.3	-	29.8	36.1
GFZ	J01			-	78.0	88.4	81.2
	J02			-	-	51.4	28.6
	J03			-	-	55.3	29.9
JAX	J01				-	39.0	18.1
TUM	J01					-	33.9
	J02					-	40.4
	J03					-	43.5

. Furthermore, the J01 RMS values of GFZ showed a clear dependency on the beta angle (the sun elevation above the orbital plane) during YS mode in 2018. The RMS value gradually increased as the absolute value of the beta angle increased. Similar dependency can also be seen in the GFZ products for J02 (top right subfigure) and J03 (bottom left subfigure) in 2018, especially in the radial component. This may have been because GFZ was the only AC applying the ECOM model during the QZSS POD process.

Comparisons of the 3D RMS values among the three components of orbit are summarized in Table 5 for individual QZSS satellites. J01, COD, ESM, JAX and WUM showed good agreement with each other at the 10–20 cm level, while GFZ and TUM exhibited comparatively large orbit differences to other ACs. For J02 and J03, COD and ESM also showed the best agreements at the 10–20 cm level, while the orbit differences among GFZ, TUM and WUM were mostly between 30 and 50 cm. The mean RMS values of J07 orbit comparisons between GFZ and WUM were 256.4, 5.7 and 6.0 cm in the along-track, cross-track and radial components, respectively.

3.2. Clock Comparisons

Satellite clocks were compared at 5-min intervals. The comparison results were aligned with a reference satellite by satellite differencing to remove systematic biases. The STD values of the aligned comparison results were used as indicators of clock quality. Differences that were three times larger than the medians of the series were regarded as outliers and excluded from the comparison results.

Figure 6 shows the clock consistency of Galileo and QZSS with respect to ESM. The mean STD of clock differences is given in

Table 6. Mean STD of Galileo clock differences among different ACs (units: ns).

	COD	ESM	GFZ	GRG	SHA	WUM
COD	-	0.09	0.24	0.08	0.51	0.09
ESM		-	0.25	0.10	0.50	0.12
GFZ			-	0.25	0.54	0.24
GRG				-	0.51	0.10
SHA					-	0.51
WUM						-

and

Table 7. Mean STD of QZSS clock differences among different ACs (units: ns).

		COD	ESM	GFZ	JAX	WUM
COD	J01	-	0.29	0.88	0.19	0.29
	J02	-	0.25	0.60	-	0.62
	J03	-	0.25	0.54	-	0.77
ESM	J01		-	1.03	0.29	0.27
	J02		-	0.52	-	0.52
	J03		-	0.50	-	0.85
GFZ	J01			-	0.84	0.91
	J02			-	-	0.62
	J03			-	-	0.71
	J07			-	-	0.40
JAX	J01				-	0.25

. The points on the horizontal axes in Figure 6 represent outliers that were excluded from the statistics. For the Galileo satellites, the mean STD values of clock differences were around 0.1 ns among COD, ESM, GRG and WUM. However, SHA showed the largest mean STD value of clock differences, which was at the 0.5 ns level. Further, it is interesting to notice that the number of outliers of GFZ was obviously larger than those of other products during the periods of March and April, which may indicate that GFZ changed strategies when estimating Galileo clocks during the periods. As for J01, the clock differences were mostly at the 0.2- to 0.3-ns level, while GFZ showed large differences of about 1 ns when compared with other products. Again, the clock differences increased dramatically during and after switching the attitude mode to ON. As a result, clock differences during the ON period were excluded when calculating the statistics shown in Table 7. For J02 and J03, COD showed the best agreement with ESM at the 0.3-ns level. For J07, the STD value of clock differences between WUM and GFZ was 0.40 ns.

Likewise, the quality of BDS-2 and BDS-3 clock products was assessed using individual satellites rather than averaging all of them. Figure 7 shows the satellite-specific clock differences with respect to WUM, and the mean STD values are given in

Table 8. Mean STD of BDS-2 and BDS-3 clock differences among different ACs (units: ns).

		COD	ESM	GFZ	SHA	WUM
COD	BDS-2 IGSO	-	0.28	0.45	1.11	0.23
	BDS-2 MEO	-	0.17	0.46	1.12	0.17
ESM	BDS-2 IGSO		-	0.47	1.03	0.34
	BDS-2 MEO		-	0.42	1.10	0.16
	BDS-3 MEO		-	-	-	0.25
GFZ	BDS-2 GEO			-	1.21	1.06
	BDS-2 IGSO			-	1.28	0.48
	BDS-2 MEO			-	1.31	0.42
SHA	BDS-2 GEO				-	1.26
	BDS-2 IGSO				-	1.14
	BDS-2 MEO				-	1.12

according to different orbit types. For BDS-2 GEOs, the STD values of clock differences reached up to 1.06 ns for GFZ and 1.26 ns for SHA with respect to WUM. Obviously, smaller STD values between 0.2 and 0.5 ns could be observed among COD, ESM, GFZ and WUM for BDS-2 IGSO/MEO clocks. SHA showed the largest STDs of BDS-2 IGSO/MEO clocks at the 1.1-ns level. For BDS-3 MEO clocks, the STD values of clock differences between ESM and WUM were mostly at the 0.2- to 0.3-ns level. One may notice that the clock differences for C35, C36 and C37 were comparatively larger than for other BDS-3 satellites. This could be attributed to the reduced amount of observations of these three satellites. C28 also showed large clock differences, which may have been because the clock of this satellite was under-testing during the test period.

3.3. SLR Validation

The optical SLR technique measures the distance between an SLR observatory and a satellite equipped with laser retroreflector arrays (LRAs). This distance can also be computed using satellite orbits and observatory coordinates. The differences between the computed distance and the distance measured by SLR, namely SLR residuals, can be used as an external criterion to validate GNSS satellite orbits independently, especially radial components [38]. In this work, SLR observations were collected from the International Laser Ranging Service (ILRS) [39] to validate MGEX orbit products. Offsets of the laser retroreflector array [40] with respect to the satellite centers of mass were used in our SLR analysis, drawing from data released by the CAO [21,22,23,24], EGSC [20] and China Satellite Navigation Office (CSNO) [41]. The mean bias and STD values of the SLR residuals were introduced as quality indicators.

3.3.1. Galileo

The mean bias and STD values of SLR residuals for Galileo satellites are summarized in

Table 9. Mean bias and STD of Galileo satellite laser ranging (SLR) residuals (units: cm).

	COD	ESM	GFZ	GRG	SHA	TUM	WUM
IOV 1	0.3 ± 3.2	0.4 ± 2.8	0.5 ± 2.9	0.8 ± 3.8	−2.0 ± 5.0	−0.6 ± 3.8	0.1 ± 3.2
FOC 2	1.0 ± 3.1	0.7 ± 2.9	1.5 ± 3.1	2.4 ± 3.2	−1.5 ± 4.4	1.1 ± 3.5	1.7 ± 3.1
FOCe 3	2.0 ± 3.0	0.5 ± 2.7	1.4 ± 2.9	1.9 ± 5.9	−0.1 ± 6.7	1.2 ± 3.9	1.4 ± 5.8

¹ IOV stands for in-orbit validation satellites. ² FOC stands for full operational capability satellites. ³ FOCe stands for FOC satellites in elliptical orbit.

. SLR residuals with absolute values larger than 0.5 m were removed as outliers for Galileo. In this section, Galileo in-orbit validation (IOV), full operational capability (FOC) and FOC in elliptical orbit (FOCe) satellites are discussed separately. Most products showed close STD values at the 3- to 4-cm level for IOV and FOC satellites, while GRG, SHA and WUM showed larger STD values at the 6- to 7-cm level for FOCe satellites. Among all the products, ESM showed the smallest STD values for IOV, FOC and FOCe satellites with mean biases less than 1 cm. The GRG solution adopted antenna thrust modeling in April 2018 [6], which resulted in an obvious increase in mean bias values. From DOY105 2018, the satellite-averaged mean bias value of GRG products increased from 0.9 to 2.4 cm. What is more, it is interesting to note that SHA was the only AC whose mean biases were negative for all of the satellites. This could be explained as a result of the antenna thrust modeling not being considered in the SHA solution.

3.3.2. BDS

Table 10. Mean bias and STD of BDS SLR residuals (units: cm).

		COD	ESM	GFZ	SHA	TUM	WUM
C01	/	-	-	−28.1 ± 22.2 3 −12.6 ± 17.6 4	−33.8 ± 21.6	-	7.1 ± 9.8
C08	YS 5	−1.0 ± 5.2	8.5 ± 9.3	0.5 ± 5.1	−2.1 ± 5.5	0.9 ± 9.7	−1.1 ± 6.2
	ON 6	−44.9 ± 37.0 1 −6.2 ± 5.0 2	1.7 ± 11.1	−2.5 ± 13.7	−0.8 ± 24.3	−18.4 ± 20.5	−3.1 ± 12.4
C10	YS	2.6 ± 5.2	9.4 ± 7.7	2.3 ± 4.9	−0.4 ± 5.6	3.8 ± 9.8	1.1 ± 5.5
	ON	−13.9 ± 24.7 1 −6.7 ± 8.3 2	4.9 ± 6.9	0.6 ± 9.3	0.2 ± 18.9	5.1 ± 15.4	1.8 ± 14.3
C13	/	1.9 ± 6.5	9.1 ± 10.4	1.4 ± 6.0	−0.6 ± 6.3	1.1 ± 10.8	0.8 ± 5.6
C11	YS	1.2 ± 3.3	4.3 ± 4.3	0.7 ± 4.0	0.6 ± 3.8	1.7 ± 6.8	1.9 ± 3.1
	ON	−17.8 ± 34.3 1 −3.6 ± 14.0 2	1.4 ± 10.2	−2.4 ± 10.8	0.9 ± 16.8	−8.7 ± 23.8	−1 ± 10.6
C20	/	-	5.5 ± 5.6	-	-	-	4.6 ± 4.8
C21	/	-	6.9 ± 5.7	-	-	-	6.4 ± 4.7
C29	/	-	−4.6 ± 4.1	-	-	-	−2.1 ± 5.1
C30	/	-	−5.5 ± 4.3	-	-	-	−2.6 ± 5.1

¹ Statistics before DOY198 2018. ² Statistics after DOY198 2018. ³ Statistics before DOY330 2018. ⁴ Statistics after DOY330 2018. ⁵ YS stands for yaw-steering period. ⁶ ON stands for orbit normal period.

summarizes the mean bias as well as STD values of SLR residuals for different BDS orbit solutions. SLR residuals with absolute values larger than 1.0, 1.0 and 0.5 m were removed as outliers for the BDS GEOs, IGSO/MEO during the ON period and IGSO/MEO during the YS period, respectively.

Note that most BDS-2 IGSO/MEO satellites use two attitude modes (YS and ON). Thus, the SLR validation results during the YS and ON periods are discussed separately for most BDS-2 IGSO/MEO satellites except C13, since C13 is reported not to enter ON mode [13,42]. In the case of the YS period, COD, GFZ, SHA and WUM presented similar STD values of 5–6 cm for BDS-2 IGSO satellites (C08, C10 and C13) and 3–4 cm for the MEO satellite C11. TUM showed the largest STD values of about 10 cm for IGSOs and 7 cm for the MEO. While ESM showed a similar STD value of 4.3 cm for the MEO satellite as it did most other products, larger STD values of 8–10 cm were observed for IGSO satellites. ESM showed large mean biases of 8–9 cm for IGSOs and 4 cm for the MEO. In addition, the SLR residuals of ESM showed a clear dependency on the satellite–sun elongation angle for BDS-2 IGSO/MEO satellites (Figure 8) that was much more obvious than it was for any other product. In the case of the ON period, orbit accuracy was degraded significantly compared to that during the YS period for most products except the IGSOs of ESM. Although a tightly constrained acceleration bias in along-track direction was added for the BDS-2 MEO and IGSO satellites in WUM solution to compensate the deficiencies of the five-parameter ECOM in ON mode [5], the STD values in ON mode were still a factor of two to three worse compared to YS mode. Whereas COD (before DOY198 2018), GFZ, SHA and TUM were all reported to use ECOM for BDS-2 satellites, GFZ performed better during the ON period than the other three products, indicating that modifications could be applied for BDS-2 IGSO/MEO satellites during the ON period in GFZ solution. COD modified the SRP model for BDS satellites since DOY198 2018 and the modified SRP model took the ON mode into consideration [36]. With the refined SRP model, the STD values of COD during the ON period were significantly improved by 87%, 67% and 59% for C08, C10 and C11, respectively, diminishing from 37.0, 24.7 and 34.3 cm to 5.0, 8.3 and 14.0 cm.

The SLR residuals of BDS-2 GEO satellite C01 are shown in Figure 9 with respect to the beta angle and the satellite–sun elongation angle, respectively. As can be seen, the SLR residuals of GFZ and SHA had a clear dependency on the satellite–sun elongation angle. The mean biases of GFZ (before DOY330 2018) and SHA residuals were about –30 cm. Compared to GFZ and SHA, C01 orbits of WUM solutions had a less obvious dependency on the satellite–sun elongation angle and showed a smaller STD of 9.8 cm. The better performance of WUM is attributable to the improved BDS-2 GEO SRP model [17]. Additionally, it is interesting to note that an improvement can be observed for the C01 orbit of GFZ beginning DOY330 2018. The STD value improved from 22.2 to 17.6 cm, and the mean bias improved from –28.1 to –12.6 cm.

For BDS-3 MEO satellites, the STD values of SLR residuals were slightly larger than those of the BDS-2 MEO C11 during the YS period. What is more, the SLR residuals of BDS-3 satellites showed a clear linear dependency on the satellite–sun elongation angle for both ESM and WUM solutions, as shown in Figure 10. The slope was negative for CAST satellites C20 and C21, but positive for SECM satellites C29 and C30. The clear dependency on the satellite–sun elongation angle is attributable to the deficiency of SRP models. In order to improve BDS-3 orbit performances, we performed BDS-3 POD experiments using an a priori box-wing model together with ECOM. For the other POD strategies, refer to Li et al. [43]. The SLR residuals of the experiment results are also presented in Figure 10. Obviously, the dependency of the SLR residuals on the satellite–sun elongation angle could be decreased with an a priori box-wing model for SECM satellites, and smaller STD values of 3–4 cm were achieved.

3.3.3. QZSS

The SLR residuals of QZSS satellites are shown in Figure 11 with respect to the beta angle. The mean bias and STD values of the SLR residuals are summarized in

Table 11. Mean bias and STD of QZSS SLR residuals (units: cm).

	COD	ESM	GFZ	JAX	TUM	WUM
J01(ON)	16.1 ± 27.3 1 −1.7 ± 13.9 2	2.8 ± 28.2	−11.8 ± 24.7 3 −8.3 ± 13.6 4	−4.7 ± 12.4	−11.9 ± 39.7 5 −4.3 ± 25.0 6	6.8 ± 18.2
J01(YS)	−5.8 ± 7.1	0.5 ± 9.0	−40.1 ± 28.2 3 −22.3 ± 22.8 4	1.4 ± 4.6	8.1 ± 19.3 5 2.2 ± 7.5 6	7.1 ± 7.1
J02	−1.6 ± 7.6	4.5 ± 6.9	−14.4 ± 26.3 3 8.1 ± 9.1 4	-	6.9 ± 21 5 10.2 ± 10.4 6	22.4 ± 16.9
J03	7.9 ± 10.3	5.8 ± 6.1	−7.4 ± 20.7 3 4.1 ± 7.5 4	-	10.7 ± 34.4 5 4.9 ± 5.8 6	33.2 ± 18.5
J07	-	-	11.7 ± 16.9	-	-	22.9 ± 19.0

¹ Statistics before DOY198 2018. ² Statistics after DOY198 2018. ³ Statistics before DOY330 2018. ⁴ Statistics after DOY330 2018. ⁵ Statistics before DOY281 2018. ⁶ Statistics after DOY281 2018.

. SLR residuals with absolute values larger than 1 m are regarded as outliers and are excluded from the statistics.

For J01, JAX showed the smallest STD values of 12.4 and 4.6 cm during the ON and YS periods, respectively. Again, the orbit accuracy degraded significantly once the attitude mode switched to ON mode. With the SRP model designed for the ON period (applied since DOY198 (July), 2018), the J01 orbit of COD significantly improved by 49% in terms of STD values [36]. WUM also established an a priori SRP model for J01 POD to enhance the ECOM model, and an additional three parameters were used to compensate for the remaining modeling deficiencies during ON periods [37]. As a result, small STD values of 18.2 and 7.1 cm were achieved in the ON and YS modes, respectively. As for GFZ, the SLR residuals of J01 showed clear half-year periodical variations, and the absolute values increased as the absolute values of the beta angle increased. Further, GFZ showed the most outliers exceeding an absolute value of 1 m during YS periods, which may have been caused by insufficient SRP modeling.

For J02 and J03, ESM showed the smallest STD values of 6–7 cm, while WUM exhibited large STD values of 16–19 cm and large mean biases of 22–33 cm. As for GFZ, the STD values improved by about 65% beginning DOY330 2018. Considering that improvements of 45% and 19% were also observed for J01 during ON and YS periods, respectively, a refined SRP model (instead of ECOM) was assumed to be used for the QZSS IGSO POD in GFZ solution beginning DOY330 2018. TUM also improved SRP modelling for QZSS by using a priori box-wing models beginning DOY281 (July) 2018, and the STD values were reduced by 50% and 83% for J02 and J03, respectively. For the GEO satellite J07, GFZ and WUM showed large STD values of 16.9 and 19.0 cm.

3.4. PPP Performance

In order to further verify the quality of precise MGEX orbit and clock products, PPP solutions were performed in the static and kinematic modes. TUM products were not evaluated by PPP, since TUM does not provide precise clock products. Observations at stations CUSV, DARW, DGAR, IISC, KARR and KAT1 on 21 October 2019 were processed for PPP solutions in the GPS-only (G), BDS-only (C), Galileo-only (E) and GPS + QZSS (GJ) modes, respectively. The distribution of the six stations are shown in Figure 12. Since GPS was established for the longest time, PPP results in GPS-only mode were computed for comparison. In the case of BDS-only PPP, ionosphere-free linear combinations of B1 and B2 were processed with COD and GFZ products applied, while the combination of B1 and B3 were processed with ESM, SHA and WUM products applied. Considering that ESA-estimated PCO/PCV values were used in the generation of ESM and GFZ BDS products, those corrections for BDS-2 satellites were used in BDS-only PPP when ESM and GFZ products were applied. However, since the ESA-estimated PCO/PCV values for BDS-3 satellites were not publicly available, MGEX PCO/PCV values were used for BDS-3 satellites instead. Note that BDS-2 and BDS-3 are processed as one system in this section. The ISBs were estimated as random walk processes in multi-GNSS PPP with GFZ products [44], while being estimated as constants with other products. To assess the positioning performance, PPP solutions using IGS final orbit and clock products were performed first to estimate the coordinates of the six stations. The estimated coordinates were introduced as a reference to compute positioning errors in this section.

3.4.1. Static PPP Solutions

Table 12. RMS of positioning errors of static PPP with different products (unit: cm).

	GPS-Only			BDS-Only			Galileo-Only			GPS + QZSS
	E	N	U	E	N	U	E	N	U	E	N	U
COD	0.4	0.3	0.7	0.7	0.6	5.6	0.5	0.4	1.1	0.5	0.3	0.7
ESM	0.5	0.2	0.7	0.4	0.3	1.4	0.3	0.2	1.1	0.4	0.2	0.6
GFZ	0.5	0.3	0.7	1.5	0.4	1.8	0.3	0.4	1.1	0.4	0.3	0.7
GRG	0.4	0.4	0.6	/	/	/	0.3	0.4	1.1	/	/	/
JAX	0.5	0.3	0.7	/	/	/	/	/	/	0.6	0.3	0.7
SHA	1.3	0.4	0.6	1.2	0.4	3.5	0.3	0.4	1.4	/	/	/
WUM	0.5	0.8	0.7	0.6	0.9	0.9	0.5	0.7	1.1	0.7	0.8	0.7

summarizes the average RMS values of positioning errors for static PPP solutions with different products. Compared to the GPS-only PPP, a similar positioning accuracy of 3–5 mm in the horizontal directions was achieved in the case of Galileo-only PPP for most products used. The positioning errors in Galileo-only mode were about 11 mm in the up direction, which is larger than in GPS-only mode. As for BDS-only PPP, the positioning accuracy in the up direction depended on the number of BDS satellites included in the products. WUM provided products for all available BDS satellites, including BDS-3. Thus, the BDS-only PPP achieved the best positioning accuracy of 9 mm in the up direction with WUM products. A worse positioning accuracy of 14–18 mm resulted with ESM and GFZ products applied. This is because ESM included BDS-3, but excluded BDS-2 GEOs, while GFZ included BDS-2 GEOs, but excluded BDS-3. Without BDS-2 GEOs and BDS-3, BDS-only PPP with COD products presented the largest positioning error of 56 mm in the up direction. The PPP performance in GPS+QZSS mode was close to that in GPS-only PPP. This is reasonable, since QZSS is designed to complement GPS in urban and mountainous areas where the number of available GPS satellites is small.

3.4.2. Kinematic PPP Solutions

The positioning errors of kinematic PPP are shown for different products at the KAT1 station in Figure 13. The average positioning error RMS values of seven stations are summarized in

Table 13. RMS of positioning errors of kinematic PPP with different products (unit: cm).

	GPS-Only			BDS-Only			Galileo-Only			GPS + QZSS
	E	N	U	E	N	U	E	N	U	E	N	U
COD	1.6	1.1	4.4	5.2	3.4	25.6	1.8	1.4	7.7	1.4	1.0	4.2
ESM	1.6	1.1	4.4	3.0	1.4	7.4	1.7	1.3	7.4	1.4	1.0	4.1
GFZ	1.6	1.2	4.4	2.5	1.8	11.0	1.8	1.4	8.1	1.6	1.1	4.2
GRG	1.5	1.2	4.3	/	/	/	1.9	1.6	7.6	/	/	/
JAX	1.6	1.3	4.3	/	/	/	/	/	/	1.6	1.2	4.2
SHA	3.0	2.1	8.3	4.1	2.3	12.8	2.2	1.8	8.5	/	/	/
WUM	1.5	1.5	4.2	1.4	1.4	4.9	3.6	3.0	9.0	1.4	1.4	4.0

. With COD, ESM, GFZ and GRG products, the positioning accuracy of Galileo-only PPP was 17–19 mm in the north and 13–16 mm in the east, which is comparable to that of GPS-only PPP. However, in the up direction, the positioning errors of Galileo-only PPP were about twice as large as those of GPS-only PPP. For BDS-only PPP, WUM products performed the best compared with other products, which is reasonable since it included all the available BDS-2 and BDS-3 satellites. The BDS-only PPP with WUM products achieved accuracies of 14, 14 and 49 mm in the east, north and up directions, respectively. As for the GPS + QZSS PPP, the addition of QZSS satellites slightly improved the positioning accuracy when compared with GPS-only solutions.

4. Conclusions

In this study, the qualities of precise orbit and clock products for Galileo, BDS and QZSS were assessed. Galileo orbit/clock comparisons showed a consistency of 3–5 cm for orbit and 0.1–0.3 ns for clock among most ACs. The orbit differences of GRG with respect to ESM significantly improved after DOY280 2018 when GRG started to calculate Galileo orbits with an undifferenced ambiguity fixing strategy. BDS-2 orbit comparisons showed agreements of 12–18 cm for IGSOs and 8–9 cm for MEOs among ESM, GFZ, SHA and WUM. COD and TUM showed large orbit differences with other ACs in the along-track and cross-track directions. The orbit differences of BDS-3 satellites manufactured by SECM were obviously larger than those of BDS-3 satellites manufactured by CAST. The RMSs of the BDS clock differences were 0.2–0.5 ns for IGSOs and MEOs, and about 1 ns for GEOs for most ACs. As for QZSS, the orbit differences of J01, J02 and J03 were within 11–16 cm among COD, ESM and JAX. However, the orbit differences of GEO satellite J07 were able to reach several meters. The clock differences of QZSS satellites varied from 0.3 to 1.1 ns.

SLR validations showed an accuracy of 3–4 cm for Galileo. For BDS-2, the STD values of MEO satellites were mostly 3–4 cm, and the STD values of IGSO satellites were mostly 5–6 cm in COD, GFZ, WUM and SHA products. The orbit accuracy of the BDS-2 GEO satellite is significantly which improved in WUM solution by applying an improved SRP model for BDS GEO satellites. For BDS-3, the STD values of SLR residuals were mostly between 4 and 6 cm. For QZSS IGSO satellites, an accuracy of 5–10 cm was achieved in COD, ESM and JAX products, while the STD values of the QZSS GEO satellite J07 were around 20 cm.

GPS-only, BDS-only, Galileo-only and GPS+QZSS PPP were performed to further evaluate the qualities of multi-GNSS products provided by different ACs. For static PPP solutions, the positioning accuracies in Galileo-only mode were 3–5 mm in horizontal directions, and 11 mm in the up direction, with different products used. As for BDS-only PPP, the positioning errors in the up direction depended on the number of available BDS satellites included in the orbit and clock products. Since WUM products included all available BDS-2 and BDS-3 satellites, static PPP in BDS-only mode achieved the best positioning accuracy of 9 mm in the up direction using WUM products. As to kinematic PPP results, positioning errors in Galileo-only mode were 17–19 mm in the north, 13–16 mm in the east and 74–81 mm in the up direction, respectively. With WUM products used, positioning accuracies of 14, 14 and 49 mm were achieved in the east, north and up directions in BDS-only mode, respectively.