ATLAS: Orbit Determination and Time Transfer for a Lunar Radio Navigation System

2.1 Orbit Geometry

We assume that the lunar constellation consists of four satellites in different ELFO orbits, with the orbital parameters reported in Table 1, in line with the parameters used for the Moonlight project (Melman et al., 2022). These orbits allow coverage of the southern polar region of the Moon (Figure 1).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Orbital geometry at the reference epoch for the LRNS constellation in a Mooncentered inertial reference frame (axes indicated by black arrows).

The colored points indicate the four spacecraft: 1 in blue, 2 in orange, 3 in green, and 4 in purple.

2.2 Ground Segment

Concerning the ground segment, we have assumed that the small antenna dishes are located close to stations of the European Space Tracking (ESTRACK) network, namely, Cebreros (Spain), New Norcia (Australia), and Malargüe (Argentina). The primary rationale behind choosing these three tracking sites is the utilization of pre-existing infrastructure, even though their distribution across the Earth’s surface is suboptimal and cannot ensure uninterrupted visibility of the constellation at a minimum elevation of 15°. The maximum visibility gap spans approximately ~6 h. This gap is clearly visible in Figure 2, which shows the constellation visibility over a period of 30 days in order to fully cover the relative Earth-Moon geometry evolution. We note that an additional fourth station located, for example, at Manua Kea (Hawaii) would completely fill the large longitudinal gap (approximately 151°) between Malargüe and New Norcia, enabling continuous tracking.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Ground station visibility of the LRNS constellation from each of the three ESTRACK sites for a 30-day period (June 2026) with a minimum elevation of 15°, taking into account the Moon occultations.

For each day, the four light-colored bars represent the visibility for the different satellites of the constellation from each ground station (according to color). The visibility windows are expressed as hours after midnight (UTC) each day.

Clearly, changing the station locations while preserving their number and relative longitude would not substantially modify the OD results, as the visibility pattern (and thus the amount of data acquired and the data acquisition times) would remain largely unchanged. A spacecraft may occasionally enter in the field of view of two ground stations simultaneously. However, each spacecraft can establish a single two-way radio link at a given time; thus, the observables are established from the ground station with the highest elevation.

3.1 Dynamical Model

The dynamical model of the spacecraft, which is used in the first step of the numerical simulation to generate a reference trajectory, may be perturbed in the estimation phase to mimic the mismodeling of poorly known quantities, as actually occurring in a realistic OD process. The dynamical model of the spacecraft in the simulation phase includes the following:

• Gravitational monopole accelerations due to the Earth, Moon, and solar system planets

• Spherical harmonic coefficients of the Moon, up to degree and order 120 (Lemoine et al., 2014)

• Spherical harmonic coefficients of the Earth, up to degree 12

• Non-gravitational acceleration due to solar radiation pressure (SRP) acting on the satellites

The relative magnitude of the considered accelerations has large variability; for example, at the orbit pericenter, the Moon’s gravitational monopole has a contribution of 3.9×10⁻¹ m/s², the degree-20 lunar spherical harmonic acceleration is approximately 2.3×10⁻¹¹ m/s², and the SRP acceleration is approximately 4×10⁻⁸ m/s². Other sources of non-gravitational accelerations, such as the lunar and terrestrial albedo and their infrared emission, as well as the spacecraft anisotropic thermal emission, have not been included in the model because of their small magnitude when compared with the SRP modeling accuracy derived from realistic assumptions. However, the simulations implicitly consider the impact of incorrect modeling of these accelerations on the final covariance matrix, by including random accelerations in the list of estimated model parameters (as discussed later).

A spacecraft shape is required to model the non-gravitational accelerations. In the absence of a consolidated satellite design, we assume typical values for the SmallSat class, both in terms of mass and volume. We assume that each satellite may be represented as two plates and a sphere, similar to the box-wing model used for the Galileo satellites (Bury et al., 2019). The first plate models solar arrays (total area of 1.5 m²) pointed toward the Sun, the second plate represents the antenna dish (area of 0.1 m²) constantly pointed toward the Earth, and the sphere (radius of ~0.34 m) models the satellite bus. In principle, a solar panel with an area of 1.5 m² could generate an onboard power of 515 W (25% efficiency), which is more than sufficient to meet power requirements during normal operations. The estimated power needs resulting from the radio frequency system (including the transponder and the clock) correspond to approximately 90 W, providing a large margin for the considered uncertainty in SRP modeling.

The thermo-optical properties of these elements have been chosen to match those of the Galileo satellite bow-wing formulation (Li et al., 2019). The mass of each spacecraft has been set to 230 kg, thus producing an area-over-mass ratio of ~0.0085 m²/kg, corresponding to an average SRP acceleration of approximately 4×10⁻⁸ m/s². This area-to-mass ratio is roughly 2.5 times smaller than that of the Galileo satellites. It is important to note that the spacecraft shape used in this paper is different from the shape described in the works of Iess et al. (2023), Iess et al. (2024), and Di Benedetto et al. (2022); in those previous works, the authors defined a larger spacecraft area-over-mass ratio, corresponding to a greater SRP acceleration, which affects the OD results.

In any OD process, inaccuracies in the dynamical model due to missing or poorly modeled accelerations unavoidably introduce biases in the estimated model parameters (Bury et al., 2020). For example, it is difficult to accurately model the SRP acceleration action on the spacecraft (generally, the errors are larger than 2% of the central value; see Park et al. (2012)). Therefore, we have assumed that the SRP acceleration is mismodeled by approximately 5%. To simulate a realistic OD process, the first-guess trajectory adopted in the estimation phase should be different from the trajectory used to generate the synthetic data, reflecting inaccuracies in the dynamical model. To this aim, we generated the first-guess trajectories by considering a 5% error in the SRP acceleration (corresponding to an error of approximately 2.0×10⁻⁹ m/s²).

In the estimation step, we included estimated empirical accelerations in the dynamical model to compensate for the erroneous representation of the SRP. Although different formulations can be used, this analysis adopts an empirical piecewise-constant acceleration model. We estimate the three components of this acceleration in the International Celestial Reference Frame (ICRF) with an update time of 4 h, for a total of 18 estimated parameters for each satellite orbit (24-h orbital period). The orbital filter processes the observables and estimates the model parameters reported in Table 2.

It is important to note that the empirical acceleration parameters have a nominal value equal to zero; thus, their estimated value should compensate for the introduced mismodeling. These accelerations are assumed to be uncorrelated between each 4-h batch interval (white-noise statistics). The a priori uncertainty of the accelerations, based on the expected mismodeling of spacecraft dynamics, is 2.0×10⁻⁹ m/s². The biases in the observable model for range and SBI have been included to account for errors in the calibration of the station and transponder delay and to mitigate any uncalibrated media effects. In our architecture, given the small dimension of the dish, the geometric delay of the antenna is smaller and more stable than that of deep space antennas. Regarding the station and transponder delay calibration system, we believe that the same approach used for BepiColombo can be implemented for LRNS as well in a straightforward manner. Thus, we estimate a single (and constant) range bias for each station (including transponder delay, station delay, and uncalibrated media contributions) during each observation arc, in agreement with the results of BepiColombo’s Mercury Orbiter Radio-science Experiment (MORE) (Cappuccio et al., 2020; Iess et al., 2021). While it is theoretically possible to resolve the phase ambiguity in the SBI observables, achieving a post-processing absolute phase measurement to within approximately 10% of a wavelength, we opted for a more conservative approach, which involves including a constant phase bias for each tracking pass and station in the list of parameters to be estimated.

3.2 Observation Model

The observation model is of paramount importance for both generating the synthetic measurements and calculating the computed observables in the estimation process. In the baseline configuration, the simulated observables are Doppler, SS ranging, and SBI measurements, assuming MSPA tracking.

The Doppler and ranging measurement models are well-known and have been presented in detail by Moyer (2005). The SBI technique offers valuable insights into the relative motion of angularly close spacecraft. The fundamental concept of SBI involves simultaneously observing two spacecraft from a single ground antenna and comparing the phases of the received carrier signals. This process typically utilizes a two-way configuration, enabling highly precise differential phase measurements with millimeter-level accuracy (Gregnanin et al., 2012). These observables provide LOS information by measuring the difference in round-trip light-time (i.e., range) between each spacecraft and the common ground station, thereby indirectly revealing the relative motion between the two spacecraft. The key advantage of this method is the extreme accuracy of phase measurements arising from the cancellation of common mode errors, which significantly reduces the media noise, with only small residual effects caused by the non-zero angular difference (in azimuth and elevation) of the signal paths.

SBI is achieved by comparing (i.e., differencing) the phases of two-way signals, which requires defining three distinct epochs: the time t₁ when the ground station transmits the signal, the time t₂ when the signal is received and retransmitted by spacecraft A and B (with a turnaround frequency ratio M), and finally, the time t₃ when the signals from both spacecraft are received at the same ground station. To develop the SBI observable model, we start by writing the two-way phase for a spacecraft (either A or B), beginning from the reception time t₃ and tracing back the signal path, as illustrated in Figure 3.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Phase definition in a two-way radio link for a single spacecraft (left) and SBI observable procedure for two generic spacecraft A and B (right)

At any given time t₃, the measured phase difference between the two received signals, which arises from differences in their respective t₁, contains information about the spacecraft state vectors. We compare the phases received from both spacecraft at a common reception time, noting that the epochs t₁(t₃) and t₂(t₃) differ for spacecraft A and B and are functions of t₃. These epochs are computed within the OD software by solving the light-time problem for a given reception epoch for each spacecraft. Similarly, the values for uplink and downlink travel times are expressed as functions of the considered reception time, denoted as τ_up(t₃) and τ_dn(t₃), respectively. For simplicity, we will omit this explicit dependence in the following notation, mentioning it only when necessary for clarity.

Referring to the signal path in Figure 3, the instantaneous phase of the received signal at the ground station for a reception time t₃ can generally be expressed as follows:

${\Phi }_{RX}^{G/S}\left(t_3\right)=M{\Phi }_{RX}^{S/C}\left(t_3-{\tau }_{dn}\right)=M{\Phi }_{TX}^{G/S}\left(t_3-{\tau }_{up}-{\tau }_{dn}\right)$3

This expression indicates that the received phase corresponds to the phase transmitted one round-trip light-time earlier, multiplied by a constant, known as the turnaround ratio M. If the ground station transmits at a nominal frequency ω₀ with an initial phase offset ϕ₀, we can rewrite Equation (3) as follows:

${\Phi }_{RX}^{G/S}\left(t_3\right)=M{\omega }_0·\left[t_3-{\tau }_{up}-{\tau }_{dn}\right]+M{\phi }_0=M·\left[{\omega }_0·t_3+{\phi }_0+{\Psi }_{RTLT}\right]$4

The phase term Ψ_RTLT = –ω₀ (τ_up+τ_dn) represents a quantity that changes very slowly over time, owing to the relative motion between the spacecraft and the ground station. This term corresponds to the two-way range ρ_RTLT of the spacecraft at the t₃ time tag: Ψ_RTLT(t₃) = –ω₀/c(ρ_up+ρ_dn) = –kρ_RTLT(t₃), where k is the wave vector and c is the speed of light. Equation (4) models the total unwrapped phase value from the beginning of the tracking pass. The measured phase is obtained at the ground station via a phase-locked loop that reconstructs the unwrapped phase (averaged over a count time), with an ambiguity corresponding to the integer number of cycles since the initial spacecraft position relative to the ground station. This ambiguity must be accounted for in the model; otherwise, a bias would appear in the OD residuals when the spacecraft trajectories are corrected. Therefore, Equation (4) must be modified as follows:

${\Phi }_{RX}^{G/S}\left(t_3\right)=M·\left[{\omega }_0·t_3+{\phi }_0+{\Psi }_{RTLT}\right]-2\pi N$5

where N is an integer, different for spacecraft A and B.

Figure 3 illustrates the concept of the SBI observable y_SBI(t₃), which is defined as the difference between Equation (5) calculated for spacecraft A and B. Assuming that both spacecraft use the same transponding ratio (enabled by CDMA) and that the ground station transmits at the same nominal frequency (enabled by CDM-M), we can express this observable as follows:

$y_{\text{SBI}}\left(t_3\right)=M\left[{\Psi }_{\text{RTLT}}^B-{\Psi }_{\text{RTLT}}^A\right]-2\pi \left(N^B-N^A\right)=\frac{M{\omega }_0}{c}\left[{\rho }_{\text{RTLT}}^A-{\rho }_{\text{RTLT}}^B\right]-2\pi K$6

Here, K represents the unknown difference in the integer number of cycles between spacecraft A and B at the start of the tracking pass. Consequently, the SBI observable corresponds to the difference in two-way range between spacecraft A and B, scaled by a constant factor. However, this observable includes an initial ambiguity that must be resolved within the formulation of the OD problem. Essentially, a phase bias, which remains constant throughout each tracking pass, must be added to the list of parameters to be determined.

Any data below a minimum elevation of 15° at each ground site are not considered. We selected a four-day observation arc, considering the trade-off between improving the positioning accuracy by using more observables and limiting the computational cost and time to process these observables. The measurement noise for each data type has been simulated based on the architecture configuration presented in Section 2 and by Di Benedetto et al. (2022), Iess et al. (2023), and Iess et al. (2024). The error budget for the radiometric observables has been computed based on recent data collected by BepiColombo (Cappuccio et al., 2020; di Stefano et al., 2023) and on the results of Iess et al. (2014).

A key factor for the observation model is the media calibration system adopted for ground observations. This system can provide either a model or measurements of the ionospheric and tropospheric (both dry and wet) path delays, with the latter approach typically outperforming empirical models, at the price of increased system complexity and larger costs. In the baseline configuration, the simulated media calibration errors include the following:

• The tropospheric path delay, considering 95% calibration of the wet component obtained with a dedicated advanced water vapor radiometer (AWVR). It is the use of the AWVR that allows an accuracy up to 0.5 cm (Linfield et al., 1996; Lasagni Manghi et al., 2023).

• The ionospheric path delay, considering 90% calibration achieved with GNSS dual-frequency data (Liu et al., 2021). Notice that the ionospheric effect is reduced when a higher frequency is adopted (the dispersive noise contribution decreases as ~l/f²).

The media calibration errors have been added to the simulated range measurements as systematic effects. In the synthetic data simulation, for each tracking pass, we approximated the tropospheric zenith path delay as a parabola whose coefficients are randomly extracted to obtain a zenith path delay profile ranging between 3.0 and 5.1 cm, in line with the work of Iess et al. (2014). Then, we mapped the zenith path delay according to the station elevation at each observation epoch and added the uncalibrated percentage of the path delay to the synthetic measurements. In contrast, the ionospheric path delay was computed via the NeQuick ionospheric model (European GNSS, 2016), and similarly, the uncalibrated portion was added to the simulated range data.

The measurement noise can be separated into three parts: media propagation effects, space segment contribution, and ground segment contribution.

The error budget for Doppler data is expressed in terms of the Allan deviation (ADEV) at a given integration time τ or σ_y(τ), a universally adopted figure of merit for the stability of a frequency standard (Riley, 2008). For an integration time of 60 s, some of the relevant contributors to the overall ADEV are the onboard transponder (6.0×10⁻¹⁵), the steerable antenna assembly (8.6×10⁻¹⁵), the spacecraft structure (4.0×10⁻¹⁵), the ground clock (4.1×10⁻¹⁵ for an H-maser), and ground electronics (1.7×10⁻¹⁵) at the tracking stations. However, the overall Doppler error budget is dominated by the media propagation, especially the variations in the wet tropospheric delay (1.2×10⁻¹⁴ after calibration with water vapor radiometers).

For the error budget of ranging measurements, the uncalibrated media effect (tropospheric and ionospheric bias) is approximately 2 cm, smaller than the ground and space segment contributions (which amount to 5.7 cm) by a factor of 3. These two last contributions are due to inaccuracies in the spacecraft and ground station group delay calibration, errors in the station location (Budnik et al., 2004), the Earth solid tide effect, and inaccuracy in the Earth orientation parameters (International Earth Rotation and Reference Systems Service, 2018). For the random portion of the error budget, we considered the ranging jitter (32.7 cm), which is the largest contributing error source in ranging measurements.

The SBI end-to-end error depends strongly on the residual media path delay, which is proportional to the angular separation between the spacecraft pair (for small angles) and to the tracking elevation. For this, to compute the SBI noise, we selected the average tracking elevation and angular separation over the analyzed time interval (30 days), which are equal to 30° and 1.25°, respectively.

At the K-band and in the baseline case, the noise levels for Doppler (see the work by Iess et al. (2014)), SS range (Consultative Committee for Space Data Systems [CCSDS], 2014), and SBI measurements used in the numerical simulation are summarized in Table 3.

3.3 Time Transfer

In the proposed architecture, the ground-to-space time transfer can be performed via two methods. The first approach is the TT-Async-Mode, which relies on two-way asynchronous/non-coherent links between an Earth tracking station and each satellite of the lunar constellation. For this method, the onboard transponder must be operating in non-coherent mode. This method is conceptually similar to the TWSTFT used on Earth in the generation of the UTC timescale (Arias et al., 2011). For Earth TWSTFT, a geostationary satellite is indeed used as a relay for the signals to close the links (in both directions) between two distant ground stations without reciprocal visibility. However, when the two terminals are mutually visible, it is possible to directly establish a two-way asynchronous link. The main advantage of this TT method is that it can be considered as virtually immune from orbital errors, at the cost of interrupting radiometric data acquisition (requiring a switch in operation mode of the onboard transponder from coherent to non-coherent) and a different signal structure with respect to the coherent mode structure. The overall uncertainty of the time transfer with the TWSTFT method can be preliminarily assumed to be below the 1-ns level (~ 0.3 ns can likely be achieved if periodic calibrations occur at both segments), in line with the uncertainties reported in Circular T by the Bureau International des Poids et Mesures (BIPM, 2023).

The second method is the TT-Sync-Mode. This method relies on two-way coherent measurements and aims at exploiting the performances of Pseudo-noise (CCSDS, 2022)/SS (CCSDS, 2011) ranging systems for time-transfer purposes. In any two-way coherent radiometric measurement, solving for the light-time solution requires the consideration of three distinct epochs:

• ° t₁: epoch of signal transmission from ground

• ° t₂: epoch of signal reception onboard satellite

• ° t₃: epoch of signal reception on ground

The onboard clock provides direct access to measurements of t₂, while the ground clock measures t₃. Then, both ${\hat{t}}_1$ and ${\hat{t}}_2$ can be computed by solving backward for the light-time solution (Moyer, 2005) by means of the OD process, where the two-way ranging observable represents a measure of the round-trip light-time (RTLT = t₃ – t₁). Therefore, the ground-to-space clock desynchronization can be inferred as follows:

$desync\left(t\right)={\hat{t}}_2\left(t\right)-t_2\left(t\right)=\left(t_3\left(t\right)-\frac{{\rho }_{23}\left(t\right)}{c}\right)-t_2\left(t\right)$7

where ${\hat{t}}_2\left(t\right)$ is the onboard time derived from the OD solution through ρ₂₃ / c, the one-way light-time solution (in the downlink leg), whereas t₂(t) is the reading of the onboard clock. The onboard time-stamping operations are triggered by a “code epoch” signal, activated by a specific chip of the SS signal. Then, the recorded epoch is sent to ground in the telemetry stream. On the ground, a similar time-tagging operation is triggered by the received code epoch signal, and the epoch t₃ provides ${\hat{t}}_2$ through the OD process, while telemetry data for t₂ are downloaded by the satellite. Finally, these data are compared to obtain the desynchronization.

The main advantage of this novel technique (Iess et al., 2023; Iess et al., 2024) is that it can be performed in parallel with nominal ground tracking operations. The desynchronization accuracy depends on both OD performance along the LOS, i.e., accuracy in the one-way light-time computation, and the precision of the time-stamping operations. The lunar orbiter laser altimeter of the LRO has achieved a time-stamping precision of approximately 0.5 ns (Bauer et al., 2017), whereas numerical simulations show an OD accuracy better than 3 cm along the LOS, corresponding to 0.1 ns. Therefore, a conservative estimate of the desynchronization accuracy attainable with this method is at the 1-ns level.

The clock error contribution to the overall SISE is due to the instability of the onboard clocks and the inherent inaccuracies of the time synchronization process. The former contribution is based on the spectral characteristics of the candidate clocks: ultra stable oscillator (USO), rubidium atomic frequency standard (RAFS), and miniaturized RAFS (miniRAFS). A stochastic noise realization of each clock is simulated based on the discretized ADEV via the method of Timmer et al. (1995). The desynchronization contribution is simulated and evaluated for the two-way synchronous time transfer, as this method is expected to provide more pessimistic results than the TWSTFT, which is not affected by the OD contribution. During the satellite tracking window, desynchronization measurements can be collected from two-way coherent ranging measurements. These observables can be used to extrapolate the clock behavior, in accordance with the following relation:

$\Delta {\tau }_s=\sum _{i=0}^D\Delta {\tau }_s^{\left(i\right)}{\left(\tau -{\tau }_0\right)}^i$8

For D = 2, this expression represents the combination of a clock offset Δτ⁽⁰⁾, a clock drift Δτ⁽¹⁾ (frequency bias), and a clock quadratic term Δτ⁽²⁾ (frequency drift). These clock correction parameters are included in the navigation message transmitted to the end user. For present purposes, we are interested in calculating the SISE clock contribution, given by the fit obtained by Δτ_s applied during the propagation interval. The onboard clock model fit is performed over a data set acquired during a single tracking pass (up to a maximum of 6 h), so that the overall process can be repeated three times per day.

4.1 Supplementary Analysis

We performed additional numerical simulations to evaluate the LRNS performance if some of the architecture hypotheses are modified as follows:

1. The SBI technique is not implemented, and only Doppler and SS ranging data are collected during the observation arcs and processed in the OD filter.

2. Water vapor radiometer calibrations of the wet tropospheric path delay are replaced by GNSS-based calibrations (Tondaś et al., 2020; Dousa & Václavovic, 2015). We assumed 80% calibration of the wet path delay (systematic effect).

3. A fourth station at Mauna Kea (Hawaii) is included, to avoid visibility gaps and enable continuous visibility from the ground.

4. An improved dynamical model is adopted, and the mismodeling of the SRP acceleration decreases from 5% to 2.5%. The empirical accelerations have the same batch interval but an a priori uncertainty of 1.0×10⁻9 m/s².

In the second scenario, the SISE increases because of the degraded media calibration level. In contrast, the third case implies a better performance for the system, given that radiometric observables are collected in a larger quantity and with more continuity. Note that for the first case, the tracking system remains the same, and the only difference lies in the data analysis procedure. Finally, the last scenario implies an improvement in ephemeris aging owing to the decreased dynamical mismodeling. These different scenarios are compared in Table 5, with the miniRAFS considered as the onboard clock for the constellation satellites.

As expected, the error in the trajectory recovery increases if the SBI observables or the AWVR are not used; in contrast, the error is reduced if an additional station is considered or a better dynamical model is used, as shown by the mean RMS values in Table 5. With regard to aging, the first three assumptions do not show significant variations, primarily because the ephemeris aging is largely driven by the inaccuracies in the dynamical model and the highly eccentric orbits. Indeed, the last scenario shows the largest variation in SISE performance.

We note that the use of SBI does not significantly change the OD performance, primarily because we are not able to solve the phase ambiguity (see Section 3.2) with the adopted ranging system accuracy. As a result, we must introduce a bias parameter for each spacecraft pair and ground station for each tracking pass, as shown in Table 2, thus degrading the performance of the SBI measurements. However, it is important to note that the SBI observables can be generated without additional hardware at the ground station, directly from an open-loop receiver using MSPA tracking, and their usage still improves the OD solution.

The improvement in the OD solution brought by the SBI is due to its capacity to cancel out a significant portion of the common noise in the observables (see Section 3.2). This feature leads to more accurate measurements, enhancing the OD accuracy, particularly along the LOS direction, and consequently reducing the OD error contribution to the time synchronization.

Note that the 99% values in Table 5 have large uncertainties, owing to the limited number of simulated arcs: the obtained ephemeris aging only excludes the worst ~6 arcs in the Monte-Carlo-like analysis (1% of 577 simulated arcs). Given this reduced number of cases for the statistic, there is a large uncertainty in the value, as shown in Figure 8.

In this work, we simulated the SISE as a function of the AOD. This definition does not consider the epoch at which the last observables are collected, thus the ephemerides update; rather, the AOD is related only to the navigation message generation time. The SISE evolution as a function of time, as shown in Figure 8, should drive the choice of frequency update of the navigation message to meet the performance requirements of the LRNS.