Integrated Design of Moon-to-Earth Transfer Trajectory Considering Re-Entry Constraints

Abstract

The exploration of the Moon has always been a hot topic. The determination of the Moon-to-Earth transfer opportunity and the design of the precise transfer trajectory play important roles in manned Moon exploration missions. It is still a difficult problem to determine the Moon-to-Earth return opportunity for accurate atmospheric re-entry and landing, through which the actual return trajectory can be easily obtained later. This paper proposes an efficient integrated design method for Moon-to-Earth window searching and precise trajectory optimization considering the constraints of Earth re-entry and landing. First, an analytical geometry-based method is proposed to determine the state of the re-entry point according to the landing field and re-entry constraints to ensure accurate landing. Next, the transfer window is determined with the perilune heights, which are acquired by inversely integrating the re-entry state under the simplified dynamics as criterion. Then, the precise Moon-to-Earth trajectory is quickly obtained by a three-impulse correction. Simulations show the accuracy and efficiency of the proposed method compared with methods such as the patched-conic method and provide an explicit reference for future Moon exploration missions.

1. Introduction

As the only natural satellite of the Earth, the Moon is an important object for studying the formation and evolution of the Earth. Exploring the Moon not only can help human beings to further understand cosmic phenomena but can also promote the rapid development of human scientific and technological progress and applications [1,2]. Human exploration of the Moon began in the 1950s [3] as an important part of deep space exploration [4,5,6,7,8,9,10,11,12]. The second round of the lunar exploration climax was revived in the 1990s after a quiet period. The ESA, China, Japan, and India have carried out lunar exploration activities one after another since then [13,14,15,16,17]. Chinese lunar exploration has developed rapidly over the recent decades. Chang’E-2, launched on 1 October 2010, obtained 3D images of the lunar surface [18], completed the expansion mission of entering the orbit of the Sun–Earth L2 libration point [19,20], and flew by asteroid 4179 Toutatis at an extremely close distance [18,21]. Chang’E-5 was launched from the Wenchang Satellite Launch Site on 24 November 2020 [22] and successfully returned to the ground with lunar samples 23 days later [23,24]. Follow-up lunar exploration schemes have also been demonstrated.

The determination of the opportunity and trajectory of the Moon-to-Earth return is the premise for the implementation of lunar sampling return missions or manned exploration missions. The manned lunar return has higher requirements for reliability and emergency response compared with unmanned missions. It should be seriously noted that the Moon-to-Earth transfer is not a simple inverse process of the Earth-to-Moon process, for the selection of initial values, model establishment, and constraints are evidently different [25,26].

For the Earth-to-Moon trajectory design, researchers have developed multiple trajectory design technologies. The traditional impulse transfer is the classic method adopted by the first Moon exploration missions [27,28], while the small thrust transfer [29,30,31,32,33,34,35] and the weak stable boundary (WSB) transfer methods [36,37,38] appeared in the 1990s. In recent years, there have been the use of Moon flybys to achieve libration point transfers and even the deployment of Earth’s high-orbit satellites [39,40].

For the design of the Moon-to-Earth return for direct atmosphere re-entry, Zhang et al. [41] obtained the transfer trajectory with matching landing field through two-level differential correction. Robinson and Geller [42] adopted the analytical method to optimize the fastest return trajectory according to the re-entry state, regardless of the fuel consumption. He et al. [43] established the solution model of landing opportunity and point return trajectory from the perspective of space geometry. Feng et al. [44] studied the return opportunity strategy departing from the lunar and provided the initial value of the transfer trajectory by applying the improved patched-conic method. Zheng et al. [45] proposed a Moon return opportunity search strategy based on the patched-conic method. The return opportunity is initially selected according to the minimum re-entry angle, and the transfer trajectory is further obtained, while it is difficult to guess the state of the probe at the sphere of influence (SOI) of the Moon. In general, what the above methods obtain is the transfer opportunity with a single impulse of alignment at the positions of the Earth and the Moon. In addition, the simplified transfer, such as the patched-conic method and space geometry strategy, causes a low accuracy of opportunity and trajectory because of the lack of consideration for third-body perturbation, and ignoring any celestial body will lead to a large deviation between the designed and the actual opportunity or trajectory.

The improvement of the calculation efficiency while ensuring the accuracy of the opportunity and the transfer trajectory are the keys to lunar return missions. Aiming at the problems above, an integrated design method for the opportunity and precise transfer trajectory is proposed in this paper. The application of geometric principles for re-entry and window and precise trajectory determination based on maneuverability are the main innovations. First, an analytical geometry method is brought up to determine the state of the probe at the re-entry point according to the re-entry constraints to ensure an accurate landing. Second, the obtained re-entry state is inversely integrated to acquire the daily transfer trajectory with a certain perilune altitude, and the Moon-to-Earth transfer opportunity is determined by delimiting the maximum perilune altitude as a boundary, which is closely related to the maneuverability of the probe. The opportunity obtained is more accurate and with a wider range compared to the patched-conic method. Then, a three-impulse correction strategy is brought forward to rapidly obtain the precise transfer trajectory within the opportunity, which enhances the emergency capability of the trajectory. Furthermore, the effects of several important design parameters on the Moon-to-Earth transfer opportunity and the performance of the transfer trajectory are also exhaustively analyzed to provide a more explicit reference for future Moon exploration missions. The proposed integrated design method obtains the transfer opportunity, which efficiently provides a reliable initial value for precise transfer trajectory.

The structure of the paper is as follows. The precise dynamical model and the definition of the re-entry coordinate frame are presented in Section 2. The Moon-to-Earth transfer opportunity determination method with strict re-entry constraints is described in Section 3. The precise fuel-saving three-impulse transfer trajectory determination is shown in Section 4. In Section 5, the numerical simulation and parameter analysis is presented. Finally, conclusions are drawn in Section 6.

2. Dynamical Model and Frame Construction

In this paper, the selection of transfer opportunity and the design of precise transfer trajectory are carried out under the simplified Earth–Moon ephemeris model and the high-precision dynamic model, respectively.

The precise transfer trajectory design should be carried out in a high-precision dynamic model, which is established as

$${\ddot{\mathit{r}}}_e=-{\mu }_E·\frac{{\mathit{r}}_e}{{∥{\mathit{r}}_e∥}^3}+{\mu }_M·\left(\frac{{\mathit{r}}_M-{\mathit{r}}_e}{{∥{\mathit{r}}_M-{\mathit{r}}_e∥}^3}-\frac{{\mathit{r}}_M}{{∥{\mathit{r}}_M∥}^3}\right)+{\mu }_S·\left(\frac{{\mathit{r}}_S-{\mathit{r}}_e}{{∥{\mathit{r}}_S-{\mathit{r}}_e∥}^3}-\frac{{\mathit{r}}_S}{{∥{\mathit{r}}_S∥}^3}\right)+{\mathit{a}}_{NE}+{\mathit{a}}_{NM}+{\mathit{a}}_{srp}$$

(1)

where ${\ddot{\mathit{r}}}_e$ is the acceleration of the probe in the Earth–Moon gravitational model; ${\mu }_E$, ${\mu }_M$, and ${\mu }_S$ are the gravitational constants of the Earth, the Moon, and the Sun, respectively; ${\mathit{r}}_e$ is the position vector of the probe in the Earth-centered inertial (ECI) coordinate frame; ${\mathit{r}}_M$ and ${\mathit{r}}_S$ are the position vectors of the Moon and the Sun in ECI, respectively; ${\mathit{a}}_{NE}$, ${\mathit{a}}_{NM}$, and ${\mathit{a}}_{srp}$ are the Earth’s non-spherical gravitational perturbation, the Moon’s non-spherical gravitational perturbation, and the solar radiation pressure perturbation, respectively. Since the gravitational effect of the probe during the Moon-to-Earth transfer process mainly comes from the Earth and the Moon, the simplified dynamics considering the gravities of the Earth and the Moon is adopted when determining the feasible return trajectory and opportunity, which greatly improves the efficiency while ensuring the reliability of the design.

To describe the state of the probe under the constraints of the re-entry condition, a re-entry coordinate frame (RF) OXYZ fixed to the center of the probe, which is shown in Figure 1, is established. The rectangle is the position of the re-entry point marked as P with an altitude of H. The Y-axis of the frame points to the zenith, and the X-axis is vertical to the Y-axis and points to the flight direction in the re-entry plane. The Z-axis forms a right-hand coordinate frame with the X-axis and Y-axis.

The re-entry point is noted as P, which is marked with a triangle in Figure 1. The longitude and latitude of the re-entry point are ${\lambda }_P,{\phi }_P$, respectively. The angle between the direction of the velocity at the re-entry point on the local sphere (noted as $\widehat{\mathit{r}}$) and the local north direction is defined as the flight azimuth Q.

The rotation matrix of the re-entry coordinate frame to the Earth-centered Earth-fixed (ECEF) coordinate frame can be obtained by Euler rotation as

$${\mathit{M}}_{EF}=\left[\begin{array}{ccc}-cosQcos{\lambda }_{P}sin{\phi }_P-sinQsin{\lambda }_P& cos{\lambda }_{P}cos{\phi }_E& sinQcos{\lambda }_{P}sin{\phi }_P-cosQsin{\lambda }_P\\ -cosQsin{\lambda }_{P}sin{\phi }_P+sinQcos{\lambda }_P& sin{\lambda }_{P}cos{\phi }_E& sinQsin{\lambda }_{P}sin{\phi }_P+cosQcos{\lambda }_P\\ cosQcos{\phi }_P& sin{\phi }_P& -sinQcos{\phi }_P\end{array}\right]$$

(2)

3. Opportunity Determination with Re-Entry Constraints

Earth re-entry is the last step in the realization of a Moon return, and it is also essential to the success of the mission. The atmospheric re-entry process generally only adjusts the altitude of the re-entry capsule without performing large-scale maneuvers, so the determination of the state of the re-entry point is the key to the design of the Moon-to-Earth transfer trajectory. Therefore, an analytical method is proposed to determine the state of the probe at the re-entry point with the necessary re-entry and landing field constraints to ensure accurate landing. Then, the high-precision Moon-to-Earth transfer opportunity is screened out by the altitude deviation of the perilune of the daily optimal transfer trajectory.

3.1. Geometry-Based Re-Entry State Determination

The state of the probe at the re-entry point can be obtained by the longitude and latitude of the landing field, the altitude of the re-entry point, the re-entry velocity, the voyage, and the inclination of the re-entry trajectory.

As shown in Figure 1, the latitude of the landing field F (marked with a five-pointed star) is chosen as ${\phi }_P$, and the inclination of the re-entry trajectory is set as $i_r$. The argument of the landing field can be determined according to the sine law of the spherical triangle as

$$U_F=arcsin\frac{sin{\phi }_F}{sin{i}_r}\left(\mathrm{Rail-up}\right),U_F=\pi -arcsin\frac{sin{\phi }_F}{sin{i}_r}\left(\mathrm{Rail-down}\right)$$

(3)

According to the mission constraints, the re-entry trajectory voyage is set as $S_r$, and the argument of the trajectory is

$$U_r=\frac{S_r}{R_e}$$

(4)

where $R_s$ represents the radius of the Earth. Then, the argument of the re-entry point can be obtained as

$$U_P=U_F-U_r$$

(5)

The longitude of the landing field is chosen as ${\lambda }_F$, and the longitude difference from the re-entry point to the ascending node of the trajectory can be obtained according to Napier’s Law of the spherical triangle as

$$\Delta {\lambda }_P=arctan\frac{cos{i}_r}{tan{U}_P}$$

(6)

Therefore, the flight azimuth of the re-entry point is

$$Q=arcsin\frac{sin\Delta {\lambda }_P}{sin{U}_P}$$

(7)

Based on the landing field argument $U_F$, the longitude difference between the landing field and the ascending node of the re-entry trajectory is calculated as

$$\Delta {\lambda }_F=arctan\frac{cos{i}_r}{tan{U}_F}$$

(8)

The longitude of the re-entry point is determined as

$${\lambda }_P={\lambda }_F-\Delta {\lambda }_F+\Delta {\lambda }_P$$

(9)

and the latitude of the re-entry point is

$${\phi }_P=arcsin\left(sin{U}_{P}sin{i}_r\right)$$

(10)

The altitude of the re-entry point is chosen as H according to the mission constraint, and the geocentric distance of the re-entry point is

$$R=H+R_e$$

(11)

Therefore, the position vector of the re-entry point in ECEF is determined as

$${\mathit{r}}_e=\left[Rcos{\phi }_{P}cos{\lambda }_P,Rcos{\phi }_{P}sin{\lambda }_P,Rsin{\phi }_P\right]$$

(12)

$\gamma$ and $v_{re}$ are denoted as the re-entry angle and the preset magnitude of the re-entry velocity, respectively; the velocity vector of the probe at the re-entry point in the re-entry coordinate frame is

$${\mathit{v}}_{re}=\left[v_{re}cos\gamma ,v_{re}sin\gamma ,0\right]$$

(13)

Then, the velocity vector of the probe at the re-entry point in ECEF is

$${\mathit{v}}_e={\left({\mathit{M}}_{EF}·{\mathit{v}}_{re}^{\mathrm{T}}\right)}^{\mathrm{T}}$$

(14)

where the superscript “T” is the transpose symbol.

Thus, the state of the probe at the re-entry point in ECEF is

$${\mathit{X}}_e=\left[{\mathit{r}}_e,{\mathit{v}}_e\right]$$

(15)

The procedures for geometry-based re-entry state determination are shown as module M1 in Figure 2.

3.2. Daily Optimal Transfer Trajectory Design

Since the Earth rotation period is one day, which is much smaller than the revolution period of the Moon, there is a daily optimal Moon-to-Earth transfer trajectory that can be obtained by inverse integration from the re-entry point to the vicinity of the moon. Due to the existence of the inclination between the Moon’s orbit and Equator, the perilune altitudes of the Moon of the optimal trajectories on some dates can reach tens of thousands of kilometers or even larger. An excessive perilune altitude means that too much fuel is consumed for adjustment, which cannot be achieved in actual missions. Therefore, the Moon-to-Earth transfer opportunity can be screened according to the altitude of the perilune every day.

The Julian date range of the default re-entry time is $D_1$ 0:00 to $D_1$ + 10:00, and the selected re-entry point time is $T_0\in \left[D_1,D_1+1\right]$. The duration of the Moon-to-Earth transfer is given as $\Delta T$; then, the time of the perilune is $T_m=T_0-\Delta T$, and the position of the Moon at that time obtained from the ephemeris is ${\mathit{r}}_M$ = $\left[x_M,y_M,z_M\right]$.

The Greenwich time corresponding to the re-entry time is GST; then, the rotation matrix from ECEF to ECI is

$${\mathit{M}}_{FI}=\left[\begin{array}{ccc}cos\left(GST\right)& -sin\left(GST\right)& 0\\ sin\left(GST\right)& cos\left(GST\right)& 0\\ 0& 0& 1\end{array}\right]$$

(16)

Then, the state of the probe at the re-entry point in ECI is

$${\mathit{X}}_i=\left[{\left({\mathit{M}}_{FI}·{\mathit{r}}_e^{\mathrm{T}}\right)}^{\mathrm{T}},{\left({\mathit{M}}_{FI}·{\mathit{v}}_e^{\mathrm{T}}\right)}^{\mathrm{T}}\right]$$

(17)

The probe is reversely integrated from the re-entry point to the perilune of the Moon to obtain a feasible transfer trajectory with an actual transfer duration of $\Delta \underline{T}$.

In order to make the actual transfer duration equal to the set value $\Delta T$, the magnitude of the re-entry velocity needs to be adjusted. If $\Delta \underline{T}<\Delta T$, this means that the re-entry velocity is set too high, and the probe will reach the perilune ahead of time; if $\Delta T>\Delta T$, this means that the re-entry velocity is set too low, and the probe lags behind the perilune.

The magnitude of the re-entry velocity is adjusted iteratively so that the transfer duration is exactly $\Delta T$. Determining the adjustment amount of the re-entry velocity in each iteration is the key, so a method is proposed to determine the adjustment step according to the deviation between the actual and the expected duration.

Assuming that the deviation between the actual transfer duration and the set transfer duration in the iteration is $\delta t=\Delta T-\Delta T$, the adjustment step magnitude of the re-entry velocity in the next iteration is

$$\Delta v_{re}=\delta t/100$$

(18)

where the unit of $\delta t$ is “day”. The magnitude of the re-entry velocity is adjusted to

$$v_{re}=v_{re}+\Delta v_{re}$$

(19)

When $\delta t<0,\Delta v_{re}<0$, the re-entry velocity is reduced, and when $\delta t>0,\Delta v_{re}>0$, the re-entry velocity is increased. Based on Equation (19), the re-entry velocity is iterated until the transfer duration accuracy requirement is met as

$$|\Delta \underline{T}-\Delta T|<{\epsilon }_T$$

(20)

where ${\epsilon }_T$ is the error limit of the accuracy of the transfer duration, which can be set according to the actual demands.

Thus far, the feasible Moon-to-Earth transfer trajectory with the transfer duration fixed as $\Delta T$ and the re-entry time as $T_0$ has been obtained.

The Earth rotates once a day. When $T_0$ changes in one day, the orbital direction of the feasible transfer will rotate once. Thus, the altitude of the perilune of the feasible Moon-to-Earth transfer trajectory obtained above changes approximately periodically in one day. However, the maximum altitude of the perilune can even be up to tens of thousands of kilometers. Therefore, it is necessary to find the re-entry time corresponding to the transfer trajectory with the lowest perilune altitude of lunar in one day. The change curve of the perilune altitude of the transfer trajectory with a fixed transfer duration of $\Delta T$ in one day shows a single valley trend, so the golden section method is applied.

For the optimal Moon-to-Earth transfer trajectory solution problem, the specific steps to realize the golden section method are as follows.

Take points $a_1,a_2$ in the time interval $[a,b]\left(a=D_1,b=D_1+1\right)$ and divide the interval into three sections, among which

$$a_1=a+0.382(b-a),a_2=a+0.618(b-a)$$

(21)

Denote $F\left(x\right)$ as the altitude of the perilune corresponding to the re-entry time of x. If $F\left(a_1\right)<F\left(a_2\right)$, set $b=a_2$, and a remains unchanged; if $F\left(a_1\right)>F\left(a_2\right)$, set $a=a_1$, and b remains unchanged. The above iterative process is implemented until the optimal re-entry time accuracy is met as

$$b-a<{\epsilon }_{rp}$$

(22)

where ${\epsilon }_{rp}$ is the error limit of the optimal re-entry time.

Therefore, when the transfer duration is $\Delta T$, the optimal re-entry time is $T_0=a$. Thus, the optimal Moon-to-Earth transfer trajectory with the minimum altitude of the perilune in $\left[D_1,D_1+1\right]$ is obtained.

The steps for optimal Moon-to-Earth transfer trajectory design are shown as module M2 in Figure 2.

3.3. Transfer Opportunity Determination

Given that the mission search date interval is $\left[D_2,D_3\right]\left(D_3-D_2>1\right)$, the optimal transfer trajectory and its corresponding perilune altitude H on each day can be solved with the above calculation process.

The fuel consumption required to achieve Moon capture after trajectory correction is highly correlated with the altitude of the perilune before correction. If the perilune altitude of the transfer trajectory is too large, the fuel consumption for correction will be too large, which will be beyond the capability of the probe. Therefore, the existence of the daily transfer opportunity can be judged according to the daily minimum perilune altitude of the Moon.

The altitude deviation limit of the perilune is set as $\chi$. If the daily minimum perilune altitude of the Moon-to-Earth transfer trajectory at date $\tau$ satisfies $H<\chi$, the date $\tau$ is judged to be included in the transfer opportunity. Then, the whole transfer opportunity $\Gamma$ within the date interval $\left[D_2,D_3\right]$ is obtained.

Module M3 in Figure 2 reveals the transfer opportunity selecting process.

4. Precise Three-Impulse Transfer Trajectory Determination

In the previous section, the opportunity for Moon-to-Earth transfer in which the daily trajectory can be reversely integrated to the vicinity of the Moon was obtained. This section aims at the design of the precise transfer trajectory whose states at the perilune and the re-entry point precisely meet the constraints on a certain day, which is included in the opportunity.

The transfer strategy proposed in this paper is to apply correction maneuvers at the boundary of the SOI of the Moon and one day before re-entry, so that the transfer trajectory can accurately meet the constraints of Moon departure and Earth re-entry under the high-precision dynamic model. This strategy also provides a more powerful guarantee for the transfer in case of an emergency.

The perilune altitudes of the daily optimal transfer trajectories vary from one hundred to several thousand kilometers, which means different magnitudes of maneuvers are acquired to ensure a high precision transfer.

The B-plane shooting method is applied in the design of the two maneuvers, while different correction parameters are chosen for the smoothness and fuel saving of the transfer trajectory. Three orbital elements of the semi-major axis, radius of ascending node, and argument of perigee are selected to adjust the position and velocity of the probe at the boundary of the SOI of the Moon. Compared to the velocity at the boundary of the SOI as the correction, this strategy may require smaller maneuvers for those cases where the feasible perilune heights are large. Usually, the positional offset of the trajectory at the boundary of the SOI after correcting is positively correlated with the initial perilune altitude deviation. Finally, a modified maneuver is applied the day before re-entry, choosing the initial velocity as the correcting parameter to precisely splice the whole transfer trajectory.

4.1. Perilune to the Boundary of SOI

Choose any $T_0\in \Gamma$ in the Moon-to-Earth transfer opportunity; the feasible transfer trajectory with a state of $X_i$ at the re-entry point that reached the vicinity of the Moon can be obtained after the calculation in the previous section.

Integrate the initial value ${\mathit{X}}_i$ inversely in the high-precision dynamic Equation (1), and the state of the perilune in the ECI is obtained. Transform it to ${\mathit{X}}_f=\left[{\mathit{r}}_f,{\mathit{v}}_f\right]$ in the Moon-centered inertial (MCI) coordinate frame. With orbital angular momentum $\mathit{h}$, semi-major axis $\mathit{a}$, eccentricity $\mathit{e}$, semilatus rectum p, true anomaly f, and eccentric anomaly H as intermediate variables, calculate the time from perilune $\tau$, altitude of perilune ${\tau }_{pm}$, and inclination of perilune $i_m$.

Note the corresponding relationship between the scalars and the vectors as

$$a=∥\mathit{a}∥,e=∥\mathit{e}∥,h=∥\mathit{h}∥,r_f=∥{\mathit{r}}_f∥,v_f=∥{\mathit{v}}_f∥$$

(23)

Then, the time from the perilune, the perilune radius, and the inclination at the perilune can be obtained as [46]

$$\tau =\sqrt{\frac{a^3}{{\mu }_M}}(esinhH-H),r_{pm}=a(1-e),i_m=arccos\left(\frac{\mathit{h}·\mathit{n}}{h}\right)$$

(24)

where $\mathit{n}=[0,0,1]$.

The design of the high-precision transfer trajectory in the SOI of the Moon is realized according to the nominal departure parameters at the perilune. Given the radius and inclination of the perilune as ${\tilde{r}}_{pm}$ and ${\tilde{i}}_m$, respectively, the differences between the quantities of the time from the perilune $\tau$, the actual perilune altitude $r_{pm}$, and the perilune inclination $i_m$ and the nominal values are chosen as the target constraint variables (the nominal time from the perilune is 0), so the solution of the Moon-to-Earth transfer trajectory can be transformed into the following nonlinear equations as

$$\left[r_{pm}-r_{pm}^{\prime },i_m-i_m^{\prime },\tau \right]=0$$

(25)

After the inverse integration of the initial value ${\mathit{X}}_i$ in the high-precision dynamic model, the state of the probe at the boundary of the SOI of the Moon in the ECI is obtained and converted into the form of orbital elements. The correction amounts are selected as the semi-major axis a, the right ascension of the ascending node $\Omega$, and the argument of the perigee w at the boundary of the SOI of the Moon, and the influence of the correction amounts on the target constraint variables are analyzed by numerical methods. Respectively giving a small disturbance $\delta$ to the initial state of each correction amount, the response of the target constraint variable to the perturbation can be obtained after inverse integration into the perilune. Then, the sensitivity matrix ${\mathbf{K}}_1$ of the target constraint variables relative to the initial states of the correction amounts can be obtained as

$${\mathbf{K}}_1=\left[\begin{array}{ccc}\frac{\partial r_{pm}}{\partial a}\hfill & \frac{\partial r_{pm}}{\partial \Omega }\hfill & \frac{\partial r_{pm}}{\partial \omega }\hfill \\ \frac{\partial i_m}{\partial a}\hfill & \frac{\partial i_m}{\partial \Omega }\hfill & \frac{\partial i_m}{\partial \omega }\hfill \\ \frac{\partial \tau }{\partial a}\hfill & \frac{\partial \tau }{\partial \Omega }\hfill & \frac{\partial \tau }{\partial \omega }\hfill \end{array}\right]$$

(26)

Then, the relationship between the differential of the correction amounts and the error of the target constraint variables is

$$\left[\begin{array}{c}\Delta a\\ \Delta \Omega \\ \Delta \omega \end{array}\right]={\mathbf{K}}_1^{-1}\left[\begin{array}{c}\delta r_{pm}\\ \delta i_m\\ \delta \tau \end{array}\right]$$

(27)

After several iterations, the initial states $a,\Omega ,w$ at the boundary of the SOI of the Moon that satisfy the nonlinear Equation (25) can be obtained. Thus, the state of the probe at the boundary of the SOI of the Moon that satisfies the given constraints at the perilune is obtained as ${\mathit{X}}_{\mathrm{SOI}}=\left[{\mathit{r}}_{\mathrm{SOI}},{\mathit{v}}_{\mathrm{SOI}}\right]$. Using this state as the initial value, perform a high-precision integration on the trajectory in the SOI of the Moon, and the perilune state ${\mathit{X}}_M=\left[{\mathit{r}}_M,{\mathit{v}}_M\right]$ is obtained.

When the probe applies the departure impulse from the lunar circular orbit with altitude $r_{pm}$ to insert the high-precision Moon-to-Earth transfer trajectory, the required velocity increment is

$$\Delta {\mathit{v}}_M=\left(∥{\mathit{v}}_M∥-\sqrt{{\mu }_M/\left(R_M+r_{pm}\right)}\right)·{\mathit{v}}_M/∥{\mathit{v}}_M∥$$

(28)

where $R_M$ is the radius of the Moon.

The steps for precise transfer trajectory design in the SOI of the Moon are shown as module M4 in Figure 3.

4.2. Boundary of SOI to the Day before Re-Entry

Compared with the state before the correction, the position and velocity of the probe at the boundary of the SOI of the Moon after the correction have small changes. Thus, the trajectory from the re-entry point to the boundary of the SOI of the Moon needs to be corrected to align with the probe status at the boundary of the SOI.

The correction impulse is applied one day before Earth’s re-entry to align the position of the probe at the boundary of the SOI of the Moon. Then, the whole transfer trajectory is divided into three segments: (1) perilune–boundary of the SOI of the Moon; (2) boundary of the SOI of the Moon–one day before re-entry; and (3) one day before re-entry–re-entry. Three impulses are applied in the whole process, that is, the connection impulse between every two segments in addition to the departure impulse at the perilune of Equation (28).

The third segment of trajectory inverse integration is based on the initial value ${\mathit{X}}_i$ of the re-entry position, and the state of the probe ${\mathit{X}}_{rr}=\left[{\mathit{r}}_{rr},{\mathit{v}}_{rr}\right]$ on the day before re-entry is obtained under the high-precision dynamic model, where ${\mathit{r}}_{rr}=\left[x_{rr},y_{rr},z_{rr}\right],$${\mathit{v}}_{rr}=\left[v_{xr},v_{yr},v_{zr}\right]$. Select the terminal constraint as the position of the probe at the boundary of the first segment in the SOI of the Moon obtained above.

The value of the probe velocity increment $\Delta {\mathit{v}}_{rr}=\left[\Delta v_{xr},\Delta v_{yr},\Delta v_{zr}\right]$ on the day before re-entry, which makes the position of the probe at the boundary of the SOI of the Moon accurately aligned, is obtained by adopting a B-plane correction. Therefore, the initial velocity on the day before re-entry in the second segment is ${\mathit{v}}_{rr}^{\prime }={\mathit{v}}_{rr}+\Delta {\mathit{v}}_{rr}$. In order to correspond to the initial state of the probe in the third segment, the velocity increment that should be applied by the probe at the boundary of the SOI of the Moon in the corresponding second segment is $\Delta {\mathit{v}}_{\mathrm{SOI}}={\mathit{v}}_{\mathrm{SOI}}-{\mathit{v}}_{\mathrm{SOI}}^{\prime }$.

The procedures for precise transfer trajectory design out of the SOI of the Moon are shown as module M5 in Figure 3.

Thus far, the high precision Moon-to-Earth transfer trajectory with the time of re-entry $T_0$ in the opportunity $\Gamma$ has been obtained, and the total maneuver required by the probe is

$$∥\Delta {\mathit{v}}_{all}∥=∥\Delta {\mathit{v}}_M∥+∥\Delta {\mathit{v}}_{\mathrm{SOI}}∥+∥\Delta {\mathit{v}}_{rr}∥$$

(29)

5. Numerical Simulation and Analysis

In this section, the specific precise transfer trajectory and the opportunity in a year are simulated to verify the reliability of the integrated determination method. Moreover, the effects of several design parameters on the opportunity and performance of the transfer trajectory are also analyzed to provide instructions for parameter selection.

5.1. Opportunity and Specific Precise Trajectory Verification

Firstly, determine the re-entry state in ECEF that meets the re-entry conditions. Select the longitude and latitude of the landing field as ${\lambda }_F=101.{45}^{\circ }{\mathrm{E}}_{,}{\phi }_F=-41.2^{\circ }\mathrm{N}$, the re-entry trajectory inclination as $i_r={45}^{\circ }$, the re-entry velocity as $v_{re}=10.7\mathrm{km}/\mathrm{s}$, and the re-entry voyage as $S_r=6456\mathrm{km}$, that is $U_r={58}^{\circ }$.

Execute the procedures of module M1 in Figure 2 to obtain the re-entry state of the probe in ECEF as ${\mathit{r}}_e=[4314.9,4783.6,851.5]\mathrm{km},{\mathit{v}}_e=[-7.033,3.535,7.248]\mathrm{km}/\mathrm{s}$.

Select a re-entry date of 3 October 2030, for optimal Moon-to-Earth transfer trajectory design. Set the transfer duration to $\Delta T=3.0$ days, and the re-entry date range is $\left[D_1,D_1+1\right]$ = [2,462,777.5, 2,462,778.5] (3 October 2030 00:00:00–4 October 2030 00:00:00).

The procedures of module M2 in Figure 2 are implemented to obtain the re-entry time corresponding to the optimal trajectory as JD = 2,462,778.43474 (3 October 2030 22:26:01).

In this case, the re-entry velocity with the transfer duration set to the preset value is adjusted to ${\mathit{v}}_{r\epsilon }=10.6541\mathrm{km}/\mathrm{s}$, and the re-entry state of the probe in ECI is ${\mathit{r}}_i$ = [5136.5, 3888.1, 851.5] $\mathrm{km},$${\mathit{v}}_i=[-6.501,5.147,7.217]\mathrm{km}/\mathrm{s}$.

The optimal transfer trajectory is shown in Figure 4, and the corresponding radius of the perilune is $2768.5\mathrm{km}$.

Analyze the opportunity of the Moon-to-Earth transfer in 2030. Set the re-entry angle as $-6^{\circ }$, and calculate the daily optimal transfer trajectory in 2030, and the corresponding perilune altitude change is shown in Figure 5.

Set the altitude deviation limit of the perilune as $\chi$ = 50,000 km. The red dotted line in Figure 5 corresponds to $\chi$; that is, the date when the perilune altitude of the optimal transfer trajectory is lower than this deviation limit is included in the Moon-to-Earth transfer opportunity. It can be seen that the transfer opportunity is about 6–8 days per month, and the opportunity width is about 74 days in 2030.

Similarly, select a re-entry date of 3 October 2030, for precise Moon-to-Earth transfer trajectory design. Choose the inclination and altitude constraints of the perilune as $i_m^{\prime }={85}^{\circ }$, $r_{pm}^{\prime }=200\mathrm{km}$, respectively, and the radius of the SOI of the Moon is selected as $r_{\mathrm{SOI}}$ = 66,200 km.

Adopt the procedures of module M4 in Figure 3 to correct the semi-major axis, argument of ascending node, and argument of perihelion corresponding to this state. It is obtained that the corrected state at the boundary of the SOI of the Moon is ${\mathit{r}}_{\mathrm{SOI}}$ = [−186,328.5, −244,850.2, −115,443.3] $\mathrm{km},{\mathit{v}}_{\mathrm{SOI}}=[0.5633,0.5225,0.1425]\mathrm{km}/\mathrm{s}$.

Table 1. Orbital elements at perilune of precise Moon-to-Earth transfer trajectory.

Orbital Elements	a/km	e	$\mathit{i}{/}^{\circ }$	$\mathbf{\Omega }{/}^{\circ }$	$\mathit{\omega }{/}^{\circ }$	$\mathit{f}{/}^{\circ }$
Value	−4808.3	1.4030	85.0	291.0	55.3	0

shows the orbital elements at the perilune of the Moon-to-Earth transfer trajectory in the MCI.

If the probe departs from a $200\mathrm{km}$ altitude circular orbit coplanar with the Moon-to-Earth transfer trajectory, the applied departure impulse is $\Delta {\mathit{v}}_M=[-217.8,687.1,496.2]\mathrm{m}/\mathrm{s}$.

Finally, corrections are applied to the second segment of the trajectory, so that the three segments are precisely connected. Perform high-precision inverse integration on the re-entry state ${\mathit{X}}_i$ obtained above, and obtain the state of the probe at the day before re-entry, ${\mathit{r}}_{rr}$ = [−96,970.8, −157,624.1, −85,389.1] $\mathrm{km},{\mathit{v}}_{rr}=[0.9457,1.0081,0.3950]\mathrm{km}/\mathrm{s}$.

Apply the differential correction of module M5 in Figure 3 to correct the velocity ${\mathit{v}}_{rr}$, so that the position at the boundary of the SOI of the Moon aligns with the position ${\mathit{r}}_{\mathrm{SOI}}$ obtained in the previous step. The velocity increment on the day before re-entry is $\Delta {\mathit{v}}_{rr}=[-5.2,-1.2,5.8]\mathrm{m}/\mathrm{s}$.

Thus, the second and third segments are precisely aligned at the boundary of the SOI of the Moon, and the required velocity increment is $\Delta {\mathit{v}}_{\mathrm{SOI}}=[4.1,1.6,-2.7]\mathrm{m}/\mathrm{s}$.

Combining all the above solution processes, the total velocity increment required to obtain the Moon-to-Earth transfer trajectory is $∥\Delta {\mathit{v}}_{all}∥=∥\Delta {\mathit{v}}_M∥+∥\Delta {\mathit{v}}_{\mathrm{SOI}}∥+∥\Delta {\mathit{v}}_{rr}∥=888.2\mathrm{m}/\mathrm{s}$.

The precise Moon-to-Earth transfer trajectory satisfying all preset constraints is shown in Figure 6.

In order to characterize some regularities of the precise transfer trajectory, Figure 7 shows the changes in the right ascension of the ascending node (RAAN) and the total velocity increment required $(\Delta V)$ for the transfer to reach the perilune in the first three months of 2030.

The RAAN and $\Delta V$ of the precise trajectories in the first three months of 2030 change in the same cycle as the perilune altitude in the opportunity screening process in Figure 5. In other words, the RAAN of the perilune of the precise trajectory based on natural transfer changes one circle synchronized with the Moon’s revolution. The $\Delta V$ required is strictly related to the perilune altitude of the optimal transfer per day, and the difference in velocity increments is mainly caused by the difference in maneuvers required to correct the perilune altitude.

In addition,

Table 2. Comparison of the obtained opportunity between this paper and the reference.

Source	$\mathit{\gamma }$	$\mathbf{\Delta }\mathit{T}$	Opportunity	Minimum $\mathbf{\Delta }\mathit{V}$
Paper	$-6^{\circ }$	3 days	1 January–4 January, 25 January–31 January	889.5 m/s
Reference [45]	$-6^{\circ }$	3 days	1 January–3 January, 26 January–31 January	890.1 m/s

represents the obtained opportunity in the first month of 2019 to compare with reference [45].

It can be seen that the opportunity obtained in this paper is approximately the same but slightly wider than that in the reference when the perilune altitude deviation limit is $\chi$ = 50,000 km. The corresponding minimum $\Delta V$ is almost equal to the reference, and the maximum $\Delta V$ required is about 1.38 km/s.

5.2. Design Parameters Analysis

There are many design parameters involved in the design process of the opportunity, such as the transfer duration, altitude of the perilune, and so on, while the location of the landing field also affects the design results. Different design parameters have distinct degrees of influence on the distribution of opportunities, the cost of transfer, etc. The analysis of the influence of several design parameters on the distribution law of opportunity and the performance of precise Moon-to-Earth transfer trajectory is as follows.

5.2.1. Position of Landing Site

On the premise that other design parameters remain unchanged, the latitude and longitude of the landing field affect the position of the re-entry point according to the geometric relationship of the spherical triangle. Figure 8 reveals the moment distribution in a day based on the single transfer trajectory example when the longitude of the landing field varies in full phase space, while Figure 9 shows the perilune altitudes of several latitudes and the distribution of opportunities when the latitude of the landing field changes within a certain range in the first three months of 2030.

It can be seen from Figure 8 that when the longitude of the landing field changes in the whole space, the daily optimal transfer moment changes with the same trend of longitude. This is easy to understand, for the direction of the velocity of the re-entry point in ECEF will change with the longitude of the landing field, while the angle of the Moon’s revolution changes very little within a day, so the date distribution of opportunities has little changes. On the same day, the directions of the optimal transfer trajectories corresponding to the landing field with different longitudes in the ECI should basically be the same. That is, when the longitude traverses a circle uniformly, the optimal transfer time in a day also changes uniformly within 24 h.

In contrast, the latitude of the landing field has a greater impact on the distribution of opportunities, which is shown in Figure 9. The velocity directions of the re-entry points in ECEF corresponding to different latitudes are different, especially in the z-direction, which leads to the fact that when the latitude changes on a feasible date, the height of the perilune, especially the distance in the z-direction, exceeds the deviation limit. Thus, this opportunity no longer exists. However, correspondingly, when the angle of the Moon can compensate for the z-direction deviation caused by latitude changes on other dates, this date is a new transfer opportunity. Thus, the difference of the latitudes of the landing field leads to the diversification of the transfer opportunity. When the re-entry inclination is set to ${45}^{\circ }$, the opportunity lengths in the first three months of 2030 are 18 days, 35 days, 18 days, 17 days, and 34 days corresponding to the latitudes of $0^{\circ },{10}^{\circ },{20}^{\circ },{30}^{\circ }$, ${40}^{\circ }$, and ${50}^{\circ }$, respectively. It should be noted that the re-entry inclination cannot be smaller than the latitude of the landing field.

5.2.2. Transfer Duration

In this paper, the method of numerical iteration is adopted to adjust the re-entry velocity, so that the duration of the Moon-to-Earth transfer can precisely meet the preset value. At the same time, the shorter the transfer time, the larger the velocity of the probe during the transfer. Given the fact that the magnitude of the relative distances and velocity of the Earth and the Moon vary little, a larger transfer velocity results in greater Moon capture maneuver. Figure 10 shows the change curves of the transfer duration and the total velocity increment $\Delta V$, as well as the atmospheric re-entry velocity. Design parameter settings are the same as in the single transfer trajectory example except for the transfer duration.

5.2.3. Altitude of Perilune

The setting of the values of the nominal altitude of the perilune is only involved in the design of the precise transfer trajectory, so the parameter does not affect the distribution of transfer opportunities. However, the altitude of the perilune closely related to the velocity increments is required to precisely transfer trajectories. Figure 11 reveals the relationship between the velocity increment and the perilune altitude. Design parameter settings are the same as in the single transfer trajectory example except for the altitude of perilune.

The curve in Figure 11 illustrates a negative correlation between the velocity increment and the perilune altitude, which can be explained by the energy difference between orbits. The higher the perilune, the larger the energy of capture orbit. When the probe enters the lunar capture orbit from the higher-energy Moon-to-Earth transfer trajectory, a higher capture orbit energy means a smaller capture velocity increment.

6. Conclusions

Considering Earth re-entry and landing field constraints, an integrated design method for Moon-to-Earth transfer opportunity and precise transfer trajectory is developed. First, this paper innovatively proposes an analytical method to obtain the re-entry state that satisfies reliable ground landing by adopting the theory of spherical geometry. Next, the Moon-to-Earth transfer opportunity is effectively determined by the permitted range of the perilune altitude according to the inverse integration of the obtained re-entry state. Then, the three-impulse Moon-to-Earth transfer strategy is adopted to satisfy the Moon departure condition from any initial inclination in the obtained opportunity. The simulation shows that there is little difference in the width of the Moon-to-Earth transfer opportunity in each month under the fuel constraint. When the deviation limit is 50,000 km, the width of the transfer opportunity in a year is about 74 days. Compared with the traditional patched-conic method, the application of the Earth and Moon gravitational model and the high-precision dynamic model can ensure higher reliability and accuracy of the resulting transfer opportunities. Moreover, the longitude of the landing field is relevant to the optimal transfer moment in a day, while the latitude influences the date distribution of the opportunity. The transfer duration and the altitude of the perilune are both closely related to the fuel consumption of the mission. The proposed design method can provide an effective Moon-to-Earth transfer opportunity and trajectory determination reference for the next stage of the Moon missions of China and other countries.