1 Introduction

1.1 Role of SLR in space geodesy

Satellite laser ranging (SLR) is a precise space-geodetic technique that provides range measurements to artificial satellites. Due to both the characteristics of geodetic satellites and precise devices for range measurements installed at SLR tracking stations, SLR plays an important role in the realization of the International Terrestrial Reference Frame (ITRF, Altamimi et al. 2016). Orbits of geodetic satellites, such as LAGEOS or Etalon, are determined with sub-centimeter accuracy using range measurements (Sośnica et al. 2014; Appleby et al. 2016). As a result, SLR contributes to the determination of geocenter coordinates, defining thus the origin of ITRF, the global scale, and station coordinates. Apart from the realization of ITRF, SLR provides the most accurate value of the standard gravitational parameter, GM, and low-degree spherical harmonics of the Earth’s gravity field (Thaller et al. 2011; Bloßfeld et al. 2015; Sośnica et al. 2015; Cheng and Ries 2017). The International Laser Ranging Service (ILRS, Pearlman et al. 2002) unifies and coordinates all activities of SLR stations that represent the ground segment of SLR. The ILRS does not only collect, archive, analyze, and distribute SLR and Lunar Laser Ranging data, but also supports different space missions by providing special SLR tracking campaigns and the priority listFootnote 1 including satellites to be tracked by the SLR stations. The priority list contains both passive geodetic satellites and various active spacecraft equipped with Laser Retroreflector Arrays (LRAs) for range measurements.

1.2 SLR tracking of GNSS

Due to the fact that all new active navigation satellites are equipped with LRAs, SLR serves as a validation tool for microwave-based GNSS orbits (Zhu et al. 1997; Appleby et al. 1999; Urschl et al. 2007; Fritsche et al. 2014; Montenbruck et al. 2015a; Steigenberger et al. 2015; Zajdel et al. 2017). The SLR validation performed by Sośnica et al. (2015) was characterized by the root mean square (RMS) at the level of 2.4 and 3.3 cm for GPS and GLONASS satellites, respectively, which coincides with the orbit accuracy declared by the International GNSS Service (IGS, Dow et al. 2009). SLR residuals play a crucial role in the quality evaluation of both operational products of the official multi-GNSS experiment (MGEX, Montenbruck et al. 2017; Prange et al. 2017) and the real-time MGEX products (Kazmierski et al. 2018). Finally, SLR residuals may serve as a validation tool for the empirical models designed for the absorption of solar radiation pressure (SRP) such as the Empirical CODE Orbit Model (ECOM, Arnold et al. 2015) or box-wing models (Rodriguez-Solano et al. 2014).

Range measurements can also be used for an independent determination of satellite orbits. Pavlis (1995) determined GPS orbits solely from SLR observations obtaining the RMS at the level of 7.7, 75.1, and 56.5 cm for GPS-35 and 9.8, 72.9, and 90.9 cm for GPS-36 in the radial, along-track, and cross-track direction, respectively. A joint adjustment of GNSS and SLR observations performed by Urschl et al. (2007) provided an improvement of the determination of the semimajor axis of GNSS orbits. Urschl et al. (2008) calculated the preliminary 9-day orbits of the very first Galileo, i.e., GIOVE-A, using solely SLR measurements and achieved the accuracy at the level of 0.1, 0.5, and 1.0 m in the radial, along-track, and cross-track direction, respectively. Montenbruck et al. (2015b) calculated 14-day orbits using solely SLR data for two NavIC satellites, IRNSS-1A and IRNSS-1B. Due to the poor geometry of SLR observations provided by 8 stations, out of which only 2 were located out of Europe, and an insufficient number of SLR observations, NavIC orbits were determined with the accuracy at the level of 2, 15, and 10 m in the radial, along-track, and cross-track direction, respectively. The combination of GNSS and SLR observations performed by Hackel et al. (2015) resulted in a mitigation of systematic errors in Galileo-IOV solutions. The degradation of the internal consistency between GNSS and SLR combination was solved by adding an offset of 5 cm to Galileo’s LRA which should rather be assigned to mismodeling of the microwave antenna thrust (Steigenberger et al. 2017), albedo (Rodriguez-Solano et al. 2012), and the existence of the satellite signature effect (Otsubo et al. 2001; Sośnica et al. 2015).

1.3 The new GNSS constellations

Due to the emergence of new navigation satellite systems, IGS initiated MGEX. Apart from continuously being modernized constellations of GPS and GLONASS, navigation satellite systems consist now of European Galileo, Chinese BeiDou, Japanese Quasi-Zenith Satellite System (QZSS), and Indian NavIC. GNSS satellites do not orbit only on Medium Earth Orbits (MEOs), but also, as in the case of BeiDou System, on Geostationary Earth Orbits (GEOs) and on Inclined Geosynchronous Orbits (IGSOs). Satellites of new constellations differ with regard to satellite bus characteristics; thus, the equipment for range measurements is also dependent on the type of the satellite. Table 1 summarizes the characteristics of the multi-GNSS constellations in terms of SLR measurements and orbit properties. LRAs mounted on each satellite type are different in shape, size, and the number of corner cubes.

Table 1 Characteristics of the MGEX constellation considered in the analysis

The GLONASS constellation consists of 24 operational MEO satellites with the majority of M-type satellites, and two experimental K-type satellites, out of which only one is operational. GLONASS satellites are equipped with the highest number of corner cubes, i.e., 112 and 123 for GLONASS-M and GLONASS-K, respectively. Retroreflectors mounted on GLONASS-K satellites surround the GNSS transmit antenna and are the only ring-shaped LRAs.

Galileo consists of 18 out of 30 planned satellites decomposed into three MEO orbits. The Galileo segment contains two types of spacecraft, i.e., four in-orbit validation (IOV) spacecraft that were launched after two, inactive now, GIOVE satellites, and 14 fully operational capability (FOC) spacecraft, out of which, two were launched into highly elliptic orbits (Montenbruck et al. 2017). Despite the fact that the satellites in elliptic orbits cannot be used in navigation, they can serve for the investigation of a gravitational redshift (Delva et al. 2015). Galileo satellites are significantly lighter than other satellites (see Table 1); thus, they are more affected by the transmit antenna thrust (Steigenberger et al. 2017) and albedo (Rodriguez-Solano et al. 2012). Galileo-IOV are equipped with arrays of the largest area (\(430.0\times 470.0\) mm) that contain one of the biggest corner cubes (23.3/33.0 mm). As a result, it is easier for stations to get good returns from the early Galileo rather than from Galileo-FOC satellites which are equipped with LRAs of a reduced number (60) of smaller (19.1/28.2 mm) corner cubes. As a result, Galileo-FOC LRAs have an area smaller (\(331.0 \times 248.7\) mm) by a factor more than 2 compared to Galileo-IOV which makes Galileo-FOC one of the hardest satellites to track by SLR stations. The current operational BeiDou-2 constellation consists of six GEO, six IGSO, and three MEO satellites. BeiDou MEO and IGSO carry fewer corner cubes (42) than Galileo-FOC satellites; however, the corner cubes are larger. As a result, the whole array has a bigger surface (316.0/280.0 mm) than Galileo-FOC. BeiDou GEO is equipped with significantly larger LRAs compared to other satellites (490.0/430.0 mm). Larger LRAs are required for geosynchronous satellites due to their high altitude. Japanese QZSS is an augmentation system for GPS and consists of three geosynchronous and one geostationary spacecraft. Similar to BeiDou GEO, all QZSS are equipped with larger LRAs than MEO satellites (400.0/400.0 mm). The other regional navigation system considered by MGEX is the Indian NavIC that consists of 3 GEO satellites and 4 IGSO satellites.

1.4 The goal of this study

This study discusses the results of precise multi-GNSS orbit determination using solely SLR observations as the first step of the combined SLR-GNSS solution. We analyze both the internal and external accuracies of the solution in order to evaluate the influence of SLR observations to GNSS satellites on the multi-GNSS solution and evaluate the role of laser ranging in combination with microwave observations. We concentrate mostly on active MEO GNSS satellites. We also evaluate the utility of SLR-derived multi-GNSS orbits according to the current requirements of Global Geodetic Observing System (GGOS, Plag and Pearlman 2009).

In existing literature on the determination of GNSS orbits using SLR, there is little information about unambiguous strategies in terms of the number and geometry of observations that are sufficient for the precise multi-GNSS orbit determination. Information about the optimal calculation strategy and the best-suited arc length is missing as well. Moreover, SLR stations are capable of tracking more than 50 GNSS satellites (Kirchner and Koidl 2015); thus, this analysis can contribute to the development of tracking strategy, in which multi-GNSS constellations would be tracked homogeneously.

This article aims at answering the following questions:

  • How many SLR observations are necessary to determine precise multi-GNSS orbits using SLR data only?

  • What is an optimal geometry of observations, and thus, how many stations should track navigation satellites to provide a homogeneous coverage with observations for the whole GNSS constellation?

  • What is an optimal length for an orbital arc, i.e., how much can we extend an orbital arc in order to both, gather the largest number of SLR observations to GNSS satellites without degradation of the orbit?

1.5 Structure of the paper

The paper is structured as follows. Section 2 describes the methodology of the solution. Section 3 summarizes GNSS special tracking campaigns held by ILRS in the period of 2014–2017. Section 4 discusses the results. After the investigation of the efficiency of solutions, we assess internal accuracy by the analysis of the mean error of the orbital semimajor axis. Then, we investigate the external accuracy of the orbit solution as a function of the number of SLR observations and the number of tracking stations. Section 5 investigates different arc lengths in order to develop an optimal processing strategy. Section 6 contains comments and summarizes the paper in terms of the improvement of the consistency of solutions based on the microwave (GNSS) and optical (SLR) observations.

2 Methodology

The calculations of multi-GNSS orbits using solely range measurements are performed in the modified version of Bernese GNSS Software 5.2 (Dach et al. 2015). In order to determine a precise orbit, a set of a priori orbit positions and Earth rotation parameters (ERPs) is needed. In our calculations, we use official ERPs and orbits from the Center for Orbit Determination in Europe (CODE, Prange et al. 2017).Footnote 2 In order to increase the consistency of orbit solutions, we use the same SRP model as CODE—the new ECOM2 (Arnold et al. 2015). The solution provided by CODE is based on the same software package and is one of the most accurate for GPS, GLONASS and Galileo constellations within MGEX (Montenbruck et al. 2017) as well as for Galileo satellites launched into highly eccentric orbits (Sośnica et al. 2018). The CODE products, however, do not contain BeiDou GEO satellites. Both SRP models developed by CODE, i.e., ECOM and ECOM2, were designed to absorb solar radiation pressure for satellites that work strictly in the yaw-steering attitude mode, whereas BeiDou GEO is maintained almost continuously in the orbit-normal mode. In our analysis, we process data from the period between 2014.0 and 2016.9, during which two changes in using orbit models occurred. Till the end of 2014, the classical ECOM was used (Beutler et al. 1994; Springer et al. 1999). At the beginning of 2015, the new ECOM2 was proposed. In August 2015, the number of estimated ECOM2 parameters in CODE products was reduced from 9 to 7 by excluding 4-times-per-revolution parameters.

The strategy of the orbit solution is similar to the official 5-system CODE solution for MGEX (Prange et al. 2017), i.e., 3-day orbital arcs are generated, and the solution for a particular day refers to the middle day of the arc. Such an approach greatly stabilizes the solution, especially for new and incomplete GNSS systems, and thus is used not only for MGEX solutions, but also for operational and reprocessed products at CODE (Lutz et al. 2016).

However, during the 3-day period, the number of SLR observations may be insufficient to provide a solution of the best possible quality. An extension of the arc length results in the increase in the number of SLR observations. On the other hand, the solution can suffer from the degradation of both the Keplerian and empirical orbit parameters as the external forces acting upon satellites may change over time. Also, with the extension of the orbital arc the risk of multiple satellite maneuvers occurs which may degrade the solution. Due to that fact, we make an attempt to arrange a strategy for the orbit determination, testing 3-, 5-, 7-, and 9-day-long arcs in order to both maximize the number of SLR observations and not allow orbit parameters to degrade.

We perform two types of final solutions: solution A with the estimation of the full set of parameters, i.e., Keplerian and empirical orbit parameters, SLR station coordinates referred to the SLRF2008,Footnote 3 ERPs, and geocenter coordinates. In solution B, we calculate only a set of Keplerian end empirical orbit parameters, whereas the remaining parameters are fixed to a priori values. The characteristics of both solutions are shown in Table 2.

Table 2 Characteristics of performed solutions: solution A with the estimation of all parameters, solution B with the estimation of orbit parameters only

After the data screening, during which observations exceeding the maximum sigma of residuals at the level of 25 cm are marked as outliers, the set of normal equations is saved. In order to minimize the influence of SLR range biases on estimated parameters, we calculate mean annual range biases for each station-satellite pair. A similar approach was used for LAGEOS satellites by Appleby et al. (2016). The computed range biases are re-substituted to a priori data for further calculations; thus, the systematic errors caused by range biases are diminished and the solution becomes more stable due to a reduced number of estimated parameters. In further calculations, re-substituted range biases are strongly constrained to a priori values. The penultimate step is the validation of SLR core stations by calculating the Helmert transformation parameters between a priori coordinates from SLRF2008 and computed coordinates. Stations with residuals from the Helmert transformation greater than 45 mm in the north, east, or up component are marked and not taken as core stations in further calculations. Coordinates of non-core stations are estimated as free parameters without any constraints imposed thereon. Having provided lists of verified core stations, we proceed to the final parameter estimation in which we test different arc length strategies by stacking 1-day normal equations into 3-, 5-, 7-, and 9-day arcs, after which we compare the middle days of orbital arcs with microwave-based orbits from CODE.

In our calculations, we consider multi-GNSS satellites both tracked by the ILRS network and satellites broadcasting signals on at least two frequencies allowing for the microwave solutions (see Table 1). We calculate orbits for the whole constellation of GLONASS satellites. Although the constellation consists of 24 satellites, 31 GLONASS spacecraft were tracked in total, due to the replacement of several satellites. We calculate orbit solutions for four Galileo-IOV and ten FOC satellites including E14 and E18 which were launched into incorrect orbits. We consider four BeiDou satellites: one MEO and three IGSO satellites. The first spacecraft of the QZSS constellation, i.e., QZS-1, is considered as well.

Fig. 1
figure 1

The number of SLR observations (right axis) to GNSS satellites (left axis: PRN number) registered to active GNSS satellites during 3-day periods between 2014 and 2016. R—denotes GLONASS, E—Galileo, C—BeiDou, J—QZSS. Asterisk (*)—satellites included in the priority list of ILRS for more than 50% of their lifetime

3 The ILRS-intensive tracking campaigns

At the 18th International Workshop on Laser Ranging in Japan in November 2013, ILRS agreed to increase efforts on tracking GNSS constellations and initiated a special study group called Laser Ranging to GNSS s/c Experiment (LARGE). The goal of the LARGE group was to define an operational strategy for tracking GNSS satellites and to improve the consistency between solutions provided by ILRS and IGS. In the frame of the LARGE project, three special tracking campaigns were announced between 2014 and 2017. Figure 1 illustrates the number of SLR observations to GNSS satellites registered for each 3-day orbit solution. The three special tracking campaigns are marked with red boxes.

The first special campaign lasted from August 1 to September 30, 2014. During the campaign, stations were asked to track all GNSS satellites or at least spacecraft from the current priority list of ILRS. The purpose of the pilot campaign was to check the capability of SLR stations to track the whole constellation of GNSS satellites, which was far more than an ordinary procedure at SLR stations. Through the first campaign the most intensively tracked GLONASS satellites were R02, R07, R18, and R21. The number of SLR observations to the Galileo and BeiDou satellites increased only slightly. For QZS-1 no changes in the number of observations were registered. One of the most important outcomes of the campaign was that the intensification of GNSS satellite tracking did not have a negative influence on the tracking of geodetic satellites, including low orbiting spacecraft (Pearlman et al. 2015).

The second campaign was held from November 22, 2014, to the end of February 2015. Instead of tracking all GNSS satellites, the scenario for the second campaign included in the first place six GLONASS satellites, i.e., R02, R07, R12 R17, R18, and R20. The tracking of Galileo-IOV and BeiDou satellites was given the second priority. The number of SLR normal points increased by 107, 154, and 107% for GLONASS R02, R12, and R18, respectively, as compared to the similar period of 100 days of the non-campaign period, i.e., between the end of April 2014 and August 1, 2014. Apart from Galileo E19, for which the number of observations increased only by 7%, the increase of SLR observations was neither registered for the other Galileo satellites, nor for the BeiDou MEO satellites. Instead, stations were focused on tracking GLONASS satellites included in the third priority list: R03, R04, R08, and R09.

The third GNSS tracking campaign was held between August 10 and October 16, 2015, and involved six GLONASS satellites: R03, R07, R09, R12, R20, and R21; Galileo-IOV and BeiDou C11. The number of observations collected in the frame of the 67-day campaign was higher by 104, 51, and 81% for GLONASS R03, R07, and R12, respectively, when compared to a similar period of 67 days before the 3rd tracking campaign. This time, the number of observations increased also for Galileo-IOV by 94, 93, and 61% for E11, E12, and E19, respectively. At the end of the third campaign, a significant increase in the number of observations at the level of 193% was registered also for BeiDou MEO C11.

Fig. 2
figure 2

The efficiency of multi-GNSS orbit solution using SLR data. The color shows a period for which a particular spacecraft was included into the priority list of ILRS. Green color means that the satellite was included in the priority list for the whole analysis period. The blue color stands for the satellites that were in the priority list for more than a half of the analysis period, whereas the yellow color implies that the satellite was included in the priority list for less than a half of the analysis period or was not included at all

4 Results

4.1 The effectiveness of multi-GNSS orbit solutions using SLR

The efficiency of the orbit solutions is evaluated as a ratio between the number of successful SLR solutions to the number of determined orbit solutions using microwave observations at CODE for the same period.

The median efficiency of the orbit solution equals 87, 86, 79, 63, and 39% for GLONASS, Galileo, BeiDou MEO, BeiDou IGSO, and QZSS, respectively. The effectiveness of GLONASS, Galileo and BeiDou MEO owes to the fact that those constellations consist only of global coverage satellites which provides an even geometry of SLR observations, whereas BeiDou IGSO and QZSS suffer from a poor geometry of observations due to their regional attitude. We distinguish satellites that were included in the priority list of ILRS during the period of analysis (see Fig. 2). Only one GLONASS satellite (R07) was included in the ILRS priority list during the whole period between 2014 and 2016. For more than a half of the GLONASS constellation, it was possible to provide precise orbit solutions for more than 80% of cases, even if the satellites were not included in the priority list. In total, 10 spacecraft of the Russian constellation were included in the priority list only for a short period or were not included at all. As a result, for some of those satellites it was impossible to determine a reliable orbit solution for about 42% cases (see Fig. 2). These satellites were tracked only by the European stations; thus, their observational geometry was inferior. Galileo satellites are placed in the priority list right after their launch. As a consequence, orbit solutions could be provided for more than 80% cases for almost all Galileo spacecraft. BeiDou and QZSS satellites were included in the priority list for the whole period of the analysis; however, all of these satellites (apart from C11) are regional, geosynchronous spacecraft. As a result, due to the insufficient number of SLR observations, the efficiency of solution is rather poor (see Sect. 4.3).

4.2 Impact of the observation number on the internal orbit quality

The impact of the number of SLR observations on the internal quality of the 3-day orbit solutions is evaluated on the basis of formal errors of the semimajor axis. Such an analysis provides information on whether the observation geometry from the normal equation systems is sufficient for obtaining high-quality orbits. The formal error of the semimajor axis \(m_{a}\) of the j-th satellite is extracted from diagonal elements of the covariance matrix:

$$\begin{aligned} m_{a,j}=m_0\sqrt{(\hbox {diag}(A^{\mathrm{T}}A)^{-1})}_{i}\ \end{aligned}$$
(1)

where A is the design matrix describing the observation geometry, and \(m_{0}\) is the a posteriori variance of the unit weight, i.e., mean error of the solution, i is the parameter index in the normal equation system corresponding to the semimajor axis of j-th satellite. Figure 3 illustrates a set of semimajor axis formal errors for particular MEO satellites for solutions A and B (see Table 2). The red line represents solution B with the estimation of orbit parameters only. The result of the solution A is slightly better, i.e., \(m_{0}\) assumes very small values. However, if considering all parameters, the number of estimated parameters approaches to the number of observations. As a result, a set of normal equations, for several satellites, is close to singular and the value of the a posteriori variance of unit weight \(m_{0}\) is underestimated. In solution A, the least squares method fits well to the observations through a large number of estimated parameters, but the solution is not as reliable as when estimating orbit parameters only, where the number of observations significantly exceeds the number of parameters. In the following analysis, we concentrate mostly on the results from the solution B which is non-singular, more stable, and thus more reliable.

Fig. 3
figure 3

Dependency of the number of SLR observations on formal errors of semimajor axes of multi-GNSS satellite orbits obtained using the 3-day solution. Solution A with the estimation of all parameters is shown in blue; solution B with the estimation of orbit parameters only is shown in red. Individual solutions are marked as faint dots, and median values are marked as large dots

The formal error of orbit semimajor axis converges quickly in the range between 13 and 50 observations (Fig. 3). In order to determine GNSS orbits with the accuracy of the semimajor axis at the level of 6 cm, we need approximately 50 SLR observations. This can be read from Fig. 3 in the case of all MEO satellites. The formal error of the semimajor axis when using 60 SLR observations equals 5.0, 4.6, 4.2, and 5.9 cm for R18, E19, E30, and C11 satellites, respectively.

4.3 Comparison to the microwave-based orbits

The external accuracy of the 3-day orbit solution is evaluated based on the analysis of differences between orbits determined using solely SLR observations and orbits from official CODE MGEX orbit products. We compare positions of satellites in 15-min intervals decomposed into three directions: R—radial—pointing from the center of the Earth to the satellite, W—cross-track—normal to the orbital plane, and S—along-track—perpendicular to other axes and approximating the direction of the satellite velocity vector. Then, the mean offset and the RMS of differences are calculated. All analyses provided in this section are calculated for the 3-day orbit solution.

Fig. 4
figure 4

A detailed 3-day solution provided from a high number of SLR observations: 110 (top), 54 (middle), and 25 (bottom) for GLONASS R18. RMS is given for the middle day (time in UTC) of the solution

4.3.1 Dependency on the number of observations

Figure 4 shows three examples of the comparison between SLR-based and microwave-based orbits of GLONASS R18 for three selected days in 2015. When the number of SLR observations is 54 (Fig. 4, middle), the differences are at the centimeter level for the radial and along-track components and at the decimeter level for cross-track. When doubling the number of SLR observations, the RMS in the radial and along-track directions does not improve significantly (Fig. 4, top). However, the cross-track component improves by the factor of 4. When the number of SLR observations is just 25, the solutions becomes unstable with differences up to several meters (Fig. 4, bottom). The along-track component exhibits a secular drift, whereas all orbit components show large periodic variations of the period corresponding to the satellite revolution and the second harmonic of the satellite revolution period.

The accuracy of an orbit based on 25 SLR observations is insufficient for geodetic purposes. However, such an orbit can be sufficient for the determination of space debris (Cordelli et al. 2016), including inactive GLONASS satellites, all of which are equipped with retroreflectors for SLR measurements. A good observational geometry, i.e., including observations collected by stations from different continents, can improve the orbit quality of inactive satellites even in a solution based on just 25 SLR observations.

Now, we provide a statistical analysis of the differences between microwave and SLR orbits for the 3-year period by estimating median values. Figures 5 and 6 show the dependency of the number of SLR observations on the RMS of orbit differences and the dependency of the number of SLR stations on the RMS values, respectively.

Fig. 5
figure 5

Dependency of the number of SLR observations on the RMS of differences between microwave and 3-day SLR solutions. The vertical scale for BeiDou IGSO (bottom-middle) and QZSS (bottom-right) is changed, due to the higher RMS values. Raw data are shown as faint dots. Median values (by five) are shown as solid lines

In all cases the radial component is characterized by the smallest value of RMS both for SLR and microwave solutions (Montenbruck et al. 2017). Due to the fact that the radial component is directly measured by SLR and is strongly subject to the a priori applied force models, it is calculated in the most reliable way, even for a relatively small number of SLR observations and a relatively small number of SLR tracking stations. In microwave solutions, the RMS of both cross-track and along-track components assumes similar values. However, the cross-track component determined using SLR is characterized by a significantly higher value of RMS, similarly to the results based on simulated data reported by Hugentobler (2017) for Galileo SLR-derived orbits. As in the case of the semimajor axis (Fig. 3), the more SLR observations contributing to the solution, the smaller the value of RMS is. In most cases, 60 SLR observations are sufficient to provide reliable GNSS orbits. For the number of 60 observations, the RMS in the radial direction equals 2.9, 2.8, 4.5, and 2.0 cm for R07, R18, E19, and C11, respectively (Fig. 5). However, the cross-track component determined from 60 observations exhibits an RMS at the level of about 20 cm with an exception of Galileo E18, for which the RMS in the along-track component reaches 10.1 cm. The solution for E18, on the other hand, is not as stable as for the other MEO spacecraft. RMS of the 3D position equals 19.3, 22.1, 22.4, and 17.6 cm, for R07, R18, E19, and C11, respectively. Characteristics of RMS values of the solutions calculated using 60 SLR observations for all satellites and constellations are shown in Table 3.

Fig. 6
figure 6

Dependency of the number of tracking stations on the RMS of differences between 3-day SLR and 3-day microwave orbit solutions decomposed into three directions: radial (red), along-track (blue), and cross-track (green). The vertical scale for BeiDou IGSO (bottom-middle) and QZSS (bottom-right) is changed, due to the higher values of RMS

The number of 100 observations collected during the 3-day period was exceeded only for several GLONASS satellites and for Galileo-IOV E19. If the number of observations exceeds the level of about 100, the RMS of the radial component does not further improve as it already equals 2.3 cm for all GLONASS satellites. However, both the along-track and cross-track components are characterized by a much lower value of RMS when the number of SLR observations reaches 100, i.e., 6.8 and 10.9 cm in the along-track and cross-track direction, respectively, which results in the 3D RMS at the level of 14.3 cm. As a result, the 3D position of satellites illustrated in Fig. 5 improves to the level of 13.3, 13.6, 15.6, and 15.1 cm for R07, R18, E19, and C11, respectively.

Table 3 Summary of RMS between microwave and SLR orbits obtained using 60 SLR observations for 3-day arcs (all values are given in cm)

More than 100 observations to one spacecraft do not significantly improve the solution. In terms of the orbit determination using SLR, stations should focus on providing a constant number of observations to the whole multi-GNSS constellation rather than focus on particular satellites.

The distribution of RMS of differences for Galileo-FOC E30 is slightly better (2.9, 7.5, and 20.7 cm in radial, along-track, and cross-track, respectively) than for GLONASS R10 (6.2, 11.1, and 20.5 cm in the radial, along-track, and cross-track, respectively), despite a similar number of observations. However, R10 is simultaneously tracked only by 5 SLR stations, whereas E30 is tracked by more than 10 stations. The high number of stations itself does not provide a recipe for an accurate orbit solution. Stations that provide range measurements have to be homogeneously and globally distributed. High quality of Galileo E30 orbit results from the fact that this satellite is tracked by worldwide distributed stations, which provide a proper geometry of observations.

Orbit solutions vary for different types of BeiDou satellites. BeiDou MEO switches from the yaw-steering mode to normal mode when the \(\left| \beta \right| \) angle, i.e., the elevation angle of the Sun above the satellite’s orbital plane, is below \(4^{\circ }\) (Montenbruck et al. 2015c). Despite the fact that the normal mode is not modeled correctly, the RMS of differences for BeiDou MEO is similar to GLONASS and Galileo constellations. The regional attitude of IGSO satellites significantly limits the number of observations and weakens the geometry of range measurements. The reference microwave orbits for geosynchronous satellites are not of the highest quality; as well, thus further investigations in terms of the geosynchronous orbit parameter estimation have to be performed.

4.3.2 Dependency on the number of SLR tracking stations

The quality of the orbit solution depends also on the number of SLR tracking stations. Our analysis confirms that the accuracy of orbit solution increases with the growth of the number of SLR stations as reported by Hugentobler (2016) who used simulated SLR data. According to Fig. 6, the solution based on fewer than 5 SLR stations provides orbits of a poor quality with the median RMS of differences at the level of 7.8, 23.0, and 48.6 cm in the radial, along-track, and cross-track, respectively, for all satellites considered in our calculations. GLONASS R10 was tracked by up to 5 SLR stations—all of which are located in Europe, which does not allow for obtaining high-quality orbits (see Fig. 6). BeiDou IGSO C08 is tracked at maximum by 6 stations (Fig. 6), whereas the QZS-1 is simultaneously tracked by only up to 4 stations and on average by 1–2 stations. Moreover, the geometry of observations is barely changeable in time. With the increase in the number of stations from 5 to 10, the value of the RMS starts to decline and reaches 3.6, 9.4, and 18.4 cm in the radial, along-track, and cross-track direction, respectively, which is by the factor of 2 better compared to the case of employing 5 SLR stations. For more than 10 SLR stations, the orbit solution stabilizes. In the case of 12 tracking stations, the accuracy of the orbit solution reaches the level of 2.9, 8.4, and 14.9 cm in the radial, along-track, and cross-track component, respectively.

The high number of SLR stations itself does not provide a reliable orbit without a sufficient observation geometry, which is a reason for the good quality of the orbit solutions for R07, R18, and E19. The median value of the RMS for GLONASS R07 tracked by 14 stations equals 2.6, 5.2, and 9.9 cm in the radial, along-track, and cross-track direction, respectively (see Fig. 6). GLONASS R18 is characterized by even smaller values of RMS, i.e., 1.7, 3.8, and 6.8 cm in the radial, along-track, and cross-track direction, respectively, while being tracked by 15 SLR stations. For Galileo-IOV E19 the best solution is provided by 14 stations and equals 2.8, 8.2, and 13.1 cm in the radial, along-track, and cross-track direction, respectively. Unexpectedly, for Galileo-FOC E30 the lowest values of RMS are obtained from the solution based on 10 tracking stations and equals 2.3, 5.7, and 6.3 cm for the radial, along-track, and cross-track direction, respectively. However, the median values are obtained from only one solution; thus, it should not be considered as representative without further investigations.

In conclusions, at least 60 SLR observations collected by 10 different SLR stations are needed to obtain 3-day SLR orbits with a quality of about 3, 9, and 18 cm for the radial, along-track, and cross-track, respectively. With the increase in the number of SLR observations and the number of stations, the solution improves and stabilizes when more than 10 homogeneously distributed stations provide at least 100 observations.

5 The impact of the arc length

5.1 Orbit overlaps

We check the internal consistency of different arc length strategies by the analysis of the orbit overlaps which are calculated as an RMS between the middle day of the arc and the corresponding day of the consecutive arc. For the longer arcs, the internal accuracy of the determined orbits increases, because for the 3-day arcs the middle (second) day is compared to the first day from the consecutive arc, whereas for the 5-day arcs the middle (third) day is compared to the second day from the consecutive arc.

For satellites for which the efficiency of solutions was higher than 80% (see Fig. 2), orbit overlaps are at the cm level in the radial direction for the 3-day solutions. With the extension of the orbital arc, the median RMS of orbit overlaps decreases by 58, 73, and 81% for the 5-, 7-, and 9-day arcs, respectively. Figure 7 shows orbit overlaps for the 5-day solutions for which the median RMS equals 4.3, 14.8, and 32.8 cm in the radial, along-track, and cross-track direction, respectively, which is better by 58% than for the 3-day solutions but is worse by 36% than for the 7-day solutions, for which almost 5% of solutions could not be calculated due to the satellites maneuvers. For the most intensively tracked satellites, i.e., GLONASS R07 the median RMS of orbit overlaps for the 5-day solution reached the level of 2.1, 6.6, and 10.8 cm in the radial, along-track, and cross-track direction, respectively. The internal accuracy for BeiDou IGSO and QZSS is of poor quality and does not significantly improve with the extension of the orbital arc. In the best case, the median RMS of the orbit overlaps for the 5-day solution equals 19, 66, and 167 cm in the radial, along-track, and cross-track direction, respectively, for BeiDou IGSO C08. As a result, the presented methodology is better suited for MEO satellites, rather than for BeiDou IGSO and QZSS.

Fig. 7
figure 7

Median RMS of the orbit overlaps calculated for the middle day of the 5-day solution B

5.2 The distribution of solutions for selected satellites

Table 4 shows the contribution of solutions with the accuracy better than 10 cm for three MEO satellites, i.e., GLONASS R02, Galileo-IOV E19, and BeiDou MEO C11. The effectiveness of the accuracy below 10 cm in the radial direction reaches 90 and 92% for 5- and 7-day solutions in the case of GLONASS R02, whereas for 3-day arcs equals only 71%. With the extension of an orbital arc, the quality of the along-track and cross-track components improves dramatically. In the case of the cross-track component, the contribution of solutions with the accuracy better than 10 cm is higher by a factor of 3 for 5- and 7-day arcs in contrary to 3-day solutions (see Table 4). However, for the 9-day solution the quality decreases even by 5% as compared to the 7-day solution A. The GLONASS constellation is the only one, for which the solution A is superior to the solution B which may be caused by the fact that GLONASS satellites are tracked by new Russian stations which do not have good a priori coordinates in SLRF2008.

Table 4 Percentage of solutions with the accuracy higher than 10 cm in the radial, along-track, and cross-track component for R02, E19, and C11 for different processing strategies. Values for solution A and B are shown

Galileo E19 is tracked by almost the same set of SLR stations as R02; however, the distribution of solutions is slightly worse as compared to the GLONASS satellites. The number of solutions below the level of 10 cm in the radial direction is higher by 16, 17 and 17% for 5-, 7-, and 9-day arcs, respectively, as compared to the 3-day solutions. In contrary to the radial component, for the cross-track component, 5-day arcs constitute the best solution (see Table 4).

BeiDou C11 is the least intensively tracked satellite in the collection (Table 4). The geometry of observations is the weakest as well, as it is tracked intensively by the European and Australian stations and Changchun. As a result, the solution improves with the extension of the orbital arc. The number of solutions below the level of 5 cm in the radial direction is greater by 20% when compared 7-day solutions to the 3-day solutions.

In conclusion, for the most intensively tracked satellites, e.g., R02 and E19, 5- or 7-day arcs constitute the optimal solution. For the 9-day solution the degradation in the along-track component can be noticed (see Table 4). For R02, the increase in solutions with accuracy below the level of 10 cm in the along-track component is more spectacular (from 38 to 57%) than for E19 (from 18 to 27%) when extending the arc from 3 to 5 days. For satellites that are tracked less intensively (mostly by the European stations), we conclude that 7- and 9-day-long arcs constitute the best solution. However, for the most intensively tracked satellites the 9-day solution causes a degradation in the along-track component. The cross-track component improves, not much enough though, to compensate for the degradation of other components. Moreover, one has to remember that the longer the arc is, the bigger the risk occurs that satellites maneuvers take place, due to which, in our case, almost 10% of solutions could not be obtained for the 9-day arcs.

5.3 Strategy development for all satellites

Figure 8 illustrates median values of the RMS of differences between microwave and SLR orbits calculated for all solutions for the whole multi-GNSS constellation in the analyzed period. Solution A with the estimation of all parameters appears to be more reliable than solution B for the whole GLONASS constellation. On the other hand, the median for GLONASS satellites contains the results for satellites that are not tracked intensively enough to determine a reliable orbit. The RMS for 7-day arcs calculated for intensively tracked spacecraft equals: 3.3, 9.0, and 14.3 cm in the radial, along-track, and cross-track component, respectively. Such satellites constitute only a half of the GLONASS constellation; thus, the median calculated for the whole constellation is not as reliable as in the case of e.g., Galileo.

Fig. 8
figure 8

The median values of the RMS of differences between microwave and SLR orbits calculated for the whole constellation of GLONASS, Galileo, BeiDou MEO, BeiDou IGSO, and QZSS. Bars contain values of the RMS decomposed into three directions, i.e., radial, along-track, and cross-track, and 3D, respectively (from left to right). Solution A corresponds to the solution with the estimation of the parameters mentioned in Table 2. Solution B considers only Keplerian and empirical orbit parameters

Table 5 Median RMS for 5- and 7-day SLR solutions compared to the internal consistency of all MGEX Analysis Centers reported by Montenbruck et al. (2017)

Due to the fact that all Galileo satellites are tracked evenly, the European constellation seems to be the most reliable one for the development of the optimal arc length strategy. According to Fig. 8, the most reliable length of the arc is 7-day. However, as in the case of 9-day solutions, the satellites maneuvers still interfere with the process of orbit determination. As a result, given the compromise between the quality of the solution and the efficiency of calculations we consider the 5-day solution B as the most reliable. The accuracy for this solution is at the level of 4.8, 15.1, and 30.9 cm in the radial, along-track, and cross-track direction, respectively.

BeiDou MEO segment is represented only by one spacecraft, i.e., C11. This satellite is not tracked as intensively as the other MEO satellites; thus, the longer the arc is the more observations supply the solution (see Fig. 8). A similar situation occurs for the BeiDou IGSO constellation. The extension of the orbital arc does not, however, improve the solution for QZSS.

Table 5 gives the median results from the 5- and 7-day orbit determination using SLR with the consistency of all Analysis Centers (AC) of MGEX provided by Montenbruck et al. (2017) which are based on microwave GNSS observations. GLONASS satellites are divided into two groups, i.e., all satellites from the Russian constellation and satellites included in the priority list of ILRS. The RMS of SLR solutions for intensively tracked GLONASS satellites (i.e., 4 cm) is only slightly higher than the bottom threshold of consistency of MGEX ACs for the radial component. A similar situation occurs for all MEO satellites, especially for Galileo and BeiDou MEO whose accuracy of the determination of the radial and along-track components of about 4 and 15 cm respectively, fits well the orbit accuracy provided by MGEX ACs. Due to the poor SLR observation geometry, the solution for BeiDou IGSO and QZSS is significantly worse than the MGEX performance. Table 5 contains also the best solutions obtained for GLONASS, Galileo, and BeiDou MEO satellites for the 5- and 7-day solutions. In all cases, the RMS is slightly lower, albeit in the threshold of the consistency of MGEX solutions.

The most consistent solution was provided for the GLONASS R07 on July 15, 2015. The 5-day solution was calculated using 129 SLR observations provided by 12 homogeneously distributed SLR stations. During the 5-day period R07 was tracked by 2 stations from North America, 3 from Asia, 2 from Australia, and 5 from Europe all of which provided an even and sufficient geometry of observations. The RMS of the best solution for GLONASS R07 equals: 0.8, 2.4, and 1.9 cm in the radial, along-track, and cross-track direction, respectively (see Table 5).

6 Discussion and conclusions

Thanks to the great effort of ILRS and SLR tracking stations, the number of SLR observations to multi-GNSS has grown significantly since late 2013. This enables us to determine GNSS orbits based on SLR observations in 87, 86, 46, and 39% cases for GLONASS, Galileo, BeiDou, and QZSS, respectively. Based on the orbit solutions, we performed analyses which answered the questions: what is the sufficient number of SLR observations and what is the optimal geometry of SLR observations in order to provide an orbit solution at the cm level. We also investigated an optimal arc length in order to both increase the considered number of observations and to avoid the orbit parameters to degrade.

The internal consistency of the orbit solutions based solely on SLR observations to multi-GNSS satellites is at a satisfying level, i.e., the median RMS of orbit overlaps reached the level of 4, 15, and 33 cm in the radial, along-track, and cross-track direction, respectively. In order to determine multi-GNSS orbits with the mean error of the semimajor axis at the level of 6 cm, one needs about 50 SLR observations. As in the case of the external accuracy, with 20–30 observations one may obtain an orbit of a meter accuracy which can be sufficient for inactive navigation satellites equipped with LRAs; thus, range measurements to GNSS satellites may contribute to tracking of space debris.

On average, 60 observations are sufficient to provide a precise multi-GNSS orbits of a scientifically useful quality at the level of: 3, 11, and 17 cm in the radial, along-track, and cross-track direction, respectively, for all MEO satellites, whereas the quality of the solution for BeiDou IGSO was worse by the factor of 4 and for QZSS the number of 60 observations has never been reached. If the number of observations increases to 100, the RMS of both the along-track and cross-track components decreases to the cm level. More than 100 observations improve the solution only marginally such that it is not worth putting effort on tracking a small subset of particular satellites in an attempt to increase the data beyond this level. It would be more efficient to track the whole MGEX constellation homogeneously.

The number of SLR stations that provide SLR observations is not without significance. More than 5 stations are necessary to provide a decent geometry of observations. In total, 10 SLR stations provide a sufficient geometry to determine multi-GNSS orbits of a useful quality. Stations that provide observations have to be, however, homogeneously distributed in order to provide a sufficient geometry of observations. BeiDou IGSO and QZSS satellites are observed by fewer than 6 stations, which is one of the reasons why we cannot calculate high-quality orbits for those satellites due to a poor geometry of observations. Moreover, the normal attitude mode, which is used by various satellites during parts of the eclipse periods, causes problems with orbit modeling for both the microwave and SLR solutions.

We also tested the impact of the orbital arc length on the orbit solution by testing 3-, 5-, 7-, and 9-day arcs. For the intensively tracked satellites 5-day arcs constitute the best solutions with the median RMS at the level of 4.8, 15.1, and 28.2 cm and 7-day arcs with the median RMS at the level of 4.1, 13.4, and 22.9 cm in the radial, along-track, and cross-track direction, respectively. As in the case of the 9-day solution, the degradation of the accuracy of the along-track component appears. Moreover, satellite maneuvers that occur in the 9-day period made it impossible to calculate 10% of solutions, whereas in the case of the 7-day arcs maneuvers deteriorated fewer than 5% of solutions.

Fig. 9
figure 9

Multi-GNSS orbit determination using SLR as a function of the number of observations, \(\beta \) angle (angular height of the Sun above the orbital plane), and the Empirical CODE Orbit Model (ECOM1/ECOM2). The number of observations and RMS values are presented as a 14-day medians. Gray areas indicate special GNSS tracking campaigns

Figure 9 summarizes the diversity of the achieved orbit quality. The accuracy of SLR-derived orbits depends on: (1) the number of SLR observations, (2) the number and distribution of SLR tracking stations, (3) the length of the orbital arc, (4) the generation of the satellite, (5) type of orbital plane, (6) \(\beta \) angle, (7) empirical models used in the calculation (ECOM1/ ECOM2).

The analysis performed in this paper provides information about the required number of SLR observations and their optimal geometry; thus, this work can serve as a recommendation for SLR stations for collecting SLR observations in terms of precise multi-GNSS orbit determination using SLR data.

Firstly, one should investigate the accuracy and capabilities of the orbit determination using solely range measurements in order to use SLR-derived multi-GNSS orbits in studies, such as the determination of global geodetic parameters. The quality of the orbit solution is at the satisfying level, with the best case of the 5-day solution for the GLONASS R07 which was calculated using 129 observations delivered by 12 homogeneously distributed SLR stations. RMS of the solution reached the level of 0.8, 2.4, and 1.9 cm in the radial, along-track, and cross-track direction, respectively, which fulfills the criteria of the precise orbit.

Multi-GNSS orbits derived using solely SLR data serve as an independent orbit solution; thus, they can contribute as an evaluation tool for the orbit modeling problems indicated by microwave solutions. Figure 9 validates the correctness of the new ECOM2 model for Galileo satellites. The influence of low \(\beta \) angle does not affect Galileo satellites to such an extent during the use of ECOM2. However, the methodology presented in this paper performs well for MEO satellites but still requires further investigation in the case of geosynchronous satellites. After that, the methodology can be used for BeiDou-3 satellites, i.e., C31, C32 (IGSO) and C34, C33 (MEO) that are not currently considered by ACs, and for the determination of Galileo-IOV E20 that broadcasts microwave signal on just one frequency since 2014, and thus, the precise MGEX orbits are not available. Finally, multi-GNSS orbit determination using solely SLR data is the first step for a combined microwave-SLR orbit solution which is crucial for the SLR-GNSS co-location in space onboard GNSS satellites independent from ground-based local ties.