Precise Point Positioning with Almost Fully Deployed BDS-3, BDS-2, GPS, GLONASS, Galileo and QZSS Using Precise Products from Different Analysis Centers

: The space segment of all the ﬁve satellite systems capable of providing precise position services, namely BeiDou Navigation Satellite System (BDS) (including BDS-3 and BDS-2), are also provided. In addition, the consistency of the precise products from different analysis centers is characterized.


Introduction
Due to the simple data processing at a single station and the high-accuracy position determination, the precise point positioning (PPP) technology based on global navigation satellite system (GNSS) has been a research focus over the past two decades. In addition to the satellite geodesy, the PPP technology has also been widely applied for the navigation of marine vehicles [1] and land-borne vehicles [2][3][4][5]. Although great efforts have been made to improve the PPP performance in recent years, the long convergence time of PPP still limits its applications in time-critical fields. Multi-system integration is a promising approach to further enhance the PPP performance, especially in terms of convergence time, as both the increased number of available satellites and the improved geometry of satellites can be beneficial to the data processing of the satellite-based technology PPP. It precise products had relatively better positioning performance. However, besides the GPS and GLONASS satellites, only 15 Galileo satellites and 14 BDS-2 satellites were included in their study. Following Zhou et al. [26], the influence of ISB handling schemes on multi-GNSS PPP performance could be distinct when using the precise products from different analysis centers, and the random walk process and white noise process were recommended in multi-GNSS PPP processing.
Although great achievements have been made in terms of multi-system integrated PPP, few studies related to the data fusion covering all available BDS-3 satellites under a fully deployed constellation, as most of the existing studies employed the datasets before 2020 [19][20][21][22]. After a further integration with BDS-3 based on GPS, GLONASS, Galileo, BDS-2 and QZSS, the number of operational satellites has been up to approximately 136 at present. As a satellite-based positioning accuracy, the joint use of as many satellites as possible will definitely benefit the data processing of PPP. Currently, the precise products from several analysis centers, such as European Space Agency (ESA), GFZ and WHU, can support all the five satellite systems. The existed studies indicated that the quality of precise products could affect the multi-GNSS PPP performance [25]. However, few studies related to the consistency of precise products among different analysis centers for BDS-3 satellites. In addition, the corresponding conclusions for GPS, GLONASS, Galileo, BDS-2 and QZSS satellites in the existed studies (using the datasets before 2020) [23,24] may be obsolete, as the quality of precise satellite products with enhanced orbit and clock determination models for these satellites has been continuously improved in recent years.
According to the above description, it is necessary to assess the latest performance of five-system integrated PPP with almost fully deployed BDS-3, BDS-2, GPS, GLONASS, Galileo and QZSS (i.e., GRECJ-PPP). In this study, the performance of GRECJ-PPP under the current constellation status using daily observations in both static and kinematic modes in terms of position accuracy and convergence time is rigorously evaluated. In addition, the consistency of the precise products from ESA, GFZ and WHU that are able to support the five satellite systems and the effects of them on the data processing of GRECJ-PPP are also investigated. For completeness, the real-time position solutions derived from four-system integrated PPP with BDS (BDS-3/BDS-2), GPS, GLONASS and Galileo (i.e., GREC-PPP) using the real-time stream from Centre National d'Etudes Spatiales (CNES) are also provided (noted that QZSS satellites cannot be supported by CNES real-time products). For comparison, the results of single-system PPP and dual-system combined PPP are also presented. The paper starts with the GRECJ-PPP model. Subsequently, the results are presented, including the description of precise products for the five satellite systems, the availability of five-system combination, the consistency of precise products from different analysis centers, and the performance of multi-system combined PPP. Then, the discussion is conducted. Finally, we summarize the main points and the conclusions.

Methods
In this section, the five-system integrated PPP model with BDS, GPS, GLONASS, Galileo and QZSS is described. In this contribution, we adopt the functional model based on ionospheric-free (IF) combined observables, which is usually employed by the analysis centers for precise orbit determination (POD) and precise clock estimation (PCE). There are uncalibrated code delays (UCDs) at both satellite and receiver ends in the code observations, and the satellite-and receiver-specific IF-based UCDs will be grouped into the satellite clock parameters and receiver clock parameters due to the strong correlation. The IF combined code and phase observations (i.e., P IF and L IF ) from GPS, GLONASS, Galileo where the superscript s and the subscript r denote a GNSS/RNSS satellite and a receiver, respectively, the superscript Q refers to G, R, E, C or J which represents the satellite system GPS, GLONASS, Galileo, BDS or QZSS, the subscript IF refers to the dual-frequency IF combination which is formed with L1 and L2 signals for GPS and QZSS satellites, G1 and G2 signals for GLONASS satellites, E1 and E5a signals for Galileo satellites, and B1 and B3 signals for BDS (BDS-2/BDS-3) satellites in this study, respectively, ρ denotes the geometric range between the receiver and the satellite, cdt r and cdt s denote the receiver and satellite clock offsets, respectively, T is the slant tropospheric delay, N IF is the IF-based phase ambiguity, ε(P IF ) and ε(L IF ) denote the code and phase measurement noises including the multipath errors, respectively, b IF,r and b s IF denote the receiver and satellite UCDs, respectively, and B IF,r and B s IF denote the receiver and the satellite uncalibrated phase delays (UPDs), respectively. Following Equation (2), the estimable phase ambiguity N s,Q IF,r absorbs the IF-based UCDs and UPDs of the receiver and the satellite (i.e., b IF,r , b s IF , B IF,r and B s IF ). When the IF combination used for PPP processing is the same as that adopted by PCE, both the satellite clock offsets and satellite UCDs (i.e., cdt s,Q ) in Equation (1) can be cancelled out. As for the estimable receiver clock cdt Q r , it is dependent on the satellite system due to the grouped receiver UCDs, and thus the setting of only a receiver clock parameter in the five-system integrated PPP model is obviously unreasonable. To solve this issue, we can estimate a respective receive clock parameter for each satellite system. Alternatively, by taking the receiver clock parameter of a selected satellite system as the reference, an ISB (i.e., inter-system bias) parameter can be introduced into the observation equations for other satellite systems. As GLONASS employs a frequency division multiple access (FDMA) technique, the receiver UCDs differ among different satellite signals. If we add many parameters to consider the inter-channel biases in the code observation equations, the strength of the functional model will be relatively weak. In this study, the GLONASS code observations are down weighted to mitigate the negative effects of inter-channel biases. After fixing the satellite orbits with precise satellite orbit products, correcting the satellite clock offsets with precise satellite clock products, and correcting the tropospheric dry delay with a priori model, the linearized observation model of five-system integrated PPP can be described as: with: where p IF and l IF denote the observed-minus-computed (OMC) IF-based code and phase observables, respectively, µ denotes the unit vector in the line-of-sight direction, x denotes the three-dimensional (3D) receiver coordinates, cdt G r denotes the estimated receiver clock offsets grouped with the receiver UCDs of GPS, m is the wet mapping function, Z is the tropospheric zenith wet delay (ZWD), and ISB G,Q r denotes the ISB between the satellite systems Q(Q = G) and G.
The estimated parameters include receiver coordinates, receiver clocks, ISB, ZWD and phase ambiguities, that is: where O is the vector of estimates.
In addition to the rigorous functional model, a proper stochastic model also plays a very important role in the data processing of PPP. It is assumed that the observations of different types (i.e., code and phase observations), from different satellites, or on different frequencies are independent. The variance of observations from a GNSS/RNSS satellite on a single frequency can be computed as follows: where ele is the satellite elevation angles, and the two terms k 1 and k 2 are constants, which are both set to 3 mm for phase observations. As for code observations, the ratio between the standard deviation (STD) of them and the phase observations is chosen as 100:1 for GPS, Galileo and QZSS satellites, 150:1 for GLONASS satellites, and 500:1 for BDS satellites [26].
It is important to notice that the observations of BDS geostationary orbit (GEO) satellites are further down weighted by a factor of ten times compared with the observations from other BDS satellites. For further analysis, we compare the positioning performance between the BDS-only static PPP (long-term performance of PPP) with a ratio of 100:1 and 500:1 based on the GFZ final precise products. The convergence time of BDS-only static PPP with a ratio of 100:1 can be increased by several minutes over the case with a ratio of 500:1, indicating that the code observations of BDS are poor, and thus the latter one is employed. According to the support of code observation types by the ground tracking stations, the code observations "C1C" and "C2W" on L1 and L2 frequencies for GPS, "C1P" and "C2P" on G1 and G2 frequencies for GLONASS, "C1C" and "C5Q" (or "C1X" and "C5X") on E1 and E5a frequencies for Galileo, "C2I" and "C6I" on B1 and B3 frequencies for BDS, and "C1C" and "C2L" (or "C1C" and "C2X") on L1 and L2 frequencies for QZSS are adopted in this study, respectively. The above code observation types are defined in Receiver Independent Exchange (RINEX) version 3.04. It should be noted that the setting for the variance of code observations shown in Equation (6) for the five satellite systems is rough, and a precise setting should be further investigated. El-Mowafy [27] proposed a method about the precise setting for the variance of code observations, which may be used to further refine the stochastic model. As for the dynamic model of the estimated parameters in the Kalman filter, the static receiver positions and the phase ambiguities are modeled as constants, while the kinematic receiver positions and the receiver clock offsets are modeled as white noise processes (10 4 m 2 ). Regarding the ISB and ZWD, they are estimated as random-walk processes. The spectral density value is set to 10 -6 and 10 -8 m 2 /s for the ISB and ZWD, respectively.

Precise Products for BDS, GPS, GLONASS, Galileo and QZSS
As GNSS moves towards multi-system integration, the MGEX project [28], aiming at the continuous tracking, collection and analysis of all GNSS (GPS, GLONASS, Galileo, BDS) and RNSS (QZSS) satellites, has been continuously expanded and improved. Currently, there are seven analysis centers capable of providing multi-system precise satellite orbit and clock products, including CNES/Collecte Localisation Satellites (CLS)/Groupe de Recherche de Géodésie Spatiale (GRGS) [29], Center for Orbit Determination in Europe (CODE) [30], GFZ [31], Information and Analysis Center (IAC) [32], Japan Aerospace Exploration Agency (JAXA) [33], Shanghai Observatory (SHAO) [34] and WHU [35]. It should be noted that the analysis center Technische Universität München (TUM) has no longer provided the precise products for MGEX project since 10 November 2019. The analysis center IAC started to synchronize the multi-system precise products to the Crustal Dynamics Data Information System (CDDIS) data server on 12 September 2020, and these products were uploaded to the internal File Transfer Protocol (FTP) (ftp.glonass-iac.ru) before this. The multi-system precise products provided by these analysis centers are all stored at the CDDIS MGEX product archive (https://cddis.nasa.gov/archive/gnss/products/ mgex/, accessed on 25 September 2021) as well as the mirror sites hosted by Institute Geographique National (IGN) (ftp://igs.ign.fr/pub/igs/products/mgex/, accessed on 25 September 2021) and École Nationale des Sciences Géographiques (ENSG) (ftp://igs. ensg.eu/pub/igs/products/mgex/, accessed on 25 September 2021). Besides, the European Space Operations Centre of the ESA [36] can also offer the multi-system precise products (http://navigation-office.esa.int/products/gnss-products/, accessed on 25 September 2021).
Because of the rapid growth in demand for real-time precise applications, several analysis centers start to transmit the real-time precise corrections for satellite orbits and clocks based on the State Space Representation (SSR) messages defined by Radio Technical Commission for Maritime Services (RTCM) through Networked Transport of RTCM via Internet Protocol (NTRIP). The analysis center CNES is the first institution to provide real-time precise products, and the transmitted real-time stream CLK93 can support all the four GNSSs. Although the CNES real-time products do not incorporate QZSS, they are used to assess the performance of real-time PPP under the current GNSS constellation status in this study. For simplicity, we employ the retrieved precise satellite orbit and clock products based on the real-time stream CLK93 and the broadcast ephemeris, which can be downloaded from the PPP-Wizard (Precise Point Positioning With Integer and Zero-difference Ambiguity Resolution Demonstrator) website (http://www.ppp-wizard. net/products/REAL_TIME/, accessed on 25 September 2021). Table 1 conducts an overview of the precise satellite orbit and clock products from multi-system analysis centers as of January 2021. Following the named rules of RINEX version three files, "COD", "ESA", "GBM", "GRG", "IAC", "JAX", "SHA" and "WUM" are used by the analysis centers to identify the CODE, ESA, GFZ, CNES/CLS/GRGS, IAC, JAXA, SHAO and WHU post-processed precise products, respectively. In this study, "CNT" denotes the retrieved CNES real-time precise products. The abbreviations GLO, GAL and QZS refer to GLONASS, Galileo and QZSS, respectively. It is seen that, in addition to the legacy GPS and GLONASS, the precise products from all the analysis centers can also support the emerging Galileo and BDS-2, except for the "GRG" and "JAX" products. As for QZSS satellites, they are absent in the "GRG", "SHA" and "CNT" products, while the "JAX" products only include QZSS J01 satellite and do not support the other three QZSS satellites (J02, J03 and J07). There are significant differences in the support of GEO satellites for the various precise products. Both the "COD" and "ESA" precise products exclude BDS-2 and QZSS GEO satellites probably due to their worse orbit determination. Besides the "GBM" products, the analysis center GFZ also releases another set of precise products, which are identified by "GFZ". The "GFZ" products do not support BDS-3 satellites, and the B1/B2 dual-frequency IF combination is used to generate the BDS-2 precise satellite orbit and clock products. Based on the compatible B1/B3 dual-frequency signals, the "GBM" products include both the BDS-2 and BDS-3 satellites. On 14 June 2020, the "GFZ" products were replaced by the "GBM" products to be uploaded to the MGEX data centers. The "GBM" and "WUM" products can also support several BDS-3 GEO satellites, making them by far the most comprehensive products in support of all the available satellite constellations. Currently, only the "COD", "GRG" and "JAX" products do not cover the BDS-3 satellites.
To sum up, only the "ESA", "GBM", "IAC" and "WUM" products can cover all the satellite systems at present, and the number of the supported GNSS/RNSS satellites in the four precise products can be up to 114, 123, 119 and 124 as of January 2021, respectively. However, the "IAC" products during the analysis period of this study are absent, and the continuous supply of "IAC" products usually cannot be guaranteed. Table 1. Overview of the precise satellite orbit and clock products from multi-system analysis centers as of January 2021. It should be noted that this paper attempts to assess the performance of PPP with all available satellites (GPS+GLO+GAL+BDS-2+BDS-3+QZS) using reliable precise products under the current constellation status (in addition to the consistency analysis of precise products), and thus the "COD", "GRG", "IAC", "JAX" and "SHA" products are excluded from our analysis. In this study, the "ESA", "GBM" and "WUM" products are adopted to carry out the post-processed five-system PPP processing, and the consistency among them is also characterized. For completeness, the "CNT" products are employed to analyze the real-time four-system (GPS+GLO+GAL+BDS-2+BDS-3) PPP, and the consistency between real-time and post-processed products is also investigated. The software used for POD and PCE is different for different analysis centers, such as NAPEOS by ESA, EPOS by GFZ, and PANDA by WHU, but all employs the undifferenced dual-frequency IF combined code and phase observables. The frequency selection is L1/L2 for GPS and QZSS, G1/G2 for GLONASS, E1/E5a for Galileo, and B1/B3 for BDS-2 and BDS-3, respectively. With the covering of all the satellite systems by the "IGS14.atx" file provided by international GNSS service (IGS) [37], all the analysis centers gradually adopt it to correct the phase center offset (PCO) and variation (PCV), except for ESA. The analysis center ESA still employs the estimated PCO and PCV corrections by itself. Currently, more detailed information about the "CNT" products is not clear.

Availability of Five-System Combination
With the number of the five-system satellites further increases, the number of visible satellites and the Position Dilution of Precision (PDOP) over the ground tracking station will also have a better performance. We evaluate the global availability of five-system (GPS+GLO+GAL+BDS-2+BDS-3+QZS) combination in terms of satellite numbers and PDOP values in this section. As the final precise orbit product from GFZ is one of the most comprehensive products in support of all the available satellite constellations, it is used to calculate the satellite positions. The sampling interval is 5 min, and the cut-off elevation angle is set to 10 • . The whole world is divided into 72 × 72 grids every 2.5 • in latitude and 5 • in longitude, respectively. The center point of each grid is regarded as a virtual ground tracking station with the geodetic height set to 25 m [38].
The global maps for the maximum, minimum and average number of visible satellites for five-system combination on day of year (DOY) 122 of 2020 are illustrated in Figure 1. The distribution of satellite numbers has two concentrated areas. One is located in the Asia-Pacific region, where more satellites can be observed since it is also the service regions of BDS-2 and QZSS. The maximum value for the maximum satellite number and the minimum satellite number within this range can be 57 and 48, respectively. The other one is located in the Americas, and the minimum value for the maximum satellite number and the minimum satellite number over this region is 33 and 21, respectively. From a global perspective, we can observe 28.6-51.5 satellites on average. This means that the continuous positioning services for whole 24 h will be hardly restricted by insufficient satellites. Figure 2 depicts the global maps for the PDOP values. Corresponding to the satellite visibility distribution, the PDOP values of five-system integration also demonstrate the same regional difference. In the Asia-Pacific region, the minimum value for the maximum and minimum PDOP is 0.8 and 0.6, respectively, while the maximum value for the maximum and minimum PDOP in the Americas is 1.5 and 0.9, respectively. The average PDOP values change from 0.7 to 1.0 on a global scale.

Consistency of Precise Products from Different Analysis Centers
In this section, the consistency (characterized by satellite systems) among the ESA, GFZ and WHU post-processed precise products as well as the CNES real-time precise products in terms of the satellite orbits in three dimensions and the satellite clocks is analyzed. For completeness, the Signal-in-Space Ranging Error (SISRE), which can comprehensively reflect the effects of satellite orbits and clocks, is also introduced. The calculation of SISRE was well described in Montenbruck et al. [39]. The analysis period here covers a time span of about a month (31 days) from DOY 122 to 152 in 2020. The 31-day data are first divided into many sub-sets (one day from GPS time 00:00:00 to 23:59:30 or 23:59:55 for each subset). The epoch-wise orbit and clock difference among different analysis centers are then obtained using the single-day data. Finally, all the single-day sub-sets of epoch-wise orbit and clock difference are used together to compute the statistics of product inconsistency. In addition, the epoch-wise consistency results of precise products between GFZ and WHU over a day are also detailed as an illustration, since they are the most comprehensive products in support of all the available satellite constellations.

GPS
The epoch-wise differences of precise satellite orbit and clock corrections among different analysis centers (ACs) for GPS satellites are first derived, and then we compute the root mean square (RMS) statistics of the epoch-wise orbit and clock differences over all the available satellites and days. The epoch-wise orbit and clock differences between GFZ and WHU on DOY 122 of 2020 are plotted in Figure 3. Different satellites are represented by different colors. The orbital differences for most of the satellites are smaller than 5 cm, except for three satellites (G18/G24/G30), which can reach approximately 0.3 m. The epoch-wise differences in clock offsets are usually less than 0.5 ns, and several satellites show very small clock differences (G01/G04/G06/G09/G25/G26/G27/G32). The obtained results about RMS statistics are provided in Table 2. To comprehensively reflect the precise product differences among different ACs, Table 2 also lists the SISRE statistics. The radial orbit differences are usually at a level of approximately 2 cm, except for those between GFZ and CNES, and WHU and CNES, which are 2.9 and 2.5 cm, respectively. As for the clock differences, the numerical values can be as small as 0.06 ns between ESA and WHU, while the corresponding differences range from 0.08 to 0.13 ns for the other five cases. The along-track and cross-track orbit differences show comparable performance for the three post-processed precise products, and are only 1.6 cm between ESA and WHU and about 2.5 cm between GFZ and ESA or WHU. The orbit differences are increased to 3.2-4.1 and 2.4-3.5 cm in the along-track and cross-track directions between CNES and the other three ACs, respectively. In terms of the SISRE, the consistency of precise satellite orbit and clock products can be up to 1.2 cm between ESA and GFZ or WHU, whereas the corresponding consistency between GFZ and WHU degrades to 3.5 cm. The SISRE differences between real-time precise products and three post-processed precise products vary in a range of 2.4 to 2.7 cm. Overall, for GPS satellites, the consistency of precise products shows best performance between ESA and WHU. In contrast, GFZ and WHU precise products exhibit a worse consistency, and it is also the case for the real-time and post-processed precise products.

GLONASS
The epoch-wise orbit and clock differences between GFZ and WHU precise products for GLONASS satellites on DOY 122 of 2020 are illustrated in Figure 4. It can be seen that the orbit and clock differences of GLONASS are larger than those of GPS, and most of the orbit and clock differences are smaller than 8 cm and 0.7 ns, respectively. The epoch-wise orbit differences of R20 satellite (marked in pink) show an anomalistic behavior, especially for the along-track and cross-track directions. The RMS values of the epoch-wise orbit differences between GFZ and WHU for R20 satellite over a day are 4.5, 11.8 and 13.1 cm in the radial, along-track and cross-track directions, respectively. For further analysis, we compute the corresponding RMS statistics of the orbit differences between ESA and GFZ, and ESA and WHU, which are 5.3, 9.3 and 13.1 cm, and 4.2, 6.0 and 3.9 cm in the three directions, respectively. Thus, the anomalistic behavior of the epoch-wise orbit differences of R20 satellite shown in Figure 4 may be attributed to the worse orbit determination by GFZ for this satellite on DOY 122 of 2020.  Table 3 lists the corresponding statistical results of GLONASS. The radial orbit differences among different ACs fall between 2.3 and 3.5 cm for GLONASS satellites. As to the orbital differences among different ACs in the other two directions, the numerical value is at a level of approximately 4 cm, except for the cases involving CNES real-time precise products. For the inconsistency of precise satellite orbit products between CNES and other ACs, the RMS differences range from 8.9 to 9.5 cm and from 5.7 to 6.6 cm in the along-track and cross-track directions, respectively. The RMS statistics of epoch-wise differences of precise satellite clock corrections over all the available epochs and satellites between ESA and GFZ, and ESA and WHU are only 0.10 ns, while the corresponding statistics vary within a range of 0.18-0.20 ns for the other four cases. Similar to the GPS satellites, the ESA and WHU precise products achieve the best consistency for GLONASS satellites in terms of the SISRE statistics, followed by the consistency between ESA and GFZ precise products, which are 2.7 and 2.8 cm, respectively. The SISRE differences are increased to 4.9 cm between GFZ and WHU precise products, and the corresponding statistics further degrade to 5.2-5.7 cm between the CNES real-time precise products and three post-processed precise products. The SISRE inconsistency of precise satellite products among different ACs of GLONASS satellites is about twice larger than that of GPS satellites, except for the results between GFZ and WHU. For further analysis, with the use of the data over 31 days, we compute the consistency statistics of precise products between GFZ and WHU for each GLONASS satellite. The single-satellite SISRE differences range from 3.5 to 10.1 cm, and the SISRE statistics are 5.2 cm for the GLONASS-K satellite R09, which are close to the results taking all GLONASS satellites into consideration (see Table 3). Thus, it seems that the consistency of precise products for the GLONASS-K satellite (R09) is comparable to that for GLONASS-M satellites.

Galileo
The epoch-wise orbit and clock differences between GFZ and WHU precise products for Galileo on DOY 122 of 2020 are provided in Figure 5. Similar to GPS, the epoch-wise orbit differences of Galileo are usually less than 5 cm, but the clock differences are smaller than 0.3 ns for most of the time, which is better than that of GPS and GLONASS. The corresponding statistical results of Galileo are shown in Table 4. Similar to GPS and GLONASS satellites, the radial orbit difference has the smallest RMS values among the three orbital directions, while the along-track and cross-track orbit differences seem to be at a similar level for the inconsistency results only involving post-processed products, and the along-track orbit difference achieves the largest RMSs for the corresponding results between the CNES real-time products and three post-processed products. Regarding Galileo satellites, the consistency of precise satellite products between ESA and GFZ is at the highest level, and the statistics are 1.9, 2.4 and 2.5 cm, 0.06 ns and 1.7 cm in terms of radial, along-track and cross-track orbit difference, clock difference and SISRE difference, respectively. The consistency between WHU and ESA, and WHU and GFZ is also at a satisfied level, with a scattering of 2.3, 3.1 and 3.7 cm, 0.06 ns and 1.9 cm, and 2.2, 3.9 and 4.6 cm, 0.08 ns and 2.3 cm for the corresponding five different statistical terms for the two cases, respectively. As for the consistency performance between the real-time and post-processed precise satellite products, the orbital difference is slightly increased by 3.4 cm at most compared with the corresponding inconsistency results among the postprocessed products, but the clock difference is enlarged by a factor of about three times, resulting in the larger SISRE difference with a value of 5.3-5.5 cm. Compared with GPS satellites, the SISRE difference of Galileo satellites is increased by approximately 0.5 cm for the inconsistency results only involving post-processed products except for those between GFZ and WHU, and is enlarged by about two times for the inconsistency between the real-time and post-processed precise satellite products.  3.3.4. BDS Figure 6 depicts the epoch-wise orbit and clock differences between GFZ and WHU precise products for BDS-2 GEO and IGSO satellites on DOY 122 of 2020. Most orbit differences of BDS-2 IGSO satellites range from −0.3 to 0.3 m for all the three orbital components. As for BDS-2 GEO satellites, the radial orbit differences also vary within 0.3 m, but the orbit differences can exceed 1 and 2 m in the cross-track and along-track components for several satellites, respectively. Except for C06 satellite, the clock differences of BDS-2 GEO and IGSO satellites have a varying range of 1.0 ns. Figure 7 depicts the corresponding results for BDS-2 MEO and BDS-3 MEO satellites. The fluctuation range of orbit differences for BDS MEO satellites is reduced to approximately 8, 15 and 15 cm in the radial, cross-track and along-track components, respectively, and the corresponding clock differences nearly always keep within 0.5 ns. The uneven distribution of tracking stations for BDS GEO and IGSO satellites will result in a weak geometric configuration between stations and satellites, which can obviously affect the accuracy of precise orbit and clock products.  The statistical inconsistency results of precise satellite products for BDS-2 GEO and IGSO satellites are exhibited in Table 5. It is noted that the precise products from ESA do not support the BDS-2 GEO satellites. The orbital difference in the along-track and cross-track directions of BDS-2 GEO satellites can be larger than 1 m, especially for the along-track orbit difference with a value of 324.1 and 355.3 cm between CNES and GFZ, and CNES and WHU, respectively. The radial orbit difference of BDS-2 GEO satellites between WHU and GFZ precise products decreases to 15.2 cm, and the RMSs of clock difference are 0.39 ns, both of which seem to be at an acceptable level, but the corresponding statistics are increased by a factor of about three and two times for the inconsistency results between the real-time and post-processed products, respectively. Regarding the SISRE statistics, the consistency of precise products of BDS-2 GEO satellites between WHU and GFZ is 22.1 cm, while the corresponding consistency between the real-time and post-processed products decreases to 54.4-55.2 cm. Compared with the BDS-2 GEO satellites, the consistency results of BDS-2 IGSO satellites are significantly improved, especially for the along-track and cross-track orbit differences. Different from BDS-2 GEO satellites, the inconsistency results between the real-time and post-processed precise products are comparable to those only involving the post-processed products in terms of the radial orbit, clock and SISRE statistics, but the along-track and cross-track orbit differences of the former ones are still larger than those of the later ones by several centimeters. The differences of radial orbits, along-track orbits, cross-track orbits, clocks and SISREs among different ACs for BDS-2 IGSO satellites vary within a range of 9.3-13.4 cm, 7.1-15.5 cm, 7.3-14.6 cm, 0.18-0.40 ns and 6.8-10.1 cm, respectively. Table 6 illustrates the corresponding statistical results for BDS-2 MEO and BDS-3 MEO satellites. It should be noted that the precise products were absent for BDS-3 GEO and IGSO satellites during the analysis period. The consistency of precise satellite products among different ACs for the BDS-2 and BDS-3 MEO satellites is further improved in comparison to the BDS-2 GEO and IGSO satellites. The precise products of BDS-3 MEO satellites achieve a slightly worse consistency results among different ACs than the BDS-2 MEO satellites. For both BDS-2 and BDS-3 MEO satellites, the WHU and GFZ precise satellite products have the best consistency, followed by that between ESA and GFZ precise products. Except for the radial orbit differences of BDS-2 MEO satellites, the inconsistency between real-time and post-processed products is larger than that among post-processed products. As to the three orbital components, the inconsistency among the post-processed products with a varying range of 2.2-5.6 and 2.6-6.5 cm for the BDS-2 and BDS-3 MEO satellites seems to be at a similar level, while the corresponding fluctuation range of the RMS difference between the real-time and post-processed products increases to 3.4-3.9, 9.8-10.9 and 7.1-8.2 cm, and 5.0-6.5, 15.3-16.4 and 10.4-11.2 cm for the radial, along-track and cross-track orbital components for the BDS-2 and BDS-3 MEO satellites, respectively. The clock differences of BDS-2 and BDS-3 MEO satellites among the post-processed products are comparable to those of GPS satellites, which range from 0.05 to 0.08 ns for BDS-2 MEO satellites, and from 0.11 to 0.14 ns for BDS-3 MEO satellites, respectively. The RMS values of clock differences between the post-processed and real-time products are enlarged by a factor of two to three times, in comparison to those among the post-processed products. The SISRE difference among the post-processed products of the BDS-2 MEO satellites is comparable to that of BDS-3 MEO satellites, and ranges from 2.6-4.1 and 3.2-3.7 cm for the two satellite constellations, respectively. As for the SISRE difference between the real-time and postprocessed products, it changes from 4.8 to 5.4 cm for BDS-2 MEO satellites and from 6.8 to 7.4 cm for BDS-3 MEO satellites, respectively. Compared with GPS satellites (MEO satellites), except for the consistency results between GFZ and WHU precise products, the SISRE difference among different ACs of BDS-2 and BDS-3 MEO satellites is expanded by a factor of two to three times. It should be noticed that there are two manufacturers for the BDS-3 satellites. One is the China Academy of Space Technology (CAST), and the other one is the Shanghai Engineering Center for Microsatellites (SECM). The CAST and SECM satellites (used for the analysis of product consistency) are equipped with two different types of atomic clocks, namely rubidium clocks and hydrogen masers, which has an impact on the quality of onboard clock stability, and thus also affects the SISRE statistics and the quality of position determination. In addition, the two types of BDS-3 satellites also have different construction and dimensions. For further analysis, the consistency results shown in Table 6 between GFZ and WHU for BDS-3 MEO satellites are re-computed by dividing them into two groups of satellites, namely CAST and SECM satellites. The results indicate that the statistics are 2.6, 5.4 and 5.1 cm, 0.11 ns and 2.8 cm for CAST satellites in terms of radial, along-track and cross-track orbit difference, clock difference and SISRE difference, respectively, while the corresponding statistics are increased to 2.5, 6.4 and 5.6 cm, 0.14 ns and 3.4 cm for SECM satellites, respectively.

QZSS
Due to the limited tracking of ground stations as the BDS-2 GEO and IGSO satellites, the epoch-wise orbit and clock differences of QZSS GEO and IGSO satellites between GFZ and WHU precise products on DOY 122 of 2020 also show large fluctuations (see Figure 8). The orbit differences are usually within 0.3 m, except for the GEO satellite J07 in the along-track direction (with an average value of 1.318 m), and the clock difference series even have obvious trends. The RMSs of epoch-wise differences of precise satellite orbit and clock corrections among different ACs as well as SISRE statistics for QZSS IGSO (J01, J02 and J03) and GEO (J07) satellites are listed in Table 7. It is noted that the CNES real-time precise products do not support QZSS satellites, and it is also the case for the ESA precise products in terms of the QZSS GEO satellite. Similar to the BDS-2 non-MEO satellites, especially for the GEO satellites, the consistency of the post-processed precise products among different ACs for the QZSS IGSO and GEO satellites is worse than that for the GPS, GLONASS, Galileo, BDS-2 and BDS-3 MEO satellites, which may be attributed to the worse accuracy of precise products of non-MEO satellites. The differences in radial orbits, along-track orbits, cross-track orbits, clocks and SISREs between GFZ and WHU precise products for QZSS GEO satellite are 5.6, 328.0 and 116.8 cm, 0.24 ns and 32.4 cm, respectively, and the corresponding five statistics between GFZ and ESA precise products for QZSS IGSO satellites are15.4, 6.0 and 6.8 cm, 0.14 ns and 14.8 cm, respectively. Compared with the differences between GFZ and ESA precise products for QZSS IGSO satellites, the corresponding differences between WHU and ESA, and WHU and GFZ precise products seem to be at a similar level in terms of radial orbits and SISRE, but are enlarged by approximately two times in terms of along-track orbits, cross-track orbits, and clocks. Regarding SISRE, the inconsistency of the post-processed precise products among different ACs for the QZSS GEO satellite is twice larger than that for the QZSS IGSO satellites. In comparison to the BDS-2 GEO and IGSO satellites, the SISRE differences among different ACs are increased by about 10 and 5 cm for the QZSS GEO and IGSO satellites, respectively.

Effects of Sun Elevations on Consistency Results
Following Kazmierski et al. [40], the quality of precise products could be affected by the sun elevations above the orbital planes. The orbit determination for low sun elevation angles is a challenge. Usually, the orbit quality is worse for low sun elevation angles due to the normal orbit mode. For further analysis, we re-compute the SISRE difference between GFZ and WHU precise products using the datasets over a month (DOY 122 to 152 in 2020) by dividing the satellites of each constellation into two groups, namely the satellites in eclipsing seasons and the satellites outside eclipsing seasons. When the sun elevations are lower than 8.6 • for BDS GEO/IGSO and QZSS GEO/IGSO, 12.3 • for Galileo, 13.1 • for BDS MEO, 13.5 • for GPS, and 14.3 • for GLONASS, the satellites enter the phase of eclipse. The eclipsing satellites for each constellation during the analysis period are listed in Table A1 of Appendix A. The time spans in the parenthesis in Table A1 refer to the eclipsing periods. There are 10, 15, 13, 2, 2, 2, and 12 eclipsing satellites for GPS, GLONASS, Galileo, QZSS IGSO, BDS-2 IGSO, BDS-2 MEO, and BDS-3 MEO during the analysis period, respectively, while there are no eclipsing satellites for QZSS GEO and BDS-2 GEO. The SISRE statistics are 3.4, 4.6, 2.0, 3.9, 2.5, 2.6, and 14.0 cm for GPS, GLONASS, Galileo, BDS-2 IGSO, BDS-2 MEO, BDS-3 MEO, and QZSS IGSO non-eclipsing satellites, respectively, while the SISRE differences are increased to 4.0, 5.4, 2.8, 5.6, 2.7, 2.5, and 15.0 cm for the corresponding eclipsing satellites, respectively. The consistency of precise products in eclipsing seasons is usually worse than that outside eclipsing seasons by several millimeters, but a significant consistency degradation of 44% can be found for BDS-2 IGSO satellites. It should be noted that the datasets on DOY 145 of 2020 (24 May 2020) for all the BDS-2 IGSO satellites as well as the datasets from 21:10:00 to 21:15:00 on DOY 135 of 2020 (14 May 2020) for all the BDS-3 MEO satellites are excluded from our analysis in this section, as the clock differences between GFZ and WHU during these periods of time show an anomalistic behavior, but it is not the case for those between GFZ and ESA. Thus, the derived SISRE statistics for the BDS-2 IGSO and BDS-3 MEO satellites in this section are better than those shown in Tables 5 and 6.

Datasets, Statistics and Processing Strategies
To assess the performance of multi-system (GNSS/RNSS) PPP with the precise products from different analysis centers, the datasets from 29 globally distributed MGEX stations spanning a week from DOY 138 to 144 in 2020 (i.e., 17-23 May 2020) are employed. It is important to notice that we do not perform a continuous PPP processing with the 7-day data in this study. The seven-day data are divided into several sub-sets (one day from GPS time 00:00:00 to 23:59:30 for each sub-set), and the cross-day data are not processed together. Actually, in this study, the 24-h observation files in RINEX format provided by MGEX are processed day-by-day, and the daily PPP solutions of each station are used for analysis (i.e., the long-term performance rather than the short-term performance for PPP). Thus, a total of 203 PPP experiments based on the daily observations are conducted for each system combination case for each precise product. The geographical distribution of the selected stations is shown in Figure 9.  Table 8. The open source program package RTKLIB is used for the data processing in this paper [41], but we have made improvements in many aspects.

Static Position Solutions
In order to show more details of the position solutions, we select the static positioning results at the station KIR0 on 17 May 2020 (DOY 138 of 2020) using ESA final products as an example, and the epoch-wise position errors for the four single-system cases and the five-system combination case are illustrated in Figure 10. The static GRECJ-PPP can achieve a centimeter-level positioning accuracy within only several minutes. The main difference of static position errors among different cases lies in the first six hours, and then the four single-system cases can almost achieve the same positioning accuracy as GRECJ-PPP. Table 9 provides the statistical positioning accuracy and convergence time of post-processed PPP with WHU products for different system combinations in the static mode, while Tables 10 and 11 show the corresponding statistics using GFZ and ESA products, respectively. It should be noted that all the following statistics are only valid for the long-term performance of PPP. The positioning accuracy does not exhibit significant difference for different system combination cases, and the difference of positioning accuracy among different combination cases is smaller than 0.4, 0. Different from the positioning accuracy, there is obvious difference in the convergence time for different combination cases. Due to the increased satellite number and improved satellite geometry, multi-system integration can shorten the convergence time. Compared with R-PPP, E-PPP and C-PPP, the improvement on the convergence time is 53%/57%/53%, 22%/47%/47% and 41%/58%/38% for the cases using WHU products, 48%/49%/47%, 28%/34%/42% and 44%/53%/39% for the cases using GFZ products, and 59%/64%/52%, 39%/37%/46% and 55%/62%/44% for the cases using ESA products, after an integration with GPS in the east/north/up directions, respectively.   As the number of available QZSS satellites is limited, the convergence time improvement of GJ-PPP over G-PPP is marginal, which is confined to tens of seconds. GRECJ-PPP achieves the shortest convergence time, the numerical value of which is 6.2/3.1/6.8, 6.9/3.1/7.1 and 5.9/2.6/6.3 min for the cases with WHU, GFZ and ESA products in the east/north/up directions, respectively. GRECJ-PPP shortens the convergence time by 59%/48%/56%, 42%/46%/54% and 43%/52%/57% over G-PPP for the three different precise products in the three directions, respectively. As for the effect of employed precise products on the convergence time, it is ignorable for C-PPP, GR-PPP, GC-PPP and GRECJ-PPP, and the difference is smaller than 1.7 min when adopting the precise products from different ACs. The convergence time of G-PPP with WHU products in the east direction is longer than that using the other two post-processed products by several minutes, and it is also the case for R-PPP with ESA products in the east direction, for E-PPP with WHU products in the up direction, for GE-PPP with WHU products in the east and up directions, and for GJ-PPP with WHU products in the east direction. The post-processed positioning accuracy of static PPP with the precise products from different ACs is found to be at a similar level, and the accuracy difference among each other is smaller than 0.3 cm. Figure 11 exhibits the epoch-wise position errors of static PPP using CNES real-time precise products for the four single-system cases and GREC-PPP at station KIR0 on 17 May 2020. It is seen that the positioning performance of E-PPP, C-PPP and R-PPP within the first five hours is significantly worse than G-PPP. The real-time static R-PPP has the longest convergence time, while the GREC-PPP only needs several minutes to obtain a centimeter-level position accuracy. In contrast to the post-processed results, the positioning performance of the single-system cases (except for G-PPP) is still inferior to GREC-PPP even after a long processing time. Table 12 illustrates the real-time PPP statistical results with CNES precise products in the static mode.
It should be noted that all the following statistics are only valid for the long-term performance of PPP. Similar to the post-processed PPP results, the positioning accuracy of G-PPP is comparable to that of the multi-system combined cases, but the R-PPP, E-PPP and C-PPP obtain worse positioning accuracy by several millimeters in the east and up directions than the other five combination cases. Similarly, the joint use of multi-system observations can reduce the convergence time. The convergence time improvement is 74%/74%/60% for GR-PPP over R-PPP, 46%/47%/47% for GE-PPP over E-PPP, 61%/68%/52% for GC-PPP over C-PPP, and 28%/22%/31% for GREC-PPP over G-PPP in the east/north/up directions, respectively. The positioning performance of real-time static PPP is only slightly worse than that of post-processed static PPP for G-PPP, E-PPP and dual-system combined PPP, but it is not the case for R-PPP and C-PPP. The performance degradation is smaller than 6 min on the convergence time, and 0.7 cm on the positioning accuracy for G-PPP, E-PPP and dual-system combined PPP. In comparison with post-processed PPP, the convergence time of real-time PPP is increased by 25.4-28.1, 14.1-15.2 and 14.4-15.6 min, and the positioning accuracy is reduced by 0.9, 0.5-0.6 and 0.6-0.7 cm for R-PPP in the three directions, respectively. Regarding C-PPP, the corresponding performance degradation in the three directions is 11.0-11.9, 9.5-10. 3 Figure 12 depicts the epoch-wise kinematic positioning errors for the four single-system cases and the five-system combination case at the station KIR0 on 17 May 2020 using ESA final products. In the kinematic mode, both the positioning accuracy and convergence time can significantly benefit from the multi-system integration. Most kinematic position errors of the single-system cases are less than 5, 5 and 8 cm after the position filters converge in the east, north and up directions, respectively, while the corresponding varying range of position errors is reduced by half for the kinematic GRECJ-PPP. Tables 13-15 present the statistical kinematic positioning accuracy and convergence time of post-processed PPP with WHU, GFZ and ESA products for different system combinations, respectively. It should be noted that all the following statistics are only valid for the long-term performance of PPP. As to the short-term performance of PPP, the relevant information could be found in the literatures [6][7][8]. Different from the static PPP solutions, both the positioning accuracy and convergence time of kinematic PPP solutions can benefit from the multi-system combination. By introducing GPS observations, with the use of WHU, GFZ and ESA products, we can achieve a convergence time improvement of 82-85%/85-90%/84-88% for GR-PPP over R-PPP, of 18-39%/53-63%/57-63% for GE-PPP over E-PPP, and of 54-68%/65-71%/59-67% for GC-PPP over C-PPP, as well as an accuracy improvement of 65-73%/65-70%/54-55% for GR-PPP over R-PPP, of 66-72%/57-75%/48-54% for GE-PPP over E-PPP, and of 59-68%/52-61%/49-52% for GC-PPP over C-PPP in the east/north/up directions, respectively. In contrast to G-PPP, GJ-PPP slightly shortens the convergence time (less than 10%), but the improvement on the positioning accuracy is relatively obvious (15-20%/14-18%/12-16% in the three directions). Compared with G-PPP, GRECJ-PPP reduces the convergence time by 72%/73%/69%, 67%/74%/69%, and 59%/73%/60% to 8.6/3.8/8.0, 8.0/4.0/8.1, and 6.8/3.3/7.8 min, and improves the positioning accuracy by 71%/53%/40%, 63%/64%/43%, and 56%/62%/40% to 1.0/0.8/2.6, 1.0/0.8/2.5, and 1.1/0.8/2.6 cm using WHU, GFZ and ESA products in the east/north/up directions, respectively.   When using the precise products from different ACs, several combination cases exhibit distinct positioning performance. G-PPP with WHU products obtains slightly worse positioning accuracy in the east direction (an accuracy degradation of 7-9 mm) than that with GFZ and ESA products, and it is also the case for GJ-PPP. The positioning accuracy of E-PPP with WHU products is improved by about 1 cm compared with that with GFZ and ESA products in all three directions. As to other combination cases, the accuracy difference is usually confined to 0.5 cm. The differences of convergence time between the post-processed PPP with GFZ and ESA products are smaller than 3 min for all combination cases, except for G-PPP, R-PPP and GJ-PPP cases. The convergence time difference is also less than 3 min for C-PPP, GR-PPP and GRECJ-PPP with the precise products from different ACs. When using ESA products, G-PPP, GJ-PPP and GC-PPP achieve the shortest convergence time in comparison with those with WHU and GFZ products. Both E-PPP and GE-PPP cases with WHU products have the longest convergence time than those with GFZ and ESA products. The difference of convergence time can be up to about 10 min for R-PPP with WHU, GFZ and ESA products in all three directions. Compared with the static post-processed PPP results shown in Tables 9-11, both the positioning errors and convergence time of the corresponding kinematic ones are increased by a factor of 3.5 times on average. Figure 13 illustrates the epoch-wise position errors of kinematic PPP for the four singlesystem cases and GREC-PPP at the station KIR0 on 17 May 2020 using CNES real-time precise products. It can be seen that the real-time kinematic positioning performance of the single-system cases is much inferior to GREC-PPP, especially for C-PPP, E-PPP and R-PPP. The epoch-wise position errors of the single-system cases vary within 8, 8 and 10 cm in the east, north and up directions, respectively, while the corresponding errors of real-time kinematic GREC-PPP have a varying range of 3, 3 and 5 cm in the three directions, respectively. Table 16 illustrates the statistical kinematic positioning accuracy and convergence time of real-time PPP with CNES products for different system combinations. It should be noted that all the following statistics are only valid for the long-term performance of PPP. Both the real-time kinematic positioning accuracy and convergence time can benefit from the multi-system combination. The convergence time improvement is 78%/86%/82% for GR-PPP over R-PPP, 39%/59%/55% for GE-PPP over E-PPP, 73%/72%/73% for GC-PPP over C-PPP, and 42%/49%/43% for GREC-PPP over G-PPP in the east/north/up directions, respectively, while the corresponding improvement on the positioning accuracy is 64%/50%/54%, 61%/55%/49%, 60%/59%/52%, and 48%/38%/31% in the three directions, respectively. The real-time kinematic GREC-PPP can achieve a convergence time of 11.5/6.9/13.0 min, and a positioning accuracy of 1.7/1.6/3.6 cm in the three directions, respectively. The real-time kinematic PPP even has slightly better positioning performance compared with some cases of post-processed kinematic PPP results shown in Tables 13-15, including the convergence time of G-PPP with WHU and GFZ products, of E-PPP with WHU products in east and up directions, and of GE-PPP with WHU products in east direction, as well as the positioning accuracy of G-PPP and E-PPP with WHU products in east direction. Similar to the static PPP, the performance degradation is obvious for real-time kinematic R-PPP and C-PPP over the post-processed kinematic ones. Compared with the post-processed kinematic R-PPP, the convergence time of real-time kinematic R-PPP is lengthened by 9.  Table 12).

Discussion
According to the above analysis, the consistency of the precise products from different ACs is usually at a centimeter level, sub-decimeter level and decimeter level for the GNSS/RNSS MEO, IGSO and GEO constellations, respectively. Except for the WHU products for GPS and GLONASS satellites, the SISRE consistency among various postprocessed precise products is usually better than 1.5, 3.0, 2.5, 4.5, 4.0, 22.5, 10.5, 32.5 and 15.5 cm for the GPS, GLONASS, Galileo, BDS-2 MEO, BDS-3 MEO, BDS-2 GEO, BDS-2 IGSO, QZSS GEO and QZSS IGSO satellites, respectively, while the increased SISRE difference between the real-time and post-processed precise products is usually smaller than 3.0, 6.0, 5.5, 5.5, 7.5, 55.5 and 9.5 cm for the GPS, GLONASS, Galileo, BDS-2 MEO, BDS-3 MEO, BDS-2 GEO and BDS-2 IGSO satellites, respectively. As to the benefits from multi-system combination for PPP processing, the convergence time in both static and kinematic modes as well as the kinematic positioning accuracy can be continually improved when including more satellite systems on the basis of GPS, while multi-system combination has little effect on the static positioning accuracy due to the enough observation information after a long processing time even for standalone GPS.
Following Chen et al. [20], with the use of the datasets in 2019, the post-processed GRECJ-PPP (only covering BDS-2 satellites for BDS) in the kinematic mode could achieve a positioning accuracy of 1.4, 1.2 and 3.3 cm in the east, north and up directions, respectively. Based on the datasets in 2018, the real-time kinematic GREC-PPP (only covering BDS-2 satellites for BDS) with CNES real-time precise products could provide a positioning accuracy of 3.4, 2.6 and 5.9 cm in the three directions, respectively [21]. In this study, the kinematic positioning accuracy is improved to 1.0-1.1, 0.8 and 2.5-2.6 cm for the post-processed GRECJ-PPP, and 1.7, 1.6 and 3.6 cm for the real-time GREC-PPP in the three directions, respectively, which can be attributed to the increased number of available satellites from BDS-3 and the improved quality of precise products recently. It should be noted that all these numbers are only valid for the long-term performance of PPP.

Conclusions
With the fully deployed BDS-3, the number of the available GNSS/RNSS satellites can be more than 130. In addition, several analysis centers have been providing the precise satellite products containing the information of all the five satellite systems, namely BDS (BDS-2/BDS-3), GPS, GLONASS, Galileo and QZSS. As a satellite-based positioning technology, the performance of PPP under the current constellation status draws an increasing attention from the GNSS/RNSS community. In this study, the post-processed five-system integrated PPP (i.e., GRECJ-PPP) with the ESA, GFZ and WHU precise products, as well as the real-time four-system combined PPP (i.e., GREC-PPP) with the CNES precise products for the purpose of completeness, is analyzed. The performance evaluation is conducted in both static and kinematic modes in terms of positioning accuracy and convergence time. For comparison, we also provide the single-system and dual-system PPP processing results.
Besides, we characterize the consistency of precise satellite products from different analysis centers. It should be noted that all the following numbers are only valid for the long-term performance of PPP based on daily observations. Based on the SISRE values, the consistency among the ESA, GFZ and WHU precise products is usually better than 1.5, 3.0, 2.5, 4.5, 4.0, 22.5, 10.5, 32.5 and 15.5 cm for the GPS, GLONASS, Galileo, BDS-2 MEO, BDS-3 MEO, BDS-2 GEO, BDS-2 IGSO, QZSS GEO and QZSS IGSO satellites (except for the WHU products for GPS and GLONASS satellites), respectively. As for the consistency between the CNES real-time products and the three post-processed precise products, the SISRE values are increased, and the corresponding statistic is usually smaller than 3.0, 6.0, 5.5, 5.5, 7.5, 55.5 and 9.5 cm for the GPS, GLONASS, Galileo, BDS-2 MEO, BDS-3 MEO, BDS-2 GEO and BDS-2 IGSO satellites, respectively. In addition, the consistency of precise products in eclipsing seasons is usually worse than that outside eclipsing seasons. The post-processed GRECJ-PPP can provide a convergence time of 5.9-6.9/2.6-3.1/6. The static positioning accuracy of post-processed PPP does not exhibit significant difference for different system combination cases, which is smaller than 0.4/0.2/0.5 cm in the three directions, respectively, and the static positioning accuracy of post-processed GRECJ-PPP with WHU, GFZ and ESA precise products is comparable to each other. The improvement (with a statistic of 18-90%) on the convergence time is significant after adding the GPS observations to the GLONASS-only, Galileo-only and BDS-only PPP processing in both static and kinematic modes in both real-time and post-processed scenarios. The post-processed GRECJ-PPP shortens the static convergence time by 42-59%/46-52%/54-57% over the post-processed GPS-only PPP in the three directions, respectively. As for the effect of employed post-processed precise products on the static convergence time, it is ignorable for GRECJ-PPP, and the difference is smaller than 1 min when adopting the precise products from different ACs. The convergence time improvement in the static mode is 28%/22%/31% for real-time GREC-PPP over real-time GPS-only PPP in the three directions, respectively. Similar to the post-processed PPP results, the static positioning accuracy of real-time GPS-only PPP is comparable to the multi-system combined cases, but the other three real-time single-system PPP cases obtain worse static positioning accuracy by several millimeters. The positioning performance of real-time static PPP is only slightly worse than that of post-processed static PPP for GPS-only, Galileo-only and dual-system combined PPP, but it is not the case for GLONASS-only and BDS-only PPP. The postprocessed static GRECJ-PPP shortens the convergence time by 21-32%/40-50%/37-44%, and improves the positioning accuracy by 50-67%/63-75%/15-23% over the real-time static GREC-PPP in the three directions, respectively. Different from the static PPP solutions, both the positioning accuracy and convergence time of post-processed and real-time kinematic PPP solutions can benefit from the multisystem combination. By introducing GPS observations, we can achieve a convergence time improvement of 18-90% and 39-86%, as well as an accuracy improvement of 48-75% and 49-64% for dual-system combined PPP over single-system PPP in the kinematic mode in the post-processed and real-time scenarios, respectively. Compared with the post-processed kinematic GPS-only PPP, the post-processed kinematic GRECJ-PPP reduces the convergence time by 59-72%/73-74%/60-69%, and improves the positioning accuracy by 56-71%/53-64%/40-43% in the east/north/up directions, respectively, while the corresponding improvement is 42%/49%/43% and 48%/38%/31% for real-time kinematic GREC-PPP over real-time kinematic GPS-only PPP in the three directions, respectively. When using the post-processed precise products from different ACs, several combination cases exhibit distinct kinematic positioning performance, and the difference can be as large as 1 cm and 10 min in terms of position accuracy and convergence time, respectively. However, the corresponding difference, which is smaller than 0.1 cm and 2 min, is not obvious for the post-processed kinematic GRECJ-PPP. Compared with the static post-processed PPP results, both the positioning errors and convergence time of the corresponding kinematic ones are increased by a factor of 3.5 times on average. The real-time kinematic PPP even has slightly better positioning performance compared with some cases of post-processed kinematic PPP results. The post-processed kinematic GRECJ-PPP shortens the convergence time by 25-41%, 42-52% and 38-40%, and improves the positioning accuracy by 35-41%, 50% and 28-31% in comparison with the real-time kinematic GREC-PPP in the three directions, respectively. The real-time kinematic PPP increases the position errors and convergence time by a factor of 2.7 times on average over the real-time static PPP.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found from MGEX.