Multi-GNSS Combined Orbit and Clock Solutions at iGMAS

Global navigation services from the quad-constellation of GPS, GLONASS, BDS, and Galileo are now available. The international GNSS monitoring and assessment system (iGMAS) aims to evaluate the navigation performance of the current quad systems under a unified framework. In order to assess impact of orbit and clock errors on the positioning accuracy, the user range error (URE) is always taken as a metric by comparison with the precise products. Compared with the solutions from a single analysis center, the combined solutions derived from multiple analysis centers are characterized with robustness and reliability and preferred to be used as references to assess the performance of broadcast ephemerides. In this paper, the combination method of iGMAS orbit and clock products is described, and the performance of the combined solutions is evaluated by various means. There are different internal precisions of the combined orbit and clock for different constellations, which indicates that consistent weights should be assigned for individual constellations and analysis centers included in the combination. For BDS-3, Galileo, and GLONASS combined orbits of iGMAS, the root-mean-square error (RMSE) of 5 cm is achieved by satellite laser ranging (SLR) observations. Meanwhile, the SLR residuals are characterized with a linear pattern with respect to the position of the sun, which indicates that the solar radiation pressure (SRP) model adopted in precise orbit determination needs further improvement. The consistency between combined orbit and clock of quad-constellation is validated by precise point positioning (PPP), and the accuracies of simulated kinematic tests are 1.4, 1.2, and 2.9 cm for east, north, and up components, respectively.


Introduction
Over the past decade, apart from the ongoing modernization of GPS and GLONASS, two new global navigation satellite systems (GNSSs) have been developed, namely the BeiDou navigation satellite system (BDS-3) of China and the Galileo system of Europe. Both BDS-3 and Galileo have already been providing positioning, navigation, timing (PNT), and other services to global users. The 'big 4' GNSSs allow users to have more choices to access reliable and accurate PNT services from multi-GNSS [1][2][3]. Precise multi-GNSS orbit and clock products are available from the international GNSS service (IGS) [4][5][6][7][8]. Combining the products from different analysis centers (ACs) improves the stability and continuity compared with the individual AC products, which enables the use of GNSS constellations for scientific research, technical tests, and various applications related to satellite navigation. Furthermore, the combined products can be used as references to objectively evaluate the service performance of the multi-constellation.
The combined solutions of GNSS orbit and clocks were firstly proposed by IGS, due to its robustness compared with the solution from a single AC. Seven analysis centers multiple ACs, and present a detailed analysis of the accuracy of iGMAS combined products in terms of internal precision and comparison with IGS products. The combined orbit products are also assessed by Satellite Laser Ranging (SLR) observations. The consistency between the combined orbit and clock is assessed by precise point positioning in this section. Finally, conclusions and a summary are provided in Section 5.

Overview of iGMAS ISC
iGMAS is a super-large and complex system, which is composed of several subsystems including the global tracking stations, the data centers, the analysis centers, the information combination and service center (ISC), the monitoring and assessment center, and the operation control center. Each sub-system supports and connects with each other in function to ensure the stable operation of iGMAS. The iGMAS structure is depicted in Figure 1. The global tracking stations of about thirty collect multi-GNSS observations and transmit them to the data centers for classification, management, and storage. The analysis centers process the observations collected by tracking stations to generate high-precision GNSS products and submit the products to ISC. Currently, twelve ACs contribute products to iGMAS ISC routinely: Beijing Aerospace Control Center (BAC) [24], Chinese Academy of Surveying and mapping (CGS), Chang'an University (CHD) [25], China University of Mining and Technology (CUM) [26], Institute of Geodesy and Geophysics (CGS), Information Engineering University (LSN) [27], National Time Service Center (NTS) [28], Shanghai Astronomical Observatory (SHA), Tongji University (TJU), Xi'an Research Institute of Surveying and Mapping (XRS) [27], Xi'an Satellite Control Center (XSC) [29], and Wuhan University (WHU) [30].
The ISC is the center of iGMAS products reprocessing, which evaluates the quality of GNSS products submitted by individual analysis centers, and then reprocesses them to generate the combined products and the corresponding accuracy indexes. The products are delivered to users through the portal website of http://www.igmas.org (accessed on 30 December 2021 for Chinese version), which is supported by monitoring stations, data centers, and analysis centers.
The combined products mainly consist of the precise orbits and clocks, the coordinates of tracking stations, the Earth rotation parameters, and atmospheric environment parameters. These combined products are stored and published by data The global tracking stations of about thirty collect multi-GNSS observations and transmit them to the data centers for classification, management, and storage. The analysis centers process the observations collected by tracking stations to generate high-precision GNSS products and submit the products to ISC. Currently, twelve ACs contribute products to iGMAS ISC routinely: Beijing Aerospace Control Center (BAC) [24], Chinese Academy of Surveying and mapping (CGS), Chang'an University (CHD) [25], China University of Mining and Technology (CUM) [26], Institute of Geodesy and Geophysics (CGS), Information Engineering University (LSN) [27], National Time Service Center (NTS) [28], Shanghai Astronomical Observatory (SHA), Tongji University (TJU), Xi'an Research Institute of Surveying and Mapping (XRS) [27], Xi'an Satellite Control Center (XSC) [29], and Wuhan University (WHU) [30].
The ISC is the center of iGMAS products reprocessing, which evaluates the quality of GNSS products submitted by individual analysis centers, and then reprocesses them to generate the combined products and the corresponding accuracy indexes. The products are delivered to users through the portal website of http://www.igmas.org (accessed on 30 December 2021 for Chinese version), which is supported by monitoring stations, data centers, and analysis centers.
The combined products mainly consist of the precise orbits and clocks, the coordinates of tracking stations, the Earth rotation parameters, and atmospheric environment parameters. These combined products are stored and published by data centers, together with the inter-frequency biases, the ionospheric scintillation indexes, and the integrity products provided by the monitoring and assessment center.
Both iGMAS and IGS provide high-precision GNSS products to global users. Although iGMAS and IGS provide the same types of products, iGMAS has achieved, for the first time, the quad-constellation, i.e., BDS/GPS/GLONASS/Galileo combination of the orbit and clock products, whereas IGS only provides GPS/GLONASS combined orbits and GPS-only combined clocks. Until now, the IGS multi-GNSS product combination experiments are still undergoing, and the experimental multi-GNSS combined orbits are available. Please be noted that the time used in iGMAS products is BDT.  Figure 2 presents the flowchart of multi-GNSS orbit and clock combination of iGMAS. Since both the ACs' orbits and the combined solutions are used for consistency correction in the procedure of clock combination, the multi-GNSS orbits are firstly combined, then it is the clock combination's turn. Apart from the orbit and clock products used as input data, the earth orientation parameters and multi-GNSS broadcast ephemerides are also utilized in the combination procedure.

Combination Strategy of Multi-GNSS Orbit and Clock
Although iGMAS and IGS provide the same types of products, iGMAS has achieved, for the first time, the quad-constellation, i.e., BDS/GPS/GLONASS/Galileo combination of the orbit and clock products, whereas IGS only provides GPS/GLONASS combined orbits and GPS-only combined clocks. Until now, the IGS multi-GNSS product combination experiments are still undergoing, and the experimental multi-GNSS combined orbits are available. Please be noted that the time used in iGMAS products is BDT. Figure 2 presents the flowchart of multi-GNSS orbit and clock combination of iGMAS. Since both the ACs' orbits and the combined solutions are used for consistency correction in the procedure of clock combination, the multi-GNSS orbits are firstly combined, then it is the clock combination's turn. Apart from the orbit and clock products used as input data, the earth orientation parameters and multi-GNSS broadcast ephemerides are also utilized in the combination procedure.

Orbit Combination
Due to different tracking observations and strategies of precise orbit determination adopted by analysis centers, there are systematic biases existing in the terrestrial reference frames (TRF) realized by precise orbits from analysis centers. Therefore, it is necessary to re-process the different orbit solutions to generate combined results, which is consistent with the international TRF (ITRF). Moreover, the combined solutions are more continuous over a large time scale, compared with the lack of some solutions from analysis centers on some days.
Generally, the multi-GNSS orbit combination consists of three steps: the orbital realized TRF calibration, the weight determination for analysis centers, and the orbit averaging.
Step (1): orbital realized TRF calibration. According to [18], the orientation of orbits, station coordinates, and earth rotation parameters (ERP) are highly correlated, it is possible to unify the orientations of the three kinds of products. Since the earth rotation parameters and orbits are simultaneously generated by iGMAS analysis centers, we

Orbit Combination
Due to different tracking observations and strategies of precise orbit determination adopted by analysis centers, there are systematic biases existing in the terrestrial reference frames (TRF) realized by precise orbits from analysis centers. Therefore, it is necessary to re-process the different orbit solutions to generate combined results, which is consistent with the international TRF (ITRF). Moreover, the combined solutions are more continuous over a large time scale, compared with the lack of some solutions from analysis centers on some days.
Generally, the multi-GNSS orbit combination consists of three steps: the orbital realized TRF calibration, the weight determination for analysis centers, and the orbit averaging.
Step (1): orbital realized TRF calibration. According to [18], the orientation of orbits, station coordinates, and earth rotation parameters (ERP) are highly correlated, it is possible to unify the orientations of the three kinds of products. Since the earth rotation parameters and orbits are simultaneously generated by iGMAS analysis centers, we decide to align the orientation of different orbits to ITRF through ERP, and the following formulas are used.
where dX p and dY p are the differences between ERPs of analysis centers and IERS standard EOP data. Due to the inconsistent model used by some analysis centers and insufficient accuracy of orbits for sometimes, the ERP residuals are checked before orientation correction. A threshold of 0.3 mas proposed by IGS is used to detect outlier analysis centers. The analysis centers with ERP residuals larger than 0.3 mas are excluded from the combination.
Step (2): weight determination. The weight of the analysis center is mainly divided into two parts, that is the initial weight and the final weight, and both are calculated by using the root mean square error (RMS) of the analysis center orbit with regard to the combined solutions. The initial combined solutions are derived from the median of all analysis centers for individuals, and the initial weights are identical for all analysis centers. Considering the differences in the accuracy of various satellite orbits, however, we compute weights depending not only on the ACs and satellites but also on the constellations, instead of one weight. Based on the orbits of combined and analysis centers, the observation equation is established as follows, and weights are iteratively updated: where O sys cmb and O sys AC are the satellite orbits of combined and individual analysis centers for system sys (i.e., GPS, GLONASS, BDS, Galileo). D AC , S AC , and R AC are the translation vector, scale, and rotation matrix, respectively, which can be expressed by and θ Z are the three rotation angles.
There are some deficiencies in the solar radiation pressure models adopted by analysis centers, especially for the newly constructed global navigation system (e.g., the ECOM-5 without a priori model is used for BDS and Galileo by most analysis centers). Moreover, the PCO and PCV adopted in data processing are not consistent among analysis centers, the values derived either IGS14.atx or calibration of analysis center are used by different analysis centers. As a result, different sets of similarity transformation parameters are designed for each satellite system in this study. Due to the much lower precision of GEO orbits, they are not included in the weight determination. The weights of analysis centers could be calculated based on the RMS derived from orbit residuals as follows: where nac indicates the number of analysis centers. RMS sys AC is the RMS of system sys, which can be computed by In Equation (4), np sys AC is the total epochs of all satellites for system sys, and P sys,Last AC indicates the weight from the last iteration with the initial value of one. Although the least square estimation method adopted by orbit combination is unbiased, it is easy to be affected by the abnormal satellites, which is always the case for the eclipse satellites. As a result, it is necessary to detect and reduce the effects induced by these satellites on the combined solutions. In this paper, the equivalent weight method is used to weaken the impacts of anomalous satellites using the following formula: Step (3): orbits averaging. The weighted averaging orbit could be obtained from the multiple ACs using the estimates of transformation parameters and AC specific weight as follows It is noticed that the orbits of some satellites are not provided by some of ACs, especially for the satellites during eclipse or attitude switch season. The number of ACs may be different for the two satellites, isat and jsat, resulting in different weights of the two satellites for the individual ACs. Fortunately, the orbit dynamics are still satisfied approximately when the sum weight of all analysis centers included in the combination for each satellite is one [10].

Clock Combination
In order to keep the consistency between combined clock and orbit solutions, a consistency correction should be applied for clock products of analysis centers by using their corresponding orbits. Meanwhile, the reference clocks and observations adopted by analysis centers are not identical, resulting in satellite-dependent bias between two analysis centers, apart from time scale difference [18,19]. Moreover, differing from the orbits, the clock products are highly susceptible to unoptimized constraints, reference clock, and other unmodelled errors. Generally, the multi-GNSS clock combination can be obtained by the following steps.
Step (1): Consistency correction. To keep the combined clock consistent with the orbit, the orbit-clock consistency correction should be added to the clock products of analysis centers. The correction of the consistency between clock and orbit (i.e., dC sys,isat AC ) can be calculated using the formula below. Considering that the interval of the orbits are always 900 s, and the sampling intervals of clock products are generally 300 s or less (e.g., 30 s for iGMAS), the Everette interpolation method is adopted in this study to interpolate the orbits [31].
where O sys,isat AC and O sys,isat cmb are the orbits of analysis centers and combined one, respectively. e is the unit vector of satellite position.
Step (2): Preprocessing of clock products. Clock preprocessing mainly includes outliers and clock jump detection. First, the clock data are converted into frequency data, and then the abnormal data in the frequency domain are detected by using the robust median detection method. It is easy to find that one frequency outlier indicates the existence of clock jump, and two consecutive outliers in the clock frequency sequence indicate the existence of phase outliers. Second, the corresponding clock phase data for the detected frequency anomaly data are identified separately, and the abnormal value is deleted and the clock phase data are segmented according to the phase jump position.
Step (3): Datum unification for clocks of analysis centers. The GPS clock combination of IGS could be used to align the time scale of clock products of analysis centers to IGST. It is noticed that there is a time delay of few days to two-week delay for the IGS products. As a result, we decided to unify the time scale of clock products to the realization of broadcast clock. Considering the insufficient accuracy of nanoseconds for broadcast clock and using one set of parameters for all satellites cannot completely remove the systematic errors between different AC clocks. It was pointed out that the clock biases among all satellites, caused by realigning the combined GLONASS clock using IGS strategy, make the analysis center unable to obtain a reliable weight, which will eventually affect the reliability of the combined products [20]. Therefore, an approach of linear fitting for individual satellites is used to reduce the satellite-dependent bias and align the clocks to the reference analysis center. Considering their relatively stable performance of precision, we select three analysis centers as the reference clock candidates in this study, that is BAC, SHA, and WHU. Moreover, the analysis center with the best data fitting accuracy and the least phase jump among the three analysis centers is used as the reference. The individual satellite clock is aligned to the reference analysis center by the following formula. Step (4): Robust estimation of parameters. The clocks after consistency correction and datum alignment are used to establish the observation equation. The median values are obtained from all analysis centers and taken as the initial combined clock. A similar observation equation as (8) is also established between clocks of individual analysis center and the combined one, and one set of offset and drift for each system is estimated. In order to reduce the impacts of undetected outliers or clock data for satellites during the eclipse, the IGGIII weight function [32] based on posterior residuals is used for the iterative estimation of parameters. The equivalent weight of the satellite j of the analysis center i at the epoch k can be computed by where k 0 generally has a range between 1.0 and 1.5, and k 1 is selected in the range of 2.5 to 3.0. This paper selects 1.5 and 3.0, respectively. q j i,k is the covariance of posterior residuals of clock products. σ i,0 is the posteriori unit weight variance factor, which could be calculated by σ i,0 = med( v j i,k )/0.6745, based on the posteriori residuals of all satellite clocks provided by analysis center i.
Step (5): Clock-weighted averaging. Considering that the product accuracy of different systems in the same analysis center may be different, the weights are calculated separately for each system and the weighted average clock combination of individual satellites is calculated by the following formula where C j i,k is the calibrated clock of analysis center with consistency correction and alignment of the time scale. w i is the weight of individual analysis centers, which is determined by the weighted RMS of clock residuals for each system by w i = 1 (C sys i,RMS ) 2 .

Assessments of iGMAS Orbits and Clocks
In order to verify the effectiveness of the orbit and clock combination strategies in this paper, the multi-GNSS orbits and clocks from 3 January 2021 to 24 April 2021 submitted by iGMAS ACs are collected and used to generate combined solutions. In this section, the internal precisions of combined solutions are first summarized. Then, the consistency between iGMAS and IGS is assessed. In addition, satellite laser ranging (SLR) normal points from ILRS are used to validate the combined orbits. Observations from five iGMAS stations are also collected to conduct precise point positioning (PPP) using iGMAS combined orbits and clocks.

Internal Precisions of Combined Orbits and Clocks
As an internal validation of combined orbits and clocks, the RMS of the differences between the AC products and corresponding combined products is statistically calculated. These statistical values can reflect the internal precision of the combined orbits and clocks. The calculation formula can be as follows: Whether for combined orbit or clock, an obvious heavy-tail distribution phenomenon exists in the internal precision, indicating that there are large differences between different AC products. It is worth mentioning that the influence of low-precision products can be partially suppressed by reducing the weight.  Whether for combined orbit or clock, an obvious heavy-tail distribution phenomenon exists in the internal precision, indicating that there are large differences between different AC products. It is worth mentioning that the influence of low-precision products can be partially suppressed by reducing the weight.          Figures 7 and 8 show the consistency between iGMAS and IGS final combined GPS orbits and clocks, respectively. The 1D errors of all GPS satellite orbits are basically between 0.6 and 3.0 cm, which are better than 2.0 cm for most satellites. Moreover, the median RMSs of all GPS satellite orbits are relatively stable, with an average value of 1.5 cm. As shown in Figure 7, the standard deviations of clocks mostly show an agreement of 20-90 ps among all GPS satellites, and the median value is 62 ps. In terms of GLONASS, the 1D errors of most satellite orbits are better than 6.0 cm with a median RMS of 4.5 cm (Figure 9). cm. As shown in Figure 7, the standard deviations of clocks mostly show an agreement of 20-90 ps among all GPS satellites, and the median value is 62 ps. In terms of GLONASS, the 1D errors of most satellite orbits are better than 6.0 cm with a median RMS of 4.5 cm (Figure 9).

SLR Validation of iGMAS Combined Orbits
As a key feature, all GLONASS, Galileo and BDS satellites are equipped with laser ranging reflector arrays, enabling high-precision two-way ranging measurements. Thus, satellite laser ranging (SLR) observations collected by the International Laser Ranging Service (ILRS) can be used as an independent validation of iGMAS satellite orbits. Figure 10 shows the quantiles of the SLR residuals for iGMAS final orbits. Table 1 summarizes the mean biases as well as the standard deviations of SLR residuals for iGMAS final orbits. For BDS-2 satellites, the standard deviations are between 4 and 6 cm for MEO and IGSO satellites. However, the SLR residuals show a bias of about 13 cm with a standard deviation of 17.5 cm for C01, which is a GEO satellite. For BDS-3, the orbits of CAST satellites (e.g., C20 and C21) have a bias of 1.2 cm, whereas the orbits of SECM satellites (e.g., C29 and C30) have a bias of −4.7 cm. Moreover, the dispersion degree of SLR residuals for BDS-3 orbits is only 3.5 cm, which is smaller than that of BDS-2. The SLR residuals of Galileo satellites generally have a negative sign of 4.0 cm, which can be explained by the insufficient solar radiation pressure (SRP) modeling, the neglection of the Earth albedo as well as the ignoring of antenna thrust of the Galileo satellites. At the same time, standard deviations of SLR residuals for Galileo and BDS-3 orbits are on a similar level of 3.3 cm. For GLONASS orbits, there is no consistent sign in the bias of SLR residuals, and a standard deviation of 4.2 cm is achieved. residuals for BDS-3 orbits is only 3.5 cm, which is smaller than that of BDS-2. The SLR residuals of Galileo satellites generally have a negative sign of 4.0 cm, which can be explained by the insufficient solar radiation pressure (SRP) modeling, the neglection of the Earth albedo as well as the ignoring of antenna thrust of the Galileo satellites. At the same time, standard deviations of SLR residuals for Galileo and BDS-3 orbits are on a similar level of 3.3 cm. For GLONASS orbits, there is no consistent sign in the bias of SLR residuals, and a standard deviation of 4.2 cm is achieved. Figure 10. Quantiles of SLR residuals for iGMAS final orbits.    In order to analyze whether a systematic error exists in iGMAS combined orbits, the SLR residuals are shown in Figure 11 as a function of the satellite elongation angle for BDS satellites manufactured by CAST and SECM. It can be seen that the SLR residuals of BDS-3 CAST and SECM satellites show an opposite trend. With the increase in the satellite elongation angle, the SLR residuals of CAST satellites are clearly reduced, whereas the SECM satellites present opposite variations. It is worth mentioning that the variation trend in SLR residuals of Galileo and GLONASS satellites can verify the existence of systematic errors, although the trends are opposite (Figure 12). The systematic errors characterized with a linear pattern in SLR residuals were also validated by [17], which should be attributed to the defective SRP model adopted in precise orbit determination. Therefore, a refined SRP model should be proposed for iGMAS analysis centers [33][34][35][36]. Moreover, the ambiguity-fixed solutions are also encouraged to improve combined solutions for the newly constructed constellation [37].
trend in SLR residuals of Galileo and GLONASS satellites can verify the existence of systematic errors, although the trends are opposite (Figure 12). The systematic errors characterized with a linear pattern in SLR residuals were also validated by [17], which should be attributed to the defective SRP model adopted in precise orbit determination. Therefore, a refined SRP model should be proposed for iGMAS analysis centers [33][34][35][36]. Moreover, the ambiguity-fixed solutions are also encouraged to improve combined solutions for the newly constructed constellation [37].

PPP Using iGMAS Combined Orbits and Clocks
In order to further verify the performance of iGMAS combined orbit and clock products, eight globally distributed stations were adopted for PPP test. The processing strategies employed in PPP are given in Table 2.

PPP Using iGMAS Combined Orbits and Clocks
In order to further verify the performance of iGMAS combined orbit and clock products, eight globally distributed stations were adopted for PPP test. The processing strategies employed in PPP are given in Table 2. Figure 13 shows the kinematic PPP results of GPS-only and GPS/GLONASS/BDS/ Galileo multi-system modes at JFNG station on April 24, the last day of the test period. Table 3 summarizes the RMS of each station's positioning errors. In statistics, the positioning results only show the accuracy after a convergence time of about 1 h (i.e., 120 epochs). It can be clearly seen that the combined GPS/BDS/GLONASS/Galileo solutions significantly shorten the convergence time and improve the position series compared with GPS-only PPP, as shown in the positioning results of the 08:00~12:00 period in Figure 13. For the positioning results of the test station, the average accuracies of GPS-only kinematic PPP are 3.2, 2.2, and 4.9 cm in east, north, and up components, respectively, and those of multi-GNSS kinematic PPP are 1.4, 1.2, and 2.9 cm in east, north, and up components, respectively. Obviously, the multi-GNSS combination significantly improves the PPP performance compared with GPS-only solutions, with improvement rates of 57%, 45%, and 41%, respectively.

Conclusions
This paper systematically introduces the strategies of orbit and clock combination

Conclusions
This paper systematically introduces the strategies of orbit and clock combination using multi-GNSS products from iGMAS ACs. In addition, the performance of the combined orbits and clocks is assessed based on half-year products. The results show that the internal precision of 5 mm is achievable for the combined orbits of GPS, whereas the internal precisions of Galileo and BDS are relatively poor at the 1.5 cm level. From the results of SLR validation, there are obvious systematic errors in the iGMAS orbits of different satellite systems, and the systematic errors are manufacturer-dependent. Moreover, there is a negative deviation of 3 cm in the Galileo satellite orbits, which may be attributed to the different non-conservative force models used in ACs' precise orbit determination processing, such as the Earth albedo radiation and the solar radiation pressure model.
Although the internal precision of iGMAS clock products is basically consistent with the orbit, they are more vulnerable to the influence of observation data quality control and processing strategies. Considering the large fluctuation among various satellite systems, it is necessary to further optimize the quality control before clock combination.
The comparison between iGMAS and IGS GPS orbit/clock shows a consistency of 1.5 cm for orbits and 60 ps for clocks. In addition, the performance of iGMAS combined orbit and clock products is further verified by PPP. As the PPP results show, the average accuracies of multi-GNSS kinematic PPP are 1.4, 1.2, and 2.9 cm in E, N, and U components, respectively; the presenting improvement rates are 57%, 45%, and 41% when compared against the GPS-only PPP solution. In the future, the iGMAS multi-GNSS combined orbit and clock products will be used as reference products for GNSS service performance monitoring and evaluation, and further validation for GNSS geoscience using iGMAS combined solutions is needed.

Data Availability Statement:
The precise orbit and clock products of iGMAS analysis centers, as well as the combined solutions conducted by this study can be made available upon request by contacting the corresponding author, G.C. by email.