Simultaneous Estimation of Lunar Ephemeris and Satellite Orbits Using a DRO-Based LiAISON Method in Cislunar Space

2.1 Concept of LiAISON

Hill proposed the principle of LiAISON (Hill & Born, 2007), demonstrating that by accumulating SST measurements over a period of time, it is feasible to determine the absolute states of satellites in inertial space when one of the satellites is located in an asymmetric gravitational field. Within such an asymmetric gravitational acceleration field, the orbit of the satellite assumes distinctiveness, that is, the shape (eccentricity) and size (semi-major axis) of the satellite orbit correspond one-to-one with the absolute orbital orientation (including ascending right ascension, perigee amplitude, and orbital inclination). The accumulated SST measurements can determine the size and shape of the orbit and the phase position of the satellite in this orbit, allowing the absolute orbital orientation to be uniquely computed.

There have been many extensive investigations into the application of LiAISON's autonomous orbit determination principle in Earth–Moon space. Hill and Born (2007) were the first to verify the feasibility of autonomous navigation by establishing an SST link between a satellite in a halo orbit about the Earth–Moon Lagrange point L1 and a lunar orbiter. Parker et al. (2012) proved that the use of the LiAISON method can increase the navigation accuracy of navigation satellites at the Earth–Moon L1 and L2 translation points. Leonard (2015) explored the application of LiAISON in manned spacecraft. Research has shown that this method can improve navigation accuracy and reduce the number of ground stations required for navigation. Hesar et al. (2015) studied the accuracy of absolute navigation on the far side of the Moon using LiAISON. Yin et al. (2024) indicated that the LiAISON navigation method can improve the accuracy of asteroid gravity field estimations. The CAPSTONE (Thompson et al., 2022) mission was launched in 2022 by NASA for the purpose of conducting verification of the autonomous orbit determination principle through SST links.

DROs, which are a specific class of three-body dynamic orbits, are characterized by trajectories in which the spacecraft rotate simultaneously around both the Earth and the Moon. These orbits are located within an asymmetric gravitational field in the Earth-Moon space, where the gravitational accelerations arising from both the Earth and the Moon are comparable. When a DRO satellite establishes SST links with other satellites, their absolute states can be accurately determined. Wang et al. (2019) and Liu et al. (2021) investigated the accuracy of the absolute state of DROs determined via LiAISON navigation. Their research demonstrated that an asymmetric gravitational field is critical for orbit determination of the satellite.

2.2 Extension to Lunar Ephemeris Estimation

Existing LiAISON studies have shown that LiAISON navigation depends on a precise lunar ephemeris, which provides orbit dynamical computations. This paper investigates the possibility of DRO-based LiAISON with simultaneous estimation of the lunar ephemeris. The theoretical foundation of this approach resides in three coupled phenomena:

1. Lunar positional inaccuracies directly perturb gravitational potential gradients throughout DRO trajectories: The position of the Moon shows a strong correlation with the gravitational field around DROs, and any deviation in the lunar position will have a significant impact on the DRO orbital dynamics.

2. Resultant orbital perturbations generate distinctive signatures in SST measurements: When biased DRO trajectories are used to compute SST measurements, substantial discrepancies emerge between actual observation and calculated measurements (referred to as prior O-C [observation minus calculation] residuals).

3. These measurement residuals can be used to estimate information about both satellite states and ephemeris errors: According to the principle of LiAISON, SST measurements from gravitationally asymmetric DROs exhibit unique sensitivity to absolute states. Thus, the prior O-C residuals enable simultaneous refinement of satellite states and lunar ephemeris parameters through their inherent dynamical correlations.

Therefore, it is feasible to use the DRO-based LiAISON method to estimate the lunar ephemeris, thereby eliminating lunar ephemeris dependencies.

To leverage SST measurements for lunar ephemeris estimation, it is crucial to identify optimal SST links. Considering the unique three-body characteristics and lunar ephemeris sensitivity of DRO satellites, one potential scenario may be to establish an SST link between two DRO satellites. Additionally, the inclusion of a LEO satellite offers distinct advantages; its short orbital period induces rapid geometric evolution in the link, thereby enhancing observability and accelerating convergence. It is worth noting that Li et al. (2022) sufficiently validated the orbit determination performance of this scenario with a known lunar ephemeris. We demonstrate that this scenario is a critical configuration for enhancing the observability and convergence speed of simultaneous lunar ephemeris estimation.

Capitalizing on the complementary strengths of DRO sensitivity and LEO geometric diversity, this paper systematically investigates the feasibility of DRO–DRO, DRO–LEO, and combined DRO–DRO–LEO links. Figure 1 illustrates these scenarios in the Earth–Moon rotating coordinate frame (Wang et al., 2019), presenting a visual representation of their respective characteristics. Crucially, Figure 1 depicts the distribution of the perturbation ratio introduced by Hill (2007), which can indicate the strength of the third-body asymmetry. It can be observed that the DRO satellites are situated in a region exhibiting significant asymmetry (approximately 20%). This high degree of asymmetry is a fundamental prerequisite for the LiAISON principle, thereby justifying the selection of DROs as ideal orbits for performing both autonomous orbit determination and lunar ephemeris estimation.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

SST link scenarios in the Earth–Moon rotating coordinate frame (DU: Earth–Moon distance)

2.3 Sensitivity of Lunar Ephemerides

This article focuses on estimating lunar ephemerides, which involves estimating the position and velocity of the Moon in the Earth-centered inertial (ECI) frame. Because deviations in lunar position directly impact computations of satellite position, we define sensitivity as the degree to which satellite position deviations are influenced by lunar position deviations. Furthermore, sensitivity plays a crucial role in assessments of the feasibility of lunar ephemeris estimation; thus, this section systematically investigates the sensitivity across various scenarios. The DE430 ephemeris is used as a true lunar position, and lunar position deviations of 10 m and 100 m are added to the orbit dynamic model to assess the DRO position integration deviation.

As illustrated in Figure 2(a), the presence of a constant lunar position deviation leads to a divergence in the DRO orbital deviations over time. Notably, the rate of divergence exhibits a periodic oscillation that aligns with the DRO orbital period, indicating that the satellite's sensitivity to lunar ephemeris varies with its orbital phase. Additionally, Figure 2(a) qualitatively demonstrates that the magnitude of the satellite's deviation scales with the lunar position deviation. This relationship is validated in Figure 2(b), where a clear linear proportionality between the lunar position deviation and the resulting DRO position deviation is shown. Specifically, over 50 d, a deviation of 10 m in lunar position results in DRO orbital deviations exceeding 300 m, while a deviation of 30 m amplifies the orbital error to the kilometer scale. The observed linear proportionality shown in Figure 2 is theoretically supported by the analytical derivation presented in the Appendix. The dynamics of the DRO, influenced by lunar gravity, can be linearized for small perturbations. Equation (A5) presents the partial derivative of the satellite's acceleration with respect to the lunar position, which is analytically defined and well-behaved. Consequently, error propagation follows a first-order linear mapping, governed by the state transition logic. Because the Jacobian matrix ∂a∂rM (derived in Equation (A5)) is independent of the error magnitude ΔrM, the resulting satellite deviation scales linearly with the lunar position deviation, matching the empirical trend shown in Figure 2.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Integration position deviation of a DRO satellite

These results clearly demonstrate the high sensitivity of DRO satellites to lunar position deviations. Consequently, it can be inferred that establishing SST links between a DRO satellite and another satellite is a feasible approach for detecting lunar position deviations.

To investigate the correlation between SST range variation and lunar ephemeris deviation in different directions, we simulated lunar ephemeris deviations of 10 m and 100 m in three orthogonal directions of the local vertical local horizontal (LVLH) coordinate system. In the LVLH coordinate system, the origin is located at the center of the satellite (or the Moon), the R-axis points in the direction of the satellite's (or lunar) position vector (i.e., radial), the N -axis points in the direction of normal to the orbit plane (i.e., normal), and the T -axis is determined by a right-handed coordinate system (i.e., tangential). Subsequently, DRO–LEO range series were computed by the perturbated orbit dynamical model based on the deviation of the lunar ephemeris. This approach enables a comprehensive examination of lunar ephemeris deviations in various directions and their influence on the computed SST ranges.

As shown in Figure 3, the deviations in lunar ephemeris in three orthogonal directions have a significant impact on the calculated DRO–LEO ranges. Because the SST range deviation is essentially the geometric projection of the satellite position deviation onto the line-of-sight vector, the linear dependency of the position deviation on the lunar ephemeris error directly propagates to the measurement domain. Specifically, a 10-m deviation in each direction will impact the ranges by several hundred meters, whereas a 100-m deviation will lead to a kilometer-level error in the ranges. If the measured SST ranges do not align with estimated satellite states, this indicates inconsistency between the measurements and the dynamic model. By adjusting the lunar position estimates, consistency between the three-body orbital dynamic model and measurements can be achieved. These adjustments enable accurate simultaneous estimation of the satellite orbits and lunar ephemeris. As illustrated in Figure 2 and Figure 3, satellite states and calculated ranges are influenced by lunar ephemerides, highlighting the effectiveness of this method in lunar ephemeris estimation.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

LEO–DRO range deviations due to deviations in lunar position for the case in which lunar position deviations of (a) 10 m or (b) 100 m are added in the R, T, and N directions

3.1 Estimated Parameters and Measurements

In the process of orbit determination and lunar ephemeris estimation, estimated parameters include the position vectors r1,r2 and velocity vectors v1,v2 of two satellites, lunar mass center position rM, and lunar mass center velocity vM. All states are expressed in an ECI frame. Thus, the estimated parameter X0 at the initial time point t0 can be expressed as in Equation (1):

X0=[r1(t0)v1(t0)r2(t0)v2(t0)rM(t0)vM(t0)]1

where r and v represent the position and velocity states of the satellites and the Moon, respectively.

In this paper, a dual one-way ranging measurement model is used for orbit determination, which can essentially eliminate any clock offset. Therefore, the effect of clock offsets is not considered in this study. A description of the measurement model and its linearization has been introduced in previous studies (Li et al., 2022). The measurement value z is the combined pseudorange received by Satellite 1 and Satellite 2 at the same time, which can be expressed as in Equation (2) (Montenbruck & Gill, 2000):

z(t)=P12(t)+P21(t)2 =12⋅[‖r1(t)−r2(t12)‖+‖r2(t)−r1(t21)‖+D+ε(t)] =12⋅(ρ12(t)+ρ21(t)+D+ε(t)) =ρ(t)+D+ε(t)22

where P12 and P21 denote the pseudorange measurement values received by the two satellites. t denotes the time at which Satellite 1 receives a signal from Satellite 2, and the signal sending time is t12 and t21.ρ12 and ρ21 denote the geometric distance between the two satellites obtained by Satellite 1 and Satellite 2, respectively. ρ denotes the combined geometric distance. D represents the influence of equipment delays in the combination. This paper assumes that the equipment delays are a pre-calibrated constant, owing to the slow variation of satellite temperature. The thermal noise is given as ε∼N(0,σε), and the thermal noise standard deviation is set to σρ=0.5m. To verify the rationality of the assumed ranging noise under cislunar distances, particularly considering the signal attenuation caused by long-distance propagation, we performed a link budget analysis.

Inter-satellite ranging techniques are critical for satellite navigation and formation flying. It is well established that ranging performance generally degrades as the inter-satellite distance increases. This degradation primarily results from a reduction in the received carrier-to-noise density ratio (C/N₀) due to path loss, as measurement noise is inversely related to C/N₀. Therefore, the validity of the 0.5-m noise assumption is critically dependent on whether a sufficiently high C/N₀ can be maintained at cislunar distances.

The received C/N₀ of the SST link can be calculated using the following link budget equation (Yang, 2025):

C/N0=EIRP+(G/T)−Ls(ρ)−Lmargin+228.6 =EIRP+(G/T)−20log10⁡(4πρλ)−Lmargin+228.63

where Ls(ρ) denotes the free-space path loss varying with distance ρ,λ denotes the wavelength and is set to 26 GHz, L_margin represents the system margin accounting for pointing errors and other implementation losses and is conservatively set to 3 dB in this study, and the constant 228.6 corresponds to the Boltzmann constant in decibels. Unlike the omnidirectional broadcast architecture typically adopted by GNSS constellations, this study assumes a point-to-point SST link employing high-gain directional antennas. Ka-band parabolic antennas are utilized to compensate for the free-space path loss associated with cislunar distances. The effective isotropic radiated power (EIRP) of the transmitter is set to 40 dBW, which is achievable using, for instance, a 10-W solid-state power amplifier combined with a 0.22-m-diameter parabolic antenna (approximately 30-dBi gain). On the receiver side, the gain-to-noise temperature ratio (G/T) is assumed to be 15 dB/K, corresponding to a representative 0.6-m-diameter parabolic antenna.

Based on these parameters, the C/N₀ values are computed using Equation (3) as a function of the SST distance.

The computed results are illustrated in Figure 4. Even when the SST distance extends to 5×10⁸ m, the link maintains a C/N₀ of approximately 45.9 dB-Hz. Recent studies by Yang (2025) have demonstrated that for Ka-band SST links, a ranging precision better than 0.06 m (1σ) can be achieved when the C/N₀ is maintained at approximately 45 dB-Hz.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Carrier-to-noise density ratio versus SST distance

In summary, the theoretical thermal noise limit at lunar distances is significantly lower than the 0.5-m noise standard deviation assumed in this study. This confirms that our assumption is reasonable and includes a substantial conservative margin.

3.2 Orbital Dynamic Model of Satellites and the Moon

Given the position and velocity states at the initial time point, the propagation of estimated states over time is determined by numerical integration of a set of ordinary first-order differential equations, as given in Equation (4):

dX(t)dt=f(X(t),u(t),t)=(v(t)a(t,r(t)))4

Here, the process noise u(t) is zero-mean Gaussian white noise, used to compensate for errors in the dynamic model. f denotes the dynamic equation of the satellite and the Moon, and r(t),v(t) denote the position and velocity of the satellite and the Moon, respectively. a(t,r(t)) denotes the acceleration of the satellite and the Moon calculated by the dynamic model.

(2) Dynamic Model of the Moon

This paper utilizes the lunar ephemeris derived from the JPL DE430 ephemeris as the true trajectory of the Moon. This ephemeris establishes a lunar dynamic model of exceptional precision, encompassing the gravitational acceleration of prominent celestial bodies, relativistic effects, non-spherical disturbances in substantial bodies, and the reciprocity between Earth's tides and the Moon. The specific calculation methods employed for each acceleration have been described in the literature (Folkner et al., 2014).

Based on the DE430 framework, we constructed a dynamic model that explicitly incorporates the primary gravitational forces acting upon the Moon. Although this model is relatively simplified in complexity compared with the full DE430 implementation, it captures the dominant dynamical environment required for the proposed navigation method. A detailed comparison between the high-fidelity DE430 model and the model implemented in this work is presented in Table 2.

As listed in Table 2, the adopted model focuses on the principal accelerations acting upon the Moon, treating other factors as secondary effects. The dynamic model in this work considers the gravitational influences of the Earth a_E and the Sun a_S (modeled as point masses), a_ME produced by the gravity of the Moon on the Earth, a_ESH arising from non-spherical perturbations of the Earth, and a₀ resulting from the gravitational pull of other planets.

Figure 5 presents the magnitude of the principal accelerations experienced by the Moon at 0:00 on January 1, 2024. Considering the orbital stability of the Moon, these values illustrate the approximate scale of forces acting on the lunar body. According to the information presented in Figure 5, it is evident that the primary accelerations experienced on the Moon stem from the gravitational forces exerted by the Earth (a_E) and the Sun (a_S), as well as the reactive force of lunar gravity on the Earth (a_ME). The accelerations resulting from other forces are comparatively insignificant, measuring approximately 1×10−9m/s2, which is four to five orders of magnitude smaller than those generated by the principal forces.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Magnitude of acceleration experienced on the Moon

As indicated in Table 2, the dynamic model of this work neglects perturbations from minor asteroids, relativistic corrections, and tidal effects. To justify this simplification, we performed a quantitative magnitude analysis. To evaluate the modeling accuracy of the lunar dynamics model, the lunar trajectory obtained by integrating the model for 50 d was compared with the DE430 ephemeris. The results of this comparison are illustrated in Figure 6.

View this table:

• View inline
• View popup

Table 2 Comparison of Lunar Dynamical Model Used in This Work and the High-Fidelity DE430 Ephemeris

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Position residuals in lunar orbits between the simplified dynamic model and the JPL DE430 ephemeris over 50 d

The results indicate that the integration error for the lunar trajectory extrapolated over a 50-d period using the dynamic model is within 700 m, equivalent to an acceleration of 7.5017×10−11m/s2. This precision can meet typical real-time deep-space autonomous navigation requirements, where operational standards for lunar missions demand positional precision at the 100-m-level to kilometer-level (Turan et al., 2022b). On the one hand, frequent measurements (60-s intervals in this study) enable the EKF to rectify modeling errors inherent in our simplified dynamics through continuous data assimilation (Montenbruck & Gill, 2000). The modeled residual acceleration accumulates a positional drift of only 1.35×10−7m per 60-s interval, substantially below the typical range noise. On the other hand, even under extreme conditions such as a 5−d measurement outage, the dynamical drift is below 10 m for lunar ephemeris. According to the analysis shown in Figure 2, DRO satellites will also experience a drift of less than 10 m. Thus, it can be concluded that the dynamic model can be adopted for lunar ephemeris estimation and real-time satellite navigation.

As a result, the acceleration model of the Moon can be expressed as in Equation (6):

aM(t,r1(t))=aS+aE−aME6

where:

aME(t)=−μM⋅rM(t)‖rM(t)‖37

3.3 Extended Kalman Filter

The EKF (Senne, 1972) used in this study represents a type of fundamental sequential estimation method and is suitable for real-time satellite navigation. Given the initial state X₀ and state covariance P₀ at the initial time t₀, the EKF processes SST measurements at consecutive measurement epochs. This process comprises two steps. The first step of the EKF is the time update, where the state and covariance evolve from the previous time t_i to the current time ti+1, as given in Equation (8):

{Xi+1∣i=X(ti+1,X(ti)=Xi)Pi+1∣i=Φi+1∣iPiΦi+1∣iT+Γi+1∣iQuΓi+1∣iT8

Here, Φi+1∣i is the transition matrix; its calculation is shown in Equation (A1). Q_u is a state noise compensation matrix, which is a diagonal matrix, and the square root of its diagonal elements is the standard deviation of the process noise. The value of the standard deviation of the process noise is set according to the magnitude of the error of the unmodeled dynamic model. Γi+1∣i is the process noise transform matrix; its calculation is shown in Equation (9):

Γi+1∣i=∂X(ti)∂u(ti−1)=[Γi+1∣i106×306×306×3Γi+1∣i206×306×306×3Γi+1∣iM]9

where:

Γi+1∣i1=Γi+1∣i2=Γi+1∣iM=[Δt22×I3×3Δt×I3×3]10

Here, Δt represents the integration step size, and I3×3 is a 3×3 unit matrix.

The second step of the EKF is the measurement update, where the residuals between the measurement data zi+1 and the calculation h(Xi+1∣i) are used to obtain the updated states Xi+1 and the covariance matrix Pi+1 at the measurement time. The measurement update is performed as shown in Equation (11):

{Ki+1=Pi+1∣iHi+1T(Ri+1+Hi+1Pi+1∣iHi+1T)−1Xi+1=Xi+1∣i+Ki+1[zi+1−h(Xi+1∣i)]Pi+1=(I−Ki+1Hi+1)Pi+1∣i(I−Ki+1Hi+1)T+Ki+1Ri+1Ki+1T11

Here, Ki+1 is the Kalman gain, and the design matrix Hi+1 denotes partial derivatives of the measurement data relative to the estimated parameters, as expressed in Equation (12). Ri+1 is the measurement noise covariance matrix, and the diagonal elements of this matrix are the variance of the measured thermal noise:

Hi+1=[∂ρ(ti)∂r1(ti)∂ρ(ti)∂v1(ti)∂ρ(ti)∂r2(ti)∂ρ(ti)∂v2(ti)01×301×3]12

The specific calculation method for Hi+1 has been reported in previous work (Li et al., 2022). h(Xi+1∣i) is calculated by the states Xi+1∣i after a time update. The calculation of h(Xi+1∣i) is given in Equation (13):

h(Xi+1∣i)=12[‖r1(t)−r2(t12)‖+‖r2(t)−r1(t21)‖+c⋅D]13

Here, r₁ and r₂ denote the position states of the two satellites after a time update. Assuming that the equipment delay D is calculated beforehand, this study will directly apply the results of the calibration.

3.4 Observability Analysis Method

The Lie derivative matrix is commonly used for assessing system observability. However, the Lie derivative matrix may not be suitable for estimating the dimension of the estimated state in systems with numerous variables. In contrast, the Gram matrix (Turan et al., 2022a) effectively captures the sensitivity of observations to the estimated states. In the context of discrete nonlinear systems, the definition of the Gram matrix W aligns with the expression given in Equation (14):

W=∑k=0MΦtk∣t0THkTHkΦtk∣t014

Here, M and k represent the number of all observations and the observation index, respectively. t₀ represents the initial time, and t_k represents the time of the k-th observation.

Generally, a reversible Gram matrix indicates system observability, with two indicators commonly used to measure observability. The first indicator is the observability index (OI), represented by the minimum eigenvalue (Butcher et al., 2017). Eigenvalues typically characterize the strength of observability in different directions. If the system's OI approaches zero, this suggests that the observability in some directions is extremely poor, implying that the system may be unobservable. A larger minimum eigenvalue indicates that the system's states can be better estimated from the observations. The second factor is the condition number (CN), which is the ratio of the largest eigenvalue to the minimum eigenvalue. This index evaluates the estimation error caused by measurement noise. When this value is large, even slight measurement noise can result in significant estimation errors, indicating an ill-conditioned system (Trefethen & Bau, 1997). A larger CN of the observability matrix indicates poorer system observability. Both indicators provide quantitative measures of observability.

4.1 Simulation Settings

This paper applies three scenarios to assess the impact of different links on autonomous orbiting and lunar ephemeris estimation, and results obtained from these analyses can provide support for future practical applications. The simulation scenarios include the LEO–DRO scenario, the DRO–DRO scenario, and the DRO–DRO–LEO scenario. Each scenario employs a combination of dual one-way pseudoranges as measurements, as shown in Equation (2). In the DRO–DRO–LEO scenario, SST links are present between all three satellites, extending to three measurement sets. The simulation time is set to 50 d with a fixed 60-s measurement interval, where measurement interruptions are considered in the case of Earth–Moon geometric occultations that block the SST link (the SST link with LEO disruptions characterized by 1.5-h operational cycles featuring 0.5-h interruptions).

Specifically, a Sun-synchronous orbit with an altitude of 5×105m and an orbital period of 1.5 h was selected to evaluate the satellite orbit determination and lunar ephemeris. Two DROs with a resonance ratio of 2:1 were selected (Wang et al., 2019); this setting allows the satellites to complete two lunar revolutions per month, resulting in an orbital period of approximately 13.5 d. This type of DRO can reach a maximum distance of 4.6×108m and a minimum distance of 2.8×108m from the Earth's center, while the distance to the Moon varies between a maximum of approximately 1.0×108m and a minimum of approximately 6.7×107m. Conventional DROs typically align parallel to the Earth–Moon plane to maintain long-term stability, but exhibit poor sensitivity to lunar ephemeris deviations in the normal direction. To enhance normal sensitivity, DRO-A was augmented with a 2×107m amplitude normal to the Earth–Moon plane, whereas DRO-B maintained near-parallel alignment (≤ 5° inclination). As shown in Figure 7, DRO-A achieves an inclination amplitude of up to 15° throughout its orbital period, which can sensitively detect lunar normal-direction errors.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

Time series of orbital inclination for DRO satellites relative to the Earth-Moon plane

4.2 Observability Analysis Results

To assess the observability for estimating lunar ephemeris, this section compares the CN and OI of the Gram matrix between different scenarios with and without lunar ephemeris estimation.

Figure 8 presents time series of the CN and OI of the Gram matrix for different scenarios without estimating the lunar ephemeris.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

Time series of the CN and OI of different scenarios without estimating lunar ephemeris for (a) the LEO–DRO scenario, (b) the DRO–DRO scenario, and (c) the DRO–DRO–LEO scenario

As demonstrated in Figures 8(a)–(c), the OI values exhibit power-law growth beyond unity across all observation scenarios, indicating significantly enhanced observability with accumulated observation data. Notably, there are considerable differences in the OI growth rates across the different scenarios. Specifically, both the LEO–DRO scenario and the DRO–DRO–LEO scenario show rapid OI increases, gaining approximately 15 orders of magnitude within hours and stabilizing (OI ≥ 1) within 10 d. In contrast, the growth rate of the OI in the DRO–DRO scenario exhibits slower growth (Figure 8(b)), requiring approximately 25 d to reach the same level. This phenomenon can be attributed to the short orbital period of LEO satellites, which results in rapid geometric changes in SST measurements, enabling rapid changes in the observation matrix. Consequently, the scenarios involving LEO satellites achieve a high level of observability in a shorter time.

In terms of CN variations, the CN values of all scenarios decrease according to a power-law trend as the observation duration increases, further corroborating the conclusion that observability improves over time. Consistent with the OI trends, the CN convergence occurs significantly more quickly for the LEO–DRO scenario than the DRO–DRO scenario. Additionally, the DRO–DRO–LEO scenario, benefiting from multi-satellite cooperative observation, accelerates the CN convergence process—reaching the same convergence level (1×1018) as the LEO–DRO link in just 1 d. Furthermore, the initial CN values differ markedly: DRO–DRO (1×1027)> LEO–DRO (1×1023) > DRO–DRO–LEO (1×1022). After final convergence, the CN values for the LEO-involved scenarios are one order of magnitude lower than that of the DRO–DRO scenario. From Figures 8(a)–(c), it is evident that incorporating LEO satellites significantly enhances the system's short-term observability, whereas the multi-satellite scenario (DRO–DRO–LEO) further improves the convergence efficiency, achieving both the lowest initial CN and the highest convergence speed.

Figure 9 presents time series of the CN and OI of the Gram matrix for different scenarios when the lunar ephemeris is estimated.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

Time series of the CN and OI of different scenarios when estimating lunar ephemeris for (a) the LEO–DRO scenario, (b) the DRO–DRO scenario, and (c) the DRO–DRO–LEO scenario

Compared with Figure 8, it is evident that the OI growth rate across all scenarios has declined and the CN convergence process is delayed for the case shown in Figure 9. This finding indicates that lunar parameter estimation will reduce system observability. This phenomenon can potentially be explained from two perspectives. On the one hand, the estimation of lunar ephemeris increases the number of estimated parameters without augmented measurements, which may require more time for parameter separation and estimation. On the other hand, when the lunar ephemeris is not estimated, the estimated parameters solely pertain to the position and velocity of the satellite, which are directly related to the measurement data. Consequently, the design matrix H does not exhibit any columns consisting entirely of zeros. However, when the lunar ephemeris is estimated, the design matrix H exhibits the presence of columns comprised entirely of zeros, as exemplified in Equation (12), increasing the time required for matrix W to achieve full rank. Despite short-term observability degradation, the results of CN and OI still indicate that these parameters are observable.

Further comparison of Figures 9(a)–(c) shows that, similar to the results in Figure 8, the OI growth rate and CN convergence speed of scenarios with a LEO satellite are still significantly better than those of the DRO–DRO scenario. The observability of the DRO–DRO scenario exhibits the greatest sensitivity to lunar estimation: the initial CN worsens to the 1×1029 level (three orders of magnitude higher than without lunar ephemeris estimation), whereas the initial CN of the LEO–DRO scenario increases to 1×1026 (still three orders below that of DRO–DRO). The DRO–DRO–LEO scenario shows optimal performance, with its initial CN only increasing to 1×1024 (only a two-order increase), highlighting the capacity of the multi-satellite configuration to mitigate poor conditioning from parameter dimensionality expansion.

These findings demonstrate that lunar ephemeris estimation reduces observability relative to the estimation-free case; however, based on the specific numerical performance indicated by the CN and OI, sufficient observability is maintained during lunar ephemeris estimation. Moreover, by appropriately extending the observation period and utilizing the dynamic characteristics of LEO satellites as well as the geometric diversity of multi-satellite observations, the system's observability can be further improved, thereby enhancing the accuracy and reliability of satellite autonomous navigation and lunar ephemeris estimation.

4.3 Simulation Results

Based on the observability analysis, this section aims to perform numerical simulation experiments to validate the performance of the three observation scenarios, as shown in Figure 1. The parameters to be estimated are the states of the satellites and the Moon, which are given in Equation (1). The initial position and velocity uncertainty of the satellite are set to 100 m and 0.01 m/s, respectively, reflecting the achievable autonomous navigation accuracy (Turan et al., 2022b). The initial states of the Moon are derived from the DE430 ephemeris, with positional discrepancies set to 300 m and 1×10−4m/s compared with contemporary ephemerides. The specific orbit determination parameters are given in Table 3.

Table 4 presents the results of 20 Monte Carlo simulation experiments for the three observation scenarios. The table quantifies the navigation accuracy through two key metrics over the final 10-d period: (1) 3σ uncertainty bounds (3σ Unc.) representing the formal error distribution and (2) actual errors (Error) computed as the root mean square (RMS) of position/velocity residuals. The convergence time is defined as the time required for the position state to converge within 50 m. Notably, in the DRO–DRO scenario, none of the estimated satellites achieved this convergence threshold, as indicated by dash markers (-).

According to Table 4, this method benefits from the DRO state's sensitivity to lunar ephemeris. By aligning the SST measurements with the dynamic model, one can estimate the lunar ephemeris with a certain level of accuracy. Moreover, Table 4 presents distinct navigation performance among the three scenarios, consistent with the prior observability analysis results.

Specifically, implementing a LEO–DRO SST link enables an orbit determination accuracy of 0.28 m (2.66 m at 3σ uncertainty) and a velocity determination accuracy of 3.05 × 10^–4 m/s (2.55 × 10^–3 m/s at 3σ) for LEO satellites. Benefiting from the short orbital period of the LEO satellite, the position residuals converge within 1 h. However, the deep-space navigation capability remains constrained. DRO satellites and lunar ephemeris estimation exhibit convergence times of 8.31 d and 8.91 d, respectively, under this scenario, with position accuracies of 27.18m(3σ:266.15m) for DRO satellites and 27.83m(3σ:220.71m) for lunar ephemeris.

When an additional DRO satellite is incorporated into the LEO–DRO scenario, significant enhancements in satellite orbit determination accuracy are observed. Whereas the LEO accuracy improves to a position error of 0.06 m, the DRO satellite position accuracy improves to 5.66m(3σ:33.74m), with the convergence time reduced to 4.67 d. The lunar ephemeris accuracy simultaneously increases to 4.04m(3σ:22.72m), with 54.61% accelerated convergence to 4.04 d.

The DRO–DRO scenario demonstrates fundamental limitations, with position errors of 225.37m(3σ:518.29m) for DRO-A and 79.30m(3σ:200.68m) for DRO-B and a lunar ephemeris error of 87.33m(3σ:203.54m)—failing to meet the 50-m convergence threshold.

As the optimal architecture is identified as the DRO–DRO–LEO scenario navigation, an analysis of the simulation results from the DRO–DRO–LEO scenario is presented below.

The measurement thermal noise for SST links is configured at 0.5 m across all scenarios. Simulation results demonstrate that the O-C residuals conform to zero-mean Gaussian distributions, with standard deviations of 0.499 m for LEO–DRO-A, 0.497 m for LEO–DRO-B, and 0.500 m for DRO–DRO links—values remarkably consistent with the predefined thermal noise level.

Furthermore, Kolmogorov-Smirnov tests (p-value > 0.75) rigorously confirm the Gaussian distribution characteristics, indicating that the selection of process noise in the filtering algorithm is well justified and appropriately models the system uncertainties.

Figure 10 displays the results of the analysis, presenting time series data for the three-dimensional (3D) position residuals of the LEO (red), DRO-A (blue), and DRO-B (green) satellites. The position accuracy of each satellite is determined by computing the difference between the true trajectory and the estimated trajectory.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 10

Time series of positioning residuals of the LEO, DRO-A, and DRO-B satellites in the DRO–DRO–LEO scenario

Based on the satellites' residual convergence in the DRO–DRO–LEO scenario, shown in Figure 10, the orbital determination process exhibits distinct convergence characteristics across different satellites. The LEO satellite exhibits the highest convergence rate, ultimately achieving the highest accuracy levels, reaching sub-meter accuracy within 5 d. The accelerated convergence for this satellite originates from its LEO orbital dynamics, where frequent geometric variations relative to background satellites generate diverse measurement geometries that rapidly constrain orbital parameters. In contrast, the residuals of the DRO satellite converge slowly, requiring approximately 10 d to stabilize at 10 m. This intermediate convergence rate reflects fundamental astrodynamic constraints of DROs, where slower geometric evolution relative to a LEO necessitates extended observation arcs to achieve commensurate precision.

Figure 11 presents the residuals of the estimated lunar states in the three orthogonal directions. The lunar position accuracy is evaluated by comparing the true lunar trajectory (from DE430) with the estimated trajectory.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 11

Time series of lunar positioning residuals in different directions: (a) positioning residuals; (b) local magnification of (a)

Throughout the extended 90-d estimation window, this methodology ultimately achieves a lunar ephemeris accuracy of 4.0 m (3D RMS) under the current configuration, which can satisfy the requirements for the dynamic model accuracy of satellites during autonomous orbit determination. Convergence is rapid in the radial and tangential directions, reaching 40 m within 2 d. The normal direction converges more slowly but still achieves the 40-m threshold within 15 d.

Furthermore, Figure 11 reveals a notable pattern in the estimated accuracy of the lunar ephemeris, characterized by periodic oscillations occurring approximately every 13.5 d, closely aligned with the orbital period of the DRO satellite. To investigate the underlying cause of this residual oscillation phenomenon, an analysis was performed to explore the relationship between the 3D position residuals of the lunar ephemeris and the orbital position of the DRO satellite. To quantify the correlation between residual behavior and orbital geometry, we isolated high-frequency trends through low-pass filtering of estimated lunar position residuals. Then, based on the filtered residuals, we classified temporal segments into convergence phases (decreasing residuals, shown in gray in Figure 12(a)) and divergence phases (increasing residuals, black). The mapping of these time series segments onto the DRO phase, presented in Figure 12(b), reveals a strong correlation between the convergence and divergence of lunar position determination residuals and the orbital phase of the DRO. When the DRO is located between the Earth and the Moon, the position residuals of lunar ephemeris are predominantly in a divergent state (black trajectory). Conversely, when the DRO is on the far side of the Moon, the residuals tend to converge (gray trajectory).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 12

Residuals in lunar position vs. DRO satellite orbit phase: (a) 3D position residuals of the Moon; (b) DRO satellite orbit phase in the Earth–Moon rotating frame (DU: Earth–Moon distance)

Further analysis reveals a significant correlation between the state of lunar ephemeris position residuals and the DRO phase during residual convergence segments, as shown in Figure 13. When the DRO is positioned between the Earth and the Moon (near side), residuals are predominantly divergent (65.4% occurrence). Conversely, during far-side DRO phases, residuals show dominant convergence, corresponding to approximately 82.3% of the time.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 13

Probability of lunar position residual variation in different DRO phases

4.4 Discussion

The current section delves into an analysis of the factors that contribute to the superior performance of the DRO–LEO scenario compared with the DRO–DRO scenario in estimating lunar ephemeris. As shown in Table 4, when two DRO satellites establish an SST link to estimate lunar ephemeris, both the trajectory accuracy and the lunar ephemeris estimation accuracy are significantly inferior to those obtained in the DRO–LEO scenario. Two primary factors contribute to this phenomenon.

Firstly, this phenomenon is attributed to the difference in the configuration of inter-satellite ranging. Figure 14 compares simulated SST ranging variation in the DRO–DRO scenario and the DRO–LEO scenario. In Figure 14(a), the range of the DRO–DRO scenario exhibits periodic fluctuations with an approximate period equivalent to the lunar orbit period of the DRO satellite (approximately 13.5 d). This slow rate of ranging variation prolongs the convergence time. Figure 14(b) illustrates the range sequence derived from satellite–satellite distance measurements in the DRO–LEO scenario. The short orbital period of the LEO satellite contributes to the short-term characteristic observed in the obtained ranges when establishing an SST link. This characteristic enables rapid convergence of the satellite's orbital state. Additionally, the dynamic characteristic of the geometric configuration of the satellite–satellite link enhances the system's robustness against model errors in the dynamic model (Wang et al., 2019). Consequently, more accurate orbit determination results are achieved.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 14

Ranging data between two satellites: (a) DRO–DRO scenario; (b) DRO–LEO scenario

Because of the significant Earth-Moon distance, current methods for estimating lunar ephemeris by utilizing LLR data impose rigorous requirements on the measurement techniques and hardware equipment of ground laser stations. Consequently, the global availability of laser stations capable of conducting lunar ranging measurements is limited. Presently, only a few prominent stations possess the capacity for a large volume of precise measurements, such as the McDonald Observatory, Observatoire de la Côte d'Azur, Apache Point Observatory, and Lunar Laser Ranging Experiment Observatory (Courde et al., 2017). The Wettzell Laser Ranging System and Matera Laser Ranging Observatory have observational capabilities, but these stations are constrained in acquiring a substantial volume of observational data. In contrast to the current lunar ephemeris estimation method, this approach is independent of ground-based or lunar-based observation instruments and does not necessitate extensive decades-long observational data. Furthermore, owing to the comparatively shorter duration of satellite orbits, SST measurements change more rapidly when compared with LLR data. Therefore, this method enables the accumulation of a substantial volume of measurement data within a relatively short time. Such an accumulation plays a pivotal role in enhancing the convergence speed of the satellite orbit and the constraints on lunar ephemeris while avoiding divergence in lunar ephemeris. However, the estimation accuracy of this method is limited, and this method still cannot replace existing lunar ephemeris estimation methods.

Another factor is the advantage of the DRO–LEO inter-satellite connection in measuring the Earth–Moon link, specifically for estimating lunar ephemeris and orbit states. As depicted in Figure 1, the SST link established between the DRO and LEO satellites simultaneously connects the Earth and the Moon. Hence, this link scenario is more sensitive for estimating lunar ephemeris than other scenarios. In contrast, the DRO–DRO scenario lacks the LEO satellite as an anchor point connecting the Earth, resulting in an overall drift of the link between the DRO satellite and the Moon. Consequently, when determining the orbit of the DRO satellite and estimating lunar ephemeris, the superior geometric configuration of the DRO–LEO inter-satellite ranging link plays a crucial role in obtaining more accurate estimation results in this scenario.