Introduction

Satellite laser ranging (SLR) over the last decades has continuously aimed to improve accuracy and precision. Most SLR stations are, however, locally organized and apply improvements and changes (hardware, software or post-processing) to the system individually, which also influences data precision and has to be carefully considered e.g. by data analysis centers. Satellites hosting different arrangements of retroreflectors while being illuminated from different angles or laser beam polarizations can introduce changes in the distribution of detected photons, which then potentially affect the precision if not considered during post processing. To make use of these signature effects for improvements is a promising approach and was already demonstrated multiple times in the SLR community e.g. for geodetic satellites but also for navigation satellites (Otsubo et al. 2001) as well as by analyzing up to decades of SLR data (Sośnica et al. 2015). This research focuses on signature effects within Galileo data, which can be utilized to improve SLR precision further.

SLR measures the range to satellites equipped with corner cube retro reflectors (CCRs) by precisely timing the time of flight of the reflected laser pulses. Since the first measurements in 1965 the precision of SLR has improved significantly. With the use of single photon avalanche diode (SPAD) detectors, picosecond pulse length lasers and precise timing and calibration, single shot precision now reaches approximately 3 mm. Measured satellites include geodetic, scientific or Global Navigation Satellite System (GNSS) satellites.

Currently there are four GNSS constellations in the world. The United States’ Global Positioning System (GPS) was the first to be established, dating back to 1978. The Russian Global’naya Navigatsionnaya Sputnikovaya Sistema (GLONASS) was launched four years later. China started its Beidou system in 2000. The first launch of the In Orbit Validation (IOV) satellites of the European Galileo system was in 2011 (European GNSS Supervisory Authority 2020). India and Japan have their own regional systems based on satellites in geostationary and geosynchronous orbits.

The current generation of Galileo, BeiDou and GLONASS satellites carry a CCR panel hosting up to 90 reflectors. The arrangement of the CCRs on the panel varies from circular, ring shaped to hexagonal or rectangular arrangements (Pearlman et al. 2002). Galileo IOV-type panels include 84 CCRs arranged hexagonally with 33 mm minimum aperture, FOC (Full Orbit Capacity) satellite have only 60 CCRs with 28.2 mm minimum aperture and hence lower optical cross section (ILRS 2023). Depending on the relative motion of the satellite with respect to the observer, the center of the reflection pattern of each individual CCR is shifted toward the direction of flight of the spacecraft. This so-called velocity aberration depends on the orbit and—especially for lower orbits—also on the elevation of the satellite in the sky and can vary between 20 µrad and 50 µrad. Due to the larger distances, the velocity aberration for GNSS orbits is approximately constant in the order of 25 µrad (Degnan 2012). Typically, larger CCR diameters on the order of 3–4 cm are used for GNSS satellites. Based on diffraction laws, CCRs are chosen so that a large fraction of the reflected energy is directed at an angle close to the velocity aberration defined by the orbit. (Nugent and Condon 1966) Further modification of the far-field diffraction pattern (FFDP) (Degnan 2012) can be achieved by intentionally misaligning the back surface—typically dihedral angles (Chandler 1960; Thomas and Wyant 1977) of 2 µrad or more shift the reflected energy from the central to the outer regions of the FFDP. The rear reflective surfaces of GNSS satellites are mostly uncoated, resulting in a hexagonal reflection pattern (Arnold 2002). The front surfaces are either broadband-coated or coated specifically for SLR wavelengths.

SLR provides an independent, purely passive technique that can be used to study the systematic effects of different techniques (Luceri et al. 2019; Thaller et al. 2015), to validate orbits (Thaller et al. 2011; Strugarek et al. 2022; Zajdel et al. 2023), to generate combined and improved orbits (Bury et al. 2018), to analyze solar radiation pressure (Kucharski et al. 2017), or to perform time transfer (Samain et al. 2015).

SLR is a technique that reliably contributes (Sośnica et al. 2014) to the International Terrestrial Reference Frame (ITRF), with new developments towards a laser ranging reference frame (SLRF) recently announced (SLRF 2020). Important parameters connected to SLR are Earth rotation, gravity field, center of mass or polar motion. Scientific missions that rely on SLR data include oceanic motions, Earth topography, sea level variations, and Earth magnetic field or ice masses. Relativistic experiments (gravitational redshift, dark matter, Lense-Thirring effect) could be performed using different satellites (Delva et al. 2019; Ciufolini et al. 2019). Future trends in the laser ranging community include further increasing the repetition rate (Dequal et al. 2021; Wang et al. 2021; Long et al. 2022), possibly performing mm-precision laser ranging of satellites and space debris (Cordelli et al. 2020; Steindorfer et al. 2020) with one system.

The precision of SLR measurements depends on several factors such as sensor type and characteristics (single vs. multi-photon detectors), laser pulse width, event timing or calibration. In addition, the satellite signature and its analysis play an important role in further improving the precision further. SLR stations capable of detecting single photons with sufficient precision are even able to identify single CCRs of low flying satellites. Geodetic satellites have a specific return signal profile that is used to filter out the reflected photons that always correspond to the same region of the satellite—the leading edge.

Galileo orbits currently reach sub-cm-level accuracies (Dilssner et al. 2022), the panels also show different reflection profiles (Degnan 2012) which depend on the panel orientation (Steindorfer et al. 2019), the angle of incidence of the laser beam (Murphy Jr. and Goodrow 2013), laser polarization (Arnold 2002), thermal effects (Dell’Agnello et al. 2012; Boni et al. 2011), return rate and detector type (Kirchner et al. 1999; Kodet et al. 2009; Prochazka et al. 2020). Gaps in the individual FFDP distributions are counteracted by rotating the individual CCRs (clocking) (Dell’Agnello et al. 2011), which can lead to additional signal strength variations when looking at the overall pattern. Satellite signatures (Kucharski et al. 2015) can, therefore, influence the mean reflectance point of the Galileo panel, as the dataset of the entire panel is currently used in the post-processing of the laser ranging data. More photons from the front of the panel will therefore also shift the measurement towards closer ranges.

We investigate alternative measurement and normal point formation techniques. The first is the peak method, which is based on panel orientations where multiple CCRs are aligned at equal ranges with respect to the observing station. These so-called symmetry conditions can be accurately predicted by adding yaw steering to nadir pointing Galileo satellites (Montenbruck et al. 2015). For these passes the satellite signature is characterized and since the panel geometry is known, individual normal points can be formed for each of these ranges. The second method is the pattern correlation method (Brunelli 2009). The expected return pattern of the residuals is correlated with the observed one to identify the center position of the panel and to form the normal points. The two techniques are compared with respect to the current method based on example passes.

Materials and methods

In the following chapter, the basics of SLR, SLR residual simulations, Galileo symmetry simulations and the measurement setup will be summarized. Simulation and observation results from four different satellites are shown and referred to by Name (Type, PRN, Norad): Galileo 101 (IOV, E11, 37846), Galileo 103 (IOV, E19, 38857), Galileo 213 (FOC, E04, 41861), Galileo 219 (FOC, E36, 43566). However, all Galileo satellites up to Galileo 219 were observed during the observation campaign.

Satellite laser ranging—measurement

Satellite laser ranging provides precise range measurements to Earth-orbiting satellites equipped with CCRs. The worldwide network of laser ranging stations is globally distributed and coordinated by the International Laser Ranging Service (ILRS). Data centers collect the data, make it publicly available and publish frequently updated orbit predictions of the satellites. Different stations follow various approaches with respect to laser and detector technology. Low repetition rate stations fire laser pulses at about 10 Hz, each pulse with an energy in the order of 100 mJ, resulting in multiple photons reflected back from the target significantly favoring the closest CCRs. High repetition rate stations use repetition rates of about 2 kHz. Pulse energies are much lower, e.g. 400 µJ for the Graz SLR station. This increases the data rate and due to the statistical distribution of the reflected photons operating single photon return rate, multiple retroreflectors can be identified in the data. To reduce the divergence of the laser beam leaving the SLR station, beam expansion optics are used to expand the laser beam to 5 cm and beyond. Some stations use a separate lens-based beam expander. Common transmit-receive stations use the same large mirror-based telescope for beam expansion and reception of reflected photons reducing the divergence to values below the atmospheric seeing, which then becomes the limiting factor in terms of return rate. Pulse durations are typically in the order of 10 ps, resulting in optimal single shot measurement precision of up to 3 mm. Laser ranging with such a precision naturally requires precise timing of the events of the photon leaving the station and returning back after reflection (event timing). Detector technology includes microchannel plates, photomultipliers or single photon avalanche diodes. Depending on the signal strength reflected back to the station, receiving telescopes with different apertures are used (typically about 50 cm diameter). With larger apertures, precise pointing and increased laser power, lunar laser ranging is also performed by a few stations. Space debris laser ranging was first demonstrated in the years of the twenty-first century and relies on the diffuse reflection of photons over the entire space object by firing higher laser powers (> 10 W). Diffuse reflections provide an indication of the size and shape of the objects. For both laser ranging and space debris laser ranging, range variations can be used for spin period and attitude determination purposes.

Satellite laser ranging—residual simulations

Earth-orbiting satellites maintain a constant orientation with respect to a given coordinate frame. Nadir-pointing GNSS satellites are orientated with their navigation antenna pointing toward the center of the Earth (radial component), the second coordinate axis pointing in the direction of flight (velocity vector, along track) and the third component out of the orbital plane (cross product of velocity and along track component). For Galileo satellites, yaw-steering is implemented to allow the solar panels to face the sun. The yaw steering frame (YSF) is also nadir pointing but the axes \({{\varvec{e}}}_{x,YS},{{\varvec{e}}}_{y,YS},{{\varvec{e}}}_{z,YS}\) are calculated using the following formulas,

$${\varvec{e}}_{y,YS} = \frac{{{\varvec{e}}_{{{\text{sun}}}} \times {\varvec{r}}_{{{\mathbf{sat}}}} }}{{|{\varvec{e}}_{{{\text{sun}}}} \times {\varvec{r}}_{{{\mathbf{sat}}}} |}}$$
(1)
$${\varvec{e}}_{z,YS} = - \frac{{{\varvec{r}}_{{{\mathbf{sat}}}} }}{{\left| {{\varvec{r}}_{{{\mathbf{sat}}}} } \right|}}$$
(2)
$${\varvec{e}}_{x,YS} = {\varvec{e}}_{y,YS} \times {\varvec{e}}_{z,YS}$$
(3)

where \({{\varvec{e}}}_{{\text{sun}}}\) is the unit vector pointing from the satellite to the sun and \({{\varvec{r}}}_{\mathbf{s}\mathbf{a}\mathbf{t}}\) is the satellite position vector. Other reference coordinate frames of importance are the inertial Geocentric Celestial Reference Frame (GRCF), which is defined with respect to the stellar background, and the geocentric International Terrestrial Reference Frame (ITRF). The True Equator Mean Equinox (TEME) is an intermediate frame—the output of the Simplified General Perturbation 4 (SGP4) propagation used to obtain satellite predictions based on Two Line Element (TLE) data sets.

To acquire simulated SLR data the following steps are implemented:

  1. (1)

    Satellite orbit predictions: sources e.g. TLE, CPF (Consolidated Prediction Format, a data format used within the SLR community); time intervals: e.g. 60 s

  2. (2)

    Lagrange interpolation: saves computation time; allows shorter time steps; e.g. 1 s

  3. (3)

    CCR parameter definition: 3D position, normal vector, field of view; defined with respect to the center of mass or the geometrical center of the panel

  4. (4)

    Optional: define rotation parameters; spin axis, spin period, initial orientation; rotate satellite

  5. (5)

    SLR residual calculation: (Rotated) CCR positions are added relative to a reference coordinate frame (RCF); e.g. YSF, GCRS, nadir, …

  6. (6)

    Observed-Minus-Calculated (O-C) residual calculation: Difference between CCR ranges and reference orbit predictions

In the following the YSF will be used as the reference frame, since the goal of this study is to analyze the Galileo SLR residuals.

Galileo symmetry conditions

As the Galileo satellites change their yaw orientation throughout the orbit, special “symmetry conditions” occur with respect to the observing ground station, where several CCRs appear at the same distance (Fig. 1). The panel orientation relative to the ground station is characterized by two parameters.

Fig. 1
figure 1

Potential symmetry conditions on a Galileo IOV panel. Relative to the x-axis, symmetry angles of 0°, ± 30°, ± 60°, and ± 90° are indicated in different colors. The insert illustrates the tilt of the panel letting individual rows appear at constant range with respect to the observing station

Full size image

\({\alpha }_{{\text{symm}}}\): The symmetry angle that defines the orientation of the symmetry axis around which the panel appears rotated with respect to the observing station. The symmetry angle is defined as the angle of the rotation axis relative to the horizontal x-axis. For our simulations and measurements the symmetry angles of interest are \({\alpha }_{{\text{symm}}}=0^\circ , \pm 30^\circ , \pm 60^\circ , \pm 90^\circ\) as defined by the hexagonal arrangement of the Galileo panel.

\({\alpha }_{{\text{laser}}}\): The tilt of the panel relative to the observer is then characterized by a rotation around the previously defined symmetry axis. This tilt also describes the angle of incidence of the laser beam on the panel.

During the residual simulations \({\alpha }_{{\text{symm}}}\) is calculated by first projecting the vector pointing from the station to the satellite \({{\varvec{r}}}_{{\text{sat}}-{\text{sta}}}={{\varvec{r}}}_{{\text{sat}}}-{{\varvec{r}}}_{{\text{sta}}}\) onto the plane given by the Galileo panel. The projection \({{\varvec{u}}}_{{\text{n}}}\) of a vector \({\varvec{u}}\) onto the plane defined by the normal vector \({\varvec{n}}\) is given by

$${\varvec{u}}_{{\text{n}}} = {\varvec{u}} - \frac{{\user2{u n}}}{{\left| {\left| {\varvec{n}} \right|} \right|^{2} }}{\varvec{n}}$$
(4)

The angle between this projected vector and the yaw steering vector \({{\varvec{e}}}_{x,ys}\) now defines the symmetry angles \({\alpha }_{{\text{symm}}}\). The incident angle \({\alpha }_{{\text{laser}}}\) of the laser beam is simply the angle between \({{\varvec{r}}}_{{\text{sat}}}\) and \({{\varvec{r}}}_{{\text{sat}}-{\text{sta}}}\).

In other words, \({\alpha }_{{\text{symm}}}\) and \({\alpha }_{{\text{laser}}}\) are the azimuth and elevation angles of the SLR station as seen from the Galileo panel. For Galileo orbits (except for the two eccentric orbits) the incident angle cannot exceed \({\alpha }_{{\text{tilt}}}=12.5^\circ\).

Applying the result from the previous chapter, it is possible to simulate the SLR residuals for half a day for Galileo 103 (Fig. 2). The individual CCRs are color coded based on the 90° symmetry condition, which means that all CCRs of the same color for this specific symmetry angle will appear in constant range (Fig. 2a). The residual simulations (Fig. 2b) clearly show that the CCRs come together at a constant range at certain times, and then move apart. A − 90° symmetry condition occurs just before the second of day 60,000. A + 30° symmetry condition occurs at second of day 83,000 which causes differently colored CCRs to join together. The overall distribution of the residuals varies between ± 20 and  ± 45 mm which directly correlates to the incident angle of the laser beam. The elevation and azimuth of the satellite (as seen from the station) show two grayed out areas where the satellite is above the station horizon (Fig. 2c). Furthermore, the symmetry angle and the laser beam incident angle correlate well with the residual simulations (Fig. 2d). It can be seen that intermediate symmetry conditions also occur throughout the day.

Fig. 2
figure 2

Simulated SLR residuals for the IOV panel satellite Galileo 103 (NORAD: 38857). a Graphical display of the distribution of CCRs on the IOV panel. b Color coded simulated SLR residuals for 2021/02/14 UTC 11:00–24:00. c Azimuth and elevation of the satellite as seen from the station. d Azimuth (symmetry angle) and elevation (laser incident angle) of the station as seen from the satellite. In gray the parts of the pass where the satellite is above the horizon are highlighted

Full size image

In-orbit measurements—measurement setup

A tool has been developed to predict symmetry passes with respect to the Graz SLR station accordingly. Based on the latest predictions, the observer is warned by a popup message shortly before a symmetry pass occurs. The measurement goal was to cover a period of ± 5 min around the central symmetry epoch with continuous observations. Preliminary tests showed that due to the small separation between the CCRs symmetry passes with \({\alpha }_{{\text{symm}}}=\pm 0^\circ\) and \({\alpha }_{{\text{symm}}}=\pm 60^\circ\) are not resolvable within the SLR data. This can be explained by the 2–3 mm single shot RMS precision limit of our SLR station measurements. Therefore, our focus in the following was to measure \({\alpha }_{{\text{symm}}}=\pm 30^\circ\) and \({\alpha }_{{\text{symm}}}=\pm 90^\circ\) symmetry conditions. Between 2020/06/11 and 2021/04/24 a dedicated campaign was performed measuring more than 100 symmetry passes in linear and circular polarization. In the results section, a method is presented to calculate the incident angle of the laser beam from raw SLR data only. Detailed data analysis sheets are calculated for each of the recorded passes.

Remote location—measurement setup

A spare IOV retroreflector panel was provided on loan from ESA to perform measurements at a remote location on a mountain at 1200 m surface elevation, 32 km outside of Graz, with direct line of sight to the SLR station (Fig. 3c, image of the green laser beam originating from Graz SLR station). The reflector panel (Fig. 3a, covered within the black housing) was mounted on an Avalon M-Zero mount, which allows its azimuthal fork arm to be moved horizontally to mount the panel in direct extension of the azimuth axis. Thus, by changing the azimuth of the mount, the incident angle of the laser beam \({\alpha }_{{\text{laser}}}\) can be varied. Changing the elevation angle of the mount is related to the symmetry angle \({\alpha }_{{\text{sym}}}\). The zero alignment of the panel with respect to the SLR station (ensuring that the panel normal vector is facing towards Graz SLR station) was done by:

  1. (1)

    Aligning the panel with the mounted add-on telescope by analyzing an image of the laser beam on a CMOS camera (Fig. 3b, small white add-on telescope) with attached camera.

  2. (2)

    Fine-tuning the alignment by covering the CCR panel and optimizing the number of return photons from the reflection of mirrors inside the panel holder.

Fig. 3
figure 3

Measurement setup at the remote location. a overall setup of the panel (with cover) mounted on an Avalon telescope mount, the white add-on telescope with attached camera was used to perform the 0° pre-alignment of the panel w.r.t. Graz station. b The setup pointing towards Graz SLR station. c Close-up of the laser firing towards the remote location (green flash)

Full size image

The mount was rotated to various combinations of incident angle and symmetry conditions using a self-written ASCOM-based Python software (Fig. 4). Each position was held for a predefined time, while storing mount parameters and timestamps for later SLR return identification. In the results section, using only the raw data sets, the incident angle was calculated and compared to the mount values.

Fig. 4
figure 4

Explanation of the orientation of the panel at the remote location: Changing the elevation angle of the mount leads to a change of the symmetry condition; a change in azimuth leads to different incident angles

Full size image

ILRS Normal point formation

Within the ILRS, raw SLR data must be post-processed in a defined way. The data cleaning pipeline starts with residual calculation, coarse outlier removal, systematic trend removal, computation of fitted residuals and outlier identification based on a 2.2 sigma criterion. The process is performed iteratively until the accepted data set (corresponding to the reflected photons from the target) does not change anymore. Normal point formation involves dividing the data into bins, calculating the mean fit residual value and the mean epoch of the accepted data, finding the observation closest to the mean epoch and then calculating the normal point and RMS from the generated data. For Galileo satellites, the routine relies on the entire data set. In the results section of this paper, alternative ways of generating the normal points based on the symmetry conditions mentioned before are discussed.

Results and discussion

In the following chapter (1) the results of the measurements the spare panel placed at a location with direct sight connection, 32 km away from the SLR station, (2) the measurement campaign to Galileo satellites and (3) the normal point analysis methods are discussed.

SLR measurements, remote location

At the remote site, 32 km from Graz Observatory, 30° and 90° symmetry conditions were set with the panel mounted on the astronomical mount. The incident angle was varied from 5° to 15° in increments of 0.5° and between 10° and 13° in increments of 0.2°. It was found that the actual position of the mount can be slightly different from the desired value after the final rotation. Therefore, at the beginning and at the end of the measurement, the current mount tilt value was read out and saved to a file along with the current GPS time. The mount value will be referred to as the predicted incident angle \({\alpha }_{{\text{pred}}}\).

Each measurement was performed for 60 s, while keeping the return rate at a single photon level by adjusting the divergence of the laser beam. In the generated data sheets, the raw data set (Fig. 5, upper left) clearly shows patterns of increased data density at constant ranges. The range values are discretized to calculate a smoothed Gaussian residual curve (Fig. 5, lower graph red curve). For each 1 mm bin containing N range measurements yi, the Gaussian smoothing is calculated according to the following equation:

$$y_{{{\text{sm}}}} = y_{{{\text{arb}}}} \mathop \sum \limits_{i = 0}^{N} \exp^{{\frac{{ - 0.5\left( {y_{i} - y_{{{\text{cen}}}} } \right)^{2} }}{{\sigma^{2} }}}}$$
(5)
Fig. 5
figure 5

Datasheet for a measurement to the remote location at 30° panel symmetry. The incident angle of the mount was set to 12.77°. The calculated incident angle derived from the range measurements resulted in 12.13°

Full size image

The parameter sigma σ corresponds to the weighting/averaging of the values in the bin and determines the smoothing of the range residuals. For the following plots, sigma is set to 0.001. The derivative is taken from the smoothed data set (Fig. 5, green curve). Both plots are compared to a histogram calculated using the same bin width (Fig. 5, green curve). A Lomb analysis of the derivative of the smoothed curve is then performed to examine the periodicity with respect to the range between the different rows. The results of this analysis shows a clear peak at certain range. Knowing the geometry of the Galileo panel, the x-axis was translated into the incident angle (Fig. 5, upper right). The maximum of the curve is detected and then compared to the predicted incident angle—the angle given by the mount. For the 30° symmetry case, where the mount was moved to a nominal incident angle of 10.8°, the value returned from the mount after slewing was \({\alpha }_{{\text{pred}}}=10.76^\circ\). Data analysis of the measurement yielded a laser beam incident angle value of \({\alpha }_{{\text{calc}}}=10.13^\circ\).

The difference between the mount-based and measured angle of incident was calculated for all measurement sets. The smaller the incident angle, the closer the individual peaks become to each other, and at some point the Lomb periodogram can no longer detect a distinct peak due to the intrinsic measurement precision of the SLR technique. For the analyzed symmetry conditions, it was possible to reliably calculate an incident angle above 9° (90° symmetry: Fig. 6a, 30° symmetry Fig. 6c).

Fig. 6
figure 6

Remote measurements to a spare Galileo IOV panel: Read out mount incident angle [°] in dependence of the calculated incident angle [°] as derived from the SLR measurements to the remote location 32 km from Graz SLR station for 30° (a) and 90° symmetry angle (c). The negative offsets between calculated (from SLR measurements) and predicted (from the read-out mount angle) result from initial 0° misalignment of the panel (b, d). Above 9° incident angle the presented method reliably determines the incident angle

Full size image

The offset between the mount angle and the calculated incident angle was also compared to the nominal mount angle (Fig. 6b, d). Furthermore, the mean (\({\Delta \alpha }_{{\text{mean}}}\)) and standard deviation (\({\Delta \alpha }_{{\text{std}}}\)) of the offsets were calculated (Table 1). Both symmetry angles showed a distribution of incident angle offsets around a mean value, and no drift was found in the data set, indicating the quality of the method used. The measurement on the second day 2021/09/14 gave better results for the 30° symmetry condition than the one on the first day 2021/07/06. For the 90° symmetry condition, the standard deviations were comparable on both days.

Table 1 Ground measurement to a spare IOV panel 32 km away from the observatory for symmetry conditions 30° and 90°. A misalignment error was detected proving an accuracy of the method below 0.1
Full size table

On the second day, the mean value of both 90° and 30° symmetry conditions was significantly better. The detected offsets are expected to come from an initial misalignment during the 0° alignment routine using CCRs and mirrors. For the measurement on the second day, an RMS analysis (\({\Delta \alpha }_{{\text{rms}}}\)) of the returns coming from the full panel moving over 0° degree resulted in comparable 0° alignment offsets for both symmetry conditions.

Taking into account alignment inaccuracies, mount positioning errors and errors due to Gaussian smoothing and period detection, it is assumed that the incident angle can be accurate to within at least 0.1° under optimal conditions (Table 1).

SLR measurements, in-orbit

A large campaign analyzing more than 100 different Galileo symmetry passes was performed in 2021. As during the remote location measurements, the symmetry conditions of 30° and 90° were analyzed. As in Fig. 5, the range residual data was analyzed by forming a histogram, together with smoothed range residuals and their first derivative (Fig. 7).

Fig. 7
figure 7

Histogram distribution (blue), smoothed range residuals (red) and first derivative (green) for three Galileo satellites (Galileo 103, Galileo 101, Galileo 219). Lower right: Lomb analysis displayed for tested incident angles between 8° and 12°, showing a clear maximum at 10.29°, 10.65° and 9.38°. Galileo 219 shows a highly asymmetric reflection pattern

Full size image

Three example passes highlight the different reflectance profiles depending on the symmetry angle and the panel type. The Lomb analysis of the derivative shows clearly identifiable peaks corresponding to certain incident angles. For the two shown examples Galileo 103 (average return rate 1.8%) and Galileo 101 (2.7%), the differences from the TLE predicted angle are less than 0.1° (Fig. 7). The measured pass of Galileo 103 corresponds to the simulated residuals as shown in Fig. 2. The larger offset for the Galileo FOC satellite 219 is connected to the low return rate (0.4%), following from the smaller FOC CCRs and the overall lower number of CCRs on the panel. This results in patterns that are harder to discriminate. Some Galileo FOC measurements show an asymmetric distribution of returns with a significant reduction in signal strength, especially around the 90° symmetry condition. For example, at an incident angle of 9.8°. Galileo 219 shows two distinct peaks at − 1.1 cm and + 2.2 cm, and a valley at + 0.7 cm. Using 2.2 sigma clipping on the entire data set of such a distribution will offset the normal point. Such an intensity distribution is expected to result from the far field diffraction pattern (FFDP) of the individual CCRs: The reflective backside of Galileo CCRs is uncoated and relies on total internal reflection. Therefore, the resulting FFDP consists of six side lobes and a central maximum. The side lobes are hexagonally distributed at angles based on Airy diffraction (depending on the diameter of the CCR). For Galileo satellites at 23,200 km altitude, the velocity aberration of 25 µrad correlates with the maximum values in the FFDP; for IOV satellites, a further distribution of energy at certain angles is achieved by intentional spoiling (tilting) of the backsides of the CCRs. Due to the minimum values in the pattern of each individual CCR, an equal alignment of all CCRs of the panel would lead to certain observer-satellite geometries where the reflected intensity is very low. For this reason, Galileo panels are “clocked”—the CCRs are rotated in 15° (FOC) or 20° (IOV) increments—resulting in a more uniform overall FFDP. However, depending on how many CCRs of which clocking orientation are within a row of the symmetry condition, the reflected intensity can still vary, leading to different intensities within each of the peaks.

Using a fixed set of analysis parameters, a total of 98 passes were analyzed and corresponding data sheets were generated. The incident angle difference (calculated minus predicted values) was plotted for each individual Galileo satellite up to Galileo-219 (Fig. 8). It can be seen that the standard deviation of the offset is larger for FOC satellites (1 sigma = 1.26°, dashed red line) than IOV satellites (1 sigma = 0.53°, dashed blue line) than. Furthermore, there is no clear dependence on individual satellites. However, Galileo satellite 213 showed a slightly larger positive offset than the other satellites.

Fig. 8
figure 8

Offset between predicted and calculated incident angle for a different Galileo satellites, b as compared to the predicted incident angle (from TLEs), c Lomb power resulting from the incident angle calculation

Full size image

The predicted incident angles were compared with the incident angle difference (Fig. 8b) for IOV type (blue) and FOC type (red) panels and plotted separately for linear (x) and circular (o) laser light polarization. It can be clearly seen that the difference to TLE values is larger for IOV panels than for FOC panels. One reason for this is that for FOC satellites, due to the smaller distance between the CCRs, the different peaks appear closer together and are therefore more difficult to separate. Second, due to the smaller aperture and fewer CCRs, the overall return rate is lower, resulting in weaker signals. For IOV, only 1 out of 31 measured passes had an incident angle difference relative to the calculated value greater than 1°. For FOC panels, 23 out of 67 calculated angles had an offset greater than 1° relative to the predicted angle. For FOC satellites, a few incident angles are much larger than 12°. These are the two Galileo satellites in eccentric orbits.

The predicted incident angles are plotted against the Lomb-Scargle power (Fig. 8c). It is clearly visible that the Lomb power correlates with the incident angle for IOV and FOC panels. Due to the larger separation and the increased number of CCRs, the pattern identification is easier for IOV and favors larger incident angles.

The calculated incident angle can be compared with the satellite telemetry as obtained from the Attitude and Orbit Control System (AOCS) of the satellite. The AOCS keeps the spacecraft pointing at the Earth mainly based on infrared Earth sensors and the Sun sensors. The infrared Earth sensors do this by detecting the contrast between the cold of deep space and the heat of the Earth. These sensors provide a direct measurement of the pitch and roll angle of the satellite. The telemetry from the satellite system has an inherent noise better than ± 0.1° in roll and pitch (Dungate et al. 2005) while the error of the calculated tilt angle is ± 0.2° assuming a 1 mm error in the sub-peak estimation while neglecting datasets with low return rate.

By taking the measurements from Fig. 8 for the ± 90° symmetry condition it is possible to get a direct measurement of the roll angle of the satellite. Figure 9 provides a comparison of the determined roll angle against the roll angle from satellite telemetry data given by ESA. Though satellites Galileo 101 and Galileo 104 present a small + 0.1° offset during the observation period, the spread and noise of the available measurements does not allow for a clear estimation of the value.

Fig. 9
figure 9

Roll angle observed from telemetry and the calculated roll angle from the ± 90° symmetry condition

Full size image

The estimation of pointing errors or intentional GNSS bias can be of added value to explain the perturbing forces on the satellite (Bury et al. 2020). It can be concluded that the SLR attitude derivation technique can provide information about the attitude pointing accuracy with ± 0.2° accuracy which can be increased with an additional number of measurements.

Additionally, it can be concluded that if the signal strength from the satellite and the separation between the CCRs is large enough, the incident angle can be accurately determined. Future ultra-high repetition rate stations operating at low pulse energies with single photon return rates should be able to identify the satellite signature more clearly.

New normal point formation techniques

Based on the residual distribution which shows increased data density under symmetry conditions, two alternative normal point formation techniques are proposed: A peak-based method and a pattern correlation-based method.

Peak-based normal point formation: After orbit cleaning and smoothing the post-processed range data set, the peaks correlating to the center of each row can be determined (compare Figs. 5, 7). For each row, the data set is then cut based on the neighboring peaks. The upper and lower data cutoffs are defined by the average of the selected peak and the neighboring peaks. The cut dataset for each peak is then treated as regular Galileo SLR data (flattened by orbit-based parameter adjustment and 2.2 sigma clipped). This method allows of up to 11 normal points (90° symmetry) and up to 12 normal points (30° symmetry) to be generated simultaneously.

Cross-correlation method: A simulated residual pattern is calculated by forming a Gaussian distribution at the average position for each CCR. The width of the Gaussian function is based on the SLR precision of 3 mm. Adding the individual Gaussian led to overall simulated distribution of the whole panel (Fig. 10a, blue). The expected return pattern of the residuals is cross-correlated (Fig. 10b, blue) with the observed one (Fig. 10a, orange) to identify the offset (Fig. 10, red vertical dots) which will then be applied to correct the normal points. It is investigated here for the same passes.

Fig. 10
figure 10

Visualization of cross-correlation method. a The histogram of the measured data (green) is plotted together with the smoothed residuals (orange) and the simulated Gaussian kernel residuals. b The maximum of the cross-correlation between simulated and measured data corresponds to the offset applied to the normal point (dashed red line)

Full size image

The calculated normal points were plotted and compared by using the available Galileo Precise Orbit Determination (POD) orbits (Fig. 11). The original normal points using the full data set are plotted in black color, the peak-based normal points are highlighted in red, and the correlation normal points in cyan. In addition, the central normal point of the peak-based method is highlighted in blue. For 90° symmetry and an odd number of 11 rows, the central column is simply the 6th normal point. For 30° symmetry, the central normal point is determined by calculating the average of the two innermost rows. The difference from the original normal point formation method is calculated to highlight the offsets between the different methods, showing an improvement of the normal points. The histogram of the fitted raw POD range residuals is generated to show the distribution of the normal point relative to the peaks in the data.

Fig. 11
figure 11

Alternative normal point formation. top row: Galileo 103, 30° symmetry condition, day 045/2021; bottom row Galileo 102, 90° symmetry condition, day 091/2021. Left column: Full rate data POD residuals with normal points calculated with different methods. Center column: Histogram and smoothed residuals. Right column: Improvement of normal point w.r.t. reference orbit

Full size image

As an example, two analyzed IOV passes are shown. A significant offset (up to 4 mm) was found compared to the original normal point formation method. Variations in the data set can be caused by a different resulting FFDP of each row, thermal changes in the FFDP, from sensor specific noise distribution based on the return rate. The analysis showed that the normal points from both alternative normal point formation techniques were 0POD reference trajectory values. Pattern correlation and peak-based methods showed similar results, differing by up to 1 mm.

The current observation strategy (connected to ILRS recommendation) in Graz is to get approx. 1000 points per normal point to achieve optimum precision. After this threshold the observers are instructed to switch to the next target. The shown pass of Galileo 103 (Fig. 11) was continuously measured for approx. 700 s covering a bit less than 3 normal point durations. The number of points per normal point are given in Table 2, indicating that a sufficient precision can be reached even for each individual normal point. Considering future developments such as MHz laser ranging with laser powers of multiple Watts, will allow for targets signature-based data analysis, also for other SLR satellites.

Table 2 For the time interval t1, t2, t3 of three normal points, the number of points within full-data normal point is compared to the individual normal points of the peak based method
Full size table

Summary and conclusion

Galileo symmetry conditions occur whenever the panel is aligned such that multiple CCRs are aligned at equal distances with respect to the observing station. Characteristic patches of dense data form within the raw SLR data sets from which a histogram or smoothed range residuals can be calculated. The range differences between the individual rows can be used to calculate the angle of incidence of the laser beam with respect to the panel.

Ground-based measurements were performed to a spare Galileo IOV panel fixed on an astronomical mount on a hill 32 km outside the Graz SLR station. Different symmetry orientations were investigated and under optimal conditions it was possible to determine the incident angle down to 9° incident angle for symmetry conditions of 30° and 90° with an accuracy of less than 0.1°. The larger the incident angle, the easier it was to identify the rows and calculate the incident angle. Constant offsets to the nominal mounting values were found and initial panel misalignments of less than 1° were identified.

Simulated SLR residuals were generated, allowing accurate prediction of the symmetry conditions of existing Galileo satellites based on yaw steering laws. More than 100 passes to Galileo IOV and FOC satellites were measured and the laser beam incident angle was calculated from the measurements. It was found that the accuracy of the measurement correlates with the separation of the CCR rows. As mentioned before IOV panels use larger CCRs than FOC and hence the distance between the CCR is also larger. Accurate incident angle measurements could play an important role in validating on-board telemetry.

In-orbit based measurements were compared with satellite telemetry for the IOV satellites. The attitude was reconstructed with ± 0.2° accuracy.

Two new normal point formation techniques based on the symmetry condition data distribution were presented, showing an offset up to 4 mm compared to the method based on the full panel data set. The application of the pattern correlation method seems promising for general non-symmetry passes. More accurate normal points could eventually help improve SLR-based orbit predictions or the consistency between the GNSS-based and the SLR-based realization of the ITRF reference frame help to In the future, MHz satellite laser ranging (Wang et al. 2021) will provide high precision results with superior temporal resolution. Due to the higher power (40 W) with pulse energies strictly in the single-photon regime, the number of returns per second can be significantly increased. This allows for improved pattern identification and thus more advanced normal point formation methods to address the issue of satellite signature in the datasets. For the peak-based normal point formation, observation strategies could be adapted to favor symmetry conditions with improved normal point precision.