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

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.


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). Table 1 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).  3 GFZ stands for GeoForschungsZentrum. 4 CNES/CLS stands for Centre National d'Etudes Spatiales and Collecte Localisation Satellites. 5 JAXA stands for Japan Aerospace Exploration Agency. 6 SHAO stands for Shanghai Observatory. 7 TUM stands for Technische Universität München. 8 WHU stands for Wuhan University. 9 The letters G, R, E, C and J denote Global Positioning System (GPS), Global Navigation Satellite System (GLONASS), Galileo, BDS and QZSS, respectively.
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. 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 main strategies of the MGEX POD and PCE processes are summarized in Table 2. 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   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].

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.

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.

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 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.
Remote Sens. 2020, 11, x FOR PEER REVIEW 6 of 24    Figure 2. RMS values of Galileo orbit differences in along-track, cross-track and radial components with respect to ESM.

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. 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.     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.

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.  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. 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. Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 24 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  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.

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 Tables 6 and 7. 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 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.   Likewise, the quality of BDS-2 and BDS-3 clock products was assessed using individual satellites    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.

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.

Galileo
The mean bias and STD values of SLR residuals for Galileo satellites are summarized in Table 9. 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.  Table 10 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.     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.

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. SLR residuals with absolute values larger than 1 m are regarded as outliers and are excluded from the statistics. 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. SLR residuals with absolute values larger than 1 m are regarded as outliers and are excluded from the statistics.
Remote Sens. 2020, 11, x FOR PEER REVIEW 17 of 24 Figure 11. SLR residuals of QZSS satellites with respect to beta angle.

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   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.

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.
Remote Sens. 2020, 11, x FOR PEER REVIEW 18 of 24 to estimate the coordinates of the six stations. The estimated coordinates were introduced as a reference to compute positioning errors in this section.  Table 12 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  Table 12 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.

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. 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. Table 13. RMS of positioning errors of kinematic PPP with different products (unit: cm).

GPS-Only
BDS-Only Galileo-Only GPS + QZSS 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.

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.

Acknowledgments:
We are grateful to IGS-MGEX for providing multi-GNSS data and products. Thanks also go to ILRS for providing SLR observation data.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The