Comparison of Autonomous Orbit Determination for Satellite Pairs in Lunar Halo and Distant Retrograde Orbits

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₁ is the Earth center and P₂ is the Moon center. μ = m₂/(m₁ + m₂) is the mass parameter, with m₁ and m₂ as masses of the Earth and the Moon, respectively. Following mass, length and time units are used throughout our study:

where P₁P₂ means the distance between P₁ and P₂.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Geometry of the CRTBP model in the Earth–Moon barycentric synodic frame. P₁ is the Earth center and P₂ is the Moon center.

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₁: particle gravity of the Earth, Moon, and Sun;

• F₂: non-spherical gravity of the Earth and the Moon, including the tidal effects of the Moon; and

• F₃: 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, μ_E, and μ_S represent the gravitational constants of the Moon, the Earth, and the Sun, respectively. r, r_E, and r_S are position vectors of the satellite, the Earth, and the Sun in the Moon-centered ICRS. Δ_E (Δ_s) is the position vector that points from the Earth (Sun) to the satellite. ∇U_E and ∇U_M 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:

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_j 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.

Equation of F_SRP is:

7

where P_SR represents the solar radiation pressure at 1 astronomical unit (AU). C_R 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 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.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^R, Satellite 2 receives that signal. τ = t^R − t^T is the time delay in this ranging process. A simple iteration process is used to obtain τ, which is shown in Figure 4. t^R is assumed to equal t^T as an initial guess. Then, R₂(t^R) 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 ρ is calculated and τ is obtained. Then, a new t^R is obtained. The iteration is finished until the difference between τ_i and τ_i−1 is less than ϵ, which is set as 10⁻⁸ seconds in this work.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

The iteration process used to generate accurate time delay τ

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_s, respectively. Usually, δt_d can be modeled as the following second-order polynomial:

9

where a₀, a₁, and a₂ 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

A time history sample of the stochastic error in the atomic clock

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_n 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; T is the prediction integration time; and T = 1s is used in this work. C/N₀ is the carrier-to-noise ratio (Hz):

13

where (N₀)_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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Variation of the 1σ of thermal noise with respect to 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7

Time slots of two satellites

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₁ and δt₂. 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.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₀ are obtained. Then, the corrections to ρ₁, and ρ₂ are expressed as:

16

Then, DOWR measurements at time t₀ 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₀ 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₀ and a, we denote:

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), ρ₋(a) and the calculated values of the estimated trajectory at the epoch t_i as:

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₀) is the state transition matrix from the epoch t₀ to the observation epoch t_i.

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⁽ⁱ⁾ 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^TH)⁻¹ 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 6N-dimensional vector and the matrix D is a 6N × 6N 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 X 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.

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₀)/p_i, where t − t₀ represents the time from the initial point (τ_i = 0) and p_i is the period of the corresponding orbit (We use one to indicate the DRO and two to indicate the halo orbit).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 8

Combinations of the two arcs of the halo orbit and the DRO

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 τ₁ ∈ [0,1] and τ₂ ∈ [0,1] to represent different combinations of the two arcs. A step size of 0.01 was chosen when we surveyed different combinations of τ₁ and τ₂.

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:

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 9

The DAGDOP factor for different combinations of initial phases of the halo orbit and for the different lengths (5 days, 10 days, 15 days, and 20 days) of arcs

• For an arc length of 5 days, the DAGDOP factor varies in a range from 10^2.5 ~ 10⁴. 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⁰ (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:

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 10

The DAGDOP factor for different combinations of initial phases of the DRO and for the different lengths (5 days, 10 days, 15 days, and 20 days) of arcs

• 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_i). The position vector at t_i obtained from the AOD process is . The 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 11

The position RMS for different combinations of initial phases of the halo orbits and for the different lengths (5 days, 10 days, 15 days, and 20 days) of arcs

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 12

The position RMS for different combinations of initial phases of the DROs and for the different lengths (5 days, 10 days, 15 days, and 20 days) of arcs

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.

View this table:

• View inline
• View popup

TABLE 4

Four Cases of Constellations

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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 13

Orbits for four cases shown in the Moon-centered synodic frame: (a) Case 1; (b) Case 2; (c) Case 3; and (d) Case 4. The orbit numbers are marked with arrows pointing to the corresponding orbit.

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 τ₁ and τ₂ in the range of [0,1]. To save the computation cost, the step size was chosen to be 0.1. For each pair of τ₁ and τ₂, 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 τ₁ and τ₂ was extremely large, which means the AOD process did not converge in those cases. We used N_n = 1, ⋯, 35 and N_s = 1, ⋯, 35 to identify the L2 north and L2 south halo orbits. For example, when N_n = N_S and τ₁ = τ₂, 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 14

The DAGDOP factor for four cases of a 10-day arc: (a) Case 1; (b) Case 2; (c) Case 3; and (d) Case 4

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_s (i.e., the north halo and the south halo have the same orbit amplitude). Actually, as we have mentioned above, in this case if τ₁ = τ₂, the AOD process does not converge. This result is understandable. For the case N_n = N_s, the two orbits are strictly symmetrical with respect to the x-y plane. We use:

42

43

to denote the states of the two satellites in the Moon-centered synodic frame. The superscript represents the satellite number. If τ₁ = τ₂, 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 τ₁ ≠ τ₂, this rank-deficiency problem can be avoided, but the DAGDOP factor is large as long as τ₁ is too close to τ₂.

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.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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 15

Two satellites moving on a halo orbit and a DRO are integrated over 20 days; the left frame shows the orbits in the Moon-centered ICRS; the right frame shows the orbits in the Moon-centered synodic frame

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_R 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.

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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 16

Position (left) and velocity (right) errors of two satellites; Satellite 1 moved in a halo orbit while Satellite 2 moved in a DRO.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 17

Clock parameter synchronization error between Satellite 1 and Satellite 2

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 18

Solutions of solar radiation reflectivity coefficient C_R

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_R value.

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.