## Summary

A comprehensive study on the autonomous orbit determination (AOD) performance of satellite pairs in halo orbits and distant retrograde orbits (DROs) is carried out. A factor called *dynamic and geometric dilution of precision* (DAGDOP) is proposed to simultaneously incorporate influences from the dynamics and geometry of satellite pairs. Based on the DAGDOP, the effect of different observation arcs on the AOD accuracy is investigated. Next, the AOD accuracy of three different types of satellite pairs—halo+halo, DRO+DRO, and halo+DRO—is systematically analyzed. The hybrid halo+DRO type shows the best overall accuracy. Finally, the AOD performance of the hybrid type is verified in a realistic model. Our studies find that the average AOD accuracy of the halo orbit is about 170 meters, and that of the DRO is about 190 meters. The relative time synchronization error of two satellites is less than 30 nanoseconds.

## 1 INTRODUCTION

Traditionally, deep space navigation for lunar satellites (e.g., Apollo [Wollenhaupt, 1970], Lunar Prospector [Beckman & Concha, 1998], Lunar Reconnaissance Orbiter [LRO; Mazarico et al., 2018], and China’s Chang’E [Duan et al., 2019]) strongly depends on radiometric tracking via ground stations, such as the Deep Space Network (DSN). For an unmanned lunar probe, the traditional tracking technique suffices. But for future crewed lunar exploration missions that are more demanding, such as NASA’s Artemis program (Smith et al., 2020), it may no longer be accurate enough due to the lack of continuous tracking. Tracking signals can be eclipsed by the Moon and Earth (if the number of ground stations is not enough).

Autonomous navigation that can be achieved on board is essential in this case. To strengthen the capability of autonomy, many autonomous navigation methods have been studied from the very beginning of the space era (Battin, 1962; Bowers Jr., 1966; Chory et al., 1984). For example, optical sensors and celestial navigation were the most common methods to obtain self-contained navigation results by tracking landmarks (Hur-Diaz et al., 2008) or celestial bodies (Christian & Lightsey, 2009). What’s more, GNSS constellations show a great capacity for the autonomous navigation of low Earth orbits (LEO) since real-time positions can be obtained immediately with a GNSS receiver.

For orbits with altitudes higher than those of GNSS, GNSS constellations also show the potential for navigation using signals from the other side of the Earth. Some experiments have been carried out in high Earth orbits (HEO; Balbach et al., 1996), geostationary Earth orbits (GEO; Kronman, 2000), and highly elliptic orbits (Winternitz et al., 2017). In the Chinese CE-5T1 mission, experiments were also conducted and the feasibility of GNSS was verified in lunar transfer orbits (Liu et al., 2017). A weak GPS receiver was developed at NASA’s Goddard Space Flight Center (GSFC; Bamford et al., 2008). Capuano et al. (2017) proposed the use of an adaptive orbital filter to aid GNSS signal acquisition and increase navigation accuracy. However, as far as the authors know, no such experiments have been conducted in lunar orbits.

Apart from the use of GNSS, a novel autonomous navigation technique known as Linked Autonomous Interplanetary Satellite Orbit Navigation (LiAISON) was proposed by Hill et al. (2005). The absolute states of the satellites can be solved on board simultaneously with the use of *satellite-to-satellite tracking* (SST), which would significantly reduce the burden of ground stations.

The utilization of SST data was originally conceived to help with the orbit determination of satellites around the Earth (Vonbun et al., 1978). However, *autonomous orbit determination* (AOD) using SST data solely cannot work in the framework of the two-body problem. Absolute states of the satellites cannot be obtained due to a rank deficiency problem. It was proven by Liu and Liu (2001) that this problem could be solved in two ways. One way is to add a ground station and combine the SST data with the ground station tracking data. Another one is to provide a-priori information about the orbit plane orientation, i.e.,(*i*, Ω) of the orbits. Namely, the (*i*, Ω) have to be fixed to determine the other elements. Furthermore, Hill and Born (2007) found that the cause of this rank deficiency problem was that the Ω_{1} and Ω_{2} only appeared in the form of ΔΩ in the two-body problem.

Tang et al. (2021) has demonstrated this idea with the use of the centralized AOD strategy by artificially fixing the Ω of at least one satellite. Further study on the observability of different satellite configurations was conducted by Qin et al. (2019). Nonetheless, cross-links have been used in GPS Block IIR satellites. However, due to the problem discussed above, GPS Block IIR satellites have to use stored navigation messages computed by ground stations to achieve autonomous navigation (Abusali et al., 1998). Significant enhancements to position and timing performance can be provided with the use of anchor stations (Rajan et al., 2003).

According to the LiAISON strategy, in the three-body problem if the satellites are deployed on some special orbits, such as halo orbits around collinear libration points (CLPs), the absolute states of the satellites will be observable. This is because of the strong perturbation of the third body to the satellites. An illustration of the strength of asymmetry in the Earth–Moon system was given (Hill & Born, 2007) by defining a parameter *α _{j}*.

A total number of *n* sources of acceleration act on the spacecraft, including the particle gravity of the Earth and the Moon. The definition of third-body depends on the satellite’s position. When the gravitational force due to the Moon is stronger, then Earth is taken as the third body. If the gravitational force due to the Earth is stronger, then the Moon is taken as the third body. The 2D and 3D maps of *α*, along with halo orbits and distant retrograde orbits (DROs) in the Earth–Moon system are shown in Figure 1. Frames (a) and (b) show the 2D contour map and frames, while (c) and (d) show the 3D contour map of *α*. Numbers on the level curves of frames (a), (b), and numbers on the nearly spherical shells in frames (c) and (d) indicate the values of *α*. It can be found that the halo orbits and DROs are mostly distributed in the area where *α* is between 30% and 40%. It is obvious that halo orbits and DROs are in the regions where the third-body perturbation is very strong, which is favorable for the AOD.

Considering the asymmetry of the dynamics, halo orbits around CLPs and DROs are good candidates and are the most used to build a conceptual navigation constellation. Plenty of studies on halo orbits have been conducted (Du et al., 2013; Gao et al., 2014; Hamera et al., 2008; Hill et al., 2006; Zhang & Xu, 2016). For lunar surface navigation, Hesar et al. (2015b) studied the feasibility of the LiAISON strategy to navigate rovers on the lunar farside surface with one satellite on a halo orbit around the Earth–Moon L2 point. According to their study, a precision of tens of meters can be achieved over the majority of the lunar farside surface.

The LiAISON strategy was also able to be extended to an Earth–Moon constellation which contains a satellite in a GEO and another in a halo orbit around Earth–Moon L1 or L2 point (Fujimoto et al., 2012; Leonard et al., 2012; Liu & Hou, 2014; Parker et al., 2012). Navigation of a crewed vehicle between the Earth and the Moon using both ground tracking and satellite-to-satellite tracking was analyzed by Leonard et al. (2013). What’s more, the LiAISON method was also employed to support interplanetary exploration (McGranaghan et al., 2013).

Apart from halo orbits, DROs can also serve as orbit candidates to achieve AOD, and they show better stability which is of vital importance for navigation constellations. Liu et al. (2014) first proposed the use of DROs to replace halo orbits in the LiAISON strategy. They investigated the AOD accuracy of a lunar satellite+a DRO configuration and compared the accuracy with that of a lunar satellite+a halo orbit configuration.

Wang et al. (2019) further exploited the idea and studied the navigation performance between the DROs and different types of cislunar orbits considering dynamic and clock model errors. Apart from navigation of satellites, Hesar et al. (2015a) presented a method to estimate the gravity field of a small body based on LiAISON strategy. Furthermore, the LiAISON strategy has progressed into actual engineering testing. Advanced Space has partnered with NASA to develop the Cislunar Autonomous Positioning System (CAPS) based on this strategy to provide navigation solutions for missions throughout cislunar space (Cheetham, 2017).

Previous studies focus on a scenario in which the AOD is based on the SST data between a navigation satellite moving in a halo orbit (or a DRO) and a user satellite. This may limit the use because the state of the navigation satellite is still unknown and it is determined together with the state of the user satellite. This means that the AOD results have to be transmitted back to the navigation satellite, which is usually not practical.

A more practical scenario is one in which a navigation constellation moves in the proximity of the Moon, and the states of the constellation have been determined in advance. The whole AOD process is self-contained. That is, the two navigation satellites would track each other and use SST data to autonomously determine their orbits, with no support from ground stations or GNSS. Then, the navigation message would be broadcast from the navigation satellites to the user satellites.

Considering the cost of building such an autonomous lunar navigation constellation, a minimum configuration of two satellites is a good starting point. Also, considering the service volume of such a constellation, the two satellites should move on orbits with high enough altitudes. Based on these two considerations, we chose the well-studied DROs and halo orbits as candidate orbits for the navigation satellite pairs. In such a case, we think it is very important to figure out the AOD performance of different configurations based on these orbits, which can serve as an impact factor when designing such a two-satellite navigation constellation. We have to admit that better AOD results can be expected if more satellites are deployed in the constellation, but as a preliminary study, we restrict our analysis to the two satellites in the current work. In the following, we call this two-satellite navigation constellation *satellite pairs*.

Even for the simplest satellite pairs, a complete AOD analysis is not trivial. The analysis can be deconstructed into two sub-problems. The first sub-problem is with the AOD analysis of the satellite pairs themselves. This is to be investigated in this work. The second sub-problem is that the AOD analysis of the user satellites uses the navigation message. In the second sub-problem, accuracy of the two satellites’ position and velocity obtained in the first sub-problem influences the AOD accuracy of the user satellites. Moreover, combinations of different user satellites and different satellite pair configurations also influence the AOD accuracy of the user satellites. The second sub-problem will be the focus of our forthcoming paper. We remark that other factors that are required to be considered when designing lunar navigation constellations include coverage ability (Gao & Hou, 2020) and maintenance (Guzzetti et al., 2017), which are beyond the scope of the current study.

Figure 2 demonstrates the idea of the AOD process for the satellite pairs discussed in this paper. Three types of satellite pairs are considered (i.e., halo+halo, DRO+DRO, and halo+DRO). The halo+halo type contains two kinds of configuration (i.e., L2 halo+L2 halo and L1 halo+L2 halo). Thus, in total, four cases of different configurations are studied. Case (a) consists of one L2 south halo orbit and one L2 north halo orbit; Case (b) consists of one L1 south halo orbit and one L2 north halo orbit; Case (c) consists of two spatial DROs; and Case (d) consists of one halo and one DRO. Of course, the satellite pair based on a halo orbit (or a DRO) and a *low-altitude lunar orbit* (LLO) is also feasible, but an LLO usually has limited navigation service capability. Also, the AOD of this type of satellite pair is exactly the same as the case of a navigation satellite and a user lunar satellite, so it will not be repeated in this work.

In frames (a), (b), and (d), the blue orbit represents a *nearly rectilinear halo orbit* (NRHO) originating from the Earth–Moon L2 point. The NRHO is a subset of the halo orbit family (Howell & Breakwell, 1984). Researchers are interested in the NRHO because of its relatively weaker instability when compared to usual halo orbits that are far away from the Moon. The red orbits in frames (c) and (d) around the Moon are the DROs, which are practically stable. There is one satellite on each orbit serving as a navigation satellite. With the SST technique, the two navigation satellites track each other, represented by a yellow dashed line. Through the LiAISON strategy previously mentioned, the absolute states of the two satellites are obtained. That is, the whole constellation is anchored. Based on this, navigation service could be provided to lunar probes by measuring the distance between them. Navigation message from the two navigation satellites to a lunar probe is illustrated in Figure 2, represented by white dashed lines.

We propose a novel factor, coined with the term *dynamic and geometric dilution of precision* (DAGDOP). It captures both the dynamic characteristics and the observation geometry of the constellation, and can be reduced to the well-known factor of *geometric dilution of precision* (GDOP) if no dynamics are involved. Our studies show that this factor agrees quite well with AOD accuracy. As a result, it can be used as an indicator of AOD accuracy without actually performing the AOD process. It is more time-efficient using this factor. This factor is especially useful for satellite constellation design when there are massive candidate orbits. We recommend its use for similar future studies.

Based on the proposed DAGDOP factor, a systematic analysis was conducted in the simple *circular restricted three-body problem* (CRTBP). This systematic analysis was aimed at analyzing the selection effects of different orbital arcs or orbits on AOD performance. To remove the influence of the measurement process, an instantaneous measurement model was adopted. The analysis was divided into two sub-problems. The first sub-problem regards the effect of different observation arcs on AOD performance. Taking the halo+DRO configuration as an example, it was found that AOD accuracy strongly depends on different observation arcs when the data length cannot cover one complete orbit whose orbital period is longer. By the use of longer arcs of data, this phenomenon is no longer obvious, and AOD accuracy can reach an accuracy level of SST measurements when the length of data covers more than one complete orbit with a longer period. For example, in the halo+DRO constellation studied in this paper, an arc length of 15 days is sufficient to achieve the accuracy of the observation data.

The second sub-problem is about the effect of different orbits on AOD performance. Three types of constellations mentioned above are investigated. Among them, the hybrid halo+DRO configuration has the smallest DAGDOP. This result confirms the advantage of hybrid halo+DRO constellations in the aspect of orbit geometry and dynamics, and it reveals the potential of combinations of halo orbits and DROs for future navigation constellations in cislunar space. Furthermore, smaller DAGDOP values can be obtained when a large-amplitude halo orbit and a small-amplitude DRO are combined. On the other hand, for the L2 halo+L2 halo and DRO+DRO configuration, it was found that a rank-deficiency problem occurs when the initial phases of two satellites are too close to each other. This situation should be avoided in actual mission designs. It is also recommended to avoid simultaneously using orbits that are symmetrical to each other with respect to the x-y plane.

Finally, take the best performing halo+DRO constellation in the simplified CRTBP model. As an example, in the second half of this study, a simulation of the real force model was carried out to check the actual AOD accuracy of the satellite pairs. The simulation results showed that the AOD process of the hybrid halo+DRO constellation works well in reality. The average AOD accuracy of the halo orbit was about 170 meters, and the average AOD accuracy of the DRO was about 190 meters. The clock synchronization accuracy can achieve an accuracy of no more than 30 nanoseconds.

The contributions of the current paper are threefold. (1) We propose a novel factor DAGDOP to describe the AOD performance between satellite pairs. It is an efficient factor to compare the AOD performance of different orbit pairs, and is recommended in future studies. This factor can be extended to multi-satellite constellations with SST measurements. It can also be applied to other types of observation data. (2) A comprehensive study of the AOD performance of halo orbits and DROs is conducted in the CRTBP. Halo orbits and DROs with different amplitudes and different configurations are compared. Such a systematic survey, to the authors’ knowledge, is seldom seen in the literature. (3) A hybrid configuration composed of a halo orbit and a DRO is proposed for the first time, and is verified in a simulated real model that includes a high-fidelity force model and clock model. Through Monte-Carlo simulations, the AOD accuracy of the satellite pair is verified. These findings can be regarded as a practical reference for the lunar navigation constellation design in the future.

This paper is organized as follows. In Section 2, the basics of two kinds of dynamic models used in this paper are introduced. One is the CRTBP, and the other is the real force model, taking into account a more realistic dynamical environment. In Section 3, the measurement system is described, including the time delay, clock model, etc. In Section 4, the principles of the AOD process are introduced. In Section 5, the factor DAGDOP is proposed and a systematic analysis of the AOD problem is carried out in the CRTBP based on the DAGDOP. In Section 6, the AOD performance of the hybrid halo+DRO configuration in the real force model is verified. Section 7 summarizes the study as a whole.

## 2 DYNAMIC MODEL

In this paper, two force models are used. First, the simple CRTBP model is used for the systematic analysis of different observations arcs and different types of satellite pairs. Then, taking the best performing halo+DRO satellite pair as an example, the AOD performance is verified in the real force model. In the following, as basics to our study, the two force models are briefly introduced.

### 2.1 CRTBP

In the synodic frame centered at the barycenter, Figure 3 depicts the geometry of the CRTBP model of the Earth–Moon system. *P*_{1} is the Earth center and *P*_{2} is the Moon center. *μ* = *m*_{2}/(*m*_{1} + *m*_{2}) is the mass parameter, with *m*_{1} and *m*_{2} as masses of the Earth and the Moon, respectively. Following mass, length and time units are used throughout our study:

where *P*_{1}*P*_{2} means the distance between *P*_{1} and *P*_{2}.

Equations of motion are as follows:

1

2

3

where:

### 2.2 Realistic Model

For the real Earth–Moon system, besides the particle gravity of the Earth and the Moon, there are many perturbations that need to be considered when studying the orbit determination problem. For libration point orbits, the most important perturbations are the real motion of the Moon with respect to the Earth and the solar gravity perturbation. The relatively minor perturbations include non-spherical terms of the Earth and the Moon, third-body perturbations from other planets besides the Sun, and solar radiation pressure. Some even weaker perturbations also exist, such as relativistic corrections, tidal effects of the Earth and the Moon, etc. It is impossible to incorporate all the perturbations into our force model and there are uncertainties for parameters of these perturbation forces. In our work, the force model adopted to approximate the real Earth–Moon system include:

*F*_{1}: particle gravity of the Earth, Moon, and Sun;*F*_{2}: non-spherical gravity of the Earth and the Moon, including the tidal effects of the Moon; and*F*_{3}: solar radiation pressure.

We point out that although the DROs and halo orbits are dynamic structures in the synodic frame of the Earth–Moon system, the AOD process is carried out in the Moon-centered *international celestial reference frame* (ICRS). That is, equations of motion are integrated in this frame and the SST data is calculated using the difference between the positions of the satellite pairs. In the Moon-centered ICRS, equations of motion follow:

4

where *μ _{M}*,

*μ*, and

_{E}*μ*represent the gravitational constants of the Moon, the Earth, and the Sun, respectively.

_{S}**,**

*r**, and*

**r**_{E}*are position vectors of the satellite, the Earth, and the Sun in the Moon-centered ICRS.*

**r**_{S}**Δ**

_{E}(

**Δ**

_{s}) is the position vector that points from the Earth (Sun) to the satellite. ∇

*U*and ∇

_{E}*U*represent the non-spherical gravity of the Earth and the Moon. [TE] and [TM] are two rotation matrices that convert vectors from the Earth’s body-fixed frame and the Moon’s body-fixed frame to the ICRS, respectively. Non-spherical gravity fields are modeled in the form of:

_{M}5

where *μ* is the gravitational parameter of the center body, (*R, λ, ϕ*) are the spherical coordinates of the satellite in the corresponding body-fixed frame, and *a* is the reference radius of the center body. is the associated Legendre function of degree *l* and order and are normalized Stokes coefficients. To make the simulations more realistic, tidal effects of the lunar geopotential are also considered. Changes caused by lunar tides are usually expressed as variations and (Konopliv et al., 2013). They are expressed as:

6

where *k _{l,m}* are the lunar love numbers,

*j*represents the disturbing body (the Earth or Sun),

*r*is the distance between the disturbing body and the Moon, and

_{j}*ϕ*and

_{j}*λ*are the latitude and longitude of the disturbing body in the Moon’s body-fixed frame. The gravity field is truncated at lower degrees and order (see Table 1) in the AOD process, and only the terms due to lunar tides is calculated.

_{j}Equation of * F_{SRP}* is:

7

where *P _{SR}* represents the solar radiation pressure at 1 astronomical unit (AU).

*C*represents the reflectivity coefficient, which is a value usually between zero and two. When we simulate SST data, we assume that the radiation term consists of a constant part plus a time-varying periodic part. The constant part is

_{R}*C*

_{R0}= 1.5. The period of the time-varying part is assumed to be one month, which is reasonable since the satellites are moving in the vicinity of the Moon. Uncertainty is included in the amplitude of the time-varying part, which is assumed as a Gaussian distribution with 1

*σ*= 10% and a mean of zero. In the AOD process, the time-varying periodic part is neglected and the constant part is corrected along with the state vectors.

*S*is the area facing the Sun.

*m*is the mass of the satellite.

**Δ**

_{S}remains the same as the one in Equation (4). The gravity model GGM05C (Ries et al., 2016) for the Earth and GRGM1200A for the Moon are used (Goossens et al., 2016; Lemoine et al., 2014). Values of the parameters of the force model are summarized in Table 1.

Due to the many perturbations in the real Earth–Moon system, the CLPs lose their meaning as dynamical equilibrium points in the synodic frame (Gómez et al., 2001; Hou & Liu, 2010). Nevertheless, dynamics around these points are proven to be similar to those of the CRTBP (Hou & Liu, 2011; Jorba et al., 2020). Numerical approaches can be taken to construct the DROs and the halo orbits (including the NRHOs) in the ephemeris model (Bezrouk & Parker, 2017; Dei Tos & Topputo, 2017; Hou & Liu, 2011; Lian et al., 2013; Qian et al., 2018; Whitley et al., 2018). Again, details are omitted here. Readers can refer to the references above for more details on how to compute these orbits in a more realistic Earth–Moon system.

## 3 MEASUREMENT MODEL

In the simple CRTBP model, it is assumed that the ranging measurements are obtained instantly. The measurement noise is modeled as Gaussian noise with a standard deviation of 1 m. We have to emphasize that the main purpose of the systematic analysis in the simple CRTBP model is to evaluate the general AOD characteristics of different constellations from the viewpoint of pure geometry and dynamics, so only this simple instantaneous measurement model is considered.

As for the simulation in the real force model, it is necessary to establish a more realistic measurement model, taking the time delay, clock errors, and other aspects into consideration. In order to improve the accuracy of the orbit determination of the satellite pair, a dual one-way ranging method is adopted. Suppose that a ranging signal is transmitted by Satellite 1 at time *t* − *τ* and received by Satellite 2 at time *t*. The general form of ranging measurement can be written as:

8

where and represent the position vector of Satellite 1 at time *t*−*τ* and the position vector of Satellite 2 at time *t*, respectively. *τ* is traveling time of signal, which is generated in an iterative way and described below. *δt*_{1} and *δt*_{2} represent the clock errors of Satellite 1 and Satellite 2, including the deterministic error and the stochastic error. The last term *ε* is the thermal noise, which is modeled as a Gaussian noise with its 1*σ* varying with the distance between the two satellites.

### 3.1 Time Delay

One of the most important effects that should be modeled is the time delay. Suppose Satellite 1 transmits a ranging signal at time *t ^{T}*. At time

*t*, Satellite 2 receives that signal.

^{R}*τ = t*is the time delay in this ranging process. A simple iteration process is used to obtain

^{R}− t^{T}*τ,*which is shown in Figure 4.

*t*is assumed to equal

^{R}*t*as an initial guess. Then,

^{T}*R*

_{2}(

*t*) can be solved using one of two methods, orbit propagation or ephemeris interpolation. The second method is more appropriate for the onboard AOD process. In this work, the first method is used. Then, the expected range

^{R}*ρ*is calculated and

*τ*is obtained. Then, a new

*t*is obtained. The iteration is finished until the difference between

^{R}*τ*and

_{i}*τ*

_{i−1}is less than

*ϵ*, which is set as 10

^{−8}seconds in this work.

### 3.2 Clock Model

One key part of a navigation satellite is the accurate atomic clock. The first satellite navigation system was TRANSIT, developed by the John Hopkins Applied Physics Laboratory. Quartz crystal oscillators were used to generate accurate time reference. In 1974, rubidium atomic clocks were carried onboard for the first time (Bhaskar et al., 1996). Deterministic and stochastic errors were involved in atomic clocks onboard, denoted as *δt _{d}* and

*δt*, respectively. Usually,

_{s}*δt*can be modeled as the following second-order polynomial:

_{d}9

where *a*_{0}, *a*_{1}, and *a*_{2} represent the clock bias, clock drift, and clock aging parameter, respectively. Only the relative difference of clock bias, drift, and aging coefficients between two clocks can be determined, which will be explained in detail in Section 4. The relative difference of clock bias, drift, and aging coefficients between two clocks onboard are:

10

where the subscripts represent clock bias, drift, and aging parameters, respectively. The superscripts indicate the satellite number. The relative difference of clock bias, drift, and aging coefficients between two clocks onboard are summarized in Table 2.

The Allan variance is usually used to describe the frequency stability of atomic clocks. Rubidium clocks on the GPS Block IIF satellites offer an Allan variance of (Vannicola et al., 2010):

11

Atomic clocks of the same stability are adopted in our work. Stochastic errors can be then generated from the Allan variance above. The conversion method is the same as Wang et al. (2019) and will not be repeated here. A time history sample of stochastic error is shown in Figure 5.

### 3.3 Thermal Noise

Thermal noise of the ranging system is related to the received signal power, which is dependent on the distance between the two satellites. The thermal noise is modeled as a Gaussian noise in this paper with 1*σ* varying with the distance. The *σ* can be approximated by:

12

where *λ _{c}* is the wavelength of each code chip (293.05 m);

*B*typically equals 0.5 Hz, representing the code loop noise bandwidth; D is the early-to-late correlator spacing (chips), with a typical value of 1 chip;

_{n}*T*is the prediction integration time; and

*T*= 1

*s*is used in this work.

*C*/

*N*

_{0}is the carrier-to-noise ratio (Hz):

13

where (*N*_{0})_{dB} is the thermal noise power component in a 1-Hz bandwidth (dBW) with a typical value of −200.9 dBW (Kaplan & Hegarty, 2005).

A free-space propagation loss model given by Kaplan and Hegarty (2005) is expressed as:

14

where (*P _{R}*)

_{dBW}and (

*P*)

_{T}_{dBW}are the receiving and the transmitting power. (

*G*)

_{T}_{dB}and (

*G*)

_{R}_{dB}are the gain of receiving and the transmitting antennas.

*λ*=

*c*/

*f*is the wavelength of the signal.

According to the International Telecommunication Union (ITU) Radio Regulations, the 300-MHz to 2-GHz range should be reserved for radio astronomy observations. Therefore, we use a frequency of 2,500 MHz in this work. *d* represents the distance between the transmitting antenna and the receiving antenna. Two antennas with the same parameters are carried on each satellite. The antenna parameters used in this work are summarized in Table 3. The standard deviation of thermal noise is shown in Figure 6, which varies linearly with the distance.

### 3.4 Dual One-Way Ranging

A *dual one-way ranging* (DOWR) method is applied in this work. Two satellites move on their respective orbits, transmitting ranging signals to each other. A time slot of 1 second is allotted to each satellite for transmitting signals. As seen in Figure 7, the first slot in blue is allotted to Satellite 1 and the other colored in red is allotted to Satellite 2. During the transmission time slot of one satellite, the other satellite is in the mode of receiving. Each of the two satellites carries one atomic clock and transmits a ranging signal at its transmitting time slot. Then, the signal is received by the other satellite.

As seen in Figure 7, Satellite 1 transmits a signal at , and then that signal is received by Satellite 2 at . The subscript indicates the satellite number and the superscript indicates the signal mode, *T* for *transmit* and *R* for *receive*. Then, the actual *time of flight* (TOF) is obtained. However, due to inevitable clock errors in the atomic clocks carried on Satellite 1 and Satellite 2, the observed TOF is slightly different from the actual TOF by *δt*_{1} and *δt*_{2}. Thus, the observed TOF is . In the time slot of Satellite 2, the ranging signal is transmitted from Satellite 2 and received by Satellite 1. Finally, dual one-way measurements are obtained:

15

One common method to deal with the DOWR measurements is to decouple state vectors with clock parameters. To achieve this, a common time in Figure 7 is used for the decoupling process. After decoupling, state vectors and clock parameters are solved separately. This is described in Section 4.

## 4 ORBIT DETERMINATION

The AOD process in the CRTBP (see Section 5) is simple. Only the absolute states of the two satellites are solved. This is a classic orbit determination problem and can be found in any textbook (Tapley et al., 2004). Thus, the AOD process in the simple CRTBP is not described here. The AOD algorithm for simulations in Section 6, which focuses on the more realistic model, is described in detail. In Section 6, apart from a more realistic force model, the dual one-way measurements are generated based on onboard atomic clocks. The absolute states and the clock parameters of the two satellites are to be estimated together with the reflectivity coefficient *C _{R}*.

### 4.1 Preprocessing the DOWR Data

Before the AOD process, preprocessing the DOWR data is necessary. Usually, the DOWR data are processed in advance in order to decouple the range measurements and clock parameters. The DOWR measurements directly obtained onboard are measurements at different times. Through orbit propagation and clock parameter prediction, the states of Satellite 1 and Satellite 2 at the common time *t*_{0} are obtained. Then, the corrections to *ρ*_{1}, and *ρ*_{2} are expressed as:

16

Then, DOWR measurements at time *t*_{0} are obtained:

17

Consequently, two virtual DOWR measurements are derived from combinations of and :

18

States and clock errors are now decoupled. Clock errors are eliminated in *ρ*_{+}. States are eliminated in *ρ*_{−}. In the following, the two virtual measurements are processed to achieve the AOD.

### 4.2 The AOD Process

In the CRTBP, only the states of two satellites are to be estimated. Thus, it is basically identical to the classical orbit determination problem. As mentioned above, this is simple and will not be repeated here.

In the simulations in the realistic model, the AOD process is more complicated and is described in detail. In our study, we use the batch filter based on the least squares principle (Tapley et al., 2004).

Since the states and the clock parameters are decoupled, the states of the two satellites and the clock parameters are solved separately. The reflectivity parameter *C _{R}* is estimated together with the states of the two satellites. We denote the integrated vector to be estimated as:

19

where the superscript denotes Satellite 1 or 2.

20

is the state vector of the spacecraft *i*, which is a 6 × 1 vector. Absolute clock parameters cannot be solved. Only relative value of two clocks can be determined. That is to say, the two clocks are synchronized. An absolute clock bias may exist, but it does not affect the SST measurement. The clock vector which contains the three relative clock error parameters is denoted as:

21

Denote the equations of motion in the real Earth–Moon system as:

22

The initial condition of the integrated state vector is denoted as:

23

In our case, the measurement data is the decoupled DOWR data. Supposing at the epoch *t _{i}*, we have:

24

The AOD process is to determine the value of *X*_{0} and ** a** through a series of SST measurements carried out at different epochs

*t*(

_{i}*i*= 1,…

*n*). Supposing we have an initial estimate and

*a*^{*}for

*X*_{0}and

**, we denote:**

*a*25

starting from and integrating the above equations of motion, we get an estimated trajectory *X*^{*}. We denote the residuals between the actual DOWR measurements *ρ _{+}* (

*),*

**X**_{i}*ρ*

_{−}(

**) and the calculated values of the estimated trajectory at the epoch**

*a**t*as:

_{i}26

We denote:

27

According to the batch estimation theory, supposing all measurements have the same weight, the following equation can be used to calculate **x** from **y**_{+} to gradually refine the estimate for the initial state vector:

28

in which the matrix **H** is defined as:

29

The observation matrix for **y**_{+} is defined as:

30

The observation matrix for **y**_{−} is defined as:

31

and **Φ**(*t _{i}*,

*t*

_{0}) is the state transition matrix from the epoch

*t*

_{0}to the observation epoch

*t*.

_{i}## 5 SYSTEMATIC ANALYSIS IN CRTBP

In this section, a factor called *dynamic and geometric dilution of precision* (DAGDOP) is introduced first. Then, based on this factor, a systematic analysis of the AOD accuracy of the satellite pair moving on the halo orbit or the DRO is carried out. The analysis is divided into two sub-problems. One regards the effects of different observation arcs, while the other focuses on different configurations of the satellite pair. In this section, the AOD problem is studied in the CRTBP.

### 5.1 The DAGDOP Factor

In the area of GNSS, a geometry factor, *geometric dilution of precision* (GDOP), is commonly used to describe the amplification of the standard deviation of the measurement errors onto the navigation solution due to the user/satellite relative geometry. Inspired by the GDOP, a similar factor (DAGDOP) is proposed.

Different from the GDOP, which only considers geometry, dynamic aspects are also included in the DAGDOP. This is understandable because it is actually an orbit determination problem. Assume that the integrated state vector of two satellites is denoted as:

32

where the definition of *X*^{(i)} is the same as Equation (20). The state deviation vector is denoted as **x**. Then, **x** is related to the range deviation vector Δ** ρ** by:

33

where:

34

The first term indicates the observation geometry. * ρ* is an m × 1 vector, where m is the total number of observations. For the

*i*-th observation:

35

The second term is the *state transition matrix* (STM), which is a 12 × 12 matrix. The STM reflects the characteristics of dynamics. The dimension of matrix **H** is *m* × 12. Δ*ρ* and **x** are assumed to be Gaussian and zero mean. The covariance of **x** is obtained:

36

Substitution from Equation (33) yields:

37

Assume each measurement is independent of each other. The standard deviation is:

38

Assuming:

Equation (37) becomes:

39

Denote (**H ^{T}H**)

^{−1}as

**D**, which is a 12 × 12 matrix. DAGDOP of position and velocity can be defined by the diagonal components of

**D**:

40

The subscript *i*, *j* of D_{i,j} indicates the row and column number, respectively. The superscript of the DAGDOP represents the satellite number. The subscript of the DAGDOP *(p* or *v*) represents position and velocity, respectively. Position accuracy is of the most concern to us. Thus, we only show computation results of DAGDOP_{p} (hereinafter referred to as DAGDOP, for simplicity) in the following.

The GDOP factor commonly used in the GNSS only considers the relative geometry between the navigation satellite and the receiver. By contrast, the DAGDOP factor considers both the relative geometry and the STM of the satellites. This factor indicates the amplification of the standard deviation of measurement errors onto the orbit determination results due to both the observation geometry and the dynamics.

Some further remarks about the DAGDOP factor are made here. First, in the case that the state of one of the two satellites (say Satellite 2) is known a priori, then we only need to estimate the state of the other satellite (Satellite 1). The state vector ** X** to be estimated is a six-dimensional vector, and the matrix in Equation (35) is then simply changed as a 1 × 6 vector:

41

We can also define the matrix **D** in the same form, but in this case, it is a 6 × 6 matrix, and we can only define this factor for the satellite to be estimated.

Second, in the case that we have *N* (≥ 3) satellites, if all of their states should be estimated using the SST data, we can define the same observation matrix as Equation (35), but in this case, the state vector ** X** is a 6

*N*-dimensional vector and the matrix

**D**is a 6

*N*× 6

*N*matrix. We can define this factor for each of the satellites. If only one of the satellites’ states is to be estimated and the states of all other satellites are known a priori, then the state vector

**is still a six-dimensional vector, and the matrix the same as Equation (41) can be defined, only with more satellites involved. As a result, the DAGDOP factor can be generalized to the case of multiple satellites.**

*X*Third, for other types of observations besides SST, we can define the same factor. The only difference is that the matrix ** ∂ρ/∂X** should be replaced with the corresponding observation matrix.

One final remark is that our following studies indicate that the DAGDOP factor agrees well with AOD accuracy. This means it is a useful indicator of AOD accuracy. Yet to compute this factor, we only need to integrate the orbits once. In other words, we don’t need to go into the iterations of the AOD process. Due to this benefit, we advocate its use in assessing the orbit determination accuracy of different constellations.

### 5.2 Sub-Problem 1—Different Observation Arcs

In this section, the effects of different observation arcs on AOD performance are studied based on the DAGDOP. When the observation arc is shorter than one orbital period, the different initial phases of the orbits will lead to different observation arcs, which will affect AOD performance. In this first sub-problem, we take the halo+DRO configuration as an example to show the results.

In order to intuitively display the different arc segments due to different initial phases, an example DRO and an example halo orbit are shown in Figure 8. The period of the halo orbit is 14 days. The period of the DRO is 6 days. The blue solid curves represent four different arcs of the halo orbit. The red solid curves represent four different arcs of the DRO. The blue stars and the red stars represent the starting points of each arc. The corresponding phase value is marked, with the unit as the orbit period. The phase is defined as *τ _{i}* = (

*t*−

*t*

_{0})/

*p*, where

_{i}*t*−

*t*

_{0}represents the time from the initial point (

*τ*= 0) and

_{i}*p*is the period of the corresponding orbit (We use one to indicate the DRO and two to indicate the halo orbit).

_{i}It can be seen from Figure 8 that the combination of different arcs of the two orbits generates different relative geometry. The DAGDOP factor for different combinations is expected to be quite different accordingly. On the other hand, if the two arcs are long enough to cover the whole of the two orbits, the relative geometry of the two orbits can be covered by the arcs. The DAGDOP factor in this case is expected to change little when different combinations of the arcs are chosen.

To demonstrate this, in our research different lengths (5 days, 10 days, 15 days, and 20 days) of the arcs are chosen. For each length, we chose different values of *τ*_{1} ∈ [0,1] and *τ*_{2} ∈ [0,1] to represent different combinations of the two arcs. A step size of 0.01 was chosen when we surveyed different combinations of *τ*_{1} and *τ*_{2}.

The DAGDOP factor of the halo orbits of different arcs is shown in Figure 9. From left to right and from top to bottom, the arc length is 5, 10, 15, and 20 days. For the halo orbit whose orbit period is about 14 days in this example, some information can be obtained:

For an arc length of 5 days, the DAGDOP factor varies in a range from 10

^{2.5}~ 10^{4}. There is a difference of 1.5 orders of magnitude between the best and the worst result.For an arc length of 10 days, the situation is similar to the case of 5 days. The DAGDOP factor varies in a range from 10

^{0.6}~ 10^{1.6}. The overall accuracy is improved, but still there is a difference of 1 order of magnitude between the best and the worst result.For an arc length of 15 days, the overall DAGDOP factor is about 10

^{0}(light blue region). This means the level of AOD accuracy equals that of the SST measurement. There is a difference of less than 1 order of magnitude between the best and the worst result.For an arc length of 20 days, when compared with the result of 15 days, the overall DAGDOP factor further improves, but the improvement is limited. The difference of the AOD accuracy between the best and the worst result is about 0.6 orders of magnitude.

For the DRO whose orbit period is about 6 days in this example, similar phenomena as Figure 9 can be observed. The difference is that the overall DAGDOP factor of the DRO is smaller than that of the halo orbit for the same time length. This is understandable, because the DRO’s orbit period is shorter and a same time length covers a longer arc of the DRO.

From the results displayed in Figure 9 and Figure 10 and other example tests not displayed, we reach the following conclusions:

For short arc lengths, the accuracy is greatly influenced by the combinations of arcs used for the AOD. Moreover, for a a 14-day halo + a 6-day DRO configuration, when the arc length is less than 3 days, the AOD usually fails. This means that for the AOD algorithm to work, a length of at least 3 days’ SST data is necessary.

For the AOD accuracy to reach the level of the SST measurement accuracy (i.e., the DAGDOP equals one), the time length should be able to cover one complete orbit whose orbit period is longer. Even longer lengths of data can further improve the AOD result, but generally the improvement is limited.

To demonstrate our argument that the DAGDOP factor can be used as an indicator of AOD accuracy, a detailed calculation of the AOD has been carried out. Position *root-mean-square* (RMS) error is used to evaluate the AOD accuracy. Suppose the reference orbit is generated, and the true position vector at epoch *t _{i}* is (

*x*,

_{i}*y*,

_{i}*z*). The position vector at

_{i}*t*obtained from the AOD process is . The

_{i}*x*,

*y*, and

*z*components of the RMS are:

Figure 11 shows the AOD accuracy of the halo orbit for different combinations of arcs the same as Figure 9, and Figure 12 shows the AOD accuracy of the DRO for different combinations of arcs the same as Figure 10. The noise on the SST measurement is modeled as a Gaussian noise with a 1*σ* of 1 m. It is obvious that the patterns of DAGDOP in Figure 9 are consistent with the AOD results in Figure 11, and the patterns of DAGDOP in Figure 10 are consistent with the AOD results in Figure 12. The edges of regions with different colors in the AOD contour map are not as smooth as those in Figure 9 and Figure 10. This is understandable since the noise in observation data and other errors are included in the AOD process.

By and large, the overall magnitude of the DAGDOP factor and the AOD accuracy are consistent with each other.

Consequently, with the use of the DAGDOP factor, the performance of different satellite pairs can be obtained in an efficient and reliable way. In the following analysis, we focus on showing the DAGDOP factor.

### 5.3 Sub-Problem 2—Different Configurations

After studying the effects of different arcs for one specific configuration, different configurations of constellations will be studied in this section. Different configurations of constellations can be formed using different kinds of orbits at different amplitudes.

#### 5.3.1 Configurations

In this sub-problem, three types of constellations will be investigated—the halo+halo type, the DRO+DRO type, and the halo+DRO type. For the halo+halo type, two cases (L2 halo+L2 halo and L1 halo+L2 halo) are studied. As a result, in total, four cases will be studied in this section as summarized in Table 4.

For each configuration, a total of 35 pairs of orbits of different amplitudes are chosen, as seen in Figure 13. The images of orbits in Cases 1, 2, 3, and 4 are shown in Figure 13(a)–(d), respectively. The orbits in are shown in the Moon-centered synodic frame. The symbol * in each figure indicates the starting point of each orbit. The orbits are numbered according to their location in the orbit family. The closer the orbit is to the Moon, the larger the number is. To be clear, the Number 1 and Number 35 orbits have also been marked in Figure 13 with the arrows pointing to the corresponding orbits.

For Case 1, as seen in Figure 13(a), the north halo orbits and the south halo orbits are symmetrical to each other with respect to the x-y plane. With the increase of orbit number, the period of the orbits decreases from 14.8 days to 8 days. The maximum out-of-plane amplitude increases from 14,674 km to 76,957 km and then decreases to 75,001 km. For Case 2, as seen in Figure 13(b), with the increase of the orbit number, the period of L1 halo orbits decreases from 12 days to 8.9 days, and the same north L2 halo orbits are used as in Case 1. As for Case 3, it should be noted that planar DROs are restricted to the Moon’s orbital plane. Thus, it is impossible to determine the out-of-plane component of the DRO if we only use planar DROs. This problem can be solved by adding a certain out-of-plane amplitudes to generate spatial DROs. Quasi-periodic orbits are usually generated in this way (Gao & Hou, 2020; Liu et al., 2014).

As can be seen in Figure 13(c), two families of spatial DROs are generated with an out-of-plane amplitude of 5,000 km. To be specific, a vertical displacement of 5,000 km is added to the planar DROs, generating the blue orbits. A vertical displacement of −5,000 km is added to the same DROs, generating the red orbits. The two families of spatial DROs are symmetrical to one another with respect to the x-y plane. With the decrease of the in-plane amplitude, the orbit number of the spatial DROs gets larger. For the hybrid configuration in Case 4, 35 L2 halo orbits and 35 DROs are shown in Figure 13(d). With the increase of the orbit number, the period of L2 halo orbits decreases from 14.8 days to 8 days, and the period of the DROs decreases from 11.8 days to 5 days.

#### 5.3.2 DAGDOP Results

Considering the fact that the DAGDOP factor is influenced by the initial phase angles of the arcs (see the previous section), to avoid the influence of different phases, the following method was used. For each combination, we surveyed the initial phase angles of two satellites *τ*_{1} and *τ*_{2} in the range of [0,1]. To save the computation cost, the step size was chosen to be 0.1. For each pair of *τ*_{1} and *τ*_{2}, we computed the average DAGDOP factor for the halo orbit and the DRO, and then we picked up the maximum as the indicator of this constellation.

We have to mention that for some combinations of the constellation, the DAGDOP factor for some combinations of *τ*_{1} and *τ*_{2} was extremely large, which means the AOD process did not converge in those cases. We used *N _{n} =* 1, ⋯, 35 and

*N*1, ⋯, 35 to identify the L2 north and L2 south halo orbits. For example, when

_{s}=*N*=

_{n}*N*and

_{S}*τ*

_{1}=

*τ*

_{2}, this phenomenon occurs. In these situations, we simply ignore them.

It is expected that longer arcs produce better results. We do find such a phenomenon in our simulations, but to save space, only the results of 10-day arcs are displayed in Figure 14.

A careful comparison between the panels in Figure 14(a)–(d) indicates that the overall DAGDOP factor of the hybrid case was smaller. For Cases 1 and 2, the DAGDOP factor was at the level of several hundreds in most areas, with the largest ones exceeding 1,000. For Case 3, the DAGDOP factor was even poorer. The largest ones even exceed 5,000. However, for the hybrid Case 4 in Figure 14, the largest DAGDOP factor was only about 300 and the DAGDOP factor in the upper right area was less than 50. Therefore, from the viewpoint of AOD accuracy, the hybrid Case 4 is recommended.

Apart from this general conclusion, some other useful insights can be obtained. An obvious feature of Figure 14(a) is that the DAGDOP values are symmetrical about the diagonal. This is because the L2 north and south halo orbit families are symmetrical to each other with respect to the x-y plane. The DAGDOP factor would be the largest when *N _{n}* =

*N*(i.e., the north halo and the south halo have the same orbit amplitude). Actually, as we have mentioned above, in this case if

_{s}*τ*

_{1}=

*τ*

_{2}, the AOD process does not converge. This result is understandable. For the case

*N*=

_{n}*N*, the two orbits are strictly symmetrical with respect to the x-y plane. We use:

_{s}42

43

to denote the states of the two satellites in the Moon-centered synodic frame. The superscript represents the satellite number. If *τ*_{1} = *τ*_{2}, then the following conditions always hold:

44

Thus, the observation matrix is written as:

45

where:

46

Therefore:

47

Meanwhile, the state transition matrix is written as:

48

Due to symmetry, many elements in the state transition matrix are the same as each other or are just the opposite number of the other. For example:

49

Therefore, the first element equals the seventh element in Equation (47). Likewise, a similar relationship exists for other elements. In summary, the *i*-th element equals the (*i* + 6)^{th} element, *i* = 1,2,4,5; the *i*-th element equals −1 × (*i* + 6)^{th} element, *i* = 3,6. Consequently, a rank-deficiency problem is caused. Of course, if *τ*_{1} ≠ *τ*_{2}, this rank-deficiency problem can be avoided, but the DAGDOP factor is large as long as *τ*_{1} is too close to *τ*_{2}.

On the other hand, for the case of a large north L2 halo and a small south L2 halo, or for the case of a small north L2 halo and a large south L2 halo, the DAGDOP result is the best. Results in Figure 14(a) suggest that if we build a navigation constellation like the one shown in Figure 13(a) from the viewpoint of AOD accuracy, the best choice is to choose halo orbits with a large difference in the orbit amplitude. For the case of L1 north halo+L1 south halo orbit, the same conclusion holds.

As for Case 2, the diagonal structure which is obvious in Figure 14(a) is no longer strictly diagonal in Figure 14(b) because of the asymmetry between L1 and L2 halo orbits. The cases in which the AOD process failed no longer appear in Case 2. That is because the rank-deficiency problem shown in Equation (44) would be impossible for Case 2.

Figure 14(c) shows the DAGDOP factor of Case 3. The best DAGDOP value for Case 3 occurs for DROs with small in-plane amplitudes. This is because spatial DROs with smaller in-plane amplitudes have shorter periods. Thus, spatial DROs with smaller in-plane amplitudes are recommended. It should be noted that the DAGDOP results of combinations in the diagonal area are also extremely large because of the same reason as in Case 1. These extremely large results are omitted in Figure 14(c).

As for Case 4, an obvious phenomenon in Figure 14(d) is that the DAGDOP factor is smallest for a combination of a large-amplitude halo orbit and a small-amplitude DRO. This suggests that we could use this combination if we were to deploy the satellite pair on such a configuration. The instability of halo orbits reduces with the increasing amplitude (Gao & Hou, 2020), which is also a property favoring this choice.

As a concluding remark to this section, we mention that the same patterns as those shown in Figure 14 are found in the contour maps of AOD accuracy. This, again, demonstrates the feasibility of the DAGDOP factor when assessing AOD accuracy.

## 6 SIMULATION IN REAL-FORCE MODEL

Based on the above analysis in the CRTBP model, it is found that the hybrid halo+DRO configuration shows better overall accuracy when compared with the halo+halo and DRO+DRO configurations, especially when short observation arcs are used. In this section, a Monte Carlo analysis of the hybrid halo+DRO constellation in a realistic Earth–Moon system is conducted.

### 6.1 Simulation Setup

Two satellites move on a halo orbit and a DRO, respectively. The initial epoch for all simulations was January 1, 2020, at 00:00:00.000 UTC. Initial states of the two satellites are shown in Table 5. The initials are in the Moon-centered ICRS. The two orbits integrated for 20 days are shown in Figure 15. The halo orbit is in blue and the DRO is in red. To validate the orbit determination results, a Monte Carlo simulation including 100 tests was conducted. The uncertainty of initial values was modeled as Gaussian noise. The 1*σ* uncertainty of the initials and other parameters to be estimated in the AOD process is summarized in Table 6.

In our work, the force model to generate the SST measurement data is different from the force model used in the AOD. A relatively more accurate dynamic model was used to generate the simulated SST data. A 100 × 100 Earth gravity field and a 100 × 100 Moon gravity field were used. To simulate dynamic model errors, a less accurate dynamic model was used for the AOD process. The lunar gravity model and the Earth gravity model were truncated at 30 × 30. Besides, statistical clones of the gravity models were used to simulate dynamic model errors.

The statistical clones were produced from the full covariance matrix and are provided along with the gravity field coefficients (Goossens et al., 2016; Lemoine et al., 2014). Lunar tides were also modeled with , and computed by Equation (6). When generating the reference orbit, *C _{R}* was a time-varying parameter. But in the AOD process,

*C*was considered a constant parameter, which was estimated together with state vectors. The force model parameter used in the AOD process can be seen in Table 1.

_{R}### 6.2 Simulation Results

100 runs of the AOD process were conducted. Results of the 100 tests are shown in Figure 16, 17, and 18. Satellite 1 moved in a halo orbit while Satellite 2 moved in a DRO. The AOD errors of Satellite 1 were more densely distributed, varying in the range of 150 *~* 220 m. The average position error of Satellite 1 was about 170 meters. The AOD errors of Satellite 2 were distributed between 150 ~ 250 m. The average position error of Satellite 2 was about 190 meters.

The velocity error of Satellite 2 varied between 1.5 mm/s and 2.5 mm/s, while that of the Satellite 1 varied in the range of 1.5 ~ 2 mm/s. Therefore, the halo orbit had a slightly better AOD accuracy than the DRO in our simulated hybrid constellation. What’s more, synchronization errors of two clocks can be corrected to less than 30 ns, as seen in Figure 17. It should be noted that only relative clock errors could be determined. The absolute clock errors could not be determined. In the estimation process, we assumed *C _{R}* as a constant. Figure 18 shows the estimated

*C*value.

_{R}In a word, simulations in this section prove that the AOD of the satellite pair composed of a halo orbit and a DRO is feasible in the real Earth–Moon system. For the case studied, the AOD accuracy of the satellite on the DRO was better than the satellite in the halo orbit. The clocks could also be synchronized, in fact, the accuracy of the clock synchronization was better than 30 ns.

## 7 CONCLUSION

This study focused on the AOD analysis of a satellite pair, which moved in either a halo orbit or a DRO. The satellite pair determined their orbits using only SST measurement data, without the support from ground stations or current GNSS. The two navigation satellites tracked each other, and the observation data was the DOWR data between the two satellites. With the LiAISON strategy, the states of the two satellites were determined simultaneously. In total, three types of configurations for the satellite pair were considered—halo+halo, halo+DRO, and DRO+DRO.

We proposed a factor named *DAGDOP.* This factor reflects the influence of AOD accuracy from both observational geometry and dynamics. Our studies indicated that the DAGDOP factor agreed quite well with the final AOD accuracy. As a result, it is a good indicator to assess the AOD accuracy. It can be easily generalized to a multi-satellite case or other types of observations, so we recommend its use in future similar studies, since we don’t need to carry out the iteration process in orbit determination. This saves a lot of computation time, especially for extensive surveys such as those carried out in the current study.

In the simplified CRTBP model, the AOD performance of combinations of different orbital arcs was studied using the DAGDOP factor. Our studies showed that the DAGDOP factor can approach one (i.e., the AOD results can achieve the level of observation accuracy) if the data length is long enough to cover one complete orbital period of the orbit, which has a longer period. Even longer data length can further decrease the DAGDOP factor, but improvement is generally limited.

In the simplified CRTBP model, the AOD performance of combinations of different configurations was studied also using the DAGDOP factor. Our studies showed that for the halo+halo type constellation, two halo orbits of the same amplitudes around the same collinear libration points should be avoided. Out of the three studied constellation types, the hybrid halo+DRO configuration performed best. The weaker instability of a nearly rectilinear halo orbit and the stability of a DRO orbit is preferred when building such a navigation constellation.

At last, the case of the hybrid halo+DRO constellation was generalized to a more realistic model of the Earth–Moon system. In our model, particle gravity of the Earth, the Moon, and the Sun, along with the Earth’s and Moon’s non-spherical gravity (lunar tides included) and solar radiation pressure were considered. The model used to generate SST measurements was different from the model used for the AOD in order to simulate the dynamic errors and uncertainties. The DOWR system was used to generate the SST measurements, with consideration of time delay, accurate clock model, signal propagation loss, and so on. Studies in such a realistic model show that the AOD process is feasible. In this case, the average position accuracy of a halo orbit was about 170 meters. The AOD accuracy of the halo orbit was slightly better than that of the DRO, which was about 190 meters on average. The two clocks were synchronized with a time accuracy of 30 nanoseconds.

The background behind the current study is the construction of an autonomous lunar navigation constellation. The basic idea is that the navigation constellation first determines its orbits autonomously with SST data and without support from ground stations or the GNSS around the Earth, and then broadcasts the navigation message to its users. Considering the cost and service volume, the simplest two-satellite constellation using the halo orbit and the DRO was considered in the current study. The proposed DAGDOP factor greatly helped us in assessing the AOD accuracy of different configurations of the navigation satellite pair.

Besides the AOD performance of the constellation itself, we have to admit that there are a lot of other factors that need to be considered when building such a constellation as well—like coverage ability, the ability to serve different users, and the possibility of including multiple satellites in the constellation. The best configuration of the criteria for AOD accuracy alone for the navigation constellation itself might not be the best choice when other factors are involved.

So, the current study serves only as a very preliminary attempt to construct such a navigation constellation. The proposed DAGDOP factor may be helpful for future studies in which more navigation satellites can be added to the constellation and the user satellites are involved. In a forthcoming paper, we will study the service ability of navigation constellations to different types of users around the Moon using this factor. This study will include a visibility analysis and a discussion of the influence on users’ position errors from the navigation satellites’ errors for different types of users around the Moon.

## HOW TO CITE THIS ARTICLE

Gao, Z.-Y., & Hou, X.-Y. (2022) Comparison of autonomous orbit determination for satellite pairs in lunar halo and distant retrograde orbits. *NAVIGATION, 69*(2). https://doi.org/10.33012/navi.522

## ACKNOWLEDGEMENTS

The authors thank the two anonymous reviewers and the editor for their valuable input, which greatly helped improve the paper.

This work is supported by national Natural Science Foundation of China (NSFC 11773017, 11703013, 11673072).

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.