Accurate ground-to-ground laser time transfer by diffuse reflections from tumbling space debris objects

Satellite laser ranging achieves precise time transfer between ground and space via the time-of-flight measurements of ultra-short laser pulses. Dedicated time transfer experiments, such as Time Transfer by Laser Link [1] and the upcoming ACES (Atomic Clock Ensemble in Space (ACES) mission [2], provide time transfer uncertainties of less than 10 ps by physically performing an Einstein synchronization, where the round-trip time-of-flight of short laser pulses is used to separate the distance between the clocks from the actual time comparison process [3]. A corner cube reflector is required for the accurate estimation of the two-way path delay in this ground-to-space time transfer, since it provides a distinct and well-defined location for the light reflection. However, it only works for a single link between ground and space because the angle of the backscattered light covers only a few seconds of arc. While several different laser stations can operate simultaneously, it is not possible to directly transfer time by laser pulses from one station to the next. In order to compare two widely spaced ground clocks in common view of the satellite, each station performs a ground-to-space clock comparison from which the offset between the two ground clocks can be estimated in post processing. In the absence of an accurate clock in space, it is still possible to compare two remote ground clocks in common view of a satellite directly via diffuse laser pulse reflections from a space-borne target, but this creates two complications. Firstly, the laser pulse is reflected somewhere from the structure of the passive target so that the offset to the center of mass is not well defined, and secondly, some larger targets, like rocket bodies, may undergo a rather undetermined tumbling motion, which changes the mean position of the laser reflection to the center of mass constantly. So, the satellite dynamics have to be established at the time of the measurements. Figure 1 illustrates the general scenario of a laser time transfer link between two ground clocks and a tumbling rocket body as the space-borne diffuse reflector. For simplicity, we ignore atmospheric and internal ranging system delays in this section here. We include these effects in section 4. The short laser pulses from station A are back reflected to A so that the round-trip delay can be established. Due to the diffuse nature of the reflection, a small portion of the optical signals are also detectable at the location of station B [4]. Figure 2 illustrates this synchronization process between the timescale of station A and station B for the case where only station A fires the laser. The round-trip laser ranging measurements are used to establish the correct orbit, the rocket tumbling motion and the time of the laser pulse arrival at the space target. With these quantities all established, we can compute the arrival time t_s = (t_A2 − t_A1)/2 + t_A1 of the laser pulse at the space target with respect to timescale A. The time offset between the two timescales then is Δt = (R_A + R_B)/c − (t_AB − t_A1). The same applies for pulses sent from station B, which can be used to provide a closure for the clock comparison and reduce the measurement error.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Rocket bodies, for example, have a cylindrical shape with a length of the order of 10 m and a diameter just under 3 m. With the laser beam illuminating the entire target, we have to establish the relative orientation and effective shape of the target structure with respect to the observing station as a function of time. Once the instantaneous orientation of the space object is constrained, we must set up the corresponding delay model, which references the mean area of reflection to the center of mass of the target. We require this step since the mean reflection points between the signals of the two time transfer stations do not necessarily coincide. This ultimately defines the achievable accuracy that we can obtain for the entire time transfer process. In this paper, we explore this scenario for a single transmitting station.

In order to reduce the complexity of the problem of finding the correct time delays between the two ends of the laser link, we restrict our discussion to cylindrical rocket body targets. Each target is characterized by its shape, namely the length l and diameter d. The surface has reflectance parameters α , which normally are not uniform across the entire surface. While the center of mass (M_c) is usually located on the main figure axis, it does not coincide with the geometrical center of the entire cylinder. We assume these investigated targets exhibit a tumbling motion with a period T_s either in the clockwise (CW) or counter-clockwise (CCW) direction. The plane in which this tumbling motion is observed is very stable over a short time and can be characterized by a normal vector, whose orientation in an inertial frame of reference is defined by angles for right ascension θ_RA and declination θ_DEC. Figure 3(a) illustrates this. Figure 3(b) illustrates the observation vectors for the mono-static laser ranging and the photometry. Only in some cases is it possible to retrieve the dimensions l and d of the rocket body from publicly available sources with sufficient accuracy in a way that is consistent with the ranging measurements. A small error in the estimation of the plane of the tumbling motion can generate a projection error that makes the rocket body shorter than in reality. The open ends of the cylinder have similar effects. As an example, we show the laser range residuals for the Soyuz rocket body (NORAD ID: 39679) taken on 29 December 2015 in figure 4(a). The orbit altitude of the object in this pass is around 450 km and the observation lasted for about 200 s. The reported diameter for the Soyuz rocket is 2.66 m [13] and the corresponding length is between (6.70–9.40) m, depending on the actual rocket variant. From figure 4(b), we obtain a peak in the histogram of the range residuals, which indicates a radius of approximately 1.34 m. The other two peaks are somewhat less obvious and represent the observed distance between each end and the center of mass. This corresponds to a total length of 3.28 m + 6.24 m = 9.52 m, which is in reasonable agreement with the expected value of length. However, we believe some small estimation errors are still present, which can be reduced by observing more passes of this target. Figure 4 also reveals that the reflectivity of the two ends is larger than that of the side. Furthermore, the reflectivity distribution of the side may not be uniform as the variable density of the echo distribution indicates. In order to reduce the error margin further, we include a second type of observation, namely photometric measurements of solar reflections in the analysis. This is shown in figure 3(b). Continuous time-of-flight measurements while the target is visible from the ground station contain both the effective albedo and the tumbling motion imprinted on the range measurements. In order to resolve the ambiguities, the additional observation of the photometric object brightness provides complementary information, because the two signal sources (i.e., the solar radiation and the laser beam) have a different angle of incidence. We combine the two methods to estimate the center of mass position of a CZ-3B rocket body (NORAD ID: 38253). Figure 5 represents the observed laser ranging residuals (a) and the corresponding light curve observation (b). The orbit altitude of 38253 is around 2000 km and the observation session lasted for about 1000 s. From the official website [14], the diameter of this rocket body is given as 3.00 m and the length is 12.375 m. The finally obtained center of mass position was estimated as 5.1 m + 7.3 m = 12.4 m. The distribution of the laser ranging echoes shows that the reflectivity of the sides is larger than that of the two ends. From the shape of the light curve, we may speculate that the reflectivity of the two ends differ from each other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The tumbling motion of the rocket body causes the most significant variation in the pathlength delay of the time transfer application. Although the orbit-related delay variation is much larger, it changes in a very predictable manner that can be modeled very well. In order to capture the variable delay of the tumbling motion, the tumbling axis orientation and the phase of the oscillatory tumbling motion have to be extracted from the available measurements. Applying a Lagrange interpolation and a fast Fourier transform suitable for unevenly spaced samples [5] yields a spectrogram from which the tumbling frequency can be extracted. Applied to the laser ranging observations shown in figure 4(a), we obtain the corresponding spectrum as shown in figure 6. This results in an apparent tumbling period of approximately 1/0.088 57 Hz ≈ 11.3 s from peak 3. Peaks 2 and 4 are the harmonics and sub-harmonics of the tumbling rate, which stand out, because the center of mass is not in the middle of the main symmetry axis of the cylindrical target body. Peak 1 is ultimately related to the length of the dataset and has no relevance here. The same method applied to the target with NORAD ID: 38253 yields an apparent tumbling period of 121.7 s for the laser ranging data. Since the observation time is much longer than the previous dataset, the observations also exhibit a considerably improved signal-to-noise ratio. The apparent tumbling motion is also encoded in the photometric data (figure 5(b)) and the spectral analysis yields 118.1 s for the apparent tumbling period derived from the light curves. In order to obtain the true tumbling period of the target, we have to account for the projection of the plane of observation onto the plane in which the tumbling motion of the target takes place. A similar geometric consideration is required for the estimation of the true tumbling period from the photometric data, where the projection of the plane containing the Sun, the point of reflection of the satellite target and the detector at the ground station is evaluated.

Zoom In Zoom Out Reset image size

We explore the relationship between the true and the observed tumbling period by adjusting a target reflectance model to the measurements. The simulation model is also required for the correct estimation of the tumbling axis orientation (θ_RA, θ_DEC). Due to some geometrical symmetry effects in the observation constellation, the simulation of the reflectance model has also to resolve some geometrical ambiguities.

Before we discuss the details of the reflection model, we need to distinguish between two types of laser ranging signal detection, namely the multi-photon and the single-photon case. The former is characterized by the fact that the returning light level generates a voltage on the detector that exceeds a predefined signal threshold. That means that the timer responds to the leading edge of the signal. The latter case applies to very weak return signals and is typical for the diffuse reflection from a rocket stage. The detector responds to the detection of a single photon in a totally statistical manner [6]. Over many repetitions of the laser ranging measurement, this mode of detection maps out the geometrical depth of the target. All our time transfer measurements were conducted in the single photon mode due to the low scatter cross sections of the targets.

In order to compute the total reflection from a cylindrical target, the surface of the rocket body is divided into many square tiles. The bidirectional reflectance distribution function (BRDF) model from Ashikhmin and Shirley [7] is applied to compute the reflection of each surface tile. Figure 7 shows the applied geometry for the light reflection for one tile. The BRDF model is composed of a specular reflection component ρ_s and a diffuse reflection component ρ_d:

$$\begin{equation}{\rho }_{\text{f}}={\rho }_{\text{s}}+{\rho }_{\mathrm{d}}.\end{equation} \tag{ 1 }$$

With L , H , N and V defined as shown in figure 7, the two components of the light reflection are given as [7]

$$\begin{equation}{\rho }_{\text{s}}=F\frac{\left({n}_{\text{u}}+1\right)}{8\pi }\frac{{\left(\boldsymbol{N}\cdot \boldsymbol{H}\right)}^{{n}_{\text{u}}}}{\left(\boldsymbol{H}\cdot \boldsymbol{V}\right)\mathrm{max}\left(\left(\boldsymbol{N}\cdot \boldsymbol{L}\right)\left(\boldsymbol{N}\cdot \boldsymbol{V}\right)\right)}\end{equation} \tag{ 2 }$$

and

$$\begin{align}\hfill {\rho }_{\mathrm{d}}& =\frac{28}{23}\frac{{c}_{\mathrm{d}}}{\pi }\left(1-{c}_{\text{s}}\right)\left(1-{\left(1-\frac{\boldsymbol{N}\cdot \boldsymbol{L}}{2}\right)}^{5}\right)\hfill \\ \hfill & \quad {\times}\enspace \left(1-{\left(1-\frac{\boldsymbol{N}\cdot \boldsymbol{V}}{2}\right)}^{5}\right).\hfill \end{align} \tag{ 3 }$$

In this representation, F describes the Fresnel fraction for which we use a simplified expression as provided in [8]:

$$\begin{equation}F={c}_{\text{s}}+\left(1-{c}_{\text{s}}\right){\left(1-\boldsymbol{V}\cdot \boldsymbol{H}\right)}^{5}.\end{equation} \tag{ 4 }$$

The parameter c_s is the material reflectance for the normal incidence. 0 < c_d < 1 is the diffuse albedo of the surface. n_u stands for the luminance caused by specular reflection. As n_u gets larger, less energy is lost by the specular term [7]. In general, these three parameters vary for different parts of a space object. They are the reflectance parameters for our work.

Zoom In Zoom Out Reset image size

For the simulation of mono-static satellite laser ranging, the light source vector and the viewer vector are antiparallel. We refer to this vector as the laser vector. Its direction is from the center of mass position of the rocket body to the reference point of the observatory. Furthermore, we assume that the whole object is illuminated by the laser beam and only back-reflections are considered. The retro-reflected pulse-width is spread in time due to the geometrical depth of the rocket body. Each echo corresponds to the location of one surface tile under ideal conditions. We assume that the detector is a highly sensitive single photon detector and the detection probability of one tile is proportional to its surface reflection R_{f_LR}:

$$\begin{equation}{R}_{\text{f}_\text{LR}}\propto {\rho }_{\text{f}}{\left(\boldsymbol{N}\cdot \boldsymbol{L}\right)}^{2}.\end{equation} \tag{ 5 }$$

The photometric illumination of the object by the Sun is not pulsed, therefore we do not have a temporal spread of the signal. Light curves are generated from the sum of the reflected photon flux by all the tiles. For all rocket body tiles, the source vector points to the Sun, whose norm is defined as R_sun and whose direction is from the center of mass position of the rocket body to the Sun. The viewer vector is the same as the laser vector in the paragraph above, whose norm is defined as R_obs. The surface reflection of each tile is

$$\begin{equation}{R}_{\text{f}_\text{LC}}\propto {\rho }_{\text{f}}\frac{\left(\boldsymbol{N}\cdot \boldsymbol{L}\right)\left(\boldsymbol{N}\cdot \boldsymbol{V}\right)}{{R}_{\text{sun}}^{2}{R}_{\text{obs}}^{2}}.\end{equation} \tag{ 6 }$$

From the satellite reflectance model, we obtain a theoretical laser echo delay distribution for the case of ranging geometry and simultaneously a temporal brightness distribution for the case of photometric observation. Please note that both ρ_f and L are related to the relative position of the rocket body and the observatory with respect to each other. By numerically adjusting the dimensions of the rocket body, the center of mass position, the reflectance parameters, the tumbling period and tumbling axis orientation in space as well as the evolution in time, we eventually obtain a physical parameter set of the observed object that provides good agreement for the measurements. This is illustrated in figure 8 for the case of the laser ranging observation. From the adjustment procedure of this measurement, we obtain two very similar solutions of the range residual distribution. This is caused by the geometrical symmetry given by the projections of the tumbling motion onto the ranging vectors, namely (θ_RA = 100°, θ_DEC = 70°) and (θ_RA = 22°, θ_DEC = −33°), basically a mirror image in the projection. In order to resolve this ambiguity, we include the simultaneously observed light curves, exploiting the different angles of illumination. Figure 9 depicts the result. Good agreement is obtained for the tumbling axis orientation of (θ_RA = 100°, θ_DEC = 70°) only, therefore it is regarded as the correct solution. While the tumbling period for each target is constant along its orbit and the orientation of the tumbling axis remains fixed in space, the actually observed tumbling period by the laser station will differ because of the orbital velocity of the satellite with respect to the ground station. If the target tumbling motion is co-rotating with the orbit motion, we observe a slightly longer period and obtain a shorter period if the orbit motion and the tumbling motion are anti-rotating with respect to each other. This property, together with the different light source incidence angle from laser ranging and photometric observations, are important to constrain the geometry for the time transfer application. The result of the modeled observation is an envelope curve, which indicates the maximum and the minimum distance between the laser station and the passive reflector target. Most of the time, only one distinct solution remains. However, in some cases, like the one shown in figure 8, a small ambiguity remains, albeit with a negligible range error. It is important to point out that this only applies to the monostatic observation case. In the desired bistatic time transfer application, there is no such ambiguity due to the additional range constraint from station B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use the reconstructed target structure envelope, as shown in figure 8, to remove detector noise from the observations and to classify the targets into two distinct categories, namely rocket bodies, which exhibit a high reflectivity from the ends of the cylindrical body and those who predominantly reflect from their sides. Both types of targets have to be treated differently for the estimation of the time transfer delay. Targets that mainly reflect from the sides show a larger scatter and the mean location of the reflection is at the ring around the center of geometry of the rocket body. The histogram has an almost Gaussian distribution. Those rocket bodies that mainly reflect from their ends are affected far more by the tumbling motion. Consequently, the histogram of the mean reflection function shows a pronounced non-normal distribution. To establish the range and tumbling motion-corrected mean clock delay to the orbiting target correctly, we apply a 2.5σ filter criterion iteratively to the histogram of the measured range residuals in order to deal with the skewed signal distribution. This is illustrated in figure 10 for both target types.

Zoom In Zoom Out Reset image size

The ground-to-ground optical time transfer facilitates the accurate synchronization of two remote clocks with no direct line of sight. The diffuse reflection of many short laser pulses from an orbiting target is used to transfer the time. Figure 2 illustrates the process in a simplified sketch. The clocks in stations A and B can be synchronized by timing the departure and arrival times of short laser pulses emitted from station A. Since the orbiting target produces diffuse reflections, their arrival times are also timed at station B in its own timescale. From the round-trip measurements of station A, we determine the true position of the target along its orbit, as well as the target delay model, to obtain the instantaneous orientation of the tumbling body in space. The distance between the reflection zone on the rocket body and stations A and B is referenced to the invariant point of each laser ranging system, which are the intersection points of the horizontal and vertical telescope axes of both ranging systems. This requires that the signal delay between the actual clock and this system reference point has to be established independently by the system calibration for the outgoing laser pulse as well as for the incoming reflection [9]. Finally, we have to apply a correction for the integral effect of the refraction in the atmosphere, which amounts to something around 33 ns for a round trip at low elevation angles and drops to about 9 ns at zenith. For this, we apply the Mendes et al refraction correction model [10], which is routinely applied to all SLR measurements. It is based on local measurements of atmospheric pressure, temperature and humidity. From the round-trip measurements of station A, we have

$$\begin{equation}\begin{aligned}\hfill {\phi }_{\text{sat}1}\left(t\right)& =\left({t}_{\text{A}2}-{t}_{\text{A}1}\right)-2{R}_{\text{A}}^{{\ast}}\left(t\right)/c-{\Delta}{R}_{1}\left(t\right)/c\hfill \\ \hfill & \quad -\left({\delta }_{\text{TX}1}+{\delta }_{\text{RX}1}+{\tau }_{\text{up}1}\left({\theta }_{\text{e}1}\left(t\right)\right)+{\tau }_{\text{down}1}\left({\theta }_{\text{e}1}\left(t\right)\right)\right).\hfill \end{aligned}\end{equation} \tag{ 7 }$$

By timing the arrival time of the laser pulses from stations A at station B, we have

$$\begin{align}\hfill {\phi }_{\text{sat}2}\left(t\right)-{\Delta}t& =\left({t}_{\text{AB}}-{t}_{\text{A}1}\right)-\left({R}_{\text{A}}^{{\ast}}\left(t\right)+{R}_{\text{B}}^{{\ast}}\left(t\right)\right)/c\hfill \\ \hfill & \quad -{\Delta}{R}_{2}\left(t\right)/c-\left({\delta }_{\text{TX}1}+{\delta }_{\text{RX}2}+{\tau }_{\text{up}1}\left({\theta }_{\text{e}1}\left(t\right)\right)\right.\hfill \\ \hfill & \quad \left.+\enspace {\tau }_{\text{down}2}\left({\theta }_{\text{e}2}\left(t\right)\right)\right)\hfill \end{align} \tag{ 8 }$$

where the parameters are as indicated in figure 2. The black dots on the respective time axes indicate the actual instants of timing, the respective δ values correspond to the transmit and receive delays (i.e., the delay correction to the system reference point) and all of these values are constant for the entire satellite pass. The outgoing laser pulse is typically timed before it passes the reference point, while the situation is the other way around for the inbound beam. All other parameters vary over the time of observation—some as a function of orbit position and tumbling motion and others as a function of telescope pointing angle θ_e(t): $2{R}_{\text{A}}^{{\ast}}\left(t\right)+{\Delta}{R}_{1}\left(t\right)=2{R}_{\text{A}}\left(t\right)$, ${R}_{\text{A}}^{{\ast}}\left(t\right)+{R}_{\text{B}}^{{\ast}}\left(t\right)+{\Delta}{R}_{2}\left(t\right)={R}_{\text{A}}\left(t\right)+{R}_{\text{B}}\left(t\right)$. ${R}_{\text{A}}^{{\ast}}\left(t\right)$ and ${R}_{\text{B}}^{{\ast}}\left(t\right)$ are calculated from the predicted orbit. ΔR₁(t) and ΔR₂(t) refer to the orbit corrections from the predicted orbit to the real orbit. The correction for the refraction delay is contained in the τ_x (θ_e(t)) parameters, separate for the uplink and the downlink. ϕ_sat1(t) and ϕ_sat2(t) ultimately correspond to the target delay model, as depicted in figure 10. With all corrections applied, we are now in the position to synchronize the two station clocks. An example with good fidelity is shown in figure 11, where the clock offset between the Graz Laser Ranging Facility and the Geodetic Observatory Wettzell is obtained as Δt = 2023 ± 3 ns, where the uncertainty is given as the 1σ error. For completeness, we mention here that the errors in meteorological correction and station position taken together are in the range of approximately 30 ps, which is too small to matter in the context of this paper. The observed target (NORAD ID: 23088), however, did not quite fit into either of our two identified target categories, so the scatter is somewhat larger than expected, owing to the limitation of the target delay model. The analysis of several mono-static rocket body observations, however, indicates that this technique of optical time transfer can achieve a time transfer accuracy well below 1 ns. When both stations operate as transmitters and receivers at the same time, as sketched out in figure 1, we expect further improvements as this adds a third independent angle of observation for the target delay model.

Zoom In Zoom Out Reset image size

The delay-compensated distribution of highly stable optical clock frequencies over fiber networks in Europe as well as in other parts of the world has considerably transformed the field of time metrology. It will allow for the comparison of two clocks with a stability of one part in 10¹⁸ over a distance of more than 1800 km [11, 12] when such high-quality clocks eventually become commonly available. However, while the frequency can be compared at a yet unprecedented level of resolution, the phase of the clock signal is lost in the process. This means that time cannot be transferred at this level of accuracy unless the exact delay of the fiber link is precisely known, including any variation of the transfer delay over time. We note, however, that time and space are not separable. Optical two-way time transfer, based on free space time-of-flight measurements, is a mature technique [1] and provides a well-defined geometrical path over the entire measurement process. Together with a high-quality clock in orbit, a global and stable distribution of time at the level of a few ps is feasible. However, with the ACES clocks [2] still waiting for launch, this technique is currently not available.

Our technique of using old rocket bodies as zero time delay reflectors allows for the transfer of the clock phase from one laser ranging station to the next under common view conditions. The only remaining uncertainty in this process is the determination of the shot-by-shot reflection point on the tumbling space target. For two different types of rocket bodies, we have established a target delay model, which statistically reduces this target-related delay uncertainty for the time transfer to values below 1 ns for the round-trip measurement in mono-static operation. For the actual transfer of time between the Geodetic Observatory Wettzell and the Graz Laser Ranging Facility, we obtained a 1σ statistical uncertainty of 3 ns, with most of this uncertainty being related to the fact that the observed target did not quite fit into the two target categories that we have modeled. For this particular time transfer experiment, we also did not have any photometric measurements of the target available, which prevents us from improving the target delay model for this experiment. For the future, we believe that a ground-to-ground time transfer over diffuse target reflections with an uncertainty of much less than 1 ns is possible on selected and well-behaved targets when both stations operate in a fully reciprocal manner, as described in the introduction.

We express our gratitude to Professor Urs Hugentobler of the Research Unit Satellite Geodesy at TUM for discussions and support. This work is supported by Sino-German (CSC-DAAD) Postdoc Scholarship Program, 2019(57460082).