Differencing strategies for SLR observations at the Wettzell observatory

The precise estimation of geodetic parameters using single- and double-differenced SLR observations is investigated. While the differencing of observables is a standard approach for the GNSS processing, double differences of simultaneous SLR observations are practically impossible to obtain due to the SLR basic principle of observing one satellite at a time. Despite this, the availability of co-located SLR telescopes and the use of the alternative concept of quasi-simultaneity allow the forming of SLR differences under certain assumptions, thus enabling the use of these processing strategies. These differences are in principle almost free of both, satellite- and station-specific error sources, and are shown to be a valuable tool to obtain relative coordinates and range biases, and to validate local ties. Tested with the two co-located SLR telescopes at the Geodetic Observatory Wettzell (Germany) using SLR observations to GLONASS and LAGEOS, the developed differencing approach shows that it is possible to obtain single- and double-difference residuals at the millimetre level, and that it is possible to estimate parameters, such as range biases at the stations and the local baseline vector with a precision at the millimetre level and an accuracy comparable to traditional terrestrial survey methods. The presented SLR differences constitute a valuable alternative for the monitoring of the local baselines and the estimation of geodetic parameters.


Introduction
Satellite laser ranging (SLR) is one of the four main geodetic techniques involved in the realisation of the International Terrestrial Reference Frame (ITRF). In particular, SLR contributes to the determination of the scale and the origin of the frame, but it is also involved in the determination of the gravity field of the Earth, and it is widely used for the validation of satellite orbits. In the operating principle of SLR, at an SLR ground station, a short laser pulse is generated and sent with an optical telescope to a single satellite. The satellite has retro-reflectors attached to its surface that reflect the laser pulse back to the ground station. With the telescope at the ground station the reflected pulse is then detected and analysed. Based on the time of emission of B Iván Herrera Pinzón Ivan.Herrera@geod.baug.ethz.ch 1 Institute of Geodesy and Photogrammetry, ETH Zurich, Zurich, Switzerland the original pulse and the time of reception of the reflected pulse, the light travel time, the range between the satellite and the telescope can be computed. Core stations of the International Laser Ranging Service (ILRS) (Pearlman et al. 2019b) are able to provide normal point measurements with a millimetre-level precision (Luceri et al. 2019). Thus, observations to SLR-dedicated satellite missions, such as LAGEOS or Etalon, and more recently to Global Navigation Satellite Systems (GNSS) satellites, turn SLR into a reliable and accurate technique for the determination of geodetic and geophysical parameters. However, systematic errors at the station level together with mismodellings of the centre of masses of the satellites, threaten the accuracy of these range observations. To this end, the ILRS performs the continuous monitoring of systematic error sources at the station level, for an enhancement of the realisation of the ITRF, searching for improving the agreement between the scale provided by SLR and the remaining geodetic techniques (Luceri et al. 2019). Systematic time biases in range observations (which amount to few μs) degrade the accuracy of the estimated satellite orbits, affecting the station coordinates by several millimetres, and have to be removed, for instance with the time synchronisation of stations using time transfer by laser links to satellites (Exertier et al. 2017).
At the Geodetic Observatory Wettzell, two SLR telescopes, increasing the observation capabilities, have been operating for quite some time. In consequence, an ideal setup for the assessment and testing of the differencing approach has become available on a very short SLR baseline. The two co-located SLR telescopes, connected by a local tie, controlled by a common timing system and affected by (almost) the same atmospheric effects, allow the study of the size and stability of instrumental biases, the quality of the local tie, and a detailed investigation of new processing strategies. Forming double differences is a standard approach in GNSS processing (Hofmann-Wellenhof et al. 2008). The advantage of forming single or double differences is the elimination (or strong reduction) of satellite-(differences between stations) and station-specific (differences between satellites) error sources and the high accuracy achieved for short baselines, as well as an adequate approach to reduce the number of parameters (especially clock corrections) during the estimation procedure. Applying these methods to SLR observations results in a tool to assess existing satelliteand station-specific biases in SLR, and to estimate accurate relative coordinates between neighbouring telescopes, an important approach to validate local ties. However, it is very difficult to obtain simultaneous SLR observations in practice. While two SLR telescopes can in principle observe the same satellite simultaneously, when given the same schedule for the observation and under optimal weather conditions, simultaneous SLR observations from one telescope to two satellites are not feasible in practice. Therefore simultaneous double differences cannot be formed, and therefore, we have to resort to quasi-simultaneous observations.
The idea of quasi-simultaneity was first introduced in Pavlis (1985) as a tool to avoid the propagation of orbital errors, implicitly included in the range observations, into the estimated parameters. To achieve simultaneous SLR observations, the ranges of two stations were interpolated using cubic splines. Based on a large set of observations to the LAGEOS satellites, as well as multiple simulations, his thorough analysis showed that the use of differenced ranges mitigates orbital errors and achieves the estimation of baseline lengths with centimetre-level accuracy, even under the presence of satellite orbit errors on the metre level. Additionally, the quality of this method was demonstrated to be related to the geometry of the network and the observed satellites, with the best results obtained from passes parallel to the baseline direction. Similarly, Dedes and Müller (1989) used the idea of quasi-simultaneous observations from pairs of stations to obtain simultaneous range differences for the estimation of baselines independently of orbital errors. Their work used SLR observations to LAGEOS satellites, and the quasi-simultaneity is achieved through the interpolation of the observed laser ranges with Chebyshev polynomials, and a careful procedure for the elimination of outlier observations. They concluded that, given enough observations, it is possible to estimate baselines up to 1000 km with centimetre-level accuracy. A concept for double differences of SLR observations was later explored by Svehla (2018) with the baseline between the telescopes GRZL at Graz (Austria) and HERL at Herstmonceux (UK), and observing only two Galileo satellites. This approach forms SLR normal points at common epochs for both satellites, using first-order polynomials fitted to the normal points of the GRZL station, so that its SLR normal points are interpolated to epochs of the normal points of the HERL station, separately for each satellite and tracking pass. By fixing the IGS orbits, his study shows that the formed double differences allow the estimation of relative coordinates for the baseline GRZL-HERL at the level of 1-8 mm. In contrast, our approach for building differences uses the advantage of a very short SLR baseline, with a custommade quasi-simultaneity strategy to build single and double differences based on the linearised observation equations derived for the original ranges, without the need to interpolate either the normal points nor the original ranges. We only require a continuous representation of the satellite orbits, i.e. dynamic orbits. Additionally we try to maximise the number of observations by focusing on more than two satellites simultaneously, including SLR observations to both, GNSS and LEO satellites. From this perspective, this paper discusses the applicability and potential of the differencing approaches, namely single and double differences, for the estimation of geodetic parameters with SLR observations. These differences, together with the original ranges (zero differences), are used to get estimates of both, satellite-and station-specific error sources, so that systematic effects common to both stations can be identified at millimetre level. This approach, therefore, may potentially improve the processing of classical SLR observations of GNSS and LEO satellites and the estimation of accurate local ties. The rest of the paper is structured as follows: Sect. 2 outlines the concept of quasisimultaneity, which constitutes the basis of our work. The available data set and the processing strategy used for the estimation of parameters are described in Sect. 3. Before summarising the paper in Sect. 5, Sect. 4 discusses the results derived from the short SLR baseline in Wettzell regarding station coordinates and range biases, and the comparison of the derived baseline vector with the terrestrial local tie.

Concept of quasi-simultaneity
The fundamental tool to build single and double differences of SLR observations is the concept of quasi- Fig. 1 Concept of quasi-simultaneity for the differencing of SLR observations, together with the error sources targeted with these approaches. Two satellites (S 1 and S 2 ) in different orbits are observed. The quasi-simultaneity of the observations to satellite S 1 is given by δ t1 , while the quasi-simultaneity of the observations to satellite S 2 is given by δ t2 simultaneous observations. Two observations are considered quasi-simultaneous if they lie within a specified time window. This is, range i from telescope 1 to satellite k observed at time t i (ρ k 1 (t i )) and range j from telescope 2 to satel- where δ t is the so-called quasi-simultaneity, a fixed value. If t j and t i satisfy this condition, they are considered quasi-simultaneous epochs. Figure 1 shows the concept of quasi-simultaneity for an SLR baseline, where the time epochs for the observations from telescope 2 with respect to telescope 1 are t 1 and t 2 for satellite 1 and 2, respectively. Using this idea, three types of differences can be built, namely 1. Single differences from two telescopes to the same satellite: By observing the same satellite from two telescopes at quasi-simultaneous epochs, orbit and retro-reflector errors are strongly reduced. For the short baseline in Wettzell the tropospheric delays are eliminated as well (apart from the effect of the height difference between the telescopes). 2. Single differences from one telescope to two satellites: This definition includes the special case of observing the same satellite at two different epochs, as long as the quasisimultaneity condition is met: If observations are made from one telescope to two satellites at quasi-simultaneous epochs, the influence of instrumental biases caused by the station setup is practically eliminated. 3. Double differences: Similarly to the single differences from one telescope to two satellites, building double differences with the concept of quasi-simultaneous observations allows the special case of having double differences to the same satellite (k = l) at different epochs. Double differences are expected to eliminate both, satellite-specific (orbits and retro-reflectors) and station-specific (instrumental biases) errors.

Single-difference observation equation from two telescopes to one satellite
From a practical perspective, the SLR ranges are processed as in a usual SLR processing to obtain the so-called zero-difference linearised observation equations. These zerodifference linearised observation equations are later differenced to obtain the single-and double-difference systems of linear equations. The use of zero difference ensures that the original ranges are processed at the exact epoch, which in turn avoids the necessity to interpolate the two ranges to a common epoch. Formally, given an observation from the telescope m to the satellite k at time t i , ρ k m (t i ), the simplified SLR observation equation reads as ρ k m (t i ) = P k m + δρ k m,atm + δρ k m,rel + δρ m,sys + δρ k sys + ε k m with P k m the geometrical distance between the satellite and the station at the observation time t i , δρ k m,atm the delay (refrac-tion) in the troposphere, δρ k m,rel the relativistic correction, δρ m,sys delays in the laser system (among others a range bias B m ), δρ k sys the retroreflector correction and ε k m the measurement error. Using standard models to account for the troposphere, such as Mendes and Pavlis (2004), applying the relativistic correction, and using precise orbits for the satellites, the linearised zero-difference observation equations are given by T is the vector of the unknown parameters containing geocentric coordinates and range biases, A Z D m contains the partials of the observation equation with respect to the unknowns and each element of the reduced observation vector (observed-computed) the result of the a-priori values applied in the observation equation. When n 1 and n 2 observations are available for telescope 1 and 2, respectively, the system of linear equations is given by with the index m = 1, 2 for each telescope. In order to build the single-difference linear equation system from two telescopes to one satellite, we adapt the approach of Beutler et al. (1986Beutler et al. ( , 1987b to the SLR case using all the available observations at the two telescopes [b ZD 1 , b ZD 2 ] T which satisfy the quasi-simultaneity condition t i − t j ≤ δ t , and define the single-difference operator C sd as where C sd is a matrix with n m rows and 2n m columns. The single-difference system is obtained as At this point, the assumption of co-located telescopes, separated by a short baseline, comes into play. For short SLR baselines, quasi-simultaneous observations from two telescopes observe the satellite in (approximately) the same relative position, provided that δ t is small enough to account for the dynamics of the orbit. Under this assumption of having the same relative geometry for the two telescopes, the partials of the two telescopes become (nearly) identical, this is A ZD 1 = A ZD 2 . Thus, to build the single-difference matrix of partials A SD 1,2 , only the partials of the first telescope are used. With the help of the single-difference operator, A SD 1,2 is calculated as where the zero in the lower part of the matrix represents a matrix with the same dimensions as A ZD 2 and whose elements are zero. The C sd operator is a function of the number of satellites and the number of epochs at which they are observed, which is ultimately defined by the quasi-simultaneity δ t used. To take into account the correlations of the single differences, we introduce the weight matrix P sd (Beutler et al. 1986(Beutler et al. , 1987b) The solution of the least squares adjustment (LSA) for the single differences, using the matrices A SD 1,2 , P sd and b SD 1,2 , provides three relative coordinates ( ) and the relative range biases B. The quasi-simultaneity δ t and, ultimately, the single-differenced linearised equations system are strongly conditioned by the quality and the dynamics of the satellite orbits. For GNSS and LEO satellites this value should not be larger than 2 h, after which it is no longer possible to ensure that satellite positions have a similar error. A detailed analysis of the optimal threshold for the quasi-simultaneity is discussed in Sect. 3.4. With this assumption in mind, single differences of two telescopes to one satellite constitute a tool for the estimation of relative coordinates and range biases differences, while mitigating errors associated to the orbits of the satellites.

Double-difference observation equations
Similarly to the single-difference case of Sect. 2.2, we build the double-difference observation equation system using all the available observations at the two telescopes [b ZD 1 , b ZD 2 ] T which satisfy the quasi-simultaneity conditions t i −t j ≤ δ t and t k − t l ≤ δ t , and define the double-difference operator C dd as where the zero in the lower part of the matrix represents a matrix with the same dimensions as A ZD 2 , and whose elements are zero. Notice that this system of linearised equations does not contain the instrumental biases any longer, as they are reduced by the differencing at the same telescope. Similarly to Sect. 2.2, to take into account the correlations between the double differences, the weight matrix P dd The solution of the least squares adjustment for the double differences, using the matrices A DD 1,2 , P dd and b DD 1,2 , provides only the three relative coordinates ( X , Y , Z ). To facilitate the understanding of the construction of the single-and double-difference systems of linear equations, Appendix A shows an example of the technical procedure to build these systems, based on real data.

The co-located SLR telescopes at the Wettzell observatory
The SLR short baseline at the Geodetic Observatory Wettzell is realised by the Wettzell Laser Ranging System (WLRS) and the Satellite Observing System Wettzell (SOS-W), two optical telescopes operating at the wavelengths 532.1 nm and 1064 nm for WLRS and 849.8 nm for SOS-W. While WLRS has been contributing for more than 30 years to the realisation of the ITRF, SOS-W was placed in operation in early 2016 to cope with the need to track more satellites and distribute the workload. The horizontal distance between systems is ca. 58 m, and the difference in height is about 2.3 m (Fig. 2). The local tie vector between the two telescopes has been determined by terrestrial measurements and is continuously compared to the SLR-derived solutions. Moreover, an external calibration target for the SLR systems is available, providing sub-mm accuracy for the definition of the SLR reference points. More details about the operation and the performance of the two systems as well as the characteristics of the local surveys can be found in Riepl et al.

SLR normal points and meteorological data
The SLR processing using the differencing strategy is based on SLR observations provided as normal point (npt) files.
For the scope of this work, we restrict our study to only LAGEOS and GLONASS observations, and we analyse the behaviour of the SLR differences using the available ILRS data (Pearlman et al. 2019b, a;Noll et al. 2019) in the 2-year time interval from 01.01.2018 to 31.12.2019. Selecting those days where WLRS and SOW-S tracked the same satellites, yielded 402 potential sessions with data, from which 172 and 255 sessions contain sufficient GLONASS and LAGEOS data, respectively, to obtain differences. The npt files are converted into range observation files and processed using a dedicated project version of the Bernese GNSS software (BSW) v5.2 (Dach et al. 2015).
In addition, the meteorological data in the npt files are used to calculate corrections for the range observations associated with the influence of the troposphere, using the standard model of Mendes and Pavlis (2004). One of the main features of the short SLR baseline is the behaviour of the relative troposphere, where for such a small distance and small height change, meteorological data are expected to vary marginally between the two stations. Analogously to GNSS, this guarantees that the dry part of the relative tropospheric delays can be calculated with the use of standard models (Beutler et al. 1987a;Dilßner et al. 2008;Saastamoinen 1972). To test this assumption, the daily values of atmospheric pressure for the two telescopes were interpolated to common epochs, using the nearest neighbour method, and subsequently subtracted, to obtain a daily mean interpolated difference. Figure 3 (top) shows the behaviour of the atmospheric pressure for 03.-04.12.2019 (days of year -DoY-337 and 338), and the daily mean differences at common interpolated points (bottom), for the 2-year interval. These estimated differences reach 0.26 mbar for the days 03.-04.12.2009 and have an average value of 0.24 mbar over the two investigated years. Additionally, using the height difference between the stations, based on the local tie, a "theoretical" pressure difference was calculated using the simplified model of Saastamoinen (1972) 5.225 ). This theoretical difference amounts to 0.25 mbar. The difference between the daily mean and the theoretical value amounts to less than 0.01 mbar, which in turn corresponds to a negligible differential hydrostatic delay in zenith direction (Beutler et al. 1987a). These results support the assumption of a similar troposphere influence on both ends of the short baseline, and they make sure that we do not introduce errors into the baseline results by using erroneous meteorological data.

Processing strategies: zero-test and baseline estimation
The data analysis is performed with two types of processing strategies. In the first approach, the so-called zero-test, coordinates of the stations are derived from the ITRF2014 coordinates and velocities, together with the local tie, and are kept fixed for each processed day. Earth rotation parameters and satellite orbits are assumed fixed and supplied by the products of the Centre for Orbit Determination in Europe (CODE) (Dach et al. 2017), for GNSS, and the ILRS orbits (Pearlman et al. 2019b, a;Noll et al. 2019), for LAGEOS. The meteorological data provided in the npt files are used to correct the ranges for the delays induced by the troposphere, with the help of the empirical model of Mendes and Pavlis (2004). Finally, the range biases were not introduced. With this in mind, each station-satellite range is processed individually at each observation epoch. Under these assumptions, namely fixed station coordinates, orbits and Earth orientation parameters, and fixed troposphere parameters, the zero-test processing of a single station-satellite range does  Saastamoinen (1972) is given next to the average over the time series not estimate any geodetic parameters, but it is used to obtain the residual of the measurements with respect to the observation model, which in turn serves as a tool to evaluate the quality and noise of the raw observations. The residuals of the two telescopes, called zero-difference residuals, are then subtracted considering the concept of quasi-simultaneity of Sect. 2.1, to build single-and double-differenced residuals.
On the other hand, for the baseline estimation, Earth rotation parameters, pole coordinates and satellite orbits are once more assumed fixed and supplied by the products (CODE & ILRS), and the meteorological data provided in the npt files are used in the Mendes and Pavlis (2004) model, and again, no range biases were used during the process. In contrast to the zero test, for the baseline estimation, the station coordinates are assumed as unknown and calculated with a weighted least squares adjustment. Once more, each stationsatellite range is processed individually, and daily station coordinates of both stations are estimated. This process deliv-ers the zero-differenced linearised equation system A ZD m and its corresponding zero-differenced reduced observation vector b ZD m , m = 1, 2, for each telescope and each processed day. These linearised equations and reduced observation vectors are then subtracted, using the quasi-simultaneity concept described in Sect. 2.2, to obtain the single-and double-difference linearised equation systems and the singleand double-difference reduced observation vectors, which constitute the main element for our estimation of geodetic parameters. Figure 4 shows a summary description of our two main processing strategies.
Aside from these two differencing strategies, "standard" SLR solutions are calculated. Within these solutions coordinates and range biases per constellation are determined to benchmark the results of the differencing methods. Their daily normal equations are stored and combined to obtain parameters with the same validity as those obtained with the differencing approaches. Fig. 4 Summary of the investigated processing approaches. While the zero test does not perform the estimation of any parameters, thus delivering only zero-difference residuals wrt. the observation model, the baseline estimation performs the calculation of station coordinates, providing the zero-difference linearised equations and reduced observation vectors Fig. 5 Concept for the determination of the best quasi-simultaneous threshold. Using only one satellite, two different orbit solutions (from different processing centres) are compared. First, perfectly simultaneous single differences between orbit solutions are formed ( r 1 and r 2 , with their corresponding radial components rad 1 and rad 2 ). Then, double differences with all possible time steps are built, and the RMS value of the double differences at a certain time step is calculated

Selection of the threshold for the quasi-simultaneity
As mentioned before, the threshold for the quasi-simultaneity δ t is conditioned by the quality and the dynamics of the satellite orbits. To determine the impact of this limit on the double differences, we designed a test using satellite orbit solutions from two different processing centres, reflecting the orbital errors to be expected, to form single and double differences. The solutions involved are COM (Center for Orbit Determination in Europe -CODE-) and GFZ (GeoForschungsZentrum Potsdam) for the GLONASS satellites, and ILRSA (International Laser Ranging Service) and BKG (Bundesamt fuer Kartographie und Geodäsie) for the LAGEOS satellites.
To simulate the effect of the orbital errors on the doubledifference observations, we form the double differences of the two orbits according to Fig. 5

and the formula
Thereby δ t = t 2 − t 1 denotes the quasi-simultaneity value of the double-difference orbit errors. Using values of 1 s, 2 s,…, for the quasi-simultaneity we can compute the RMS of the radial orbit error double differences as a function of δ t . This procedure is done individually for each satellite, and we use "equal" instead of "smaller than" for the quasi-simultaneity. Of special interest is the radial component of the orbit, as it is closely related to the SLR ranges. As we perform the differences with exactly certain time step (QS = δ t ), instead of less than or equal to (QS ≤ δ t ), the Fig. 6 RMS of the double differences to a single satellite, in relation to the quasi-simultaneity, for the radial component. Top: Satellite GLONASS 13. Bottom: Satellite LAGEOS 2. These differences have been built using two different orbital products. Moreover, the bins for the RMS calculation are considered with quasi-simultaneity QS = d t instead of QS ≤ d t number of differences decreases when increasing the lapse for the difference. From each of these groups of differences (with quasi-simultaneity 1 s, 2 s, …) we take the RMS and use it as a metric for our analysis. Figure 6 shows these RMS values, for the radial component, in relation to the quasi-simultaneity, where the top plot shows the RMS values at each quasi-simultaneity for satellite GLONASS 13 and the bottom plot for satellite LAGEOS 2. In both cases the quasi-simultaneity was truncated at about the revolution period of each satellite: 223 min for LAGEOS 2 and 11 h 15 min for GLONASS 13. As expected, the RMS increases with the quasi-simultaneity. Moreover, it is evident that if all differences smaller or equal to a certain quasisimultaneity were considered, the RMS values would grow in a similar way. It is also visible that a threshold for the quasi-simultaneity of 0.5 h ensures a relatively small RMS value. In this way we guarantee that possible errors of the orbits are not propagating into the differences.

Analysis of the residuals from the zero test
With the zero-test strategy described in Sect. 3.3, the SLR data corresponding to the two telescopes in the time interval between 01.01.2018 and 31.12.2019 were processed in daily sessions, using the BSW software. Figure 7 shows the behaviour of the GLONASS and LAGEOS residuals as a function of the azimuth and elevation angle, for two of these days, namely 07.03.2018 (2018, DoY 184) and 03.12.2019 (2019, DoY 337), displayed by telescope. These residuals represent the level of noise of the observations plus orbit errors and possible instrumental biases. A summary of the statistics of these residuals is shown in Table 2. Aside from the difference in the number of observations and the presence of few outliers, the residuals of both days are among nominal values for the SLR technique, between −10 cm and 10 cm. However, there is a large contrast between the mean values obtained for the data set on the left wrt. the data set on the right in Table 2, mainly caused by the range bias of the station.
To further investigate the presence of instrumental biases associated with the stations, the single difference of the zero-difference residuals is built. For this, we considered a quasi-simultaneity of 2 h, and built the differences from two telescopes to one satellite. For the 03.07.2018, 513 singledifference residuals are formed, from which 276 correspond to differences to GLONASS satellites, and the remaining 237 are to LAGEOS. The same procedure produced 200 singledifference residuals for day 03.12.2019 with 98 differences for GLONASS and 102 for LAGEOS. These difference residuals are displayed in Fig. 8. As the range bias for the WLRS station was not applied during the processing, the singledifference residuals of day 03.07.2018 show a distinctive bias with a mean of −30.2 mm and −16.8 mm for GLONASS and  LAGEOS, respectively. However, the difference residuals of the second data set are almost unbiased, with mean bias values of 0.5 mm and −3.4 mm for GLONASS and LAGEOS, respectively. This difference in bias results might be due to the different target response of the two SLR systems wrt. the detectors.
To determine the behaviour of the instrumental bias and establish the time of possible bias changes, the time series of single-difference residuals is analysed. As seen in Fig. 8, when using a quasi-simultaneity of 2 h, it is possible to build multiple single differences per day. In light of the results of Sect. 3.4, where a threshold of quasi-simultaneity of 0.5 h was found as the best option to minimise errors from the orbits, we consider as unique "representative" for each day the differences with quasi-simultaneity ≤ 0.5 h. These daily values are displayed in the top plots of Figs. 9 and 10, where the mean single-difference residuals, together with the daily standard deviation, are displayed. While the standard deviation of the daily mean difference residuals is considerably larger for the GLONASS satellites than for LAGEOS, both time series display a distinctive break ca. 01.06.2019, confirming the results seen in Fig. 8. Therefore, we calculated the weighted mean value of these daily results imposing a break on DoY 156 of 2019, and the corresponding error for the mean is calculated. The middle plots of Figs. 9 and 10, show these mean values over the entire time series, where the vertical black line indicates the date of the break imposed for the calculation. Notice that the imposed break ensures the calculation of different values for the instrumental biases, 23.19 mm and 5.03 mm for GLONASS, and 12.56 mm and 0.98 mm for LAGEOS, with sub-mm errors for these mean values. If these instrumental biases are removed, the (unbi-ased) time series show a zero-mean behaviour (bottom plots of Figs. 9 and 10). These results demonstrate that the use of SLR single differences can be a useful tool for the accurate estimation of instrumental biases.
When now forming the double-difference residuals, the effects of these instrumental biases, and possible errors in the orbits of the satellites, are eliminated or mitigated. Figure 11 shows the double-difference residuals of days 03.07.2018 and 03.12.2019 allowing a quasi-simultaneity of 2 h. In these figures, the x-axis indicates the quasi-simultaneity allowed building the single-differenced residuals from the SOS-W telescope to two satellites, while the y-axis shows the quasisimultaneity allowed to build the single differences from the WLRS telescope to two satellites. Finally, the colour bar indicates the value of the residuals. With this approach we built 4596 and 3374 double differences, for day 03.07.2018 Fig. 9 Time series of the daily mean single-difference residuals of the zero test, at quasi-simultaneity of 0.5 h, for the GLONASS satellites. Top: Mean daily single difference at quasi-simultaneity 0.5 h, with its corresponding daily standard deviation. Middle: Mean daily single dif-ference at quasi-simultaneity 0.5 h, averaged over the two years to obtain the (relative) range biases. Bottom: Unbiased time series of daily mean single differences, after removing the range biases  and 03.12.2019, respectively. These residuals are no longer biased, with mean values of 5.1 mm for the former and 3.3 mm for the latter day. These small mean daily values, together with mm-level standard deviations, are an indication that SLR double-difference observations are potential candidates for the estimation of geodetic parameters.

Baseline estimation based on single differences
The first step of the baseline estimation process performs the determination of coordinates and instrumental biases based on single differences from two telescopes to one satellite, for each day individually. All the linearised equations satisfying the quasi-simultaneity condition t i − t j ≤ δ t with δ t = 0.5 hours, are stacked, as described in Sect. 2.2, and a weighted least squares adjustment is used to derive the corresponding 5 parameters, three relative coordinates ( X , Y , Z ) and two range biases, one for GLONASS and one for LAGEOS. For instance, the behaviour of the aggregated instrumental bias at the day 03.07.2018, shows the trade-off between a larger quasi-simultaneity value and While large values of quasi-simultaneity guarantee a larger number of observations, which in turn favour the estimation process and improve the formal errors, the errors caused by the orbits (intrinsic to our method) grow larger. Due to the large amount of data available on day 03.07.2018 (Fig. 12), the values at quasi-simultaneity 0.5 h, 21.9 ± 2.3 mm for GLONASS and 14.8 ± 2.1 for LAGEOS, are in agreement with the preliminary values observed for the same day in Sect. 4.1, during the entire session. Similarly, the difference of the estimates for the east, north, and up (ENU) components of the baseline with respect to the local tie, shows that the daily estimated values are in a close agreement with the local tie, with mean differences per component of [4.7, 4.5, −10.1] ± [1.2, 0.9, 2.9] mm, for east, north, and up, respectively. As a compromise between the number of observations and the influence of the orbital errors, we have selected those values obtained for a quasi-simultaneity of 0.5 h for each day's weighted least squares estimation. With this, we form the time series of daily estimates, three coordinate components, and two biases (GLONASS and LAGEOS), over the period of the two years of interest. Figure 13 shows these time series, with the time (year) of observation in the x-axis. This daily estimation addition- ally shows the need to differentiate between the instrumental biases before and after the observed break at day 01.06.2019. These results support the results of Sect. 4.1 regarding the changes in the instrumental biases and the need for a break in the data during the estimation process. For the coordinate estimates, the differences between the ENU estimates and the local ties are displayed. This daily analysis also facilitates the assessment of the performance of the differencing strategies regarding the number of formed differences and possible outliers.
Despite the high level of agreement between the daily estimates and the local tie, for a more rigorous and accurate solution, a complete (unified) weighted least squares solution is preferred. This least squares solution is computed by stack-ing all the daily linearised equations with quasi-simultaneity ≤ 0.5 h over the two years and producing a unique weighted least squares solution. When considering this approach with the available data, three cases are possible (Table 3): 1. Estimation of five parameters, three coordinates, and two range biases, where jump in the biases is not considered 2. Estimation of ten parameters, three coordinates, and two range biases before and the same 5 parameters after the observed bias jump 3. Estimation of seven parameters, three coordinates, and four range biases. These four biases correspond to two biases for GLONASS and two for LAGEOS, when considering the bias jump in the observation time series Figure 14 shows the residuals of the weighted least squares adjustment for these three cases. The time series with respect to the day of observation, shows an evident break for the residuals when estimating only five parameters (plot on the top-left), while the cases where breaks in the instrumental biases were considered, show unbiased zero-mean residuals in time. The right plot of Fig. 14 also shows that for larger values of quasi-simultaneity the residuals of the estimation begin to slowly increase, justifying the selection of a relatively small quasi-simultaneity value of 0.5 h.   Table 4 summarises the estimated values, together with the formal errors, of these three cases. As done before, the ENU estimates were subtracted from the local tie, so the values depict the difference of the SLR estimates with respect to the terrestrial measurements. While all the coordinate estimates are in a millimetre to sub-millimetre agreement with the local ties, there is an improvement in the formal errors for those solutions considering a break for the instrumental biases in the data. A tendency of lower values of the formal errors in the north component is seen in the three approaches, as a result of the south-north orientation of the baseline. As in the case of the GNSS-based solutions, the formal errors of the up component show the largest values. When comparing the values for the estimates of the instrumental biases, we see a clear difference between the approach with only five parameters and the remaining two cases. Looking at the solutions considering breaks for the range biases, we see that the biases for the GLONASS satellites change from 24.3 mm to 0.1 mm and 24.0 mm to 0.4 mm, for the 10 and 7 parameters solutions, respectively. These changes are statistically significant, when considering the formal errors of these values: 0.49 mm and 0.23 mm for the 10 parameters, and 0.50 mm and 0.23 mm for the 7 parameters. Similarly, with the formal errors found for the range biases, the change in the range biases is statistically significant and cannot be ignored. These results support the preference for the solution including a break for the estimation of the instrumental biases. As the coordinates of the stations are not expected to change due to the change in the bias, the solution with 7 parameters, where the components of the baseline vector are estimated once, is preferred. The correlations among the estimates for this approach are shown in the left part of Fig. 16. The small correlations ( ρ ≤ .25) among the coordinate components, zero correlation among the biases, and the expected large correlations between the height component and the range biases, support this selection.

Baseline estimation based on double differences
Although the estimated coordinates based on single differences from Sect. 4.2 show a millimetre agreement with the local tie, with relatively low formal errors, possible temporal variations of the instrumental biases are still present in the adjustment and may have an impact on the final solution.
To avoid these issues with the determination of instrumental biases and to reduce and mitigate the influence of other error sources, a double-difference weighted least squares estimation has been performed. Moreover, as discussed in Sect. 3.2, due to the short distance between the two telescopes, the influence of the troposphere on the SLR signal is expected to be the same (apart from the height difference). Therefore, with the use of double-difference observations, the tropospheric delay affecting the original observations is mitigated or strongly reduced. To perform the estimation process, we used the idea of Sect. 2.3 to select the daily linearised equations which satisfy the quasi-simultaneity condition δ t ≤ 0.5 h. This set of linearised equations is free of the influence of the troposphere, instrumental biases, and with a reduced influence of orbital errors. Instead of performing a daily estimation of parameters and then a calculation of the corresponding mean, we stack the linearised equations for one rigorous 2-year (unified) weighted least squares adjustment, where the only unknowns are the components of the baseline between the telescopes. Figure 15 shows the behaviour and distribution of the residuals of this weighted least squares adjustment, where the absence of any systematic influence of the instrumental biases is noticed. Moreover, for relatively short periods of quasi-simultaneity (≤ 600 s) the majority of the residuals of the double differences are between -10 mm and 10 mm. As for the single-difference case, we calculate the differences between the resulting coordinates and the local tie. The ENU components of this difference amount to: East : 0.7 ± 0.2 mm, North : − 0.9 ± 0.2 mm, Up : −0.6 ± 0.2 mm Furthermore, the correlations among these estimates (right part of Figure 16) show small values. The high level of agreement with the local tie and favourable correlations among the estimates, support the use of SLR double differences for the estimation of the local short baseline.

Conclusions
The single-and double differences of SLR observations for the short baseline in Wettzell have been investigated, and a novel approach to build differences of SLR observations has been developed, based on the concept of quasi-simultaneity. These differences are built with the linearised equations of the zero-difference processing. Therefore, the interpolation of the SLR ranges is no longer required, and system- atic errors common to the baseline can be targeted. The experiments over two years of SLR observations with the co-located telescopes in Wettzell, showed the advantages of the proposed method, namely the estimation of relative coordinates, suitable for the validation of local ties, and the estimation/elimination of instrumental range biases. Table 5 summarises the estimated values obtained using all the investigated methods, in addition to the "standard" results obtain for the same time interval (Bernese Zero Estimates). In particular, the single-difference rigorous solution with 7 estimates (for 3 relative coordinates and 2 relative range biases for LAGEOS and GLONASS) shows a much higher level of agreement with respect to the local ties, especially in the up component, with a difference of -2.1 mm against the 10.5 mm of the standard zero-difference solution. Furthermore, the formal errors of the ENU components for the single-difference method are significantly better than those of the zero differences. While there is no benchmark for the values of the range biases, the analysis of the formal errors of the single differences in relation to the standard method of the zero differences, shows a much higher performance of the former, with a sub-mm level for the errors. However, in contrast to the zero-difference approach, the single-difference strategy is only able to provide range bias differences between the stations. Nevertheless, to avoid the estimation of range biases, the double differences were studied. With the elimination of these biases, the estimation of relative coordinates, and therefore the validation of the local ties, is performed clearly more accurately. We found that the agreement of the relative coordinates and the local tie is within 1 mm for each of the ENU components, with corresponding formal errors in the sub-mm domain. Moreover, with a σ apost of 5.1 mm of the double-difference solution, against 7.1 mm for the solution based on single differences, the double-difference approach shows an additional and more effective reduction of systematic biases. The main improvement is observed in height component, which corresponds to the expectations, as the error sources mitigated with our approach (range biases, orbit errors, troposphere) are influencing the height the most. These characteristics support the usability of the proposed differencing methods for the estimation of geodetic parameters with a high degree of accuracy and their usability for the validation of local ties. Future activities will include the study of longer baselines, where we expect that the principle still works, however with a less pronounced reduction of orbit and tropospheric errors. org), CODE (http://www.aiub.unibe.ch/CODE), GFZ (https://www. gfz-potsdam.de/), BKG (https://www.bkg.bund.de), and the ILRS for providing the necessary orbital products and reference coordinates required for the realisation of this work.
Author Contributions I.H.P. and M.R. designed the study and developed the methodology. S.R provided inputs to refine the analysis. I.H.P developed the processing tools and analysed the data with the help of M.R. I.H.P wrote the manuscript with support of M.R. and S.R, and all authors read and approved the final manuscript.
Funding Open access funding provided by Swiss Federal Institute of Technology Zurich.

Data Availability
The SLR data sets used for this study are available from ftp://edc.dgfi.tum.de/pub/slr/data. The GNSS orbits used can be found at ftp.aiub.unibe.ch/CODE_MGEX/CODE/, and the LAGEOS orbits can be found at ftp://edc.dgfi.tum.de/slr/products/orbits/. The data sets generated during this study are available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy-right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Table 6 shows a sample of observations collected at the two telescopes in Wettzell for DoY 184, 2018. Considering a quasi-simultaneity condition of δ t ≤ 1 h, the single differences are built with the help of the C sd operator