Adaptation of One-Way Radiometric Range and Range-Rate Errors to the Lunar Environment

2.2.1 Ephemeris and Timing Accuracy

From Equations (1) and (2), it is clear that the user must obtain the satellite positions and velocities at the time of signal transmission, calculated from ephemerides, to estimate their own state. This information is traditionally broadcast from GPS satellites to the receiver in the form of time-dependent equation parameters encoded in the navigation message. The Control Segment tracks each GPS satellite and relays ephemeris updates at least daily (Anthony & Kerns, 2022). In addition, each satellite’s onboard clocks drift independently with respect to GPS Time. The Control Segment monitors this drift and transmits clock correction parameters along with the ephemerides. However, both the clock and ephemeris corrections include residual errors that propagate to the user. These errors appear as biases over short time intervals when the ephemeris and clock models differ from each satellite’s true state. The equations used to compute the position and velocity of the satellite antenna phase center are given in Table 30-II of the GPS Interface Specifications (Anthony & Kerns, 2022); the computation process involves calculating orbital elements at the time of measurement, calculating the GPS satellite’s position and velocity in the perifocal frame, and then transforming these into an Earth-centered Earth-fixed (ECEF) reference frame.

Satellite clocks are referenced to GPS Time, which is a fixed whole number of seconds ahead of Coordinated Universal Time (UTC). The satellite clock bias is approximated by the Control Segment using a quadratic polynomial equation, given by Equation (3):

$\begin{array}{ccc}{\delta t}_i& =& t_i-t_{\mathrm{GPST}}\\ & =& a_{f0}+a_{f1}(t_{\mathrm{GPST}}-t_{\mathrm{oc}})+a_{f2}{(t_{\mathrm{GPST}}-t_{\mathrm{oc}})}^2+{\Delta t}_r,{\Delta t}_r=\frac{-2\sqrt{{\mu }_{\oplus }}}{c^2}{e}^{\sqrt{a}}\end{array}$3

where a_f0, a_f1, and a_f2 are the coefficients for clock bias, drift, and drift rate correction, respectively; t_oc is the reference time of the clock model; and Δt_r is a relativistic correction. These polynomial coefficients are broadcast to the user as part of the satellite’s navigation message.

The U.S. government provides estimates of expected errors in its GPS Performance Specifications (PS). For example, Table 1 shows the expected performance for a civilian user equipped with a single-frequency receiver (U.S. Department of Defense, 2020). As shown in this table, the root-sum-squared error contribution from the ephemeris and clock uncertainties (including stability, estimation, prediction, and curve fitting) reaches a standard deviation of $3\sigma {\rho }_{,E/C}=14.0$ meters under the maximum Age of Data (AOD) in normal operating conditions (i.e., 24 hours). AOD refers to the time since the Control Segment last updated the satellite ephemerides and clock correction parameters.

2.2.2 Transmission Delays (Ionospheric, Tropospheric, Group)

Terrestrial GNSS systems experience ionospheric and tropospheric delays due to refraction of RF signals in the atmosphere. Although these effects are not relevant in cislunar space, group or instrument delays caused by antennas, hardware, and software in both the user and satellite electronics are still a concern in the cislunar environment. In GPS, group delays on the user side are absorbed into the receiver clock bias estimates, while group delays on the satellite side result in signal transmission delays. These delays vary with broadcast frequency and PRN code. Because GPS broadcasts across multiple frequencies, users with a dual-frequency receiver can combine measurements to effectively eliminate this error. However, a user with a single-frequency receiver on L1 C/A can expect an uncertainty contribution of approximately 3σ_GRP = 4.2 meters after applying broadcast corrections. This uncertainty represents the root-sum-square of the group delay’s stability and the uncertainty in the correction terms (U.S. Department of Defense, 2020).

2.2.3 Multipath

Multipath refers to the phenomenon in which a signal is reflected off some part of the environment near the user’s antenna and thus arrives at the receiver both later and weaker than the direct-path signal. Due to the autocorrelation properties of PRN codes, multipath effects diminish significantly for reflected signals delayed by more than 1.5 chip wavelengths; higher chipping rates therefore provide greater immunity to multipath effects. Additional mitigation strategies include using a highly directional antenna and/or elevation masks. Although multipath effects are highly dependent on the user’s environment, the GPS PS reports potential error values ranging from 3σ_Multi = 3.6 meters for traditional single-frequency receivers down to 0.3 meters for modern single-frequency receivers, as shown in Table 1 (U.S. Department of Defense, 2020).

2.2.4 Receiver Noise

In general, there are two ways to measure a GPS signal at a given frequency: tracking the PRN code embedded in the signal (code tracking), or tracking the phase or frequency of the carrier signal on which the data is modulated (carrier tracking). For a given GPS signal and frequency, carrier tracking typically yields errors that are two orders of magnitude smaller than errors that arise from code tracking. However, carrier tracking suffers from an additional accuracy limitation known as an integer ambiguity. To remove this ambiguity, most modern receivers use a hybrid approach known as carrier-aided code tracking. In these architectures, the output of the PLL or FLL (used for carrier tracking) is scaled and applied to the DLL output as a correction factor. This approach yields an exceptionally accurate pseudorange measurement compared to unaided code tracking loops (Kaplan & Hegarty, 2017; Misra & Enge, 2006). In addition, the Doppler shift or differenced carrier-phase measurements can be used to estimate the pseudorange-rate. Kaplan and Hegarty provide a comprehensive discussion of GNSS receiver design and measurement errors (2017); the errors vary widely over different receiver designs and are therefore closely tied to the mission design process. This paper includes both a description of different receiver-related error sources and a sample receiver design for comparative analysis.

The dominant sources of error in each receiver tracking loop are typically thermal noise and dynamic stress. For a given code-tracking DLL, the associated measurement uncertainty can be approximated as:

$\begin{array}{cc}3{\sigma }_{DLL}\left(t\right)=3{\sigma }_{tDLL}\left(t\right)+R_e\left(t\right)\left(m\right),& R_e\left(t\right)=\frac{1}{{\omega }_0^n}\left(\frac{d^n}{d{t}^n}\delta r\left(t\right)\right)\end{array}$4

where σ_tDLL is the thermal noise, R_e is the dynamic stress, δr is the line-of-sight range, n is the code loop order, and ω₀ is the loop natural frequency. The maximum loop error is proportional to its tracking threshold; as such, receiver design often involves a tradeoff between accuracy and dynamic stress tolerance. Architectures that use carrier-aided code tracking mitigate the dynamic stress of code tracking, allowing R_e to be ignored. Equations for computing σ_tDLL and guidance on selecting loop parameters can be found in Chapter 8 of Kaplan and Hegarty (2017).

For a carrier-tracking FLL, the errors are dominated by the time derivatives of the same underlying sources as the code-tracking DLL, as follows:

$\begin{array}{cc}3{\sigma }_{\mathrm{FLL}}\left(t\right)=3{\sigma }_{\mathrm{tFLL}}\left(t\right)+{\dot{R}}_e\left(t\right)\left(m/s\right),& {\dot{R}}_e\left(t\right)=\frac{1}{{\omega }_0^n}\left(\frac{d^{n+1}}{d{t}^{n+1}}\delta r\left(t\right)\right)\end{array}$5

Here, σ_tFLL is the frequency thermal noise, ${\dot{R}}_e$ is the dynamic stress, and n and ω₀ are the loop order and natural frequency, respectively. Because the FLL tracks the carrier frequency, its errors are typically expressed in either in m/s or Hz. FLLs have a wider pull-in range than PLLs but are also more noisy, making them more suitable in environments with higher dynamic stress or more challenging RF.

Carrier-tracking PLLs are affected by similar error sources as DLLs and FLLs (thermal noise and dynamic stress) but with additional contributions from oscillator imperfections. PLL errors can be formulated as:

$3{\sigma }_{\mathrm{PLL}}\left(t\right)=3\sqrt{{\sigma }_{\mathrm{tPLL}}{\left(t\right)}^2+{\sigma }_v{\left(t\right)}^2+{\sigma }_A{\left(t\right)}^2}+R_e\left(t\right)\left(m\right)$6

where σ_tPLL describes the thermal noise, σ_v is the vibration-induced oscillator phase noise, and σ_A is the Allan deviation oscillator phase noise. In this paper, we consider a receiver using FLL-assisted PLL carrier tracking loops and DLL code tracking loops. Under nominal operations with continuous phase lock, the pseudorange measurement uncertainty for this receiver is ${\sigma }_{\rho ,\hspace{0.17em}rec}\left(t\right)={\sigma }_{tDLL}\left(t\right)$. For velocity estimation using differenced carrier-phase measurements, which are assumed to be statistically independent, the associated uncertainty is:

${\sigma }_{\dot{\rho },\mathrm{rec}}\left(t\right)=\sqrt{2}{\sigma }_{\mathrm{PLL}}\left(t\right)/T_m$7

where T_m is the time between measurements. The choice for T_m represents a design trade-off, as larger T_m reduce ${\sigma }_{\dot{\rho },\hspace{0.17em}rec}$ but increase quantization error in the velocity measurement. The appropriate design choice ultimately depends on user dynamics. According to the GPS PS, measurement uncertainties range from 3σ_{ρ, rec} = 4.5 to 0.6 meters depending on receiver quality (U.S. Department of Defense, 2020).

3.3.1 State Propagation in Lunar Orbits

Orbital dynamics in cislunar space are significantly more perturbed than in Earth orbit. The nonspherical effects of the Moon’s gravity on satellite trajectories are proportionally larger, and third-body perturbations from the Earth and Sun play a major role in shaping the types of orbits currently being considered for lunar navigation satellites. For the purpose of this paper, we define the acceleration of a satellite around the Moon ${\ddot{r}}_s$ using Equation (9):

${\ddot{r}}_s={}_{\mathrm{ICRF}}{}^{\mathrm{PA}}T\nabla \frac{{\mu }_{⦆}}{r_s}\sum _{n=0}^N\sum _{m=0}^n\frac{R_{⦆}^n}{r_s^n}{P}_{nm}\left(\sin \varphi \right)\left(C_{nm}\cos m\lambda +S_{nm}\sin m\lambda \right)+\sum _{i=1}^k{\mu }_i\left(\frac{r_{si}}{r_{si}^3}-\frac{r_i}{r_i^3}\right)-v{P}_{\odot }{C}_R\frac{A}{m}\frac{r_{s\odot }}{r_{s\odot }^3}{\left(1\mathrm{AU}\right)}^2$9

The first term on the right-hand side represents the contribution of the Moon’s gravity, including nonspherical effects (Montenbruck & Gill, 2000). In this term, μ_⦆ is the Moon’s gravitational parameter, R_⦆ is its reference radius, P_nm is the associated Legendre polynomial of degree n and order m,ϕ and λ are the satellite’s latitude and longitude, and C_nm and S_nm are the Moon’s spherical harmonic coefficients (Konopliv et al., 1998). This gravity expression is normally given in the Moon Principal Axis (PA) frame for lunar gravity fields (Williams et al., 2013), so for our purposes, it must be rotated to the Moon-centered inertial frame, known as the International Celestial Reference Frame (ICRF), using the rotation matrix ${}_{\mathrm{ICRF}}{}^{\mathrm{PA}}T$. The second term on the right-hand side accounts for third-body perturbations from up to k planets. The final term is the cannonball model for the effects of solar radiation pressure (SRP). Here, v is the shadow function, P_⊙ is the SRP at 1 AU, C_R is the coefficient of reflectivity, and A/m is the spacecraft area-to-mass ratio. These perturbations are discussed in Chapter 3 of Montenbruck and Gill (2000).

For computational feasibility, we truncate all spherical harmonics calculations at a maximum degree and order N, noting that the algorithm scales with complexity following $O\left(N^2\right)$). For a case study, we take the initial position and velocity of the Lunar Reconnaissance Orbiter (LRO) on October 15, 2009 at midnight UTC, and propagate its trajectory over 7 days using models with varying degrees of fidelity. The propagated results are then differenced from the true mission data (Burtnick et al., 2010) to compare the accuracy of different dynamical models. The results from this case study are plotted in Figure 2(a). For low lunar orbits, trajectory propagation using spherical harmonics up to degree and order 165 performs nearly three orders of magnitude better than models using Keplerian dynamics alone or models that only include third-body perturbations from the Earth, Sun, and Jupiter. The N = 165 model is also two orders of magnitude more accurate than the N = 10 model, though at a higher computational cost. Including the effects of SRP further improved prediction accuracy, though this effect was relatively small. However, even with the highest-fidelity model that we considered, the projected trajectory drifts by 11–12 m per orbital period, resulting in a cumulative position error of approximately 1 km after one Earth week. This error is most likely due to uncertainties in the provided LRO trajectory, as NASA only guarantees positional accuracy under 500 meters (Burtnick et al., 2010). Other unmodeled perturbative forces, like lunar tides or more complex SRP interactions, may also have contributed to this error.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Accuracy of various dynamical models for low lunar and elliptical lunar frozen orbits. (a) Models compared with truth data from the LRO mission. (b) Models compared with simulated data from GMAT.

Figure 2(b) shows a corresponding plot for elliptical lunar frozen orbits (ELFOs). For this analysis, the reference trajectory was generated using a full-order force model in NASA’s General Mission Analysis Tool (GMAT) (Hughes et al., 2017), including additional planets and relativistic effects beyond those used in our simulation. A full summary of the propagation conditions for the ELFO analysis is provided in Table 2.

In the ELFO scenario, third-body perturbations play a much larger role, improving the propagated trajectory by two orders of magnitude over purely Keplerian dynamics. An additional two orders of improvement are gained by including spherical harmonic terms; however, improving model fidelity from N = 10 to N = 165 showed limited improvement. A final two orders of magnitude in accuracy are gained by including SRP, which reduced the integration tolerance to the meter level. The remaining error can be attributed to unmodeled third-body effects from other planets, reference frame misalignment, or integration tolerances. We note that the lunar gravitational parameter used in GMAT differs from that used in SPICE and elsewhere, so we modified this value to ensure consistency across the propagation analyses presented here.

Letting Equation (9) define the acceleration model F(r_s, t), we can represent the satellite dynamics as:

$\begin{array}{ccccccc}\dot{x}& =& f\left(x,t\right)+{\epsilon }_f,& & x& =& {\left[r_s\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\dot{r}}_s\right]}^T\\ \left[\begin{array}{c}{\dot{r}}_s\\ {\ddot{r}}_s\end{array}\right]& =& \left[\begin{array}{c}{\dot{r}}_s\\ F\left(r_s,t\right)\end{array}\right]+\left[\begin{array}{c}0\\ {\epsilon }_{\text{accel}}\end{array}\right]& & & & \end{array}$10

where ε_accel represents the unmodeled errors inherent in Equation (9). The model derived in (9) is therefore an estimate $\overline{x}\left(t\right)$ of the true spacecraft state, x(t).

Given an initial state estimate $\overline{x}\left(t_0\right)={\overline{x}}_0$ from a navigation solution provided by the spacecraft computer, and a corresponding covariance $E\left(\left({\overline{x}}_0-x_0\right){\left({\overline{x}}_0-x_0\right)}^T\right)=P_0$, the propagation of uncertainty over time can be approximated via linearization using the state transition matrix:

$P\left(t_{k+1}\right)=\Phi \left(t_{k+1},t_k\right)P\left(t_k\right){\Phi }^T\left(t_{k+1},t_k\right)$11

where $\Phi (t_{k+1},t_k)$ is the state transition matrix of x from t_k to t_k+1. Assuming perfect knowledge of the dynamics, we can ignore process noise. This iteration can then be used to find the state covariance at any time t from a starting time t₀. The state transition matrix itself is defined by:

$\begin{array}{cc}\dot{\Phi }\left(t,t_0\right)=A\left(t\right)\Phi \left(t,t_0\right),& \Phi \left(t_0,t_0\right)=I_6\end{array}$12

where A(t) is the Jacobian of the nonlinear dynamics f with respect to the state vector x. Numerical methods for approximating (12), including details on step size and computational accuracy, are outlined by Carpenter and D’Souza (2018).

Once A is obtained, Equation (11) can be used to forecast the uncertainty at some future time. A is given by:

$\begin{array}{ccc}A\left(t\right)={\left(\frac{df}{dx}\right)}_{x\left(t\right)}& ={\left(\left(\begin{array}{cc}\frac{\partial {\dot{r}}_s}{\partial r_s}& \frac{\partial {\dot{r}}_s}{\partial {\dot{r}}_s}\\ \frac{\partial F}{\partial r_s}& \frac{\partial F}{\partial {\dot{r}}_s}\end{array}\right)\right)}_{x\left(t\right)}& ={\left(\left(\begin{array}{cc}{0}_{3\times 3}& I_3\\ \frac{\partial F}{\partial r_s}& 0_{3\times 3}\end{array}\right)\right)}_{x\left(t\right)}\end{array}$13

The measurement uncertainty at any given time is then obtained by projecting the covariance P(t) along the line-of-sight direction ${\hat{u}}_i$, as described in Equation (2). To reduce complexity for the user, the line-of-sight direction may be approximated as the vector connecting the Moon’s centroid to the satellite.

The covariance matrix can be decomposed into four 3 × 3 submatrices:

$P\left(t\right)=E\left(\left(\begin{array}{c}\overline{r}\left(t\right)-r\left(t\right)\\ \overline{\dot{r}}\left(t\right)-\dot{r}\left(t\right)\end{array}\right){\left(\begin{array}{c}\overline{r}\left(t\right)-r\left(t\right)\\ \overline{\dot{r}}\left(t\right)-\dot{r}\left(t\right)\end{array}\right)}^T\right)=\left(\begin{array}{cc}{P}_r\left(t\right)& P_{\mathrm{rv}}\left(t\right)\\ P_{\mathrm{rv}}\left(t\right)& P_v\left(t\right)\end{array}\right)\cdot$

We can then write the uncertainty in the range and range-rate due to propagation effects as

$\begin{array}{cc}{\sigma }_{\rho ,\mathrm{prop}}\left(t\right)=\sqrt{{\hat{u}}_i^T{P}_r\left(t\right){\hat{u}}_i},& {\sigma }_{\dot{\rho },\mathrm{prop}}\left(t\right)=\sqrt{{\hat{u}}_i^T{P}_v\left(t\right){\hat{u}}_i}\end{array}$14

Due to the complexity of the dynamics in (9), we approximate the spherical harmonics of the Moon as $J_2=-C_{2,0}$. Under this approximation, the acceleration on the satellite due to the J₂ effect in the Moon PA frame is then described as follows (note that, for brevity, state variables not explicitly written as a function of time):

$\begin{array}{c}\begin{array}{ccc}{}_{\mathrm{PA}}\stackrel{\cdot \cdot }{r}{}_{J_2}& =& \frac{3}{2}{J}_2{\mu }_{⦆}{R}_{⦆}^2\left(\frac{5z^2}{r_s^7}{}_{\mathrm{PA}}r{}_s-\frac{1}{r_s^5}S{}_{\mathrm{PA}}r{}_s\right),\end{array}\\ \begin{array}{ccc}{}_{\mathrm{PA}}r{}_s& =& {\left(\begin{array}{ccc}x& y& z\end{array}\right)}^T,\end{array}\hspace{0.17em}S=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 3\end{array}\right).\end{array}$15

We then approximate F(r_s, t) as:

$F\left(r_s,t\right)=-\frac{{\mu }_{⦆}}{r_s^3}{r}_s+{}_{\mathrm{ICRF}}{}^{\mathrm{PA}}T{}_{\mathrm{PA}}{\ddot{r}}_{J_2}+{\mu }_{\oplus }\left(\frac{r_{s\oplus }}{r_{s\oplus }^3}-\frac{r_{\oplus }}{r_{\oplus }^3}\right)$16

The partial derivatives provided by Vallado and McClain (2007) can be extrapolated to provide a concise representation of the generalized derivative:

$\frac{\partial }{\partial r}\frac{-\mu }{{\left(r\right)}^n}r=n\mu \frac{r{r}^T}{{\left(r\right)}^{n+2}}-\mu \frac{1}{{\left(r\right)}^n}{I}_3\cdot$17

Applying Equation (17) to (16), where possible, yields:

$\frac{\partial F}{\partial r_s}=3\left({\mu }_{⦆}\frac{r_s{r}_s^T}{r_s^5}+{\mu }_{\oplus }\frac{r_{s\oplus }{r}_{s\oplus }^T}{r_{s\oplus }^5}\right)-\left(\frac{{\mu }_{⦆}}{r_s^3}+\frac{{\mu }_{\oplus }}{r_{s\oplus }^3}\right)I_3+\dots +\frac{3}{2}{J}_2{\mu }_{⦆}{R}_{⦆}^2{}_{\mathrm{ICRF}}{}^{\mathrm{PA}}T\left(\frac{5z^2}{r_s^7}S-\frac{35z^2{}_{\mathrm{PA}}r{}_{s\mathrm{PA}}{r}_s^T}{r^9}-\frac{1}{r_s^5}S+\left(\begin{array}{ccc}1& 1& 3\\ 1& 1& 3\\ 3& 3& 3\end{array}\right)\odot \frac{5\left({}_{\mathrm{PA}}r{}_s\right){\left({}_{\mathrm{PA}}r{}_s\right)}^T}{r_s^7}\right){}_{\mathrm{PA}}{}^{\mathrm{ICRF}}T$18

Through these methods, we arrive at an analytic time-varying expression for A using Equation (13). This expression can be used in Equation (11) to propagate the satellite’s state covariance through time, and it will be used later to develop a broader time-dependent expression for user measurement uncertainty. Importantly, this analytic expression is reasonably compact and computationally efficient—significantly more so than numerically differentiating the full system dynamics, which took an order of magnitude longer to evaluate during testing using both central difference and complex step methods. The dynamics are simplified by neglecting third-body perturbations beyond Earth, nonspherical gravity higher than degree and order (2, 0), and SRP. These forces are either lower in magnitude or change minimally with perturbations to the state.

To validate this model, we use various Monte-Carlo simulations to compare the modeled and computed covariances for several lunar frozen orbits. We propagated eight distinct orbits 1000 times each for eight hours, with varying initial states selected according to (19):

$\begin{array}{cc}E\left({\hat{x}}_0\right)=x_0,& E\left(\left({\hat{x}}_0-x_0\right){\left({\hat{x}}_0-x_0\right)}^T\right)=\mathrm{diag}\left(\left({\sigma }_0^2,{\sigma }_0^2,{\sigma }_0^2,{\dot{\sigma }}_0^2,{\dot{\sigma }}_0^2,{\dot{\sigma }}_0^2\right)/3\right)\end{array}$19

where σ₀ = 2.24 m and ${\dot{\sigma }}_0=0.20$ mm/s are half of the maximum allowable position and velocity SISE, respectively, under NASA’s Lunar Relay Services Requirements Document (SRD): 13.43 m 3σ for position and 1.2mm/s 3σ with T_m = 10 s for velocity (Speciale et al., 2022). We used a 100 × 100 lunar gravity field for truth modeling and accounted for third-body perturbations from the Earth, Sun, and Jupiter as well as SRP using the coefficients listed in Table 2.

Table 3 lists the coefficients of determination for the fit between (i) the simplified model given in (11) and (13) and (ii) the sample position and velocity variance from the integration of (10). Overall, R² ≈ 1 for each test case. Slightly lower performance was observed during orbit passes over perilune (0°), likely due to the higher influence of unmodeled nonspherical gravity terms at lower altitudes. According to Folta and Quinn (2006), these unmodeled terms terms begin to influence the dynamics at altitudes below 750 km and dominate at altitudes below 100 km. Figure 3 illustrates a sample simulation run. In general, the maximum uncertainty gain occurs around perilune for each orbit, while the uncertainty evolves more slowly around apolune. This behavior is advantageous for LNSS service to the lunar south pole (LSP), as satellites will be near apolune when in view of surface users and will therefore require fewer navigation message updates to remain within accuracy requirements.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Sample Monte-Carlo analysis with 100 runs. (a) ELFO passing through perilune. (b) Position and velocity errors for ELFO runs.

3.3.2 Lunar Trajectory Model Fitting

To relay satellite states to users at future times, the navigation message must include a model of the satellite ephemeris. While high-fidelity dynamics can be used to propagate the current satellite state forward in time, this result must be approximated such that the ephemeris can be relayed using minimal data. Cortinovis et al. (2024) provide a comprehensive treatment of ephemeris approximation via surrogate modeling, including the message length requirements for various models. Some of the most accurate results were obtained using a Chebyshev polynomial basis surrogate model up to order 14, with Chebyshev nodes used as interpolation points (Atkinson, 1989). For this approach, the spacecraft state is first expressed as a function of time since some epoch $x(t-t_0)=\left[r\right(t-t_0\left)v\right(t-t_0\left)\right]T$, and the interval of validity [t₀, t_f] is normalized to [–1, 1]. The spacecraft state is then approximated using the Chebyshev polynomial basis.

Here, we modify the approach of Cortinovis et al. (2024) by first computing the Keplerian state of the satellite and then providing correction terms as necessary. Kepler’s equations can be solved using either p-iteration or a universal variable formulation to find the satellite state at some future time given its state at a starting epoch. This calculation is described in Chapter 4 of Bate, Mueller, and White (1971). From there, we use a Chebyshev basis to model the difference between the true and Keplerian solutions:

$\begin{array}{cc}r\left(\tau \right)=r_{\mathrm{Kep}}\left(\tau \right)+\sum _{i=0}^{n_{\mathrm{Kep}}}{\alpha }_i{b}_i\left(\tau \right)+{\epsilon }_{\mathrm{mdl}}\left(\tau \right),& v\left(\tau \right)=v_{\mathrm{Kep}}\left(\tau \right)+\sum _{i=1}^{n_{\mathrm{Kep}}}{\alpha }_i{\dot{b}}_i\left(\tau \right)+{\dot{\epsilon }}_{\mathrm{mdl}}\left(\tau \right)\end{array}$20

where r_Kep and v_Kep represent the solution of Kepler’s equations for a time-of-flight τ, b_i(τ) are the Chebyshev basis functions, α_i are the fitted coefficients, and ε_mdl represents the modeling errors. In the approach described by Cortinovis et al., a polynomial of order n requires that 3n + 3 coefficients be transmitted. Our approach, outlined in (20), requires that 6 additional coefficients be transmitted to convey the initial state of the spacecraft, either as Cartesian position and velocity or as Keplerian orbital elements. This requirement results in a total of 3n_Kep + 9 transmitted elements. To maintain the same number of transmitted coefficients as Cortinovis et al., we therefore use n_Kep = n – 2.

Figure 4 compares the maximum feasible approximation interval—defined as the longest propagation interval over which model-fitting errors remain less than one-sixth of the NASA SRD requirements discussed in Section 3.3.1—between the baseline Chebyshev basis approach described by Cortinovis et al. and our modified Keplerian model with polynomial corrections. The test case shown in this figure uses the ELFO defined in Table 4. In these results, the largest feasible approximation interval is found for various starting true anomalies spaced 12° apart. We note that our approach differs slightly from Cortinovis et al. as the coefficients α_i are obtained using a simple least-squares solution rather than constrained optimization.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Comparison of feasible approximation intervals for the polynomial and Keplerian models using the ELFO described in Table 4

As shown in Figure 4, the Keplerian model outperforms the Chebyshev basis for both low and high numbers of coefficients for this orbit. The two models perform comparably near perilune, but the Keplerian model better predicts the spacecraft trajectory throughout the remainder of the orbit. Maximum performance was achieved around 100° true anomaly, corresponding to approximately 3/4 of an orbital period. This represents an ideal placement for the LNSS satellites considered in this analysis, as they could update their ephemerides shortly before becoming visible to the LSP, those ephemerides would remain valid (in terms of model error) for its entire visibility period to the user service volume. The post-fit uncertainty for the model over its validity window can be computed onboard as σ_ρ,mdl and ${\sigma }_{\dot{\rho },\hspace{0.17em}mdl}$.

3.3.3 Clock Prediction and Uncertainty Modeling

Every navigation satellite carries an onboard clock, or oscillator, that provides the reference frequency for signal generation and other time-based decisions. Any error in this oscillator will shift transmission times and carrier frequencies from their specified values, thereby introducing error into the user’s measurement. Zucca and Tavella (2005) describe these oscillator behaviors using a random walk- and run-type process model, given by the following equations:

$\begin{array}{c}{x}_{\mathrm{clk}}\left(t_{k+1}\right)=\Phi \left(t_{k+1},t_k\right)x_{\mathrm{clk}}\left(t_k\right)+c\cdot w_k\\ \begin{array}{cc}\left(\begin{array}{c}b\left(t_{k+1}\right)\\ \dot{b}\left(t_{k+1}\right)\\ a\left(t_{k+1}\right)\end{array}\right)=\left(\begin{array}{ccc}1& \tau & \frac{{\tau }^2}{2}\\ 0& 1& \tau \\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}b\left(t_k\right)\\ \dot{b}\left(t_k\right)\\ a\left(t_k\right)\end{array}\right)+c\cdot w_k,& \tau =t_{k+1}-t_k\end{array}\end{array}$21

$E\left(w_k\right)=0_{3\times 1},E\left(w_k{w}_k^T\right)=S_k\left(t_{k+1},t_k\right)=\left(\begin{array}{ccc}{\sigma }_{\mathrm{WFM}}^2\tau +{\sigma }_{\mathrm{RWFM}}^2\frac{{\tau }^3}{3}& {\sigma }_{\mathrm{RWFM}}^2\frac{{\tau }^2}{2}& 0\\ {\sigma }_{\mathrm{RWFM}}^2\frac{{\tau }^2}{2}& {\sigma }_{\mathrm{RWFM}}^2\tau & 0\\ 0& 0& 0\end{array}\right)\cdot$22

Here, $b=c\cdot \delta t$ is the oscillator phase offset from some reference time, $\dot{b}=c\cdot \delta \dot{t}$ is the phase drift or frequency offset, and a is the frequency drift or aging rate. σ_WFM and σ_RWFM correspond to white frequency modulation (WFM) and random walk frequency modulation (RWFM) processes, respectively. This model ignores both random run frequency modulation (RRFM) and flicker frequency modulation (FFM), the latter of which is not directly translatable to a finite-state linear system. Depending on the oscillator being used, a more detailed model that accounts for FFM may be required to accurately approximate the uncertainty over time. For example, Van Dierendonck et al. (1984) assume that FFM exists only within a certain frequency range and incorporate it into a two-state model, while Davis et al. (2005) model FFM as the linear combination of several Markov noise processes.

Both the user and the broadcasting satellite carry oscillators that drift from some reference time system—whether that reference is GPS time, UTC, or some other standard. When the user makes pseudorange or pseudorange-rate measurements, they estimate their own clock offset and drift, shown as b_u (t) and ${\dot{b}}_u\left(t\right)$ in Equations (1) and (2). Similarly, the satellite’s clock offset from the externally provided reference time can be described using Equation 3. These coefficients are derived from an onboard estimate of the clock state, given by Equation (21). When broadcasting new navigation message parameters at t_oc, these coefficients are assigned as $a_{f0}=b\left(t_{oc}\right)$,$a_{f1}=\dot{b}\left(t_{oc}\right)$, and $a_{f2}=a\left(t_{oc}\right)$. At t_oc, the spacecraft will have some uncertainty in the clock terms, represented by P_clk (t_oc). This uncertainty is conferred directly to the user and, like the state uncertainty discussed in Section 3.3.1, will accumulate over time. The clock state covariance can be propagated forwards in time in a process similar to Equation (11), but in this case process noise is modeled explicitly:

$P_{\mathrm{clk}}\left(t_{k+1}\right)=\Phi \left(t_{k+1},t_k\right)P_{\mathrm{clk}}\left(t_k\right){\Phi }^T\left(t_{k+1},t_k\right)+S_k\left(t_{k+1},t_k\right)\cdot$23

Based on this equation, making range and Doppler measurements from the satellite’s transmitted signal implicitly involves measuring the satellite oscillator’s state and including it as error in the user’s final measurements. Pseudorange measurements directly observe the phase bias b(t), so the clock error contribution to the pseudorange error can be written as ${\sigma }_{\rho ,\hspace{0.17em}clk}\left(t\right)=\sqrt{P_{11}\left(t\right)}$ (where t is the signal transmission time). However, it is not possible to measure an oscillator’s instantaneous frequency (Kaplan & Hegarty, 2017; Lombardi, 2017). For example, the carrier tracking loops discussed in Section 2.2.4 determine Doppler shift by differencing successive phase measurements. In an FLL, the phase itself is not continuously tracked, so Doppler is estimated by differencing measurements separated by the loop’s predetection integration time. In contrast, a PLL tracks the phase continuously, so measurements can be differenced over longer intervals.

When measuring frequency over short time intervals, the short-term stability of the oscillator, which refers to oscillator’s ability to produce a consistent frequency over a specified interval (Lombardi, 2017), becomes increasingly important. Notably, frequency stability is independent of the reference frequency and is typically quantified using either the Allan or Hadamard deviation (σ_y(τ) or σ_H(τ)). Given a measurement interval T_m, the oscillator-induced error in the velocity measurements becomes:

${\sigma }_{\dot{\rho },\mathrm{clk}}\left(t\right)=\sqrt{P_{22}\left(t\right)+c^2\cdot {\sigma }_y^2\left(T_m\right)\cdot }$24

Like the GNSS systems discussed in Section 2.2.1, the satellite clock for LNSS systems will be impacted by both general and special relativity. Fortunately, these relativistic effects are well understood and highly predictable. To account for these effects, appropriate correction parameters should be provided to users through the navigation message. Our simulations in the following sections assume this correction is implemented appropriately and do not model it explicitly. Relativity also introduces timing discrepancies between the Earth and Moon; these discrepancies will need to be corrected in any system implementation that involves time transfer between the two bodies.