Abstract
Plans to establish a satellite network around the Moon to support communication, position, navigation, and timing services are rapidly evolving. Satellites that are part of this system broadcast their ephemeris as finite parameters to lunar users for user state estimation. In this work, we investigate lunar satellite ephemeris design to identify the optimal parameterization to broadcast to a lunar user. The proposed framework directly approximates the lunar satellite position and velocity in the inertial frame and obtains the conversion parameters necessary for state representation in the lunar fixed frame. The framework leverages signal-in-space-error requirements as constraints in the parameterization process to guide the search for the best ephemeris parameter set. We evaluate the performance of our proposed framework for satellites in a low lunar orbit and an elliptical lunar frozen orbit. The performance of different methods is assessed based on the precision of the ephemeris prediction, fit interval, and message size. We showcase the ability of the developed framework to approximate satellite ephemeris for both orbits to the desired precision by adjusting the fit interval and the number of parameters to broadcast. In particular, we demonstrate that formulations with a standard polynomial basis and a Chebyshev polynomial basis produce feasible solutions for ephemeris approximation at varying epochs in orbits, abiding by signal-in-space-error requirements.
1 INTRODUCTION
Several space agencies, including the National Aeronautics and Space Administration (NASA), the European Space Agency (ESA), and the Japanese Aerospace Exploration Agency, plan to establish a network of satellites orbiting the Moon to provide reliable communication and position, navigation, and timing (PNT) services to lunar users (Schönfeldt et al., 2020). Recently, NASA and ESA introduced the LunaNet architecture concept to promote interoperability amongst industry, academia, and international partners to establish a dedicated lunar satellite-based navigation and communication system (NASA Goddard Space Flight Center, 2022a). This type of navigation system would allow users to perform autonomous lunar operations without relying on continuous interaction with terrestrial ground stations, as well as enhance overall real-time position accuracy. Satellites orbiting the Moon as part of this network are required to deliver satellite ephemeris to users over an augmented forward signal (NASA Goddard Space Flight Center, 2022b). Ephemeris data define the navigation satellite’s location and on-board clock parameters, which are required for user positioning (Misra & Enge, 2012). Ephemeris parameters that are broadcast by different navigation satellites are processed by the user receiver to compute satellite position, velocity, and acceleration at the time point that each signal was transmitted. Users combine this information with other signal-derived observables to calculate their state and clock offset.
Satellite ephemeris parameters must provide a precise representation of the true orbital trajectory and timing to provide PNT services: the error in ephemeris parameterization contributes to the error in user PNT. Therefore, the choice of parameterization must produce precise parameters that enable the desired user PNT performance. These parameters are regularly updated with high-fidelity orbit determination data for the satellite’s orbital state. Moreover, the ephemeris parameters must be provided in a common lunar reference frame or broadcast sufficient data such that the user can perform coordinate transformation to the frame of interest (NASA Goddard Space Flight Center, 2022b).
1.1 Related Works on Satellite Ephemeris Parameterization for PNT
The choice of ephemeris parameterization is a subject of research for terrestrial satellite-based navigation systems. The models discussed in this section are summarized in Table 1. Time-related parameters are omitted in counting the total number of parameters broadcast by a navigation satellite.
Ephemeris Parameterization Methods for Terrestrial Satellite Navigation
ICD: interface control document, LEO: low Earth orbit, OE: orbital elements
On Earth, global navigation satellite systems (GNSSs) provide PNT services and may choose to broadcast satellite ephemeris as a parameterization of a satellite’s orbital elements. GNSS constellations such as the Global Positioning System (GPS) and Galileo constellations employ this model (European Global Navigation Satellite Systems, 2021; Navstar GPS Directorate, 2022). The ephemeris portion of the navigation message of these constellations includes six classical Keplerian orbital elements, nine perturbation elements to capture medium Earth orbit (MEO) perturbations, and timing-related parameters (e.g., clock offset from timing reference, relativistic effects, ionospheric delays), which are used to derive a satellite state and time of broadcast in an Earth-centered Earth-fixed (ECEF) frame of reference (Navstar GPS Directorate, 2022). A set of GPS ephemeris parameters is valid for a 2-h interval and is determined by the master control station via a least-squares curve fit. In contrast, a set of Galileo satellite parameters is valid for a 3-h interval and is determined by the Galileo ground segment. Furthermore, the BeiDou GNSS constellation employs a similar parameterization methodology, broadcasting three additional elements to capture MEO, geosynchronous orbit (GEO), and inclined geosynchronous orbit (IGSO) perturbations (China Satellite Navigation Office, 2020).
Researchers have proposed modifications to the default GNSS parameterization algorithms to capture a larger variety of satellite orbits supporting Earth PNT. Reid et al. (2013) introduced a methodology that decreases the number of orbital parameters from 16 to 10; however, the fit interval, or the interval of time over which an ephemeris parameter set is valid, was reduced to the order of 10 min. Additionally, Xiaogang and Mingquan (2017) proposed a modification to the BeiDou broadcast ephemeris model based on non-singular orbit elements. This approach addresses failures of the classical broadcast framework under conditions of low eccentricity and low inclination angle. In reducing the number of parameters to 14, the authors found no loss in precision for BeiDou’s inclined geosynchronous satellite orbits, while decreasing communication resources dedicated to ephemeris.
Alternatively, satellite ephemeris may be broadcast as a direct approximation of the Cartesian satellite position, velocity, and acceleration in a set reference frame. This model includes propagation, interpolation, and extrapolation of the Cartesian satellite states. The navigation message broadcast by the Russian GNSS (GLONASS) provides position, velocity, and acceleration values in the PZ-90 coordinate frame of reference at a reference epoch time. To retrieve the satellite state at any valid time, the user numerically integrates the satellite equations of motion using the fourth-order Runge–Kutta technique to propagate the satellite orbit from the reference epoch time to the desired time (Russian Institute of Space Device Engineering, 2008). The propagation interval for which a set of ephemeris parameters is valid with this parameterization is 15 min within an epoch. The ephemeris message for the minimum operational performance standards of the L1 satellite-based augmentation system (SBAS) consists of nine parameters (the satellite’s three-dimensional ECEF position, velocity, and acceleration) and the reference epoch time. The satellite position can then be approximated by using a second-order polynomial extrapolation, which is valid for 6 min (Reid, 2017). Dobbin and Axelrad (2023) investigated numerical interpolation of orbit position and velocity with B-splines to create compact and accurate on-orbit transmitter ephemeris. The authors showed that this parameterization can provide ephemeris at a precision comparable to that of the precise orbit data distributed by the International GPS Service. However, to match the accuracy of the GPS ephemeris parameters over the same update period, this method requires approximately twice the bit volume of the standard GPS ephemeris broadcast.
Finally, orbital state approximation via numerical interpolation is applied for purposes beyond PNT. At the Jet Propulsion Laboratory, planetary and satellite ephemeris parameters are derived using Chebyshev polynomial coefficients obtained via interpolation of high-precision planetary orbit propagation (Newhall, 1989).
1.2 Challenges in Lunar Ephemeris Parameterization
The models discussed above for terrestrial satellite ephemeris parameterization are not directly applicable to the lunar navigation case. First, the methodology developed for lunar ephemeris parameterization must be flexible for the varying orbit geometries of interest to the lunar navigation community. Moreover, the differences in dynamic environments between lunar and terrestrial conditions present additional challenges in the design of satellite ephemeris, including the following:
Relatively stronger perturbations from both third-body gravitational influences and solar radiation pressure compared with Earth orbits, primarily due to the weak central gravity of the Moon (Pereira & Selva, 2020)
Lunar mass concentrations that alter the local gravitational field to the extent that most low-altitude satellite orbits are unstable or will be significantly altered on the timescale of months if no corrections are applied (Meyer et al., 1994)
Stable elliptical lunar orbits for which the orbiting satellite experiences vastly different dynamic environments between the apolune and perilune (Vallado, 1997)
To tackle these challenges, Stallo et al. (2023) proposed extrapolating position and velocity via second-order polynomial approximation to estimate the orbital positions of satellites in an elliptical lunar frozen orbit (ELFO) as part of a lunar radio navigation system. The authors found satisfactory position and velocity accuracy with this method. However, their work did not include an analysis of the ephemeris message size and examined only one type of lunar orbit.
1.3 Proposed Approach
Our work contributes to the effort of lunar satellite ephemeris parameterization. Here, we holistically investigate the design space of ephemeris modeling, considering multiple lunar orbit geometries and different parameterization models. Figure 1 exemplifies the approach of obtaining and analyzing satellite ephemeris parameterization as implemented in this paper.
In the proposed framework, we generate orbit data for an ELFO and LLO to test different orbit parameterization methods and obtain a parameterized predicted orbit; we then evaluate the performance of the methods with regard to three criteria.
This study follows three core steps: orbital data generation, orbital data fitting, and performance evaluation. First, the lunar navigation satellite trajectory is simulated through an in-house orbit propagator with integrated multi-body dynamics, solar pressure perturbations, and high-order central gravity modeling (Iiyama et al., 2023). We use these data as “truth data” to design surrogate models that directly approximate the six Cartesian components of satellite position and velocity over a bounded time interval. Orbital-element-based parameterization methods assume small variations between predicted and true orbital elements throughout the fit interval, or parameterization interval. This assumption is challenging to uphold in lunar orbits for large intervals because of strong orbital perturbations. Therefore, for this work, we consider the direct parameterization of satellite position and velocity. Note that we do not consider orbit determination errors in this study, and we assume that the orbital data to be used for parameterization fitting are known without error. Integrating orbit determination simulations (either on-board or from ground stations) in the error analysis is left as an extension for future work.
Next, we propose fitting methods to generate a parameterized predicted orbit (shown in cyan in Figure 1) from the true orbit data (shown in red in Figure 1). Because of nonlinear perturbations in the orbiter dynamics, the parameterization process introduces errors between the predicted orbit and the true orbit (shown in purple in Figure 1). These errors arise from limitations of the surrogate model to represent the full design space: the lower the error in parameterization over the entire domain, the better the model performs in providing estimated orbital data to the user. The signal-in-space-error (SISE) is the chosen metric for quantifying the ephemeris parameterization error. Following lunar relay service requirements (NASA Goddard Space Flight Center, 2022b), we define the SISE as the scalar instantaneous difference between the position, velocity, and time (PVT) of the predicted orbit and the true satellite PVT, as expressed in the lunar reference frame. Our approach leverages known lunar PNT service requirements to construct ephemeris surrogate models by finding the optimal parameter set that solves the SISE-constrained regression problem between the predicted and true orbital states. Finally, we analyze the trade-off between the precision of ephemeris prediction, parameterization interval, and message length. Different components of the broadcast message are expected to be constrained in length; therefore, a smaller message length is favorable. We can improve the precision of ephemeris by fitting the orbit in higher-order models; however, this approach increases the number of parameters sent to the user, thus increasing the message size. In contrast, decreasing the number of orbital parameters sent to a user has been shown to decrease the interval over which the parameterization is valid, requiring more parameter updates in orbit. The goal of this research is to map different ephemeris parameterization methods into the trade-off space of the three criteria and to identify the method that best supports lunar PNT system design.
1.4 Key Contributions
This work is based on our prior conference paper presented at the 37th International Technical Meeting of the Satellite Division of the Institute of Navigation (Cortinovis et al., 2023). Our key contributions are listed below.
To the best of the authors’ knowledge, this work is the first of its kind to study the parameterization of satellite ephemeris for lunar PNT activities, tackling challenges unique to the parameterization of lunar satellite orbits.
We formulate the task of ephemeris parameterization as a constrained convex optimization problem, leveraging known SISE requirements to find an accurate representation of the true satellite state in orbit.
We consider two types of basis for surrogate models (polynomial and Chebyshev) and evaluate the parameterization accuracy of the designed models by testing initialization at a variety of reference epoch times. We compare the performance of these models based on precision as well as orbital coverage.
We quantify the fit interval and size of the message broadcast by the satellites (in a 2’s-complement binary representation) to guide the selection of parameterization methodology.
We investigate ephemeris parameterization methods for two orbits of interest to the lunar navigation community: the ELFO and the low lunar orbit (LLO). These orbits vary significantly in geometry, influencing the stability and contributions of orbital perturbations throughout a satellite’s trajectory. We show that our methodology is flexible to either orbital geometry.
The remainder of this paper is organized as follows: Section 2 details the satellite dynamics and relevant reference frames for this work. Section 3 outlines the formulation of the parameterization models as constrained convex optimization problems, along with message length characterization. Section 3 also provides instructions on how to extract a satellite’s state with this approach. Section 4 addresses the setup for experimental evaluation and provides a discussion of the results. Section 5 summarizes the work accomplished and outlines future directions.
2 LUNAR DYNAMICS AND REFERENCE FRAMES
We simulate two orbits of interest to LunaNet (Schönfeldt et al., 2020), which encompass various perturbation influences because of their differing geometries. We utilize an in-house orbit propagator to retrieve true orbital data at any desired time instance. We also describe the forces acting on a satellite in orbit considered for simulation, as well as the lunar reference frames relevant to this work.
2.1 Lunar Satellite Dynamics
In orbit, a lunar satellite is subject to gravitational acceleration, calculated as the sum of higher-order gravity terms from the primary body of influence, the Moon, and third-body perturbations from the Earth, Sun, and Jupiter as follows (Montenbruck & Gill, 2013):
1
2
where Ug is the gravitational potential of the primary body (the Moon) with ∇Ug representing its gradient, G = 6.6743 × 10−11 m3/kg s2 is the universal gravitational constant, M☾ and R☾ are the Moon’s mass and radius, respectively, ru is the satellite’s position vector in orbit with respect to the Moon’s center of mass, rp is the position vector of the third-body contributors (Earth, Sun, and Jupiter) with respect to the Moon, Mp is their corresponding mass, λ is the longitude of the spacecraft position in the Moon-fixed frame, 
is the normalized Legendre function, and 
and 
are the normalized geopotential coefficients, respectively. Finally, Nsph is the maximum degree of the higher-order gravity terms. The satellite also experiences acceleration due to solar radiation pressure in lunar orbit, calculated via the cannonball model (Kenneally, 2016):
3
where ϕ⊙ = 1360 W/m2 is the solar flux at 1 AU, Au and Mu are the satellite’s surface area and mass, respectively, CR is the solar radiation pressure coefficient, c = 299792458 m/s is the speed of light, d is the heliocentric distance to the satellite, and ê⊙ is the Sun’s unit direction vector from the satellite. Further details on lunar orbital dynamics modeling are outlined in the work by Iiyama et al. (2023).
2.2 Lunar Reference Frames
There are two Moon-centered Moon-fixed reference frames commonly used by the lunar exploration community: the mean Earth/polar (ME) and principal axes (PA) frames. The ME frame has axes aligned with the Earth’s mean direction and the Moon’s axis of rotation, whereas the axes of the PA frame are aligned with the Moon’s principal axes.
To generate ephemeris for lunar surface users, we are interested in the conversion from a Moon-centered inertial (MI) frame to these Moon-fixed reference frames. In this study, we define the MI frame of reference as the J2000-oriented Moon-centered frame. Figure 2 shows how one can convert between the MI and ME frames with a 3-1-3 Euler angle transformation. Three angles are needed for rotation at time t: ϕt, θt, and ψt. These angles vary with time primarily because of the Moon’s motion about its axis of rotation, in addition to precession and nutation of the Moon’s rotational axis.
3-1-3 Euler angle transformation from an inertial (red) to rotational (blue) frame of reference
We employ a Euler angle frame rotation to enable ephemeris representation in both inertial and rotational Moon frames of reference (figure adapted from Folta et al. (2022)).
3 EPHEMERIS PARAMETERIZATION METHODOLOGY
We formulate the problem of ephemeris parameterization as a regression problem. We find the optimal basis parameter set that minimizes the difference in position and velocity between the predicted and true ephemeris subject to SISE constraints.
3.1 Unconstrained Regression for Cartesian Coordinate Parameterization
We extract the satellite orbital position (x, y, and z) and velocity (vx, vy, and vz) Cartesian components at time t from orbital data. The relationship between a position coordinate and its derivative in time (i.e., the respective velocity coordinate) is used to create three time-varying objective functions (fi(t) for i = {x, y, z}):
4
such that:
5
The domain for parameterization is t ∈ [t0, tf], where t0 is the initial time and tf is the final time for parameterization. This domain may also be characterized in terms of the satellite true anomaly, ν(t), where ν0 and νf are the true anomalies at the initial and final time, respectively. We design a surrogate model 
to represent an objective function over this domain: for a particular set of parameters (αi), the surrogate model predicts fi(t). Fitting a model to a set of design points requires optimizing for the parameters αi such that the difference between the true evaluations and those predicted by the model is minimized (Kochenderfer & Wheeler, 2019). We accomplish this by minimizing the L2 norm squared of the residual at each sampled design point. For the i-th coordinate, we minimize both the position and velocity residuals:
6
where 
is the surrogate model evaluation at each sampled point, fi is the true model evaluation at each sampled point, and 
and 
are their respective time derivatives. The optimization problem in Equation (6) is widely known as regression, which outputs the optimal set of ephemeris parameters to broadcast for the i-th coordinate, αi. Overall, we seek to find three sets of parameters (αx, αy, αz) with satellite data expressed in an inertial reference frame. We include velocity explicitly in the minimization problem to reduce potential noise that would otherwise appear if the surrogate model were differentiated alone (Newhall, 1989).
3.2 Constraint Definition Based on SISE
The SISE is the chosen metric for assessing the precision at which a surrogate model represents the true satellite state in orbit. The SISE includes contributions from relay on-board state knowledge, broadcast signal integrity, on-board path delays, and state and time representations, as well as ephemeris parameterization. These contributions are limited to the orbiter system and do not include delays and phase noise introduced by the receiver system. Although the SISE is influenced by a satellite’s clock bias and drifting with respect to the network’s timing standard, modeling time-related errors is beyond the scope of this study. We assume that these time-related components are perfectly modeled and that any corresponding parameters are implicitly broadcast to the user. Given this assumption, the SISE for position and velocity is mathematically defined as follows:
7
8
for true position (x, y, z), estimated position 
, true velocity 
, and estimated velocity 
.
SISE variations for lunar navigation satellites have been defined in the Lunar Relay Services Requirements: the SISE in position shall not exceed 13.34 m (3σ), and the SISE in velocity shall not exceed 1.2 mm/s (3σ) over any 10-s interval (NASA Goddard Space Flight Center, 2022b). We reformulate the optimization problem in Equation (6) using these requirements to constrain the magnitude of the residual at the sampled points via L∞ norms:
9
where rtol and ytol are upper bounds on the position and velocity residuals defined to ensure that the SISE contributions from the parameterization of each coordinate and frame transformation do not exceed the requirements. Ephemeris parameterization is not the only contribution to the SISE; therefore, rtol = 13.34/γ m and vtol = 1.2/γ mm/s. The scale factor γ is chosen by the designer based on the desired contribution magnitude of the the ephemeris parameterization to the SISE over the entire fit interval. For the cases investigated in this article, γ = 6, such that the residuals at the sampled points must remain below 1/6th of the total SISE budget. The optimization problem in Equation (9) is solved for each coordinate expressed in the inertial reference frame. Next, we discuss the surrogate model choice and sampling plan.
3.3 Surrogate Models
A surrogate model can be formulated as a linear combination of basis functions. In general, the surrogate model 
evaluated at a set of design points t is given as a linear combination of q-basis functions b(t) (Kochenderfer & Wheeler, 2019):
10
The objective of the least-square minimization problem in Equation (9) can be rewritten with this surrogate model formulation as follows:
11
where B is the basis matrix constructed from m data points:
Without constraints, the regression problem in Equation (6) with basis functions has an optimal parameter solution αi = B+Fi, where B+ is the Moore—Penrose inverse of B. Three bases are considered for optimization, such that the constrained optimization problem in Equation (9) over the time interval [t0, tf] has the following form:
12
Given this formulation, the parameters αi are broadcast to the user to obtain the satellite state in the inertial frame of reference corresponding to the set of basis coefficients obtained by solving the optimization problem shown in Equation (9).
3.3.1 Polynomial Basis
Polynomial basis functions consist of the product of the design vector components such that the i-th polynomial basis is as follows:
13
By setting the number of bases, n, solving the regression problem in Equation (12) provides a feasible solution of n + 1 coefficients for coordinate parameterization, totaling in 3n + 3 coefficients for position and velocity parameterization. We explore the polynomial basis because, for time periods smaller than an orbit’s orbital period, the position and velocity of the satellite follow a polynomial-like trend with respect to time. Furthermore, the computation of polynomial basis coefficients involves only exponentiation of the time value, which is relatively efficient compared with the other bases explored.
Note that polynomial bases are defined over the domain t ∈ [−1, 1]. We normalize time to ease the implementation of the chosen sampling plan (Section 3.4) such that the surrogate model evaluated at t ∈ [t0, tf] is as follows:
14
3.3.2 Chebyshev Polynomial Basis
An arbitrary, smooth, continuous function can be represented in terms of a linear combination of orthogonal polynomials. These bases provide high resolution near the domain boundaries and well-defined convergence properties. Chebyshev polynomials are part of this class of orthogonal polynomials. Chebyshev bases are strictly defined over the domain t ∈ [−1, 1], calculated recursively as follows (Rivlin, 2020):
15
such that the surrogate model evaluated at t ∈[t0, tf] is also defined by Equation (14).
Although the Chebyshev polynomial is part of the polynomial basis family, it is not redundant to explore this model along with the regular polynomial model. Apart from other orthogonal bases, Chebyshev polynomials provide near-minimax representation: the linear combination of the n + 1 Chebyshev bases will be close to the value of the function parameterization that guarantees minimization of the maximum possible parameterization error (Fraser, 1965). An n-order surrogate model provides n + 1 coefficients to solve for with regressions. This approach results in a total of 3n + 3 coefficients, which is identical to that of the standard polynomial model described above. However, it is more computationally expensive to obtain Chebyshev basis values compared with regular polynomial values because of the recursive nature of the Chebyshev basis.
3.4 Conversion to Lunar Reference Frames
The optimal parameter set is found by solving Equation (9) using data from the in-house propagator expressed in the J2000 MI reference frame. However, it is useful for a user on the Moon to obtain a satellite state in a common Moon-centered Moon-fixed reference frame, which is not inertial. With the current framework, it is not possible to generate a set of parameters for a position and velocity coordinate already in this common frame, as Equation (5) would not hold. Therefore, the user must obtain additional parameters to convert from inertial ephemeris to Moon-centered Moon-fixed reference frame ephemeris.
For the time intervals investigated, a linear approximation of the three angles ϕt, θt, ψt introduced in Section 2.2 is chosen to convert from the MI to ME frame. Separately from optimizing the ephemeris surrogate model, we extract the three Euler angles and perform linear regression, constraining the parameterization to coincide at the starting point. Mathematically, the following optimization problem is solved:
16
where 
for angle v = {ϕ, θ, ψ} and X is the basis matrix for a polynomial basis of order n = 1 for each angle. For further details on frame rotation and creating rotation matrices, we refer the reader to the work by Folta et al. (2022). Note that, in this framework, we only need to broadcast six additional parameters to the user to express ephemeris in the ME frame. If we were to directly approximate the ephemeris in the ME frame, we would require three additional parameter sets, as Equation (5) would be violated. For example, when using a polynomial basis, direct parameterization in the ME frame requires 6n + 6 parameters, which, for n > 1, produces more parameters than the proposed framework. Beyond providing a lower number of ephemeris parameters, this approach allows the user to choose either PA or ME ephemeris representation, as the conversion between the PA and ME frames is time-invariant and well characterized (Folta et al., 2022).
3.5 Sampling Plan
To reduce the computation time needed to generate ephemeris parameters, a sampling plan is used to identify which sampled positions and velocities provide sufficient coverage to precisely parameterize a satellite state for the entire parameterization period. We explore uniform sampling to obtain equally spaced samples over the parameterization interval, defined in either the time domain or the true anomaly domain. We also explore sampling at Chebyshev—Lobatto nodes: for j = 0, …, m, we sample m + 1 points (Rivlin, 2020):
17
These nodes are often used to handle Runge’s phenomenon, which refers to the presence of high-frequency oscillations in the surrogate model near the domain endpoints. To ensure that the least-squares problem is never under-determined while avoiding over-fitting, we sample n + 1 points for each n-order polynomial and Chebyshev basis tested.
To decrease the effect of Runge’s phenomenon, the samples obtained in this work are Chebyshev–Lobatto nodes. We provide an example to further motivate this choice in the context of lunar orbits. Consider the x-coordinate of the position of a satellite in an ELFO near perilune, plotted on the left in Figure 3. We attempt to find a feasible parameterization with a standard polynomial basis by solving the constrained optimization problem outlined in the previous sections. The parameterization interval is 2 h, and the basis order is 14. We solve the constrained optimization problem using the three kinds of sampling plans addressed previously and analyze the absolute truncation error in position for the three sampling plans. This error is plotted on the right in Figure 3, which also includes a dashed line at 13.34 m representative of the total SISE requirements.
(Left) x-position of a satellite in an ELFO near perilune approximated with a standard polynomial model obtained using three different sampling strategies and (right) the absolute truncation error resulting from parameterization
The truncation error from both uniform sampling strategies is larger in magnitude compared than for the Chebyshev–Lobatto sampling plan near the endpoints. Chebyshev–Lobatto sampling attenuates the higher-frequency error near the domain endpoints, such that the parameterization outside of the sampled points is sufficiently precise at all times.
3.6 Message Length Determination
We quantify the length of the ephemeris component of the navigation message broadcast by a navigation satellite in terms of bits. With our formulation, the satellite broadcasts the parameter set to represent the satellite orbital state in the inertial frame of reference as well as the parameter set to convert the state representation from the MI to the ME frame of reference. We assume that these parameters will be broadcast as 2’s-complement binary numbers.
For the inertial position and velocity model, the number of bits associated with each parameter αk is influenced by the precision necessary to meet the SISE requirements, as well as the largest expected range of that parameter. Mathematically, we determine the number of bits (Nb) associated with each parameter αk as follows:
18
where m refers to the smallest integer order of magnitude greater than or equal to the largest parameter’s order of magnitude and LSB precision is the precision of the least significant bit (LSB). The additional bit is required to represent the sign of that number. A scale factor is provided to the user to obtain the appropriate parameter unit. The overall scaling and quantization procedure is illustrated in Figure 4.
Overview of the scaling and bit quantization process
The left figure shows the original curve, and the right figure shows the scaled curve. The curve is first scaled with the magnitude of the curve and then quantified using Nb + 1 bits (1 is for the sign). The maximum quantization error is 
, which should be smaller than the LSB. Therefore, the required number of bits is calculated as Equation (18).
The metric m represents the expected magnitude range of a parameter. We seek to represent different orbits with the same message format. However, the range of αk values varies based on orbital altitude, eccentricity, and true anomaly. To ensure that we always capture every possible integer digit of a parameter with the chosen bit representation, we consider the next closest magnitude order. Furthermore, the LSB precision is set by the required SISE accuracy: we must obtain the position at the sub-meter level (= 10–4 km) and the velocity at the sub-millimeter level (= 10–7 km/s). Thus, we consider the worst-case truncation error by broadcasting αk parameters in units of km/sk and find that the LSB precision for position and velocity is 10–8.
Concerning the angular parameters, {ϕt, θt, ψt}∈[0, 2π], we consider the maximum value 2π explicitly instead of 10m to determine the message length of the initial angular displacement. For the angular rates, the Moon’s rotation about its axis dominates other rotational effects in magnitude, which is on the order of 10–6 rad/s. The most significant bit is then of order 10–5; therefore, we set m = −5 in Equation (18) to guarantee appropriate representation. The LSB precision for angular parameters is also influenced by SISE requirements: by performing the experiments outlined in Section 4, we find empirically that 10–7 rad and 10–15 rad/s are the highest allowable angular displacement and rate precision values.
Table 2 outlines the scaling factors, LSB precision, and units for the parameter sets obtained with our framework.
Scaling Factors, LSB Precision, and Units for the Parameter Sets
3.7 Optimization Pipeline
The ephemeris regression problem in Equation (12) and the linear regression problem for finding the Euler angle parameters in Equation (16) are constrained quadratic programs. Additionally, all of the objective functions, the L2 norms squared, are convex. Leveraging this property of the problem, we utilize the Python library CVXPY (Diamond & Boyd, 2016) to find the best-fitting parameters for each surrogate model. CVXPY is an open-source Python-embedded modeling language for convex optimization problems that supports the optimization of various programs. Furthermore, this library employs open-source solvers, such as the embedded conic solver (ECOS). We chose one of the default solvers for this work, ECOS, because of its ability to handle constrained quadratic programs. We tested various time intervals, n-basis orders, and orbits. Figure 5 shows the optimization pipeline implemented in this work. We note that the same basis order is set for all three coordinate parameterizations.
Optimization pipeline for ephemeris parameters
After propagating the orbit of choice in MATLAB, we set up the surrogate model and perform optimization to find all necessary parameters, storing the parameterization metrics. These metrics include starting time, parameterization interval, and message length (in bits and number of parameters).
The optimized ephemeris parameters (α, β) are uploaded to the navigation satellites and are broadcast from the navigation satellite. We assume that the reference epoch time, toe, is also broadcast or that the means to obtain the reference epoch time are delivered. The reference epoch time is the initial time at which the parameterization is valid. In Appendix A, we outline the algorithm that would be utilized by the user (receiver) to obtain the navigation satellite position and velocity in the lunar fixed frame from these parameters at time t.
4 EXPERIMENTAL SETUP AND RESULTS
Here, we evaluate the performance of our methodology in simulation through a trade study, where we vary the basis order and the initial state of the lunar satellite. We analyze the resulting feasible parameterizations in terms of accuracy, orbital coverage, and message length.
4.1 Experimental Setup
4.1.1 Simulated Lunar Orbits
We simulate two lunar orbits of interest to the lunar navigation community, LLO and ELFO, with an in-house MATLAB-based orbit propagator (Iiyama et al., 2023). Each satellite orbit is propagated according to the dynamics model outlined in Section 2. Figure 6 showcases the simulated orbits as viewed in the MI frame of reference. Table 3 outlines the dynamics model parameters used in the simulation for this research. Table 4 includes the initial conditions for each orbit in terms of the orbital elements and the orbital period instantiated at the starting epoch. The simulator leverages NASA’s information system SPICE to generate satellite ephemeris data in the lunar frames of reference described in Section 2 at the specified time.
(Left) LLO and (right) ELFO as viewed from the MI frame
Dynamics Model Parameters for Simulation
LLO and ELFO Initial Conditions (in a J2000 Oriented Moon-Centered Frame) and Orbital Period
RAAN: right ascension of the ascending node
4.1.2 Trade-Off Study Details
We evaluate the performance of the surrogate models at different initial times of parameterization t0 by varying the starting true anomaly v0 ∈ [0, 360°]. The initializations are equally spaced by 12° in the true anomaly, providing 30 different starting points for parameterization throughout the entire orbit. At initialization, we change the parameterization interval and basis order. For this study, parameterization intervals vary from 4% to more than 100% of the orbital period for both the ELFO and LLO. We present results for which the basis orders vary from n = 6 to n = 14, as no feasible solutions were found for n < 6 for the parameterization intervals investigated.
We find feasible solutions to the regression problems (Equations (12) and (16)), evaluate the SISE between the approximated and true ephemeris at every available time instance from the propagator (every 1 s), and verify whether the model produces SISE-compliant parameterization as required by the NASA Goddard Space Flight Center (2022b). First, we discuss differences in the number of feasible solutions between the ELFO and LLO. We then address the performance of the surrogate models in terms of the SISE and orbital coverage. Finally, we analyze the message length produced by SISE-compliant models.
4.2 Feasible Parameterizations for Different Orbits
Figure 7 shows the number of feasible solutions found in the ELFO and LLO as the parameterization interval relative to the orbital period changes. If a feasible SISE-compliant solution is obtained for any of the surrogate models in question for a variety of basis orders, we record that a solution was found for that initial condition and time interval. Finding 30 feasible solutions for a specific parameterization interval implies that we are guaranteed to provide precise satellite ephemeris with our proposed approach.
Feasible SISE-compliant solutions found in orbit with respect to the parameterization interval (agnostic to the choice of surrogate model)
The number of solutions found for the LLO is consistently lower than the number found for the ELFO for the majority of the parameterization intervals. We want to achieve complete orbital coverage of the ephemeris broadcast from a satellite. Therefore, in the current framework, LLO parameters must be updated at a higher frequency throughout orbit than ELFO parameters. The satellite’s Cartesian position and velocity change at a higher rate in the LLO because of the characteristic low altitude of the satellite. These variations could be captured by large basis orders, generating a high parameter number. However, this would severely increase the number of bits to broadcast to the user.
4.3 Evaluation of Surrogate Model Performance
We evaluate the performance of the surrogate models investigated in terms of the resulting SISE variation. We obtain 3σ variations in SISE for different parameterization intervals by varying the starting point in orbit and basis order, as detailed in Section 4.1.2. The highest observed variations are tabulated in Tables 5 and 6 for the ELFO and LLO, respectively.
SISE Performance of Different Surrogate Models for an ELFO for Select Parameterization Intervals
SISE Performance of Different Surrogate Models for an LLO for Select Parameterization Intervals
For the parameterization intervals investigated, we find feasible solutions with the polynomial and Chebyshev surrogate models. These two models exhibit near-identical performance, which matches our expectations given that both models are generally classifiable as polynomial surrogate models. These models employ similar function combinations to approximate true values. The SISE variations in position tend to increase with parameterization interval for the ELFO, whereas they remain consistent for the LLO. SISE velocity variations remain consistent with each other for both orbits.
For a solution to be feasible, all constraints enforced at the sampled points must be satisfied. Outside of the sampled points, the residual between the true and predicted Cartesian satellite state in orbit does not need to satisfy the prescribed constraints. Because of how these constraints are enforced, the contribution to the velocity SISE from ephemeris parameterization exceeded the allocated budget (1/6th for γ = 6) for all data reported in Tables 5 and 6. However, the contribution to the position SISE exceeded the allocated budget only for the ELFO case at a parameterization interval of 8 h. All solutions remained below the total SISE requirement for lunar navigation satellites.
We can further increase the precision of parameterization by enforcing tighter constraints, increasing the scale factor γ. We can also consider reformulating the constraints as chance constraints, ensuring that the probability that the residual exceeds the allocated SISE budget remains below 99.73% (3σ).
4.4 Orbital Coverage
We record the parameterization interval associated with each SISE-compliant solution found from a specific initialization in orbit at a particular basis order. For each initialization, we then obtain the highest parameterization interval to assess how the starting parameterization location impacts the update frequency and number of parameters.
Figures 8 and 9 show the results of this analysis for the polynomial and Chebyshev surrogate models for both orbits. Only the smallest and largest basis orders investigated are shown in the figures. As noted in the previous section, both models exhibit identical performance, as they belong to the same surrogate model family.
Maximum parameterization interval at different true anomaly initializations for (left) polynomial and (right) Chebyshev surrogate models in an ELFO
Maximum parameterization interval at different true anomaly initializations for (left) polynomial and (right) Chebyshev surrogate models in an LLO
As the time interval for parameterization increases, the solution requires a higher order of basis functions, thus increasing the number of parameters to transmit. No feasible solutions are found for the LLO for intervals larger than a quarter of the orbital period (approximately 30 min). This phenomenon is observed for the ELFO for intervals above six-tenths of the orbital period (approximately 8 h). LLOs have small orbital periods; thus, for larger time intervals, positions and velocities evolve in a periodic trend. These two models are unable to simultaneously capture this trend and meet SISE requirements.
As shown in Figures 8 and 9, the lowest basis order investigated is not able to provide solutions at every starting point in orbit. In particular, for the ELFO, the lowest basis order provides solutions concentrated around apolune, at 180° in the true anomaly. We further verify the orbital coverage provided by a specific basis order for these two models: we record the instances at which a basis order was capable of providing a solution for the lowest parameterization interval. Table 7 summarizes these results for the ELFO and LLO, demonstrating that as the basis order increases, the orbital coverage increases as well.
Orbital Coverage of Polynomial and Chebyshev Surrogate Models at Different Basis Orders for LLO and ELFO
Considering the histograms above, near apolune (where v = 180°), lower-basis-order models produce feasible and SISE-compliant solutions. In contrast, near perilune (where v = 0°), we find that we need higher-basis-order models. We often cannot find feasible solutions with low-basis-order models. This behavior is explained by the differences in dynamic environments between apolune and perilune. The satellite experiences relatively stronger gravitational attraction to the Moon when close to its surface, i.e., near perilune. Specifically, the acceleration term due to high-order lunar gravity characterization increases in magnitude closer to the Moon’s surface, as shown in Equation (2). The satellite moves relatively faster in orbit at perilune than it does at apolune, which implies that position and velocity trends vary faster over the same time period at perilune than at apolune. Therefore, a higher basis order is required to capture this effect.
4.5 Message Length Evaluation
Along with model accuracy, we evaluate the required length of the message of SISE-compliant parameterizations found for varying initial parameterization times, basis orders, and parameterization time intervals. The message length is calculated following the procedure from Section 3.6, and the results are plotted in Figure 10.
Surrogate model message length for (left) ELFO and (right) LLO
The Chebyshev basis surrogate model consistently outputs a smaller message bit length than the polynomial basis surrogate model. We also observe increasing message length with increasing basis order and fit interval, although the latter does not always occur (note the ELFO orbit for higher n).
As expected from Section 3, the order of the maximum basis (n) is the primary driver determining the parameter number and total bits of message length. Furthermore, as the parameterization interval increases, the bit message length generally increases for the same n. We observe that for most cases, the parameter magnitude increased with longer parameterization periods for the same n, leading to an increase in message bit length.
We previously observed that both models offer comparable SISE variations for the time intervals and basis orders investigated. However, the Chebyshev basis model consistently requires fewer bits for message length for the same n and parameterization interval than the standard polynomial basis model. The recursive basis construction characteristic of Chebyshev polynomials tends to reduce the parameter magnitude. Because of this property, there are inherent coefficients for each time component. Multiplying these coefficients by α parameters decreases the maximum parameter magnitude, thus decreasing the overall message bit length.
Overall, given the comparable orbital coverage and ephemeris precision performance, the different in message length sets the two models apart. We prefer a lower message length. The smaller the message length, the faster a user can extract the necessary parameters to determine satellite ephemeris. Therefore, the Chebyshev polynomial model is favored for lunar PNT applications because of its consistently smaller message length compared with the standard polynomial model. To facilitate computation, the user could store Chebyshev coefficients of the greatest expected basis order and avoid a recursive coefficient calculation each time satellite states are obtained.
5 CONCLUSIONS
Precise ephemeris parameterization is necessary to enable lunar PNT and communication. In this work, we investigated ephemeris parameterization with a surrogate-model-based approach to tackle challenges unique to the parameterization of lunar orbits. We investigated the trade-off between the precision of ephemeris prediction, fit interval, and message length. In this study, we found optimal parameterization for the ELFO and LLO with polynomial and Chebyshev polynomial bases that satisfy NASA requirements across different parameterization intervals. Polynomial and Chebyshev polynomial models offer comparable precision and orbital coverage performance. However, the Chebyshev model produces a broadcast message with fewer bits to broadcast. This framework provides sufficient parameters for a lunar user to find the satellite position in the inertial reference frame or a Moon-fixed Moon-centered frame of reference.
Future work includes combining orbit determination and time synchronization with ephemeris parameterization, evaluating the effect of unmodeled perturbations on the performance of the parameterization methods via covariance analysis. Furthermore, we may consider modeling time-related parameters to transmit to the user along with position and velocity parameters, accounting for light-time transferring delays and relativistic effects to extract the time of ephemeris transmission. Finally, we may investigate reformulating the constraints in our proposed framework as chance constraints, ensuring that the SISE budget allocated to the satellite ephemeris remains below a given threshold on a probabilistic basis.
HOW TO CITE THIS ARTICLE
Cortinovis, M., Iiyama, K., & Gao, G. (2024). Satellite ephemeris parameterization methods to support lunar positioning, navigation, and timing services. NAVIGATION, 71(4). https://doi.org/10.33012/navi.664
ACKNOWLEDGMENTS
The authors thank the members of the Stanford Navigation and Autonomous Vehicles (NAV) Lab for insightful discussions that contributed to this work. We also thank Dr. Tara Mina and Kaila Coimbra for reviewing this paper.
A: USER ALGORITHM FOR STATE DETERMINATION IN THE ME FRAME
A satellite providing navigation services to lunar users broadcasts the parameter sets (α, β) and reference epoch time toe such that any user can determine the satellite state in orbit in either the MI or ME frame of reference. Table 8 specifies the necessary operations for obtaining the satellite position and velocity in the ME frame of reference.
Receiver-Side Algorithm for Obtaining the Navigation Satellite Position and Velocity in the Lunar
Fixed Frame from the Coordinate-Based Parameter Set at Time t
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.