A prototype for a Multi-GNSS orbit combination

Since 1994, the International GNSS Service (IGS) provides a combination of orbit and clock offset products from its different Analysis Centers. These products are used as input by many software for countless scientific and/or operational applications. They can also be used as independent reference for benchmark experiments. Nevertheless, those products include GPS/GLONASS-only data, despite the fact that nowadays several GNSS constellations are available. We performed modifications on the existing IGS combination software and made it compatible with multi-GNSS products. We present here the results of a combination offive Multi-GNSS Experiment Analysis Centers (ACs) for five constellations (GPS + GLONASS + Galileo + Beidou + QZSS) from GPS week 1800 to 2000 (July 2014 to May 2018). The RMS of combined orbit w.r.t. input ACs ones are ~ 30 ±15 mm for recent weeks of the test period.


Introduction
The International GNSS Service (IGS) is a worldwide scientific organization gathering since the early 1990s numerous Geosciences institutes, universities, along with mapping and space agencies [17]. Some IGS members host so-called Analysis Centers (ACs), having the operational task to provide to the public on a continuous basis high accuracy products, i.e. station positions, satellite orbits & clocks offsets, and Earth orientation parameters (EOPs) for instance.
The IGS Analysis Center Coordinator (ACC) uses all individual AC orbit products as input, to determine an combined orbit solution provided to the community. This combination consists basically of an aligned and weighted mean of all individual AC products. It will be called hereafter the legacy combination [19].
tracking network with stations able to record the maximum of new signals associated with the new systems or satellite generations, along with the enhancement of the different AC processing software. More and more new MGEX ACs (including some established IGS members) start to provide multi-GNSS products, i.e. containing namely Galileo, BeiDou, and QZSS 2 .
Nowadays, the multi-GNSS environment can be considered as enabled since Galileo, BeiDou and QZSS have been officially declared in their operational phases by their respective operating agencies [6,3,16]. Nevertheless, no official IGS combination including the new GNSS is provided so far, despite the fact that some preliminary studies have been performed by [7]. This constitutes an important limitation for the exploitation of the new GNSS for the scientific and engineering communities which usually use IGS products for its applications, in geophysics, oceanography, surveying, etc.. This article exposes the results for a preliminary multi-GNSS orbit combination based on an extension of the legacy combination software. We present hereafter the modifications performed on the software, and the different multi-GNSS AC products used for a test period of 200 weeks between the GPS weeks 1800 and 2000 (July 2014 to May 2018). The orbit combination results are compared with each individual AC, for all satellites and for each separated constellation.

Modifications to the Combination Software
The combination strategy follows the method developed by [36] and described in more details by [20]. The general workflow and can be summarized as follows: 1. A preliminary alignment step, where each AC orbit solution is transferred to a common reference. Historically, for the legacy combination, this alignment was performed with a rotation of the residual values between the Earth Orientation Parameters (EOP) determined by each AC and the ones provided by IERS Bulletin B [20]. Since 2000, this alignment is done directly based on the IGS realization of the International Terrestrial Reference Frame using the transformation parameters provided by the IGN [33]. Those two strategies are equivalent: the underlying idea is to remove potential misalignments at the ACs level. Nevertheless, the latter can be considered as better, since the alignment is done with respect to the IGS reference frame. Unfortunately, this second approach cannot be used, since the transformation parameter are not available yet. Thus, for the current work, this step is done on the EOP approach.
2. A first equally weighted orbit mean is estimated.
3. An Helmert transformation is estimated between each AC solution and the step 2 mean.
4. Weights for each AC are computed based on the formula: where S cent is the number of satellites per AC, E sat cent is the number of positions per AC per satellite. The position of the transformed orbits and the mean orbit is represented by P sat cent (step 3) andP sat (step 2), respectively.
5. The weighted mean is estimated based on: Where N sat cent is the number of Centers submitting a solution for that satellite, and P sat cent is an AC solution realigned on the step 2 mean but using a satellite-weighted Helmert transformation (see [20] for details).
6. Satellite outliers are determined: if the RMS of a satellite computed by an AC exceeds three times the overall satellite RMS mean used as a threshold, then this AC's satellite is excluded for the next iteration. Note that an exclusion does not mean that the satellite will be completely eliminated. It means only that this satellite will not participate in the determination of Helmert transformation parameters and weights, but a combined orbit for an excluded satellite is still provided using the AC weight computed based on the valid satellites. Note also that the software can manage satellite provided by only one AC: if a satellite is included in only one input, it will be included in the combined solution as it is (but still corrected of the pole alignment, step 1).
7. While there is still a satellite exclusion, the combination restarts at step 2.
This method can be considered as reliable, since only few modifications have been performed on it since the mid-1990, and it is still this software which is used on a continuous basis by the ACC.
We started from a modified version of the original ACC combination software designed for a GPS and Galileo orbit combination within the GGSP consortium (Galileo Geodetic Service Provider) [8,34], and performed the necessary modifications to make it capable of handling the new constellations. The improvements carried out are: • A proper handling of the new Galileo, BeiDou and QZSS constellations for the Pole Alignment, Satellite & AC weight determination and the final statistics computations.
• The handling of the SP3-d format, i.e. the most suited one for a multi-GNSS environment, since this version can support an unlimited number of satellites [13]. The previously implemented SP3-c format could only handle 85 satellites, and this limit has been reached for the last weeks of the test period (GFZ products include 86 satellites starting from the day 1947-6, 2017/05/06).
Nevertheless, no modification on the algorithm itself has been made. Thus, no specific weighting is applied for each constellation. Then, internally, all satellites are assumed as belonging to one and the same "big" constellation. We used this strategy instead of a constellationindependent combination, to maintain consistency in the AC weighting for the different constellations, in the continuity of the scheme designed for the GGSP [8].
TUM products containing only orbits for Galileo, QZSS and BeiDou were not used in this study. Indeed, for historical reasons, the combination software necessarily requires GPS data to work properly, and adequate modifications would have been beyond the scope of the enhancements we wished. SHAO products covered a too short time span of only few weeks at the end of our study period (from week 1990 on CDDIS server [28]). Therefore, we did not use them. Lastly, we precise that the ESOC, Darmstadt also generates multi-GNSS products [26], but they are not submitted to the MGEX project, thus they were not used in our study.
Hereafter, we use the RINEX format conventional abbreviations for each GNSS: G for GPS, R for GLONASS, E for Galileo, C for BeiDou, J for QZSS. The different AC products per constellation processed during the test period are represented in Figure 1(a), and the number of satellites used in the combination in Figure 1(b). We note that during the test period, GRGS and JAXA started to provide their products (from week 1826 and 1945 respectively on CDDIS server, GFZ started to process GLONASS and Galileo (from day 1820-4, 2014/11/27) along with QZSS (from day 1846-3, 2015/05/27), and Wuhan University started to process QZSS from day 1826-3, 2015/01/07.   As mentioned above, the first step of the combination is an alignment to a common reference frame. This is done using the EOP approach described in section 1. Nevertheless, some ACs do not provide EOP data in the dedicated .erp format (namely TUM and GRGS, the latter providing their EOP estimation in the SINEX file). Thus, no preliminary alignment is performed and they are assumed to be correctly pre-aligned to the common reference, This step aims to reduce at the maximum the standard deviation of the first mean combination.

Orbit combination results
The RMS differences between the combined solution and the individual ACs solution for all satellites are computed using the formulas described in [20]. They are represented in Figure  2(a). We use a similar representation as the ones usually provided by the ACC: the daily RMS are represented by dots, and a smoothing curve based on a Gaussian filter (with a 7-day window in our case) is used to represent the long-term tendency. We remark that important differences are dominated by WU solution, due to the BeiDou geostationary satellites (hereafter designated as GEO) which degrade their general orbit RMS. The GEO orbit determination is known to be more difficult, since the satellite position variations with respect to the ground tracking network is small, especially on the along-track component [14]. The attitude law for BeiDou-2 GEO (orbit-normal) is also different from the regular yaw-steering law normally used for GNSS [21] and require a proper implementation in the AC processing software. If we exclude the BeiDou GEO and IGSO (inclined geosynchronous satellites) from the RMS difference computation (Figure 2(b)), the RMS difference for all ACs is around 30±15 mm. We decided to exclude IGSO too in Figure 2(b), even if their consistency is one order of magnitude better than the geostationary ones [37], in order to have a MEO-only visualization more consistent with the other constellations.
Figures 2(a) to 2(h) represent the RMS differences for each individual constellation. For GPS, differences of individual AC with respect to the combination is between 8 mm and 18 mm with an average value around 10 mm for the second part of the test period. The differences with the IGS legacy combination (in dark blue for GPS, in dark red for GLONASS) are also shown, with a RMS value below the centimeter level for the second part of the test period for GPS and around 20 mm for GLONASS. RMS differences with respect to the individual AC solutions for GPS is at the centimeter level, which is consistent with the average difference obtained for the legacy combination [15] (between 9 mm and 20 mm for 10 ACs when these lines are written i.e. for a representative period between weeks 1990 and 2055) or for the repro2 campaign combination [10] (between 8 mm and 13 mm for 10 ACs for a representative period between weeks 1250 and 1831). For GLONASS (Figure 2(d)), the differences are between 20 mm and 60 mm with an average around 25 mm for the second part of the test period. We remark that when GRGS starts to provide GLONASS (on week 1826) the CODE solution RMS increase brutally, meaning that a consistent solution between the other ACs is then reinforced with the new inputs. This behavior is corrected around week 1900, when CODE implemented new attitude models in their processing according to IGS technical report [40]. Furthermore, this improvement is also visible for GPS (Figure 2(c)).
Regarding the new constellations, a noticeable improvement of Galileo orbit quality can be remarked (Figure 2 (Figure 2(h)), the RMS difference, although stable, are significantly higher, centered around ∼ 150 mm, but with variation from ∼ 80 mm up to ∼ 420 mm.

SLR residuals
We used Satellite Laser Ranging (SLR) observations to perform an independent quality assessment of the combination. All active Beidou, Galileo, GLONASS, and QZSS satellites are suitable for such analysis since they are equipped with Laser Retroreflector Arrays (LRA, [4]). We used as input the normal points provided by the International Laser Ranging Service [30].
The processing is performed with GFZ's EPOS8 software [39], designed for GNSS precise orbit determination and SLR residual estimation. SLR station coordinates are fixed to SLRF2008 handling for satellites with new kind of orbits, such as geostationary orbit and inclined geosynchronous orbits (some BeiDou satellites and QZSS) or even the two elliptical orbit Galileo satellites. Since it is more challenging to achieve a highly accurate orbit determination for those satellites because of different orbit parameterization and software capabilities, they penalize the weight of the AC which process them, even if the determined orbits for more well-modeled satellites (i.e. circular MEO such as GPS vehicles) are precise enough. Moreover, the necessity of a combined product for those "unconventional" satellites is maybe the most relevant, since the differences between individual solutions appear larger. Lastly, the alignment to a common reference, which is the prerequisite step before the weight determination, need Earth Orientation Parameters data. Unfortunately, some ACs does not provide such products. The preliminary alignment is then skipped for some ACs, which biases the combination. Then, we would like to encourage MGEX Analysis Centers to provide a full set of Multi-GNSS products, including especially the EOP.
Thus, and even if the current IGS software could actually handle the new GNSS constellation and achieve a relevant combination quite easily, it would be beneficial to start the development of a new algorithm. Several ideas for this new strategy are under study: we first developed a workaround method to align the different MGEX solution on a common pole reference [25]. We also started to develop an optimized weighting scheme, considering the different constellations.