Determination of SLR station coordinates based on LEO, LARES, LAGEOS, and Galileo satellites

The number of satellites equipped with retroreflectors dedicated to Satellite Laser Ranging (SLR) increases simultaneously with the development and invention of the spherical geodetic satellites, low Earth orbiters (LEOs), Galileo and other components of the Global Navigational Satellite System (GNSS). SLR and GNSS techniques onboard LEO and GNSS satellites create the possibility of widening the use of SLR observations for deriving SLR station coordinates, which up to now have been typically based on spherical geodetic satellites. We determine SLR station coordinates based on integrated SLR observations to LEOs, spherical geodetic, and GNSS satellites orbiting the Earth at different altitudes, from 330 to 26,210 km. The combination of eight LEOs, LAGEOS-1/2, LARES, and 13 Galileo satellites increased the number of 7-day SLR solutions from 10–20% to even 50%. We discuss the issues of handling of range biases in multi-satellite combinations and the proper solution constraining and weighting. Weighted combination is characterized by a reduction of formal error medians of estimated station coordinates up to 50%, and the reduction of station coordinate residuals. The combination of all satellites with optimum weighting increases the consistency of station coordinates in terms of interquartile ranges by 10% of horizontal components for non-core stations w.r.t LAGEOS-only solutions.


Introduction
The International Association of Geodesy (IAG) established services, which by using space geodetic techniques, provide valuable infrastructure for observations and products, used for description and understanding of the dynamic Earth system (Plag and Pearlman 2009). Four IAG services and their basic observation techniques are used for deriving global geodetic parameters and the realization of the International Terrestrial Reference Frames (ITRF; Altamimi et al. 2016), including Satellite Laser Ranging (SLR; Pearlman et al. 2019a) and Lunar Laser Ranging (LLR, used in previous ITRF realizations; Hofmann et al. 2018), Very Long Baseline Interferometry (VLBI; Schlüter and Behrend 2007), Global Navigation Satellite Systems (GNSS; Johnston et al. 2017), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS; Willis et al. 2010). The IAG founded the Global Geodetic Observing System (GGOS; Plag and Pearlman 2009;Gross et al. 2009) to provide the geodetic expertise and infrastructure necessary for the monitoring of the Earth system and global change research by integrating geodetic, geodynamic and geophysical communities. The GGOS requires 0.1 mm/year stable-in-time station velocities and the accuracy of 1 mm for the station coordinates. The realization of the ITRF requires the position of geocenter (GCC), scale, and Earth orientation parameters determined with the accuracy of 1-mm or better to appropriately realize the origin, metrics, and the orientation of the reference frame (Gross et al. 2009).
These requirements cannot be fulfilled without a network of high-quality observing, co-located GGOS sites, which provide the necessary infrastructure to detect and reduce sources of systematic errors of observation systems, and to fully exploit the current constellation of satellite missions. Nowadays, the combination of space geodetic techniques is based on co-located GGOS sites equipped with different techniques, however, the future co-location onboard active satellites is foreseen, as well. Co-location of space geodetic techniques onboard GNSS satellites and active low Earth orbiters (LEOs) pave the way for new opportunities to future inter-technique connections, for the ITRF realizations through the integration of various types of geodetic observation systems and their products.

SLR observations to spherical geodetic, LEO, and GNSS satellites
The SLR technique employs optical range measurements to various types of satellites equipped with laser retroreflector arrays (LRA). The main contribution of the SLR technique to GGOS originates from measurements to spherical geodetic satellites, such as LAGEOS-1/2, Etalon-1/2, which are used for the determination of GCC, global scale, station coordinates (Thaller et al. 2011;Glaser et al. 2019;Zajdel et al. 2019), standard gravitational parameter GM, gravity field parameters (Bloßfeld et al. 2015;Sośnica et al. 2015a), and for the verification of general relativity effects (Ciufolini et al. 1998). The distribution of all normal points (NPs) collected by SLR stations in 2016-2018 reveals that in total 81% of all SLR observations were collected to LEOs, over 10% of the observations to GNSS satellites, and the remaining observations comprising less than 9% of all SLR observations to spherical geodetic satellites 1 . Only these latter were used so far for deriving global geodetic parameters, and for the computation of geodetic reference frames using the SLR technique. The spacecraft and orbit characteristics of different satellite types with the SLR retroreflectors are summarized in Table 1 and Fig. 1.
SLR spherical geodetic satellites, such as LAGEOS-1/2 and LARES, are covered with corner cube retroreflectors, with reduced area-to-mass ratios (Pearlman et al. 2019b). Those satellites are used to estimate geodetic parameters because out of all artificial satellites, their orbits are not very much affected by non-gravitational perturbations. The mission main objectives include determining the terrestrial reference frame and Earth rotation parameters (ERPs) and improving the gravity field models. The main objective of the LARES satellite are the measurement of relativistic effects, including the Lense-Thirring effect, space geodesy and geodynamics (Ciufolini et al. 2015). SLR measurements to active LEOs, such as GRACE-A/B (Tapley et al. 2004), Sentinel-3A/B (Donlon et al. 2012), and Jason-2 (Lambin et al. 2010), or GNSS constellations are typically employed for the validation of microwave-based or DORIS-based orbits of active LEOs, and microwave-based orbits of GNSS satellites. The number of active LEO satellites, which were equipped with at least two different positioning techniques, i.e., SLR+GNSS, SLR+GNSS+DORIS, or SLR+DORIS has been increasing for the last 20 years. The co-location of SLR and GNSS can be found on satellites, such as GOCE (Drinkwater et al. 2007), GRACE-A/B, CHAMP (Reigber et al. 2000), SWARM-A/B/C (Friis-Christensen et al. 2008), TerraSAR-X (Buckreuss et al. 2003), TanDEM-X (Krieger et al. 2007), Sentinel-3A/B, Sentinel-6/Jason-CS, and Jason-1/2/3. Depending on the mission requirements and objectives, the LEOs have different shapes, constructions, orbit characteristics, and specific scientific instruments or payload onboard, such as GNSS receivers, accelerometers, communication devices, solar panels, propulsion systems, and LRAs. The LRAs onboard LEOs have either four cubes arranged on a regular 45° pyramid-shaped array, seven cubes arranged on a spherical surface, or eight corner cubes symmetrically mounted on a hemispherical surface with one nadir-looking corner cube in the center (Table 1; Arnold et al. 2019). To achieve the mission objectives, the requirements for LEOs orbit accuracy reaches the level of 1-2 cm for altimetry and gravity field recovery missions (e.g., Sentinel-3, GRACE), thus the determined orbits must be of the superior accuracy.
The worldwide GNSS constellations are: the American GPS, the Russian GLONASS, the European Galileo, and the Chinese BeiDou. The Galileo system is expected to complete constellation deployment in 2022. In terms of GNSS satellites, the SLR observations already proved to be useful for the determination of the GNSS orbits using combined GNSS and SLR observations (Bury et al. 2021). The International GNSS Service (IGS; Johnston et al. 2017) initiated the Multi-GNSS Experiment (MGEX; Montenbruck et al. 2017) to track, collate, and analyze all available GNSS signals for providing highprecision products accessible to a wide range of engineering and science communities. All active Galileo satellites are equipped with LRA for SLR. Between 2014 and 2020 the International Laser Ranging Service (ILRS; Pearlman et al. 2019a) initiated several intensive tracking campaigns of GNSS in order to improve the scheduling, increase the data volume especially in the daytime, and to find the best possible tracking strategy . The Galileo constellation achieved partial operational capability in late 2018 with 24 active satellites. In this analysis, we use two types of Galileo satellites: In-Orbit Validation (IOV), launched in 2011, and Full Operational Capability (FOC), launched in 2014-2018 . The IOV LRAs consist of 84 corner cubes, whereas FOC LRAs consist of 60 corner cubes, both types are arranged as flat panels 2 ( Table 1).

Determination of geodetic parameters based on SLR observations to LEO and GNSS
The realization of terrestrial reference frames based on active LEO was previously studied by, e.g., Haines  Männel and Rothacher (2017), Koenig (2018), who focused on selected missions and the possibility of deriving origin, scale, and the orientation of the reference frames using LEO. SLR measurements to single active LEO satellites were used by Guo et al. (2018) together with the GPS-based kinematic orbits of GRACE-A for the determination of SLR station coordinates in 2012. The authors achieved a consistency at the level of 20-30 mm, w.r.t. SLR-specific ITRF2014 realization, SLRF2014. Couhert et al. (2018) used DORISonly and SLR-only observations to Jason-2. Arnold et al. (2019) used SLR observations and GNSS-based orbits of LEO satellites for the determination of corrections to SLR station coordinates, range biases, and timing offsets and achieved 5-10 mm precision of SLR residuals. Strugarek et al. (2019a) combined SLR measurements to Sentinel-3A/B and LAGEOS-1/2 and achieved a consistency of estimated station coordinates, characterized by interquartile range, better than 13 mm and GCC with an RMS of 6 mm. Considering SLR to GNSS-only, and combinations with spherical geodetic satellites, Thaller et al. (2011) combined SLR observations to GPS and GLO-NASS satellites for the determination of the terrestrial reference frame, the GNSS satellite antenna offsets, and the SLR range biases and achieved 1-2 cm agreement of estimated station coordinates with a priori terrestrial reference frame. Sośnica et al. (2018) combined SLR measurements to Galileo, GLONASS, BeiDou and QZSS with LAGEOS-1/2 and achieved an improvement of the SLR station coordinate repeatability w.r.t. LAGEOS-only solutions at the level of 6 to 15%, whereas in Sośnica et al. (2019) the authors achieved an improvement in terms of non-core SLR station coordinate repeatability from about 20-30 mm to the level of 15-20 mm. Furthermore, the determination of Galileo orbits based solely on SLR observations can be found in Urschl et al. (2008), Hackel et al. (2015, Bury et al. (2019).

Objectives of this study
The goal of this research is to derive station coordinates based on SLR observations to eight active LEO satellites: Sentinel-3A, GRACE-A/B, TerraSAR-X, Jason-2, SWARM-A/B/C; 13 high-orbiting Galileo satellites (three IOV and ten FOC); and 3 passive spherical geodetic satellites: LAGEOS-1/2 in medium orbits and LARES in a low orbit. We aim to investigate SLR data to as large number of satellites as possible, especially in the case of LEOs for improving the determination of SLR station coordinates and to assess the consistency level between SLR observations collected to different missions. The Sentinel-3A mission started on February, 2016, GRACE-A/B mission ended on October, 2017, and also other LEOs, such as SWARMs orbited the Earth in 2016. In addition, the Galileo constellation have been being developed since the end of 2014. Thus, we conduct the analysis for the 2016.0-2017.0 period. We derive the best combination strategy of SLR observations to passive and active satellites at different heights, from 330 to 26,200 km, different revolution periods, and equipped with different onboard devices and LRA (Table 1, Fig. 1).
We check whether SLR measurements to other than spherical geodetic satellites either can be used as an alternative or as an enhancement of standard SLR solutions, or whether they are insufficient for the determination of SLR station coordinates of better quality because of the systematic errors embedded in GNSS-based orbits and SLR measurements to non-spherical LRAs. We investigate the quality of determined SLR station coordinates and the quantity of possible solutions. In the past, such an analysis was not possible because of the difficulties in precise orbit determination (POD) of active satellites equipped with large-size solar panels and resulting in substantial non-gravitational forces perturbing the orbits. Modeling of precise orbits of heavy, spherical, geodetic satellites is much easier and hence has been employed for the determination of geodetic parameters. In this study, we use state-of-the-art modeling strategies for LEO and GNSS orbits employing precise GNSS-based LEO and Galileo orbits with satellite attitude corrections measured down to 1-s intervals for LEOs (Table 1).
We calculate the solutions based on LEOs-only, Galileo-only, and various LEOs, Galileo, and spherical geodetic satellite combinations, also including different weighting of Normal Equation Systems (NEQs), and different constraining of estimating parameters. In this study, we assume the high accuracy of the GNSS-based orbits for the Galileo constellation provided as a part of the MGEX project, and for the LEOs determined by the Astronomical Institute of the University of Bern (AIUB), and by German Aerospace Center (Deutsches Zentrum für Luft-und Raumfahrt; DLR).

SLR normal points and data modeling
The ILRS provides time observations measured by SLR stations in form of the normal points (NPs) in the consolidated data format 3 . The NPs are formed using different temporal size of data bin, which depend on the mission type, and satellite orbit altitude, which, as shown in Table 1, may vary from 5 to 300 s. The modeling of SLR observations requires introducing various corrections. Corrections are related to signal propagation effects or LRA models, such as: the general relativistic correction due to the Earth gravitational bending of the light path (Shapiro effect; Petit and Luzum 2010); the tropospheric range correction for optical wavelength (Mendes and Pavlis 2004); range biases which are station-specific and satellite-specific corrections; the LRA reference point correction, which considers the change of optical signal propagation in the LRA prisms, and offsets between LRA position to the satellite center-of-mass, and the center of mass of trajectory obtained from the satellite orientation data and POD products.
For spherical geodetic satellites, one center-of-mass correction is sufficient to correct the measured range and refer it to the satellite center-of-mass. These corrections are typically station-dependent because of the different detector types and screening procedures employed at different SLR stations. For LEO satellites, two corrections are employed: one referring the center of LRA to the satellite center-of-mass and the second one which depends on the elevation and nadir angle of the laser incident angle. This allows for considering the specific construction of different LRA types consisting of different numbers of corner cubes. For GNSS, flat LRAs are used with only one correction between the center of LRA and the center-of-mass of the satellites. For GNSS, the dependencies between the elevation and nadir angle of the laser incident angle are thus neglected, which directly maps into the quality of registered SLR observations as the satellite signature effect.
The spacecraft attitude is irrelevant for POD and data processing in the case of spherical geodetic satellites. However, for active LEO and GNSS satellites, information about the spacecraft orientation is fundamental to properly refer the LRA position to the satellite center-ofmass. Most of the GNSS satellites follow the yaw-steering mode, for which the spacecraft orientation and the LRA vector orientation can be calculated using the positions of the Sun, Earth, and the satellite. Active LEO satellites follow different orientation modes, which can only approximately be calculated using analytical formulas, whereas the exact position has to be measured directly onboard satellites, using, e.g., star trackers and accelerometers. Orientation of the LEO spacecrafts is then stored in the attitude files and used for the calculation of the LRA vector orientation with respect to the satellite center-of-mass.

POD of satellites
Precise orbits are needed for deriving the high-quality SLR station coordinates. In this study, the microwave data were used for POD of LEO and Galileo satellites, whereas the laser observations are used for POD of LAGEOS-1/2, and LARES satellites. The Sentinel-3A, GRACE-A/B, SWARM-A/B/C orbits were generated at AIUB 4 (orbit products: the ESA Copernicus Sentinel Data Hub 5 ), whereas the TerraSAR-X and Jason-2 orbits were provided by the DLR (Hackel et al. 2017). LEO orbits are generated in the 1-day solutions and are based on undifferentiated GPS phase measurements, the ionosphere-free linear combination of two frequencies L1 and L2, using a reduced-dynamic approach (Wu et al. 1991) with float ambiguities and with applying the GPS antenna phase-center offsets and variations ). The empirical accelerations were constrained using piecewise constant parametrization and estimated in the radial, along-track, and cross-track directions. The non-gravitational perturbing forces were accounted for by empirical accelerations without any a priori satellite macro-models. The LEO POD is analogous to that described in the study by Arnold et al. (2019).
In the case of the Galileo orbits, we use the 1-day orbits from the Center for Orbit Determination in Europe (CODE) MGEX (Prange et al. 2018), which are based on network processing employing carrier phase double differences, and the ionosphere-free linear combination of two frequencies E1 and E5a. For Galileo, the extended empirical CODE orbit model ECOM2 (Arnold et al. 2015) is used assuming the yaw-steering attitude model ).
In the case of both LAGEOS-1/2 and LARES, we generate 7-day orbit solutions with the estimation of the six Keplerian orbit parameters, five empirical parameters, i.e., constant acceleration and two periodic accelerations in the along-track direction, two periodic accelerations in the cross-track direction, and stochastic pulses in the along-track directions (Sośnica et al. 2018. We use the ILRS realization of the ITRF2014, SLRF2014, for stations, and the IGS14 for satellite microwave orbits (Rebischung and Schmid 2016) as a priori station coordinate frames and the reference in station coordinate comparisons. All three reference frame realizations: ITRF2014, IGS14, and SLRF2014 share the same origin, scale, and the orientation and their rates over time, thus, are fully consistent.

SLR observations
The screening of SLR measurements is performed for each spherical geodetic solution, each LEO satelliteonly solution, and Galileo-only solution. The screening is based on residuals between range observations and orbits. As a priori values, we use the IERS 14 C04 ERPs (Bizouard et al. 2019), and SLRF2014 for station coordinates. We remove SLR observations with absolute residuals larger than 150 mm, or if the standard deviation of residuals is larger than 50 mm, the observations with largest residuals are being iteratively removed. Figure 2 shows the number of SLR observations (NPs) to LARES (LAR), LEOs, Galileo (GAL), and LAGEOS-1/2 (LAG) and the sum of all (ALL) in 7-day batches used in the analysis after data screening. In total, 358,767 NPs are used for LEOs, 121,801 NPs for LAG, 73,551 NPs for LARES, and 19,856 NPs for Galileo. The sum of all observations is 573,975. The number of LEO observations is thus almost three times of that of LAGEOS NPs, whereas the Galileo constellation constitutes only 4% of all used NPs, which can be explained by difficulties in tracking high-orbiting satellites by ILRS stations. For GAL, the number of weekly SLR NPs is stable at the mean level of 430 NPs, whereas for the LAR, LAG, and LEOs, we can observe an increase of the NPs number near the 115th and 225th day of the year, which may be explained by the beginning of Sentinel-3A tracking in April 2016 and the position of selected LEOs and spherical geodetic satellites on ILRS tracking priority list in 2016 6 . Including SLR data to LARES, Galileo, and in particular to LEOs, can increase the number of observations from the level of 2,400 NPs per week to the level of 11,500 NPs per week.
Next, we investigate the global distribution of SLR NPs. Figure 3 shows the global distribution of NPs collected by SLR stations to LAG, all LEOs, Galileo and LARES satellites used in this study for the example period of June 1-7, 2016. The number of NPs to other than LAG satellites has increased in the Europe, Australia, and Asia regions.
However, also other regions, such as South America and Africa, benefit in terms of the number of NPs. Moreover, we investigate the distribution of SLR NPs in more local scale, i.e., for azimuth and zenith angles of NPs for the example three stations in the period of June 1-7, 2016 (Fig. 4). The group of core stations, i.e., Yarragadee, Greenbelt, Matera, Hartebeesthoek, Haleakala, Zimmerwald, Mt Stromlo, Graz, Herstmonceux, and Potsdam, which are considered as top-performing, are characterized by the largest number of NPs to LAGEOS satellites, in all directions [e.g., Yarragadee (7090) station, Fig. 4, left]. Furthermore, most of core stations provide a large number of NPs to LEOs, especially at high zenith angles, and a few NPs to Galileo. Monument Peak (7110) and Arequipa (7403) (Fig. 4) show few or no NPs to Galileo, a moderate number of NPs to LAGEOS or LARES satellites, but relatively large number of NPs to LEO targets when compared to the core stations. The NPs in 2016 to Galileo satellites constitute the lowest percentage of all collected NPs, and are generated, in the most of cases, based on observations collected at high elevation angles. The reason for more NPs to Galileo at high elevation angles is that the return link is much stronger thanks to reduced atmospheric refraction and turbulences, lower noise (during daytime), and closer satellite-station distance. The NPs to LEOs, are mostly collected at low elevation angles, due to the low probability of the satellite pass through the station zenith, and a short period of flyby above the SLR station. Thus, the LEO NPs should mostly influence the horizontal components of station coordinates, whereas the Galileo NPs should mostly influence the Up component of station coordinates, when considering combined solutions based on SLR data to different satellites.

Generation of NEQs
After screening of SLR measurements, we generate the NEQs for particular satellites. We use precise, daily, microwave-based orbits of each LEO and Galileo, as well as SLR-based 7-day orbits of LAGEOS-1/2 and LARES. The generation of NEQs is based on linearized observation equations in the Gauss-Markoff model with considering the differences between SLR observations and modeled station-satellite ranges reduced by a series of corrections. We generate 1-day NEQs for each LEO separately, 1-day NEQs for group of all used Galileo, 7-day NEQs for both LAGEOS-1/2, and 7-day NEQs for LARES. NEQs include the parameters: station coordinates, range biases, GCC, the X and Y pole coordinates, and the Length of Day (LoD). In the case of all LEOs and Galileo, the orbital parameters in NEQs are fixed to the microwave-based orbits, whereas in the case of spherical geodetic satellites the orbit parameters are estimated. For LEO NEQs, we additionally introduce attitude orientation files with quaternions. We used information about the characteristics of the LRA reference point 7 and the center-of-mass corrections 8 provided by the particular mission supplier that are distributed by the ILRS. Data screening and generation of NEQs were conducted using the modified version of the Bernese GNSS Software .

Range bias determination
Range biases constitute one of the major error sources when processing SLR data to LAGEOS (Drinkwater et al. 2007;Appleby et al. 2016), GNSS Sośnica et al. 2015b), or LEOs (Arnold et al. 2015;Strugarek et al. 2019a). The omission of introducing station-satellite-specific range bias may deteriorate the solutions. Biases originate from erroneous sensor offsets, such as satellite center-of-mass corrections and GNSS antenna offset vectors, detector-specific and satellite reflector-specific signature effects (Strugarek et al. 2019b), uncalibrated biases occurring in station device circuits (Otsubo et al. 2001;Sośnica et al. 2015b), or deficiencies in SLR data modeling including the tropospheric delays ). However, the estimation of range biases increases the number of parameters. Thus, to stabilize the solutions and to reduce the number of parameters, we calculate range biases in separate processing as mean values for the entire year and we introduce them to combined solutions as known quantities. Only the station coordinate corrections and annual average range biases for each SLR station to a particular satellite are estimated in the process of deriving annual biases. Annual range biases are treated as station-specific, in analogy to Appleby et al. (2016), and satellite-specific parameters. The orbit parameters of spherical geodetic satellites are pre-eliminated in the NEQs, whereas the ERPs, GCC, are fixed to a priori values to retrieve the range bias embedded in the SLR observations. In range bias processing, we use 1-day NEQs for LEOs, and Galileo data, and 7-day NEQs for LAGEOS-1/2 and LARES, all together to generate annual solutions referenced to the middle day of the year. We also tested different approaches, where annual range biases and station coordinates were calculated using LAGEOS-1/2 only NEQs, then separately from LEO groups, and from Galileo NEQs, etc. However, these approaches resulted in cm-level biases for stations and a degradation of station coordinate results and thus are not shown here. Figure 5 shows the annual range biases for spherical geodetic, LEO, and Galileo satellites, calculated for each station-satellite pair. For Galileo satellites and most of stations the range biases are negative within − 10 to − 35 mm. For some stations, such as Komsomolsk (1868), Arkhyz (1886), Katzively (1893) the annual range biases reach over − 50 mm. Also Wettzell (8834), Monument Peak (7110), Brasilia (7407), and Yarragadee (7090) show large range biases, up to − 40 mm to Galileo satellites. In the case of LEOs and the most of stations, the range biases are within ± 15 mm (Fig. 5). The largest, negative corrections are derived for Arkhyz (1886) and Katzively (1893) station (up to − 50 mm), and also to Wettzell (8834) (up to − 30 mm). For some of the stations, e.g., Simeiz (1873), Altay (1879), and Arequipa (7403), range biases have large positive corrections at the level of 30 mm. Spherical geodetic LARES and LAGEOS-1/2 satellites are characterized by the positive range biases at the range of 5 to 15 mm for most of stations, with few The large annual range biases for Galileo are related with POD products of satellites. In 2016, the CODE Galileo orbits did not consider albedo, Earth infrared radiation, and antenna thrust modeling, which in total result in a shift of GNSS-derived satellite orbits of 36 to 40 mm in the radial direction (Bury et al. 2020). The value of estimated range bias depends on antenna offsets, phase ambiguity resolution, and POD product quality, and require a good SLR observation geometry to decorrelate station height and range bias estimates, which are calculated in annual processing as a mean annual value. Figure 6 shows the influence of considering the annual range bias corrections by comparing the solutions with and without annual station-satellite range biases. We compare the interquartile range (IQR, that is the first quartile subtracted from the third quartile) of coordinate residuals, derived as a difference of calculated SLR station coordinates to the SLRF2014 coordinates. In the case of the Up component, the reduction of IQR is 4 mm for LAGEOS, LEO, and Galileo solutions, whereas the North component is characterized by the reduction of IQR at the level of 2 mm for all tested solutions. The East component is characterized by IQR reduction of 5, 3, 9 mm for the LAGEOS, LEOs, and Galileo solutions, respectively. Thus, introducing annual range biases is somewhat beneficial for processing SLR observations to different types of satellites, and enhances the process of accounting for systematic errors, which lead to the more accurate solutions.

Combined SLR solutions
Most of the parameters, i.e., SLR station coordinates, GCC, and ERPs are common in the Galileo, LEO, LAGEOS, and LARES processing and can be stacked at the NEQ combination level. The combined solutions require consistent models, including, e.g., models describing relativistic effects, gravity field, solid Earth, ocean, and pole tides, and others. For the generation of NEQs, determination of annual range biases, processing of individual, and combined solutions, we use the same models listed in Table 2. The issue of possible reference frame differences from SLRF2014 in SLR processing with IGS14 from LEO and Galileo orbits is discussed in the further part of article, that is 'Network constraining tests' section. Tests conducted by Strugarek et al. (2019a) showed that 7-day LEO satellite solutions based on Sentinel-3A/B are needed for a good and global distribution of SLR observing stations. Depending on the case study (see Table 3), we generate LEO-only (LEO), Galileo-only (GAL), LAGEOS-1/2-only (LAG), LARES-only (LAR) solutions and combinations, i.e., LEO+GAL, LAG+LAR, LEO+GAL+LAG, ALL (all used satellites). Thus, an example of combined 7-day solution includes one 7-day NEQ for LAG, seven 1-day NEQs for Galileo, and seven 1-day NEQs of particular LEO satellites.
In the individual and combined solutions, we impose no-net-rotation and the no-net-translation constraints for the SLR core station network group. The core station list is provided by the ILRS Discontinuities File 9 . We estimate SLR station coordinates, GCC, LoD, and pole coordinates. Orbit parameters and annual range biases are fixed to annual values. Station coordinates and GCC are estimated as one set of parameters per 7-day solution. The pole coordinates and the LoD parameter are estimated with a 1-day resolution, parameterized as piecewise linear (Table 2). For outlier detection of core stations, the resulted coordinates of core stations were verified in every 7-day solution using a Helmert transformation with imposing a threshold of 20 mm for each coordinate component, which follows the strategy described by Zajdel et al. (2019). After outlier detection, we updated the SLR core station list and recalculated the solutions. Individual and combined solutions were calculated using the modified version of the Bernese GNSS software  for the period of 2016.0 to 2017.0.
We test excluding some satellites (e.g., LEOs, LARES) from combinations or employing the weighting of NEQs in combined solutions to reduce the negative impact of the lowest-performing satellites. Excluding the satellites from combinations led to a decrease of the number of SLR observations and a decrease of the percentage of successful solutions for SLR stations, which resulted in no or rather minor improvement of determined parameters with respect to the LAGEOS solutions. Thus, we test combinations with NEQs weighting, in which weighting scale factors for the particular satellites are employed. The scale factor can be expressed as ratio of the square roots of variances. In this study, we modify weighting scale factors based on pre-analysis results. The remaining processing scheme and estimated parameters are the same as in the case of combined solutions, except for NEQ weighting.
To select weighting scale factors and investigate the quality of individual SLR solutions, we used the parameter estimation approach. We analyze the estimated SLR station coordinate repeatability by means of IQR 3D statistics of individual solutions based only on particular LEO, LAGEOS-1/2, LARES, and 13 Galileo. The 3D IQR value is calculated as a square root of the sum of the North, East, and Up IQRs squared. Figure 7 shows that solutions based only on LAG, Sentinel-3A, Jason-2 are characterized by the lowest 3D IQR of SLR station coordinates at the level of 21 mm. The quality of LEO satellite orbits is better for the satellites at altitudes higher than 500 km (see Fig. 1, Table 1), due to the lower influence of non-gravitational forces. Moreover, the shorter, 1 s attitude data intervals, were used in processing for Sentinel-3A, which may affect the particular LEO solution (Table 1). However, the Jason-2 solution despite 30 s attitude interval is also characterized by low IQR 3D statistics. Based on pre-analysis results (Fig. 7), we calculated the weighting scale factor for each satellite NEQs, as a ratio between 3D IQR squared for LAGEOS to 3D IQR squared of particular satellite solution. Therefore, the weighting scheme resulted in scaling factors: 1.00 for LAG, Sentinel-3A, and Jason-2, 0.44 for SWARM-B/C, 0.25 for TerraSAR-X and Swarm-A, 0.16 for GRACE-A and LARES, and 0.11 for GRACE-B and Galileo. The solutions based on weighted NEQs of all analyzed satellites are named ALL+W (Table 3). We tested also other weighing strategies based on observation residuals and the a posteriori sigma of a unit weight of individual solutions, and standard deviations of estimated station coordinates. However, the strategy based on weights adjusted by using IQR values turned out to give the best results of the combinations out of all tested strategies, therefore has been adopted for further analyses.

Station network
We analyze calculated coordinates for all used SLR stations, core stations, and the non-core group. First, we compare specific satellite-type solutions, combinations of different satellite types, and the combinations with the weighting of observations in terms of IQR. Then, we investigate the formal errors of calculated station coordinates, as well as number of solutions with improved coordinate residuals, and time series of station coordinates. Finally, we conduct the tests with different constraining approaches to discuss the possible issue of reference frame inconsistencies. Figure 8 illustrates the SLR station coordinate IQR based on LAG, LEO, LAR, and GAL solutions. Considering all stations, the lowest IQR values of the Up component are at the level of 17 mm for LAG and LEO solutions, whereas the East and North components are at the level of 10 and 11 mm for LAG and LEO, respectively. The LAR and GAL solutions are characterized by an IQR worse by a few to over a dozen mm in the case of all components (Fig. 8, left). Considering the core stations (Fig. 8, right), the lowest IQR values for the Up component are at the level of 9-10 mm for LAG and LEO. Both horizontal components are at the lowest level of less than 5 and 7 mm of IQR for LAG and LEO, respectively. LEO-only solutions obtain similar results w.r.t LAG solutions, the LAR solutions result   in a few-mm worse results, whereas the GAL solutions result in station coordinate statistics that are worse by a dozen of mm. LAGEOS satellites allow for deriving high-quality station coordinates because of their mission characteristics, such as orbit altitude, low area-to-mass ratio and thus minimized orbit perturbing forces, and the homogeneous arrangement of corner cubes, which enables to compute an accurate correction for the satellite centerof-mass. The LARES satellite has also low area-to-mass ratio, but its orbital altitude, lower than that of LAGEOS, makes LARES more sensitive to the residual atmospheric drag, time-variable geopotential changes and errors in background models, such as ocean tidal model that are not compensated by additional parameters, all of which may generate orbit modeling errors and cause the inferior quality of LAR solutions. Moreover, LAR solutions are based on measurements to only one satellite, which results in a decreased ratio of the number of observations to estimated parameters. However, we employ the same orbit modeling in terms of the number of empirical orbit parameters for LARES and LAGEOS-1/2 to keep a consistency of both solutions and to investigate their impact on further combinations.
The high-quality of SLR station coordinates based on altimetry LEO satellites at altitudes above 500 km (Sentinel-3A, Jason-2) confirms the high quality of microwave LEO orbits. Microwave GNSS-based orbits benefit from: pseudo-stochastic accelerations, which absorb non-gravitational orbit perturbations and geopotential variations; frequent attitude information data; high-quality LRA vectors w.r.t. satellite center-of-mass; and satellite azimuthnadir-dependent range correction of the LRA reference point (Table 1). Galileo solutions are insufficient for the determination of the highest-quality SLR station coordinates for the year 2016 due to the low number of observations to single satellites, microwave orbit errors related to non-gravitational orbit perturbations (Bury et al. 2020). For Galileo, the LRA correction does not consider satellite azimuth-nadir angle corrections; therefore, the SLR observations at high-zenith angles are possibly affected by errors for stations which are not designed to operate at a low signal level from medium Earth orbit altitudes. Finally, Galileo satellites are typically observed only in the near-zenith direction, therefore the horizontal station coordinates cannot be well determined as opposed to the LEO satellites, for which the majority of observations are collected at high zenith angles.
Before the final combinations, we investigate the influence of SLR data to other than spherical satellites. Figure 9 shows how many station coordinates in % can be successfully calculated from 7-day solutions for SLR stations when using LAG, LEO, and ALL observations. In total 52 weekly solutions were generated in 2016.0-2017.0. Numerous stations observe LAGEOS, LARES, LEO and Galileo satellites on a regular basis, such as Yarragadee (7090), Changchun (7237), Mount Stromlo (7825), or Greenbelt (7105). However, some stations, such as Baikonur (1887), Komsomolsk (1868), Altay (1879), do not frequently observe LEO satellites, which is related with their main purpose to support GLO-NASS, rather than tracking LEO targets, which have much shorter passes. Most of stations increased the percentage of obtained solutions, at the level of 10-20% for, e.g., Monument Peak (7110), Wettzell (8834), Potsdam (7841) or even doubled the percentage of solutions, such as Areqiupa (7403), Badary (1890), while comparing the ALL with LAG solutions. LEO targets caused the largest increase of possible solutions to obtain for stations in the ALL approach. Figure 10 illustrates the repeatability of estimated station coordinates by means of IQR for combined solutions. When considering all stations (Fig. 10, left), IQR values for the Up component are at level of 17, 16.5 mm for LAG, and ALL+W, respectively, and at the level of 18 mm IQR for other tested solutions. The horizontal components provide IQR values at the level of 9-10 mm, for the ALL+W combinations. The LAG solutions are characterized by IQR values worse by 2 mm for the both, East and North components. When considering the core stations (Fig. 10, right) and the Up component, IQR results are less than 10 mm for ALL+W, and LAG, at the level of 10-11 mm for remaining solutions.
The results of IQR statistics for both ALL+W and LAG are at a similar level with differences of 0.5 mm, which is less than a change of 5% in the case of all and core stations. The horizontal components of non-core stations in ALL+W show an improvement of 2-mm IQR, which corresponds to 10% change, w.r.t LAG-only. The LEO+GAL and LAG+LAR show the worst IQR values for all tested components and SLR station groups. The improvement of station repeatability for non-core stations might be caused by introducing a large number of SLR observations, especially to high-quality LEO orbits. The LEO targets, characterized by precise microwave orbits, altitude higher than 500 km with frequent attitude information, have a positive impact on the combination results. The LAG+LAR solutions are characterized by a deteriorated IQR statistics, which suggests that LARES orbit modeling may require more frequent determination of dynamic orbit parameters than one set per 7-day to compensate for the residual atmospheric drag, time-variable gravity, and errors in the background models.
The proper weighting of SLR data to other than spherical geodetic satellite targets improves the quality of the horizontal coordinate components for non-core stations.
Core SLR station coordinates are neither improved nor deteriorated, however, the number of solutions increases. The determination of core station coordinates may be improved when using more refined POD products, especially for LARES, Galileo, and LEOs at altitudes lower than 500 km.

SLR station coordinate quality
We check the observation geometry and the parameter observability for station coordinates in NEQs in terms of station coordinate formal errors (Fig. 11). The median values of formal errors of all components are reduced for ALL, and especially for ALL+W solutions, w.r.t. LAG solutions. While considering core stations, the median values of formal error for all components are less than 3 mm for LAG, and ALL solutions, and less than 1.5 mm for ALL+W solutions. Also IQR values and the maxmin ranges of formal errors of all components are about 1-3 mm lower for ALL+W than for LAG solutions, e.g., Hartebeesthoek (7501), Haleakala (7119), Graz (7839). In the case of most of the non-core stations, median values are higher by a few mm than for core stations, and vary from 1 to 5 mm, 1 to 7 mm, 1 to 6 mm, for the Up, North, and East components, respectively. Median values of non-core station formal errors show an improvement of 0.5 and more than 1 mm for the Up and horizontal components, respectively, for ALL, ALL+W, w.r.t. the LAG solutions. Some of the stations, e.g., Arequipa (7403), Borowiec (7811), Simeiz (1873), and Katzively (1893) are characterized by an increased reduction of 3-4 mm of the median formal error for horizontal components, when comparing ALL+W with LAG solutions. A similar reduction of a few mm occurs in terms of IQR, and maxmin range of formal errors for horizontal components.
Almost all stations in ALL+W solution are characterized by a substantial reduction of formal error statistics, w.r.t. LAG solutions. The level of reduction is from 1 mm to even several mm, especially for horizontal components. The reduction of station coordinate formal errors is caused by the improvement of SLR observation geometry when introducing observations to satellites at different zenith angles (LEO for high angles and Galileo for low angles) with a good sky distribution as seen from the SLR stations and the reduction of correlations between estimated parameters.
A comparison of coordinate statistics for each SLR station (see Table 4) shows that in the case of most of the stations, ALL+W solutions have IQR values lower by up to a few mm than the LAG solutions. The reduction of IQR occurs in the case of the core stations, e.g., Hartebeesthoek (7501), Haleakala (7119), Matera (7941), and especially for non-core stations, such as, Simosato  (7080), are characterized by the coordinate statistics at the cm-level, which may be caused by a low number of SLR observations, or station-specific errors. The use of SLR measurements to other than spherical geodetic satellites increases precision for most of the non-core station coordinates for all components. The core stations show neither a degradation nor a large improvement of the Up and horizontal components when compared to LAG results. Although, introducing LEO, LARES, and Galileo SLR data can rather be treated as an enlargement of the number of observations, successful solutions, and the improvement of observations geometry, which results in lower formal errors of station coordinates, rather than the substitution of SLR data to LAGEOS.
Next, we investigate whether the ALL+W combinations increase the number of solutions with lower coordinate residuals and reduce variability of coordinate time series w.r.t. the LAG solutions. Figure 12 shows two examples of stations: Potsdam and Arequipa. For the Potsdam (7481) station, the number of solutions with residuals accumulated within the range − 4 to 4 mm for ALL+W is greater by 14, 10, and 8 than in LAG, for the Up, North, and East component, respectively (Fig. 12, left). In the case of Arequipa (7403) (Fig. 12, right), ALL+W and LAG solutions are less concentrated around zero residual values than for Potsdam (7841). However, ALL+W shows over 7, 8, and 10 more solutions concentrated in the range between − 4 and 4 mm than LAG, for the Up, North and East components, respectively. Thus, for Potsdam, and Arequipa there are even two times more ALL+W solutions with lower coordinate residuals than in LAG solutions. Figure 13 shows an example of the Changchun (7237) station time series, decomposed into the Up, North, and East components for LAG, ALL, and ALL+W solutions. In the case of the Up component, we can observe a 20 mm offset, which is related to the inferior a priori station coordinates in SLRF2014 affected by a range bias and not fully valid non-linear post-seismic deformations after 2014. In the beginning of 2016, all tested solutions show a variability of residuals in the range of − 15 to 15 mm for all components. The ALL and ALL+W solutions from April 2016 show more stable results with variability in the range of less than ±10 mm, in the case of the North and Up components, when compared to LAG. LAG solution is affected by the larger variability and the larger formal errors than ALL and ALL+W.

Network constraining and SLR-PPP solutions
The combination of SLR observations to satellites, whose orbits are calculated from SLR and GNSS techniques, implies the possibility of reference frame inconsistencies. We analyze different parameter constraining methods of weighted combinations to investigate the influence of fixed, microwave-based LEO and GNSS orbits on the realization of the reference frame. We generate combined solutions (ALL+W+NoPar) with no-net-rotation and no-net-translation network constraining, as in the case of ALL+W, LAG, and previous solutions. However, remaining parameters, such as GCC and ERPs are strongly constrained to a priori values or to values calculated in the separate processing (annual range biases). ALL+W+NoPar should be characterized as very stable solutions, because only station coordinates are determined.

Table 4 Median and interquartile ranges (IQR) of SLR station coordinates (in mm)
Bolded station numbers and location denote the core stations, IQR are given in parenthesis In other tested combinations (ALL+W+PPP), we estimate SLR station coordinates as free parameters, with fixing all remaining parameters (except for LAG and LAR orbits). ALL+W+PPP solutions can be considered as SLR Precise Point Positioning (PPP) solutions in analogy to GNSS-PPP solutions, due to fixing microwave-based orbits and clocks, and not imposing any constraints on the station network. The station coordinates are calculated independently for each station in PPP, because there are no common parameters that could establish the network solution. In SLR-PPP we do not need to estimate station clocks or troposphere delays as it is typically done in GNSS-PPP, because SLR measurement principle employs only one clock (the station clock), and modeling of troposphere delays in SLR is more accurate than in GNSS for the wet component of the atmosphere. We also tested SLR+PPP solutions based on SLR to LAGEOS-1/2 only (not shown), however the results of station coordinates showed deficiencies in the proper recovery of the East station coordinate components. Large differences of the East station coordinate are caused by the singularity of solution, due to unobservability of the ascending node of LAGEOS orbits causing a singularity of the rotation of reference frame around Z-axis. This test showed that the SLR-PPP solution is possible only on fixed orbits. Figure 14 shows the repeatability of coordinates for the group of all stations for LAG, ALL+W, ALL+W+NoPar, and ALL+W+PPP solutions. ALL+W+NoPar solutions are characterized by a similar repeatability of the Up and North components, w.r.t. LAG and ALL+W solutions. In the case of the East component, the IQR are 3 mm higher for all tested groups of stations, w.r.t LAG and ALL+W. Nevertheless, the repeatability of estimated SLR station coordinates does not substantially change when changing the network constraining or when even estimating all station coordinates independently in SLR-PPP. Network constraining in, e.g., LAG, ALL+W, ALL+W+NoPar, enforces the determined reference frame origin and orientation to those from the SLRF2014. The origin of the frame is then in the center-of-figure of the Earth or, more precisely, the center-of-network because the origin is realized by the network of the stations. ALL+W+NoPar solutions characterized by the lowest number of estimated parameters show similar or worse by 3-mm results of station repeatability, which may be related to small inconsistencies in the reference frames and over-constraining of the solution. For ALL+W+PPP solutions, the realized frame is delivered by microwave orbits of LEO and GNSS in the IGS14 reference frame supported by SLR-based orbits of LAGEOS and LARES. Moreover, ALL+W+PPP solution results enforces the origin of the reference frame to be in the Earth's centerof-mass because of no network constraining. The reference frame origin in the PPP solutions is thus provided by the satellite orbits. Estimating stochastic parameters in the POD of LEOs allows them to orbit around the center-of-mass of the Earth. The vector of center-ofmass referenced to the center-of-figure represents the GCC. Thus, the relatively low 1-5 mm increase of IQR in ALL+W+PPP solution is mostly caused by GCC motion, and to a small extent by possible inconsistencies of reference frames delivered by different types of orbits. Moreover, introducing SLR measurements to other than LAGEOS satellites in SLR-PPP solutions removes the singularities related to the drift of ascending node parameter, and allows for the determination of SLR station coordinates.
To confirm this interpretation, we calculate the Helmert transformation between ALL+W and ALL+W+NoPar solutions, and between ALL+W and  Strugarek et al. Earth, Planets and Space (2021) 73:87 ALL+W+PPP solutions. Figure 15 shows the X, Y, and Z translation parameters compared with GCC estimated from ALL+W and LAG solutions. The X and Y GCC show a high consistency of results, when comparing LAG with ALL+W solutions. In the case of Z GCC, ALL+W solutions are less scattered, but represent similar variations as LAG solutions. The X, Y, and Z translation parameters of ALL+W and ALL+W+PPP show similar results to GCC from ALL+W and LAG, and confirm that station coordinates from ALL+W+PPP include the GCC motion. Thus, ALL+W solutions are not affected to a large extent by differences between SLRF2014 and IGS14 frames. Also, near-zero translation parameters of ALL+W and ALL+W+NoPar solutions show that the influence of the GCC motion is compensated by GCC that is estimated in ALL+W.

Conclusions and discussion
Expansion of so far used SLR measurements to LAGEOS-1/2 satellites, by considering also LARES, LEO, and Galileo data for the realization of the terrestrial reference frame and determination of station coordinates requires proper handling of range biases, high-quality satellite orbits, attitude information, and correct weighting of observations. In this study, we used SLR observations from 33 laser stations to 24 satellites of different characteristics, i.e., LAGEOS-1/2, LARES, Jason-2, Sentinel-3A, SWARM-A/B/C, TerraSAR-X, GRACE-A/B and 13 Galileo to determine station coordinates based on the methodology which combines laser data and microwavebased orbits. The processing considers GNSS-based and SLR-based POD products of satellites, generating a series of 7-day solutions with no-net-rotation and nonet-translation network constraining, with verification of the list of core stations, weighting tests, estimated parameter constraining tests, and re-substitution of annual SLR range biases. The processing of SLR observations to different types of satellites for the estimation of parameters needs introducing range bias corrections. Range biases are used for accounting of SLR systematic errors, which may be related to the orbit modeling, solution processing deficiencies, delays in station circuits, troposphere refraction biases, and others. We propose using the station-satellite-specific annual range biases. Introducing annual range bias correction increases the repeatability of station coordinates by 2 to 9 mm, depending on station coordinate component and solution. In the case of LEO satellites and most of stations, the annual the range biases are within the range of − 15 to + 15 mm, whereas for the Galileo satellites the annual range biases are within the range of − 10 to even − 50 mm. SLR observations to spherical geodetic, LEO, and Galileo satellites increased the number of used NPs by almost three times, and increased the percentage of successfully determined solutions for most of the SLR stations, w.r.t. LAG-only solutions. Most of stations showed an increase of the number of combined solutions at the level of 10-20%, w.r.t LAG-only solutions. Moreover, some stations, such as Arequipa or Badary, even doubled the number of successfully determined station coordinates.  The determination of SLR station coordinates based on SLR data to LEO-only satellites results in an IQR at the level of 17 mm for the Up component and 11 mm for the horizontal components for all stations. For the core stations, the IQR is 10 mm for the Up, and 7 mm for the horizontal components. The IQR values of all stations for multi-satellite combination are increased by about 2 mm for the horizontal components, and at a similar level for the Up component, w.r.t LAG solutions. However, the combined solution with weighting of observations (ALL+W) does not deteriorate the statistics. Weighting of observations is based on parameter estimation approach in terms of station coordinate results for a particular satellite. Here, the IQR values of all stations from ALL+W and LAG solutions are at the similar level of 17, 10, and 9 mm for the Up, North, and East components, respectively. The non-core SLR stations benefit the most from the ALL+W approach. The horizontal components of non-core stations show IQR lower by about 2-mm, which corresponds to the improvement of 10% w.r.t. LAG solutions.
Also, the ALL+W solutions are characterized by better statistics than LAG solutions in terms of geometry of observations. The reduction of the formal errors in terms of median and IQR reaches even 50%, i.e., one to a few-mm level in the case of almost all stations. The improvement of formal errors is largest for horizontal components. The spread of coordinate residuals is also better in the multi-satellite solutions which means that the coordinate time series become more stable. Some stations, such as Arequipa and Potsdam, have even two times more solutions characterized by lower coordinate residuals in ALL+W than in LAG solutions, whereas the Changchun station shows more stable time series of coordinates that are reduced from about ± 15 to ± 10 mm. Moreover, combining SLR measurements to different types of satellites allows for the determination of SLR station coordinates without network constraining and independently for each station. The SLR-PPP solutions (ALL+W+PPP) show repeatability of all tested station coordinates at the level of less than 20, 15, and 12 mm for the Up, East, and North components, respectively. The SLR-PPP solution has its origin the the Earth's center-of-mass because the network origin is transferred by satellite orbits.
This paper shows that SLR observations to spherical geodetic, LEO and Galileo satellites together with precise microwave orbits of LEO and Galileo satellites, allow for the determination of high-quality SLR station coordinates. The differences of reference frames delivered by SLR-and GNSS-based orbits in SLR combination are minor and do not affect the solution. The determination of station coordinates based on weighted combination resulted in an improvement, especially for horizontal coordinates of the non-core SLR stations, and no degradation of coordinates for the remaining stations. Findings correspond well with the research by Sośnica et al. (2014) based on the combination of different types of spherical geodetic satellites, which also resulted in the most prominent improvement for the horizontal SLR station components. Also, the combinations results confirms the previous findings of Strugarek et al. (2019b), where LAGEOS-1/2 and Sentinel-3A/B combinations increased the consistency of station coordinates and allowed for the SLR-PPP approach.
The multi-satellite combination benefits the most from typically used spherical geodetic LAGEOS-1/2 satellites, and two LEOs, i.e., Sentinel-3A, Jason-2, characterized by the altitude higher than 500 km, the high quality of attitude and spacecraft orientation information. However, the SLR data to the European GNSS, Galileo, lead to inferior quality of station coordinates with respect to SLR tracking of LEO and spherical geodetic satellites. The worse performance of Galileo solutions is caused by a small number of observations in the analyzed period, microwave orbit errors, and the signature effect of flat LRA installed onboard GNSS, which increases the spread of measured ranges, when a satellite is observed at high zenith angles. Therefore, better modeling of the GNSS signature effect and nadirdependent corrections for GNSS LRAs are indispensable for deriving high-quality geodetic parameters in the near future. LARES-only solutions also are of inferior quality in terms of estimated station coordinates, due to orbit quality, which require frequent determination of dynamic orbit parameters.
The study allowed for deriving the following highlights: • The vast majority of all SLR observations constitute measurements to LEO satellites. When considering them for other purposes than just the orbit validation, i.e., determination of station coordinates, the number of observations is increased, even by factor of five, the geometry of observations improves, and the formal errors of station coordinates are reduced. • The SLR-GNSS combination onboard LEO and Galileo satellites allows for the inter-technique combination through the co-location in space. • Mitigating of station-satellite-specific range biases is beneficial for the processing of SLR observations to different types of satellites, as these parameters can absorb systematic errors caused, e.g., by deficient orbit modeling. • The solutions based on multi-satellite constellation of geodetic, LEO, and GNSS satellites show a promising perspective for the estimation of SLR station coordinates. The SLR-PPP solution is plausible when using multiple satellites. • The multi-satellite combination with the proper weighting of observations shows better statistics of SLR station coordinates than the LAGEOS-only solution. The determined coordinates are characterized by decreased formal errors, decreased variability of residuals and IQR values, and increased stability over time, especially for the non-core stations. However, when the weighting is omitted, the estimated parameters may be deteriorated with respect to the LAGEOS-only solutions.
Methods of improving the quality of station coordinates and other parameters based on SLR technique consider increasing the number SLR stations and their global distribution (e.g., Otsubo et al. 2016;Kehm et al. 2018;Glaser et al. 2019), or as in this study, exploiting the full potential of SLR stations by using observations to various satellite types at different altitudes that are currently being tracked by SLR stations. The processing in our approach requires considering the features of two observation techniques, as well as specific characteristics of observed satellites for the determination of parameters. The combined SLR data processing must employ: precise orbits of considered satellites, attitude information for LEOs, consistent reference frames, compatible models describing phenomena of Earth's and space systems, individual processing scenarios for different types of satellites, and a reduction of systematic effects, related to the observation techniques. The quality of determined parameters may be limited due to time-variable number of satellites equipped with laser retroreflectors, sparse distribution of SLR stations, relatively low number of SLR observations, w.r.t GNSS, and the performance of SLR, and GNSS techniques. Recent and ongoing developments conducted by various scientific groups and space agencies led to improvement, e.g., in the terms of POD products, especially for LEOs (e.g., CPOD Quality Working Group, Fernández et al. 2015) and GNSS satellites (e.g., MGEX project, Montenbruck et al. 2017), or upgrading processing models in spherical geodetic SLR satellites. Our expanded methodology may improve the realization of the SLR reference frame, and therefore, should be taken into consideration in the further realizations of terrestrial reference frames. We can expect a further improvement of results, based on the derived approach when considering the development of GNSS-based POD of LEOs and GNSS satellites and including other spherical geodetic satellites, e.g., Starlette, Stella, Ajisai or other observation techniques, such as DORIS for the co-location in space onboard LEO. The improvement can likely be reached by using precise correction of center-of-mass position for geodetic spherical satellites (Rodríguez et al. 2019), ambiguity-fixed LEOs orbits, with better consideration of orbit perturbations, such as non-gravitational force modeling and satellite macromodels (Mao et al. 2020) or by using more precise Galileo orbits based on solutions which consider the boxwing model with the estimation of the small set of empirical accelerations (Bury et al. 2020). Our approach can also benefit from the modernization and the emergence of other global GNSS systems, as well as the improvement in SLR processing, such as range bias handling, upgrading models describing geopotential, time-variable gravity field, atmosphere and ocean dealiasing, updating center-of-mass correction, and modeling of satellite signature effect and troposphere refraction biases.
The approach described in this paper employs the GNSS-based LEO orbits. Thus, the estimated parameters of interest derived from the combination rely on observations to spherical geodetic satellites, and consequently, to the performance of SLR technique. However, in the case of observations to LEOs and Galileo, both SLR and GNSS performance impact the combinations. As shown in our study, the quality of derived station coordinates confirms advancement in the consistency of SLR and GNSS as a step towards the full integration of space geodetic techniques. This can be enhanced in the future by the simultaneous POD of LEO using GNSS and SLR data and co-location of both techniques onboard LEO satellites, which may diminish technique-specific systematic errors and provide high-quality parameters. Nevertheless, this paper shows a great potential and the excellent consistency between GNSS-based LEO orbits, GNSS-based Galileo orbits, and SLR-based LAGEOS/LARES orbits in the combined multi-GNSS satellite solutions, which is of crucial importance in terms of the GGOS goals aiming at the Earth monitoring with the accuracy of 1 mm.