Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
The ASTRI Mini-Array: A New Pathfinder for Imaging Cherenkov Telescope Arrays
Previous Article in Journal
Origins and Natures of Inflation, Dark Matter and Dark Energy
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multiple SSO Space Debris Flyby Trajectory Design Based on Cislunar Orbit

1
Technology and Engineering Center for Space Utilization, Chinese Academy of Sciences, Beijing 100094, China
2
University of Chinese Academy of Sciences, Beijing 101408, China
*
Author to whom correspondence should be addressed.
Universe 2024, 10(3), 145; https://doi.org/10.3390/universe10030145
Submission received: 12 January 2024 / Revised: 1 March 2024 / Accepted: 4 March 2024 / Published: 16 March 2024

Abstract

:
This paper investigates the trajectory design problem in the scenario of a multiple Sun-synchronous Orbit (SSO) space debris flyby mission from a DRO space station. At first, the characteristics of non-planar transfer from DRO to SSO in the Earth–Moon system are analyzed. The methods of large-scale ergodicity and pruning are utilized to investigate single-impulse and two-impulse DRO–Earth transfers. Using a powered lunar flyby, the two-impulse DRO–Earth transfer is able to fly by SSO debris while satisfying the requirements of the mission. After the local optimization, the optimal result of two-impulse DRO–Earth transfer and flyby is obtained. A multi-objective evolutionary algorithm is used to design the Pareto-optimal trajectories of multiple flybys. The semi-analytical optimization method is developed to provide the estimations of the transfer parameters in order to reduce the computations caused by the evolutionary algorithm. Simulations show that transferring from the 3:2 resonant DRO to a near-coplanar flyby of a SSO target debris using a powered lunar gravity assist needs a 0.47 km/s velocity increment. The mission’s total velocity increment is 1.39 km/s, and the total transfer time is 2.23 years.

1. Introduction

A distant retrograde orbit (DRO) is a special kind of plane-symmetric orbit in the circular restricted three-body problem. A DRO in the Earth–Moon system is a geocentric orbit in the cislunar space which is subject to the gravity of the Earth and the Moon. It has high stability and good coverage of both the Earth and the Moon. The application of DRO has gradually become a research hotspot. In 2016, NASA’s Asteroid Redirect Mission planned to tow tons of asteroids to DRO for mooring [1]. In November 2022, NASA launched an unmanned Orion spacecraft to the Moon as a part of the Artemis I mission, and it is planned to fly to DRO through a lunar gravity assist [2]. DRO is well suited for use for a space station in deep space exploration missions and near-Earth observation missions, and it has important research value for on-orbit service and maintenance, Earth coverage, and observation.
In 1968, Broucke revealed the existence of DROs in the Earth–Moon system [3]. The DRO dynamics of three-body systems such as the Earth–Moon system, the Sun–Earth system, and the Ganymede system have been studied [4,5,6]. Existing studies of DRO transfers in the cislunar space focus on how to transfer from Earth or other locations in the cislunar space to DRO [7,8,9,10,11,12,13]. The transfer orbit from DRO to the vicinity of the Earth can be seen as a negative integration from the vicinity of the Earth to DRO [14,15,16,17]. These studies show that the spacecraft requires only small velocity increments to transfer from DRO to the vicinity of the Moon and then to other types of near-Earth orbits such as Low Earth Orbit (LEO) and Geosynchronous Earth Orbit (GEO) via lunar gravity assist, which has significant fuel-saving advantages for space debris removal missions.
Space debris has seriously affected the safety of space missions. Between the years 2000 and 2013, only 12% of spacecraft performed reentry or deorbit maneuvers at their end-of-life, and only 10% of spacecraft successfully performed 25-year deorbit maneuvers [18]. NASA research shows that even if 90% of constellations successfully perform 25-year deorbit maneuvers, there will still be 290% more space debris in the next 200 years than if there were no giant constellations [19]. Even if no new man-made objects were put in space, the space debris in the LEO region would still grow exponentially [20].
Active removal is a direct approach to reduce the amount of space debris. Space robotic arm capture [21] and flying net capture [22] have been the subjects of long-term research in various space agencies. The above “one-to-one” methods need complex techniques for grappling an uncooperative object, which is costly and inefficient.
Apart from the deorbit removal methods for spacecraft at the end of mission life, i.e., drag sail [23], inflate-sail [24], and electrodynamic tether [25], contactless active removal methods such as Coulomb thrusting [26], laser beam [27], artificial atmosphere [28], and expanding foam [29] have been proposed in the application of space debris mitigation. Coulomb thrusting removal is applicable only when the Debye shielding is weak [30,31]. Laser beam removal can only remove very small debris [32,33]. The method of artificial atmosphere is a technique in which the spacecraft releases an artificial atmosphere in a region near the path of the space debris to decelerate the debris and accelerate natural orbital decay [28,34,35]. The approach of expanding foam is similar to the artificial atmosphere method, except that the material released by the spacecraft is expanding foam. The artificial atmosphere and expanding foam removal methods are suitable for removing LEO debris. Moreover, the expanding foam removal method can deorbit large LEO debris (1.4–3.5 tons per year) [29,36]. Both of these two methods require a spacecraft to arrive in the vicinity of the debris and release the deorbit material.
Regarding the analysis of space debris removal mission, one author [37] designs a space mission in which, when the spacecraft rendezvous with a piece of debris, a deorbit device is attached to the object by a robotic arm and the device is activated, and then the spacecraft travels to the next piece of debris. However, the results show that this removal mission would take more than 30 km/s of total impulse to deorbit 35 large objects in 7 years. Another author [38] considers the removal mission of two Cosmos-3M second stages. The spacecraft performs maneuvers to approach the object, and then the spacecraft captures the object and attaches the thruster deorbit kit. A further author [39] considers a mission to deorbit five large pieces of debris each year. They suggest that the object in the orbit close to coplanar should be deorbited first because changing the orbital plane is costly. One author [40] assesses the effectiveness of three heuristic algorithms in the optimal design of the active debris removal mission by comparing the solutions and running times. Another author [41] designs an on-orbit serving mission to deorbit a non-cooperative object by using space robotic systems to manipulate and deorbit the object. These studies require spacecraft rendezvous with debris to allow the robotic arm or other removal kit to perform deorbit, which requires large velocity and faces complex technical issues.
Inspired by the contactless removal methods of artificial atmosphere and expanding foam, this paper proposes a concept that uses one spacecraft to fly by much Chinese space debris in SSO orbit. When the spacecraft flies by one piece of space debris, it can release a gaseous cloud to deorbit the space debris. This “one-to-many” method has high efficiency and low cost, and is easier to realize in space debris active removal missions. This paper focuses on the multiple SSO space debris flyby trajectory design, and the detail of the deorbit by the gaseous cloud is beyond the scope of this paper.
This paper is organized as follows. Section 2 presents the problem statement and the equations of motion. Section 3 selects the target debris in SSO. Section 4 designs the orbit by which the spacecraft departs from DRO and performs a near-coplanar flyby of a SSO target piece of debris using two impulses. Section 5 uses the multi-objective evolutionary optimization algorithm and a semi-analytical optimization method to design the multiple highly elliptical orbits. Section 6 presents the design results of the flyby mission. Section 7 concludes this paper.

2. Problem Statement and Equations of Motion

2.1. General Scenario

Envisioning the active removal of space debris, this paper investigates the mission design problem of a multiple-object flyby such that the detailed information of the object can be obtained by the on-board inspection payload. Debris in SSO are chosen as the target database because SSO is one of the busiest orbits around the Earth.
Figure 1 shows the mission scenario. Initially, the spacecraft is in the parking orbit. At point P1, the spacecraft undertakes the departure maneuver to leave the parking orbit. At point P2, the spacecraft flies by one piece of target debris. After the first flyby, the spacecraft flies to point P3. When the spacecraft reaches P3, it applies the second maneuver to adjust its orbit. As a result, the spacecraft flies by another piece of target debris at point P4. In this way, the spacecraft performs a multiple debris flyby mission through several maneuvers. Points P1 and P3 are calculated and optimized in Section 4 and Section 5.
The motion of the spacecraft is in cislunar space and described in the circular restricted three-body system. The orbit of the spacecraft is presented in the Earth–Moon rotating coordinate, i.e., the center of mass rotating frame { C x ^ , y ^ , z ^ } , where C is the center of mass of the Earth–Moon system, the x ^ y ^ plane is the orbital plane of the Moon, x ^ points from the Earth to the Moon, and z ^ is the normal of the x ^ y ^ plane. The spacecraft’s equation of motion in the Earth–Moon rotating coordinate is as follows [42]:
r ¨ + 2 [ y ˙   x ˙   0 ] = ( W r ) T
where r = [ x , y , z ] T is the spacecraft’s position in the Earth–Moon system. r ˙ = [ x ˙ , y ˙ , z ˙ ] T is the spacecraft’s velocity and r ¨ is the spacecraft’s acceleration.
W = 1 2 ( x 2 + y 2 ) + U ( r 1 , r 2 ) U ( r 1 , r 2 ) = ( 1 μ ) / r 1 + μ / r 2 μ = M 2 M 1 + M 2 J t = 2 W | | v | | 2
where r 1 and r 2 is the distance from the spacecraft to the Earth and the Moon, respectively. J t is the Jacobi constant. M1 is the mass of the Earth and M2 is mass of the Moon.
r 1 = ( x + μ ) 2 + y 2 + z 2 r 2 = ( x + μ 1 ) 2 + y 2 + z 2
In the circular restricted three-body problem (CRTBP), the relevant common parameters are as follows [43]:
μ = 0.0121506683 TU = 104.8585095   h DU = 388405   km VU = 3665.892272   km / h
where TU, DU, and VU are the normalized units of time, distance, and the modulus of velocity in the Earth–Moon rotating coordinate. In addition, the ephemeris model should be used if the eccentric orbit of the Moon and the perturbations due to the gravity of the Sun and other celestial bodies are considered. The orbits in CR3BP cannot be directly used in the ephemeris model, since the CRTBP is a low-fidelity dynamical model. It is necessary to correct and optimize the orbits obtained from CRTBP in the ephemeris model [44]. This could be analyzed in future research.
In this paper, the SSO space debris is affected by the J 2 perturbation. The orbit of the SSO space debris is presented in the J2000 ECI coordinate { O X ^ , Y ^ , Z ^ } , where O is the center of the mass of the Earth, the X ^ Y ^ plane is the equatorial plane, X ^ points from O to the location of the equinox in 2000 years, and Z ^ is the normal of the X ^ Y ^ plane. The equation of motion of the SSO space debris in the J2000 ECI coordinate is as follows:
a d = μ e r d r d + a J 2
where a d is the acceleration vector of space debris. r d is the inertial position vector of space debris. r d = | | r d | | . a J 2 is the acceleration vector caused by the J 2 perturbation. μ e is Earth’s gravity constant.

2.2. Parking Orbit

In designing the multiple-object flyby trajectory, several prerequisites should be clarified. First of all, let us determine the type of the parking orbit in this multiple-SSO-object flyby mission. This paper identifies the following requirements in selecting the parking orbit:
(1)
The parking orbit should be long-term stable.
(2)
A small velocity increment is demanded for the spacecraft to transfer from the parking orbit to flyby multiple SSO targets.
(3)
The apogee of the parking orbit should be large so that it would be easier for the spacecraft to adjust the orbital inclination, which is more favorable for multiple-object near-coplanar flyby missions.
Within the cislunar space, typical potential candidates include Geosynchronous Earth Orbit (GEO), Distance Retrograde Orbit (DRO), Polar Lunar Orbit (PLO), Near Rectilinear Halo Orbit (NRHO), Lagrange point 1 (L1) halo orbit, Lagrange point 2 (L2) halo orbit, and so on.
Table 1 compares the velocity increments for transferring from LEO to different target orbits. Note that without considering the perturbations, the velocity increment transferring from LEO to a designated target orbit is nearly identical to that transferring back. As shown in Table 1, the minimum velocity increment for transferring from LEO to the 2:1 DRO using lunar gravity assist and weak stability boundary is 3.1934 km/s. This is a rather small number compared to the velocity increment budgets transferring to GEO, PLO, NRHO, L1 halo orbit, and L2 halo orbit. Therefore, this paper chooses DRO as the parking orbit for the spacecraft.
This section does not specify specific parameters of DRO in this multiple-debris flyby mission. This problem will be discussed in the following sections.

2.3. Flyby Mission Description

The mission scenario is shown in Figure 2. The steps of the flyby mission are illustrated in Table 2 and Figure 3. As shown in Figure 2, the spacecraft is located in DRO at the initial moment τ ¯ . At τ 0 , by imposing ΔV0, the spacecraft leaves DRO and approaches the Moon. At τ 1 , the spacecraft applies ΔV1 in lunar sphere of influence (SOI) to take a powered lunar gravity assist. Then, the spacecraft enters the Moon–Earth transfer orbit. At τ 2 , the spacecraft near-coplanar flies by a piece of space debris near the perigee. After that, the spacecraft flies towards the apogee and applies an impulse near the apogee. By adjusting the orbit and the phase, the spacecraft can near-coplanar fly by the next piece of space debris in the highly elliptical orbit (HEO). By iterating Step 4–Step 5 in Table 2, the multiple-object flyby mission is accomplished.
This paper divides the multiple-object near-coplanar flyby mission into two phases. As shown in Figure 2b, Phase 1 is the segment in which the spacecraft transfers from DRO to Earth and implements the first flyby, i.e., Step 1–Step 3 in Table 2. Phase 2 is the segment in which the spacecraft near-coplanar flies by all the target debris, i.e., Step 4–Step 5 in Table 2.
The requirements of the flyby mission are shown in Table 3. In Table 3, ζ 1 = 31   days , ζ 2 = 0.4   km / s , ζ 3 = 0.2   km / s , ζ 4 = 0.1   km / s , ζ 5 = 50   km , ζ 6 = 100   km , ζ 7 = 4   km / s , h ˜ 1 = 100   km , and h ˜ 2 = 400   km . Note that the errors of measurement and controlling are not considered. The two impulses of Phase 1 are applied on DRO (ΔV0) and the perilune (ΔV1), respectively. The impulse ΔV of Phase 2 is applied in HEO. h p m o o n h ˜ 1 = 100   km means that the PLF cannot impact the Moon. h p e a r t h h ˜ 2 = 400   km means that the spacecraft cannot enter the Earth’s atmosphere. When the spacecraft with an apogee radius close to the Moon coplanar flies by a space debris with an orbital height of 800 km, the relative velocity is approximately 3 km/s. Thus, this paper requires the modulus of the relative velocity when the spacecraft flies by the target debris to be less than or equal to 4 km/s.

3. Target Selection of SSO Space Debris

3.1. Analysis of SSO Space Debris

SSO has great and scientific research value. Polar-orbiting meteorological satellites and imaging satellites are often located in SSO. The main characteristics of SSO include the fact that the angle between the normal orbital plane and the projection of the solar direction on the equatorial plane remains constant. Due to the J2 perturbation, the mean rate of change of the right ascending node is as follows:
Ω ¯ ˙ = 3 n J 2 R e 2 2 a 2 ( 1 e 2 ) 2 cos i
where n = μ e / a 3 , a is the semi-major axis, e is the eccentricity, and i is the inclination. J 2 = 0.001082 . The rate of change of the ascending node of SSO is the same as the Earth’s mean motion around the Sun, which is as follows:
Ω ¯ ˙ = 360 365.24 days = 1.991 × 10 7   rad / s
Substituting Equation (7) into Equation (6) yields the following:
1.991 × 10 7   rad / s + 3 n J 2 R e 2 2 a 2 ( 1 e 2 ) 2 cos i = 0
Because of the computational error, this paper defines the orbit satisfying the following equation as SSO:
| 1.991 × 10 7   rad / s + 3 n J 2 R e 2 2 a 2 ( 1 e 2 ) 2 cos i | 0.001
Selecting objects from the Space Track Database according to Equation (9), it can be seen that the region with perigee altitude of [600, 900] km, apogee altitude of [650, 1050] km, and orbital inclination of [96.3, 99.3]° is the most densely distributed region of Chinese SSO space debris. The space debris in this region have mainly come from the Fengyun satellites and Changzheng rockets.
Reference [49] proposes that states and international intergovernmental organizations should remove spacecraft and launch vehicle orbital stages at the end of the mission life, and also encourage the development and use of relevant technologies for the measurement, monitoring, and characterization of the orbital and physical properties of space debris. All countries should uphold the concept of building a community of human destiny, share the responsibility of maintaining outer space security, and address threats to all aspects of outer space security.
For the purpose of maintaining the long-term sustainable peaceful use of outer space, and for better measurement and detection of space debris, this paper offers a solution to obtain detailed information on our space debris so that active removal can be better implemented.

3.2. Selection of SSO Target Debris

Next, let us build up the target debris database for this multiple-SSO-object flyby mission. Based on the distribution characteristics of the orbital location, type, source, and size of SSO space cataloged objects, this paper selects space cataloged objects that satisfy the following conditions to construct the target debris collection for the flyby mission:
(1)
The types of target debris are debris or rocket body.
(2)
The perigee altitude has not yet begun to decrease.
(3)
The target debris is generated by Changzheng rockets.
(4)
The radar cross-section area of the target is as follows: S min > 1   m 2 .
(5)
The mass of the target is as follows: m > 1   kg .
(6)
h p [ 600 , 900 ]   km , h a [ 650 , 1050 ]   km , and i [ 96.3 , 99.3 ] .
The second condition means that the orbit of space debris remains stable. In the fourth condition, S min > 1   m 2 means that the target debris is large space debris. The last condition means that the target debris is in the area with the highest density of debris distribution.
A total of 22 pieces of debris were selected that satisfied the above conditions. Their orbits in the J2000 ECI coordinate are shown in Figure 4. The physical and orbital parameters of these target debris are shown in Table A1 and Table A2 in Appendix A.

4. Phase 1: Two Impulses DRO–Earth Transfer Orbit Design

As shown in Figure 5a, this section discusses how to use a two-impulse maneuver [7] to transfer from DRO to the Earth and fly by a target debris near the perigee. The objective of this section is to find an orbit that minimizes the total impulse while satisfying the flyby requirements. As a rule of thumb, in the typical Earth–Moon three-body problem, a key factor in this transfer orbit design problem is to fully utilize the lunar gravity assist.
This section proposes the spacecraft’s flight procedure as follows: at τ 0 , the spacecraft applies a departure impulse Δ V 0 in DRO to transfer to the Moon. At τ 1 , the spacecraft applies the impulse Δ V 1 when it arrives at the perilune. At τ 2 , the spacecraft flies to the Earth and flies by the SSO target debris near the perigee. The transfer time of the spacecraft from DRO to the perilune is Δ t 0 = τ 1 τ 0 . The transfer time from the perilune to flyby is Δ t 1 = τ 2 τ 1 .
Typical DRO–Earth transfer methods include the internal transfer, which is mainly accomplished based on the CRTBP system, and the external transfer, in which the spacecraft flies to the Sun–Earth gravitational equilibrium region and achieves a low-energy Moon–Earth transfer using the weak stability boundary theory. Note that the external transfer method requires a much longer transfer time (100 days to 130 days), which does not conform to the requirements shown in Table 3. So, this paper chooses the internal transfer approach.
Taking the perilune as a node, the complete DRO–Earth transfer and flyby orbit is composed of two segments: (1) The segment from DRO to the perilune, denoted by G 0 , as shown in Figure 5a. (2) The segment from the perilune to the perigee, denoted by G 1 . The two maneuvers are denoted as Γ 0 ( X 0 ) and Γ 1 ( X 1 ) . The notation { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } means the trajectory of the spacecraft is composed of G 0 and G 1 , with two orbital maneuvers Γ 0 ( X 0 ) and Γ 1 ( X 1 ) . The notation { G 0 ( Γ 0 ) + G 1 ( 0 ) } means the trajectory of the spacecraft is composed of G 0 and G 1 , with only one maneuver Γ 0 ( X 0 ) .
The plan is constructed as follows: First, the DRO–Earth transfer orbit with a single impulse { G 0 ( Γ 0 ) + G 1 ( 0 ) } is to be investigated. Then, based on the results of { G 0 ( Γ 0 ) + G 1 ( 0 ) } , the DRO–Earth transfer orbit with two impulses { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } is designed. The parameters of { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } are to be optimized to solve the local optimization problem of the first flyby.

4.1. Single-Impulse DRO–Earth Transfer

This subsection intends to investigate the feasibility of using a single-impulse strategy to perform the SSO object flyby mission. The DRO–Earth transfer orbit with a single impulse is denoted as { G 0 ( Γ 0 ) + G 1 ( 0 ) } . The spacecraft applies Γ 0 ( X 0 ) in DRO and flies to the perilune. Then, the spacecraft takes the thrustless lunar flyby, i.e., Δ V 1 = 0 , and flies towards the Earth. The requirements of the single-impulse DRO–Earth transfer orbit are as follows:
Δ τ ζ 1 ,   Δ t 0 ζ 1 | | Δ V 0 | | ζ 2 ,   i 1 i ˜ 1 h ˜ 1 h p m o o n h S O I h ˜ 2 h p e a r t h h ˜ 3
where ζ 1 , ζ 2 , h ˜ 1 , and h ˜ 2 are identical to those of Table 3, h S O I is the height of the lunar SOI, and the radius of lunar SOI is set to 61,525 km. h ˜ 3 = 2000   km . i 1 is the orbital inclination of the spacecraft when it arrives at the perigee and i ˜ 1 = 80 .
In order to investigate the feasibility of the single-impulse approach, this section studies whether the requirements in Equation (10) can be satisfied under different parameters determining Γ 0 ( X 0 ) . The parameters determining Γ 0 ( X 0 ) in the CRTBP system are as follows:
X 0 = [ Δ τ , | | Δ V 0 | | , α 0 , β 0 ]
As shown in Figure 5b, α 0 is the angle between Δ V 0 and the projection of Δ V 0 in the x ^ y ^ plane, and β 0 is the angle between the projection of Δ V 0 in the x ^ y ^ plane and x ^ . The equation of Δ V 0 with respect to ( α 0 , β 0 ) is shown in Equation (12).
Δ V 0 = | | Δ V 0 | | ( cos α 0 cos β 0 , cos α 0 sin β 0 , sin α 0 ) T
The calculation procedure includes large-scale ergodicity according to Table 4 and verification of the satisfaction of Equation (10) for each Γ 0 ( X 0 ) , followed by pruning. In Table 4, this paper chooses the stepsizes as δ 1 = 1   day , δ 2 = 0.005   km / s , and δ 3 = 10 . Let us traverse the parameters according to Table 4 and calculate the values of the following variables in { G 0 ( Γ 0 ) + G 1 ( 0 ) } for each set of parameters: Δ t 0 , the perilune height h p m o o n , Jacobi constant J t for G 1 after one thrustless lunar gravity assist, and the perigee height h p e a r t h . From Equation (10), we can see that the maximum flight time of the spacecraft from the perilune to the perigee is 31 days. Since there may be more than one perigee, G 1 is the orbit of the spacecraft from the perilune to the first perigee that conforms to the pruning condition of h p e a r t h . The pruning conditions are shown in Equation (13).
J t C 1 h ˜ 1 h p m o o n h S O I h ˜ 2 h p e a r t h 36,000   km
where C 1 = 3.024150064 is the Jacobi constant of point L 1 in CRTBP. Equation (13) includes three inequalities. The first one means that the cislunar space is connected. The second one means that the spacecraft could use a lunar gravity assist without colliding with the Moon. The third inequality means that the spacecraft could reach the Earth without entering the Earth’s atmosphere.
Note that, in this subsection, the second maneuver Δ V 1 is assumed to be zero. After a thrustless lunar gravity assist, the resulting height of the perigee in the ergodicity varies in a large range. Therefore, the upper limit of h p e a r t h in the third condition is set to be 36,000 km, which is a rather large value. In addition, because the i 1 of all { G 0 ( Γ 0 ) + G 1 ( 0 ) } are less than i ˜ 1 , the requirement of i 1 in Table 4 is not added to the pruning conditions.
After traversing all parameters shown in Table 4, the h p m o o n , J t , and h p e a r t h of each { G 0 ( Γ 0 ) + G 1 ( 0 ) } are obtained. The traverse results are pruned according to Equation (13). There are 2872 pairs of { G 0 ( Γ 0 ) + G 1 ( 0 ) } after pruning.
Let S 0 denote the group composed of all orbits { G 0 ( Γ 0 ) + G 1 ( 0 ) } after being pruned by Equation (12). [ a * , e * , i * , Ω * , ω * , M * ] are the orbital elements of the spacecraft when it exits the Lunar SOI. Figure 6 shows the distribution of α 0 and β 0 about i * in set S 0 . Each data point represents a { G 0 ( Γ 0 ) + G 1 ( 0 ) } in S 0 . In Figure 6a, the orbit of the data indicated by the red box and blue box are shown in Table 5 and Figure 7a–d. It can be seen that { G 0 ( Γ 0 ) + G 1 ( 0 ) } is symmetrically distributed about α 0 = 90 and α 0 = 270 when a * [ 0 , 180 ] and α 0 [ 180 , 360 ] , respectively. In Figure 6b, the parameters of the data indicated by the red box and blue box are shown in Table 5. It can be known that the period of variation of { G 0 ( Γ 0 ) + G 1 ( 0 ) } about β 0 is 180°.
Figure 6 also shows that the orbital inclination of the spacecraft can be increased from 20.98° to 62.70° by the thrustless lunar flyby. The orbit when i * = 62.70 is shown in Figure 7e,f. The parameters of Γ 0 ( X 0 ) when i * = 62.70 are Δ τ = 31   days , | | Δ V 0 | | = 285   m / s , and ( α 0 , β 0 ) = ( 10 , 210 ) . The orbital elements when the spacecraft enters and exits the lunar SOI are shown in Table 6.
From Figure 6, we can see that lunar gravity can reduce the orbit altitude and raise the orbital inclination of the spacecraft.
Figure 7g,h show the distribution of i 1 and h p e a r t h about | | Δ V 0 | | . Figure 7g shows that orbital inclinations of all { G 0 ( Γ 0 ) + G 1 ( 0 ) } in S 0 are i 1 < 80 . Orbits of spacecraft are symmetrically distributed in the plane of | | Δ V 0 | | i 1 about i 1 = 24.88°, and 24.88° is the orbit inclination of the Earth–Moon plane when the spacecraft departs from DRO. This is because of the fact that i 1 is related to the departure impulse. When α 0 = 0 , the spacecraft’s orbit is in the Earth–Moon plane both before and after the thrustless lunar gravity assist.
Figure 7h shows that the minimum height of the perigee is 3559 km, which does not conform to the orbital height requirement of the flyby. It can be seen from Figure 7g,h that it is impossible for the DRO spacecraft to fly by the SSO target by applying a single impulse while satisfying the requirements shown in Equation (10). Thus, it is necessary to investigate strategies with at least two impulses.

4.2. Two-Impulse DRO–Earth Transfer

This subsection investigates the two-impulse strategy { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } . The spacecraft applies Γ 0 ( X 0 ) in DRO and flies to the perilune. Then, the spacecraft takes a PLF Γ 1 ( X 1 ) and flies to the Earth. The requirements of the two-impulse DRO–Earth transfer orbit are as follows:
Δ τ ζ 1 ,   Δ t 0 ζ 1 ,   | | Δ V 0 | | ζ 2 | | Δ V 1 | | ζ 3 ,   Δ t 1 ζ 1 ,   i 1 i ˜ 1 h ˜ 1 h p m o o n h S O I ,   h ˜ 2 h p e a r t h h ˜ 3
In order to investigate the feasibility of the two-impulse approach, this section studies whether the requirements in Equation (14) can be satisfied. The calculation procedures are composed of two-dimensional ergodicity: traverse Γ 0 ( X 0 ) in the set S 0 and Γ 1 ( X 1 ) according to Table 7. The parameters determining Γ 1 ( X 1 ) in the CRTBP system are as Equation (15) shows:
X 1 = [ | | Δ V 1 | | , α 1 , β 1 ]
where α 1 is the angle between Δ V 1 and the projection of Δ V 1 in the x ^ y ^ plane, and β 1 is the angle between the projection of Δ V 1 in the x ^ y ^ plane and x ^ . The equation of Δ V 1 with respect to ( α 1 , β 1 ) is shown in Equation (16).
Δ V 1 = | | Δ V 1 | | ( cos α 1 cos β 1 , cos α 1 sin β 1 , sin α 1 ) T
Let us traverse the parameters according to Table 7, and calculate the following parameters of { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } : Δ t 1 , the perilune height h p m o o n , Jacobi constant J t for G 1 ( Γ 1 ) after a powered lunar flyby (PLF), and the perigee height h p e a r t h . For each { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } , we take these parameters into the pruning conditions shown in Equation (17).
J t C 1 h ˜ 1 h p m o o n h S O I h ˜ 2 h p e a r t h h ˜ 3
The third condition means that the height of the perigee after the PLF is similar to the orbital height of LEO, and the spacecraft will not enter the Earth’s atmosphere. After pruning, there are 334,841 pairs of { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } . The group composed of all orbits { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } is denoted as set S 1 .
When the parameters of Γ 0 ( X 0 ) are Δ τ = 25   days , | | Δ V 0 | | = 0.4   km / s , and ( α 0 , β 0 ) = ( 0 , 0 ) , Figure 8a,b show the distribution of α 1 and β 1 about i 1 . Each data point represents a G 1 ( Γ 1 ) orbit. [ a 1 , e 1 , i 1 , Ω 1 , ω 1 , M 1 ] are the orbital elements of the spacecraft when it arrives at the perigee. The orbital elements of the data of the red and blue box in Figure 8a,b are shown in Table 8. Their orbits are shown in Figure 8c–f. From Figure 8a,b, it can be seen that G 1 ( Γ 1 ) is symmetrically distributed about α 1 = 90 and α 1 = 270 , and that the period of variation of G 1 ( Γ 1 ) about β 1 is 180°.
Therefore, when calculating the parameters of Γ 0 ( X 0 ) and Γ 1 ( X 1 ) , in order to simplify the computation, the ranges of the ergodicity of α 0 and α 1 can be set as [ 0 , 90 ] [ 180 , 270 ] and the ranges of the ergodicity of β 0 and β 1 can be set as [ 0 , 180 ] .
Figure 8g,h show the distribution of i 1 and h p e a r t h about ( | | Δ V 0 | | + | | Δ V 1 | | ) . It is known that there are some orbits that satisfy i 1 > 80 . Comparing Figure 7g,h and Figure 8g,h shows that Γ 1 ( X 1 ) greatly increases the number of orbits that meet the requirement of DRO–Earth transfer. Within two maneuvers Γ 0 ( X 0 ) and Γ 1 ( X 1 ) , the spacecraft is able to transfer from DRO to Earth and fly by the SSO debris.

4.3. Local Optimization

The last subsection obtained the pruning results of the two-impulse DRO–Earth transfer orbit. In this subsection, { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } in set S 1 with i 1 > 80 are used as the initial values for local optimization to obtain the optimal two-impulse DRO–Earth transfer and flyby orbit. At first, the target debris of the flyby must be selected. The group composed of { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } with i 1 > 80 is denoted as set S 2 . In order to choose the target debris of the flyby for each { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } in set S 2 , let us traverse every { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } in set S 2 and all target debris, and then calculate the angle between the orbital plane of the debris and the spacecraft when the spacecraft at the perigee. For each { G 0 ( Γ 0 ) + G 1 ( Γ 1 ) } in set S 2 , the debris whose orbital plane is closest to the spacecraft is selected as the target debris.
After choosing the target debris of the flyby, the following optimization problem is to be solved: find the departure time τ 0 on DRO and two maneuvers Γ 0 ( X 0 ) and Γ 1 ( X 1 ) such that the total impulse | | Δ V 0 | | + | | Δ V 1 | | is minimized and the flyby condition of one piece of SSO target debris is satisfied. The framework of the optimization problem is shown in Equation (18). The minimization strategy is obtained via the built-in fmincon tool.
Minimize   f ( x )   subject   to : x min x x max c ( x ) 0
where f ( x ) = | | Δ V 0 | | + | | Δ V 1 | | is the cost function. x denotes the optimization parameters. x min and x max are the lower and upper boundaries of x . c ( x ) denotes the inequality constraints. The optimization parameters are as follows:
x = [ Δ τ , | | Δ V 0 | | , α 0 , β 0 , Δ t 0 , | | Δ V 1 | | , α 1 , β 0 , Δ t 1 ]
The lower and upper bounds of the optimization parameters are as follows:
x min = [ Δ τ δ 1 , | | Δ V 0 | | δ 2 , α 0 δ 3 , β 0 δ 3 , 0 ,   | | Δ V 1 | | δ 2 , α 1 δ 3 , β 1 δ 3 , 0 ] x max = [ Δ τ + δ 1 , | | Δ V 0 | | + δ 2 , α 0 + δ 3 , β 0 + δ 3 , ζ 1 ,   | | Δ V 1 | | + δ 2 , α 1 + δ 3 , β 1 + δ 3 , ζ 1 ]
where δ 1 , δ 2 , and δ 3 are defined in Table 4, and ζ 1 is defined in Table 3. The inequality constraints are as follows:
c 1 = ( h p m o o n h ˜ 1 ) 0 c 2 = h p m o o n h S O I 0 c 3 = J t C 1 0 c 4 = ( h p e a r t h h ˜ 2 ) 0 c 5 = h p e a r t h h ˜ 3 0 c 6 = ( i 1 i ˜ 1 ) 0 c 7 = ( | | δ r | | ζ 5 ) 0 c 8 = | | δ r | | ζ 6 0 c 9 = | | δ v | | ζ 7 0
After optimization, the optimization result with the minimum | | Δ V 0 | | + | | Δ V 1 | | is chosen as the optimal two-impulse DRO–Earth transfer and flyby orbit.

4.4. Optimal Result of Two-Impulse DRO–Earth Transfer

In this subsection, the optimal two-impulse DRO–Earth transfer orbits when the spacecraft departs from different DROs are to be discussed. Note that there are numerous DRO orbits. The most commonly discussed types of cislunar DRO include 3:2 resonance (i.e., the initial position in the Earth–Moon rotating frame is x0 ≈ 0.73 DU), 2:1 resonance with x0 ≈ 0.81 DU, and 3:1 resonance with x0 ≈ 0.86 DU. A comparison of optimal transfer orbits starting from different DROs is shown in Table 9.
In Table 9, Δ V ˜ = | | Δ V 0 | | + | | Δ V 1 | | , and Δ t ˜ = | | Δ t 0 | | + | | Δ t 1 | | . As Table 9 shows, when the spacecraft departs from the 3:2 resonant DRO, Δ V ˜ is the smallest, being 15.6% less than the Δ V ˜ of the 2:1 resonant DRO and 11.3% less than the Δ V ˜ of the 3:1 resonant DRO. Therefore, the subsequent design is based on the case that the spacecraft departs from the 3:2 resonant DRO. The orbital maneuver parameters of Γ 0 ( X 0 ) and Γ 1 ( X 1 ) when the spacecraft departs from the 3:2 resonant DRO are shown in Table 10. The orbits of the spacecraft in different coordinates are shown in Figure 9.
From Table 10, it can be seen that the spacecraft departs from the 3:2 resonant DRO at 2025/10/31/23:58:43, and the phase of departure is 236.70°. The impulse of departure | | Δ V 0 | | is 285 m/s. After 27.46 days, the spacecraft reaches the perilune. Near the perilune, the impulse of the PLF | | Δ V 1 | | is 185 m/s. Then, after 13.89 days, the spacecraft flies by the first selected piece of target debris, 25733. The relative distance of the flyby is 67.75 km, and the relative velocity is 3.99 km/s. The total velocity increment is 0.47 km/s, and the total transfer time of the DRO–Earth transfer and flyby is 41.35 days.
In Phase 1, the spacecraft finishes a DRO–Earth transfer and flies by one piece of SSO debris. On a computer with an AMD 3900X 64-core processor and 128 GB RAM, it takes about 3.5 h and 3.3 days to complete the traversal and pruning calculations of the single-impulse DRO–Earth transfer and the two-impulse DRO–Earth transfer, respectively. Next, the orbits needed for the spacecraft to fly by all the SSO debris in Phase 2 are to be designed.

5. Phase 2: HEO Based SSO Objects Flyby

Starting with the two-impulse DRO–Earth transfer orbit shown in Table 10, this section investigates the approach needed to fly by all of the SSO targets shown in Table A1. As shown in Figure 10a, the flyby procedure in this section is as follows: the spacecraft applies an impulse Δ V when it arrives at the apogee, and it flies by a piece of target debris near the perigee. Multiple transfer orbits of the spacecraft are designed to fly by the debris until all of the pieces of space debris in the database have been flown by at least once.
The transfer orbit from the perilune to the first flyby shown in Table 10 is the highly elliptical orbit (HEO). The apogee of HEO is near the altitude of the Moon, and its perigee is near the altitude of SSO, which allows the spacecraft in HEO to make better use of the lunar gravity assist and to quickly fly by SSO debris. Moreover, it is expected to maximize the use of the limited propellent to fly by debris. Lowering the altitude of the apogee will reduce the flight time but increase the fuel consumption. Thus, in this paper, the orbit of the multiple-object flyby is designated as HEO. Other types of orbits with lower apogees are left for later studies.
In order to make better use of the lunar gravity assist, the orbit of the multiple-object flyby is required to be a lunar-resonant orbit. The types of lunar-resonant orbits with their apogee near the lunar orbital altitude and their perigee near the SSO debris orbital altitude are 8:3, 5:2, 7:3, and 9:4, as shown in Figure 10b. Because the range of the apogee of the 8:3 resonant orbit is closer to the orbital radius of the Moon (which means that it is easier to adjust the orbit using the lunar gravity), in designing the flyby orbits, this paper tends to choose the 8:3 resonance HEO.
For the debris flyby mission, this section expects that the fuel consumption is small and the time of flight is short. These two optimization objectives cannot be optimal at the same time. This section uses a multi-objective evolutionary optimization algorithm to design the spacecraft’s orbit. By coordinating and compromising between each objective, the optimal solution, i.e., the Pareto-optimal solution, can be found among all the compromise solutions.
There exist several multi-objective evolutionary optimization algorithms, such as the non-dominated sorting genetic algorithm (NSGA-II), the strength Pareto evolutionary algorithm (SPEA2), and the multi-objective particle swarm optimization algorithm (MOPSO). Among these algorithms, the NSGA-II has lower computational complexity and running time [50]. This paper chooses the NSGA-II algorithm to solve the Pareto-optimal orbit with minimum impulse and minimum transfer time.
The calculation procedure of Phase 2 is shown in Figure 11. At first, the debris whose orbital plane is closest to the spacecraft after one orbital period of the spacecraft is selected as the flyby target. Then, the NSGA-II algorithm is used to design the Pareto-optimal flyby orbit of the spacecraft. In the Pareto frontier, the optimal solution with the minimum impulse is chosen as the flyby orbit because of the fact that the differences in the time of flight of all solutions in the Pareto frontier are not significant.
Because the multi-objective evolutionary algorithms have limitations on the diversity and convergence of solutions, the NSGA-II algorithm is not guaranteed to find the Pareto frontier when the parameters of the maneuver are not within the ranges for generating the initial population. One may think of trying to compensate for this shortcoming by expanding the number of populations and the ranges of generating the initial populations, but this approach will lead to a tremendous increase in computation time. This paper proposes a semi-analytical optimization method to design the optimal transfer orbit when the NSGA-II algorithm cannot obtain the Pareto-optimal orbit. The semi-analytical optimization method is to find the parameters of maneuver that can fly by a target debris in the two-body system and then substitute these parameters into the CRTBP system for optimization. In the case that the spacecraft’s orbital plane is far from all the target debris, the spacecraft will orbit the Earth without flying by any debris.

5.1. Multi-Objective Evolutionary Optimization Algorithm

The parameters of the NSGA-II algorithm in the CRTBP system are shown in Table 11.
The optimization parameters in the CRTBP system as shown in Table 12. The optimization objectives are as follows:
f 1 ( x ) = | | Δ V | | f 2 ( x ) = d t 1 + d t 2
The lower and upper bounds of the optimization parameters are as follows:
x min = [ 0 , 0 ° , 0 ° , 0.5 TU , 0.5 TU ] x max = [ ζ 4 , 2 π , 2 π , 1.35 TU , 1.35 TU ]
The inequality constraints are as follows:
c 1 = ( h p m o o n h ˜ 1 ) 0 c 2 = ( h p e a r t h h ˜ 2 ) 0 c 3 = ( | | δ r | | ζ 5 ) 0 c 4 = | | δ r | | ζ 6 0 c 5 = | | δ v | | ζ 7 0
If the spacecraft would enter the Moon’s SOI, the inequality constraints should contain the constraint of lunar gravity so that the inequality constraints would be as follows:
c 1 = ( h p m o o n h ˜ 1 ) 0 c 2 = h p m o o n h S O I 0 c 3 = ( h p e a r t h h ˜ 2 ) 0 c 4 = ( | | δ r | | ζ 5 ) 0 c 5 = | | δ r | | ζ 6 0 c 6 = | | δ v | | ζ 7 0
As discussed before, the orbit with the smallest | | Δ V | | in the two-objective Pareto-optimal frontier is chosen as the transfer orbit because of the fact that the differences in the time of flight of all solutions in the Pareto-optimal frontier are not significant. Next, the semi-analytical optimization method is proposed to solve the transfer orbit with the minimum impulse.

5.2. Semi-Analytical Optimization Method

This section presents the semi-analytical optimization method to solve the transfer orbit with the minimum velocity increment.
Assume that, in the CRTBP system, τ 3 is the previous flyby moment. The spacecraft’s position and the velocity at τ 3 in the center of mass rotating frame { C x ^ , y ^ , z ^ } are ( r 3 , v 3 ) . At τ 3 m , the spacecraft applies an impulse Δ V to change its orbit. Before changing the spacecraft’s orbit at τ 3 m , the position and the velocity in the { C x ^ , y ^ , z ^ } frame are ( r 3 m , v 3 m ) , which are denoted as ( R 3 m , V 3 m ) in the ECI coordinate { O X ^ , Y ^ , Z ^ } . The corresponding orbital elements are denoted as ( a 3 m , e 3 m , i 3 m , Ω 3 m , ω 3 m , f 3 m ) . After changing the spacecraft’s orbit at τ 3 m , the position and the velocity in the ECI coordinate { O X ^ , Y ^ , Z ^ } are denoted as ( R 3 m + , V 3 m + ) , and the corresponding orbital elements are ( a 3 m + , e 3 m + , i 3 m + , Ω 3 m + , ω 3 m + , f 3 m + ) .
At τ 3 m , the states should satisfy the following:
R 3 m = R 3 m + ,   Δ V = V 3 m + V 3 m
h 3 m + R 3 m ,   h 3 m + R 3 m +
h 3 m + R 3 m = 0 ,   h 3 m + R 3 m + = 0
where h 3 m + is the unit angular momentum after changing the spacecraft’s orbit at τ 3 m .
The flyby target is denoted as Target-1. The orbital elements of Target-1 at τ 3 are ( a t 3 , e t 3 , i t 3 , Ω t 3 , ω t 3 , f t 3 ) . At τ 4 , the spacecraft flies by Target-1. The orbital elements of Target-1 at τ 4 are ( a t 4 , e t 4 , i t 4 , Ω t 4 , ω t 4 , f t 4 ) , and the orbital elements of the spacecraft at τ 4 are ( a 4 , e 4 , i 4 , Ω 4 , ω 4 , f 4 ) .
It is required that the spacecraft will fly by the target debris after changing its orbit. At τ 4 , the orbital plane of the spacecraft should satisfy the following:
i 4 i τ 4 ,   Ω 4 Ω τ 4
The specific calculation procedure for the optimization method is composed of five steps as follows:
Step 1: Input the spacecraft’s states ( r 3 , v 3 ) at τ 3 in the { C x ^ , y ^ , z ^ } frame, and the states of Target-1 ( R t 3 , V t 3 ) at τ 3 in the ECI coordinate { O X ^ , Y ^ , Z ^ } . Integrate the spacecraft’s orbit an orbital period away from τ 3 in the two-body system. The moment that satisfies Equation (28) is denoted as τ e . The spacecraft’s states r e , v e at τ e in the CRTBP system are recorded.
In Equation (28), the equations of h 3 m + = [ h 3 m + x ; h 3 m + y ; h 3 m + z ] are calculated through the following:
{ h 3 m + x = sin ( Ω 3 m + ) 1 h 3 m + z 2 h 3 m + y = cos ( Ω 3 m + ) 1 h 3 m + z 2 h 3 m + z = cos i 3 m +
where i 3 m + = i τ 4 ,   Ω 3 m + = Ω τ 4 .
Step 2: τ e , r e , v e are converted into the ECI coordinate { O X ^ , Y ^ , Z ^ } and then are denoted as τ e , R e , V e . We calculate the spacecraft’s orbital elements ( a e , e e , i e , Ω e , ω e , f e ) after it changes its orbit at τ e . Then, let us implement the ergodic calculation as Equation (31) shows.
r p [ 6771   km : 20   km : 7171   km ] r a [ 380000   km : 500   km : 390000   km ] e e [ 0.8 : 0.01 : 0.99 ] i e = i 3 m + ,   Ω e = Ω 3 m +
The argument of latitude u e is calculated by Equation (32). From Equations (33) and (34), ω e , f e can be obtained.
R e = [ R e x R e y R e z ] = a e ( 1 e e 2 ) 1 + e e cos f e [ cos Ω 3 m + cos u sin Ω 3 m + sin u e cos i 3 m + sin Ω 3 m + cos u + cos Ω 3 m + sin u e cos i 3 m + sin u e sin i 3 m + ]
f e = a e cos ( ( a e ( 1 e e 2 ) | | R e | | 1 ) / e e )
ω e = u e f e
Step 3: Each set of ( a e , e e , i e , Ω e , ω e , f e ) is converted into the spacecraft’s position and velocity R e + , V e + . The impulse required to change the orbit is Δ V = V e + V e .
When | | Δ V | | 0.3   km / s , R e + , V e + is integrated in the two-body system. The time span of the integration is [ τ e , τ e + T e ] and the stepsize is set to 50 s. T e = 2 π a e 3 / μ e . We calculate the following parameters: the integral time τ 4 , and the spacecraft’s orbit R 4 , V 4 at τ 4 .
In the J 2 perturbed two-body system, the orbit of Target-1 is R t e , V t e at τ e , and R t e , V t e are integrated in the J 2 perturbed two-body system. The time span of the integration is [ τ e , τ e + T e ] and the stepsize is set to 50 s. Multiple sets of Target-1′s orbit R t 4 , V t 4 are obtained. The relative distance and the relative velocity between the spacecraft and Target-1 are as follows:
R r e l = R 4 R t 4 , V r e l = V 4 V t 4
For each set of [ τ 4 , R 4 , V 4 , R t e , V t e ] and τ e , the spacecraft’s orbital elements ( a e , e e , i e , Ω e , ω e , f e ) after orbit changes and d t = τ 4 τ e are recorded if | | R r e l | | 2000   km ,       |   | V r e l | | 6   km / s . After the pruning, the orbital parameters of the spacecraft after the maneuver and the moment of the maneuver that is able to fly by a piece of target debris in the two-body system are obtained.
Step 4: Each set of ( a e , e e , i e , Ω e , d t ) is used as the initial values of the optimization parameters in the CRTBP system. The objective function is defined as follows:
f ( x ) = | | Δ V | |
The optimization parameters are as follows:
x = [ a e , e e , i e , Ω e , d t ] T
The lower and upper boundaries of the optimization parameters are as follows:
x min = [ 1.9339 × 10 5   km ,   0.8 ,   80 ° ,   0 ,   0 ] x max = [ 2.2879 × 10 5   km ,   0.99 ,   110 ° ,   360 ° ,   13   days ]
The inequality constraints are as follows:
c 1 = ( h p m o o n h ˜ 1 ) 0 c 2 = ( h p e a r t h h ˜ 2 ) 0 c 3 = ( h p e a r t h h ˜ 3 ) 0 c 4 = ( | | δ r | | ζ 5 ) 0 c 5 = | | δ r | | ζ 6 0 c 6 = | | δ v | | ζ 7 0
The optimization frame in the CRTBP system is shown in Equations (35)–(38). From Step 4, the orbital parameters of the spacecraft and the moment of the maneuver obtained in Step 3 are optimized in the CRTBP system.
Step 5: The optimized solution with the smallest | | Δ V | | is selected as the design result for the spacecraft’s orbit. We output the following variables: the parameters of the maneuver [ | | Δ V | | , α , β , d t 1 , d t 2 ] T , where d t 1 is the time from the previous flyby moment to the application of the impulse in the current orbit, and d t 1 = τ e τ 3 and d t 2 = τ 4 τ e .
The calculation procedure of the semi-analytical optimization method is summarized as follows: in Step 1, the states of the spacecraft in the CRTBP system are converted into the two-body system. Then, in Step 2–Step 3, the orbital parameters of the spacecraft and the moment of the maneuver that are able to fly by a target piece of debris in the two-body system are obtained. In Step 4, these parameters are used as the initial values and are optimized in the CRTBP system. In Step 5, the optimized solution with the smallest | | Δ V | | is chosen to be the flyby orbit.
The calculation procedure of Phase 2 is shown in Figure 12. τ 3 is the previous flyby moment. D i is the serial number of the debris. angle i is the angle between the spacecraft and D i debris at τ 4 = τ 3 + 10.5   day . flyby i is the number of times that D i debris was flown by. The colors of the block diagram in Figure 12 correspond to Figure 11. The green parts indicate that the NSGA-II algorithm is utilized for optimization. For each flyby target, the transfer orbit is optimized using the NSGA-II algorithm up to two times when the Pareto-optimal results are not obtained. When a piece of debris cannot be flown by using the NSGA-II algorithm, another piece of debris will be chosen as the flyby target until all three pieces of debris have been used as flyby targets and Pareto optimization results have still not been obtained. The yellow part means that the semi-analytical optimization method is used and the purple diagram means that the spacecraft flies without flying by any pieces of debris.
Figure 13 shows the calculation procedure when the spacecraft orbits the Earth without flying by any target. T e is the orbit period of the spacecraft. As shown in Figure 13, when the spacecraft flies without flying by any target and the height of perigee and the inclination of the spacecraft satisfy the conditions in Equation (40), the spacecraft does not perform orbit maneuver, i.e., | | Δ V | | = 0 .
h ˜ 2 h p e a r t h h ˜ 3   i ˜ 1 i 1 i ˜ 2
where i 1 is the orbital inclination of the spacecraft; i ˜ 1 = 80 and i ˜ 2 = 110 . The reason that Equation (40) is not satisfied is that the spacecraft is influenced by the lunar gravity. There are two general cases. One is that an impulse should be applied to prevent the spacecraft from entering the Earth’s atmosphere when the height of the perigee is less than 400 km because of the fact that the lunar flyby is in front of the Moon. The other one is that when the height of the perigee is increased because of the fact that the lunar flyby is behind the Moon and the spacecraft is far away from the nearest debris, the spacecraft will continue to fly several revolutions without maneuver rather than changing the orbit with a large impulse to comply with the conditions of debris flyby. When an impulse should be applied, the optimization parameters are as Table 12 shows. The objective function is as follows:
f ( x ) = | | Δ V | |
The lower and upper boundaries of the optimization parameters are as follows:
x min = [ 0 , 0 ° , 0 ° , 0 , 0 ] x max = [ ζ 4 , 2 π , 2 π , 3   TU , 3   TU ]
The inequality constraints are as follows:
c 1 = ( h p m o o n h ˜ 1 ) 0 c 2 = ( h p e a r t h h ˜ 2 ) 0 c 3 = ( h p e a r t h h ˜ 3 ) 0 c 4 = ( i 1 i ˜ 1 ) 0 c 5 = ( i 1 i ˜ 2 ) 0
Equations (41)–(43) indicate that if the lunar gravity did not help the spacecraft to fly by a piece of debris, the impulse could be applied near the perigee to adjust the semi-major axis or the phase of the spacecraft to prevent the spacecraft from entering the lunar SOI or to reduce the effect of the lunar gravity.

6. Numerical Simulations and Analysis

The two-impulse DRO–Earth transfer orbit shown in Table 10 is defined as the 0th revolution orbit of the spacecraft. After the spacecraft first flies by the target debris in the 0th revolution, the orbit from the position of the first flyby to the next perigee is defined as the 1st revolution orbit. Similarly, the subsequent orbits are defined as the 2nd revolution orbit, 3rd revolution orbit … and Nth revolution orbit.
The spacecraft departs the 3:2 resonant DRO at 2025/10/31/23:58:43. The target debris 25733 is flown by in the 0th revolution orbit. The subsequent transfers and flybys are shown in Table 13. The total impulse of the flyby mission is 1.39 km/s and the total transfer time is 2.23 years. All of the spacecraft’s orbits are shown in Figure 14.
As shown in Table 13, it takes 83 revolutions for the spacecraft to flyby all the debris shown in Table A1. In 62 of the 83 revolutions, the spacecraft cruises around the Earth without flying by any debris.
The relative velocity versus distance of the 21 flybys are shown in Figure 15a. It shows that the minimum relative velocity of a flyby is about 3.23 km/s. The statistics of the impulses and the flight times are shown in Figure 15b,c. From Figure 15b, it can be seen that all impulses are less than 70 m/s. The average impulse of 83 revolutions is 11.05 m/s. The smallest impulse is 0.16 m/s (6th revolution) and the largest impulse is 60.70 m/s (50th revolution). In the 50th revolution orbit, the spacecraft applies a 60.70 m/s impulse to avoid the spacecraft entering the Earth’s atmosphere due to the lunar gravity. In addition, the time of flight of the 50th revolution orbit is reduced to about 5 days, as shown in Figure 15c. Figure 15c shows that before the 63rd revolution, the flight times of the most of HEO transfer orbits are between 10–12 days. This proves that the period of the lunar-resonant orbit is not constant due to the influence of the Moon in the CRTBP system. After the 63rd revolution, the flight times are all about 8 days. It is because of the fact that in the 63rd revolution, the spacecraft applies 23.31 m/s impulse near the perigee to reduce the height of the apogee; otherwise, the spacecraft would enter the Earth’s atmosphere due to the lunar gravity.
A key point in the HEO transfer orbit design problem is to appropriately utilize the lunar gravity. In the 2nd revolution, the lunar flyby is behind the Moon; therefore, the height of the perigee of the spacecraft is increased. However, in the 2nd revolution, the spacecraft is far away from the nearest debris. Thus, the spacecraft cruises around the Earth without taking any maneuvers in the 2nd and 3rd revolutions. In the 4th revolution, the lunar flyby is in front of the Moon and the spacecraft takes an impulse of only 6.61 m/s to decrease the height of the perigee with the help of lunar gravity. In the 30th revolution, the lunar flyby is behind the Moon, while in the 32nd revolution, the lunar flyby is in front of the Moon. Therefore, the spacecraft cruises around the Earth without taking any maneuvers in the 30th and 31st revolutions. Using the lunar gravity appropriately, the spacecraft only needs a 43.86 m/s impulse to fly by a piece of debris in the 32nd revolution. In the 42nd revolution, the spacecraft takes only a 16.52 m/s impulse to fly by a piece of debris with the help of the lunar gravity. In the 17th revolution, the spacecraft takes a 17.16 m/s impulse to lower the height of the perigee as the lunar flyby is behind the Moon. Otherwise, the spacecraft would not be able to fly by two pieces of debris in the 20th and 21st revolutions. The same occurs in the 45th revolution. The spacecraft takes 51.41 m/s impulse to increase the orbit height as the lunar flyby is in front of the Moon so that the spacecraft can fly by three pieces of debris in the 47th, 48th, and 49th revolutions, respectively.
The time that the spacecraft cruises around the Earth without taking any maneuver is about 0.89 years. This is because of the fact that in order to reduce fuel consumption, except for necessary maneuvers such as adjusting the orbital altitude to prevent the spacecraft from entering the Earth’s atmosphere, the spacecraft does not apply any of the impulses until the orbital plane is close to a piece of space debris.
On a computer with an AMD 3900X 64-core processor and 128 GB RAM, it takes about 3 days to complete the design of Phase 2.
The optimized result obtained by the semi-analytical optimization method in the 1st revolution is as follows:
| | Δ V | | = 8.60   m / s ,   α = 6.24   rad ,   β = 0.62   rad d t 1 = 3.99   days ,   d t 2 = 6.40   days
Next, Equation (44) is used as an example to analyze the optimality of the maneuver solved by the semi-analytical optimization method. The two-body analytical optimal single-impulse transfer theory is [51]: as shown in Figure 16a, the spacecraft applies an impulse at P 1 and then flies to P 2 . v 0 is the initial velocity of P 1 , and the velocity increment is as follows:
Δ v = v 1 v 0
Decomposing v 1 in the chord direction and in the radial direction yields the following:
v 1 = v ρ 1 + v c 1 v ρ 1 = | | v ρ 1 | | e r 1 , v c 1 = | | v c 1 | | e c | | v ρ 1 | | = c μ p | | r 1 | | | | r 2 | | sin θ , | v c 1 | = μ p 1 cos θ sin θ
where
e r 1 = r 1 | | r 1 | | , e c = r 2 r 1 | | r 2 r 1 | | c = | | r 2 r 1 | | , p = a 1 ( 1 e 1 2 )
The optimal v 1 should satisfy the following:
( v ρ 1 v c 1 ) ( v ρ 1 + v c 1 v 0 ) = 0
which means that the optimal Δ v should satisfy the following:
( v ρ 1 v c 1 ) Δ v
Equations (45)–(49) present the conditions that need to be satisfied for the optimal single-impulse transfer in the field of the two-body system. At τ 2 , the spacecraft completes the 0th revolution. In the 1st revolution, τ 2 m = τ 2 + d t 1 , τ 3 = τ 2 m + d t 2 . The spacecraft’s position at τ 2 m on the 1st revolution orbit is denoted as P 1 , and it applies an 8.7 m/s impulse at P 1 . Then, the spacecraft’s orbit is integrated from τ 2 m to τ 3 in the two-body system. The position at τ 3 is denoted as P 2 . Substituting these parameters into Equations (45)–(49) yields that the angle between Δ v and ( v ρ 1 v c 1 ) is 91.24° or 271.24°. Therefore, in the two-body system, the impulse in the two-body system is near optimal.
Figure 16b shows the spacecraft’s orbit from τ 2 to τ 3 in the CRTBP system and the two-body system, respectively. The spacecraft applies an impulse at τ 2 m and flies by the target debris 25732 at τ 3 . During this period, the Moon is far away from the spacecraft; as a result, the Moon has little influence on the spacecraft. Thus, the impulse calculated by the semi-analytical optimization method is near optimal in the CRTBP system.
It is worthwhile to compare the transfer strategy from DRO to the debris orbit in this paper to the other low-energy transfer. For the problem of LEO constellation design, reference [52] proposes a strategy in which all satellites carried by several vehicles transfer from a single launch to the lunar L1 halo orbit; each vehicle needs 0.6 km/s impulse for halo orbit insertion. Each satellite can return to the desired LEO with small propulsive burns (minimal ΔV in Earth return trajectory and 0.12 km/s to adjust the perigee). Compared with reference [52], the DRO–Earth transfer strategy presented in this paper is more fuel-efficient because the DRO parking orbit is not only long-term stable but also requires less than 0.3 km/s impulse to depart.
Based on reference [52], reference [53] proposes a similar transfer strategy that uses the manifold, with six vehicles carrying 36 satellites transferring to L1 halo orbit. Then, the satellites perform station-keeping and transfer back to LEO in different epochs when the orbital planes of the desired orbits are appropriate. The average total time of flight of the return trajectories is 15.98 days (12.34 days for station-keeping and 3.66 days for transfer). Differing from reference [53], the transfer strategy in this paper takes the departure time as one of the optimization parameters. Therefore, once the spacecraft leaves DRO, it does not need to perform station-keeping while entering directly into the Moon–Earth transfer orbit. Compared with the total time of flight of the return trajectory in reference [53], it can be seen that the flight time in this paper is lower (12.89 days).
Reference [54] presents a strategy to build a lunar constellation from the same launch trajectory in the Sun–Earth–Moon system. First, the launch vehicle approaches the region of L1 equilibrium via a low-energy orbit. Then, the satellites are released and transferred to low-energy lunar capture orbits. Although this method is designed for lunar constellations, the LEO–Moon transfer is nearly identical to that transferring back when no perturbations are considered. Comparing the Moon–LEO transfer in this paper with the LEO–perilune transfer in Reference [54], it can be seen that although the transfer time of this paper is longer than that of Reference [54], the impulse near the perilune of this paper is 0.043 km/s less than that of Reference [54] (0.228 km/s).

7. Conclusions

The paper investigates a one-to-many Chinese SSO space debris flyby mission. The debris in this multiple-SSO-object flyby mission are chosen by analyzing the distribution characteristics of the orbital location, type, source, and size of the SSO space cataloged objects. Based on the ergodic searching method, pruning function, and local optimization, the optimal orbit of the optimal two-impulse transfer and flyby when the spacecraft departs from three types of resonant DRO are obtained. The results show that transferring from the 3:2 resonant DRO has the smallest velocity increment. A semi-analytical optimization method is proposed to design the optimal large elliptical transfer orbits when the NSGA-II algorithm cannot obtain the Pareto-optimal orbit. The simulation results show that all target SSO space debris in the database are near-coplanar flown by at least once. The total velocity increment of the mission is 1.39 km/s, and the total time of flight is 2.23 years. The results of this paper can provide the basis for subsequent multiple space debris removal missions.

Author Contributions

Conceptualization, S.Z. and S.W.; methodology, S.Z.; validation, S.Z. and S.W.; formal analysis, S.Z.; investigation, S.Z.; writing—original draft preparation, S.Z.; writing—review and editing, S.W. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDA30010500, and the Lab funding of the Key Laboratory of Space Utilization, Chinese Academy of Sciences, Grant No. CXJJ-22S021.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

The physical parameters of 22 pieces of target debris are shown in Table A1. The orbital parameters are shown in Table A2. The catalog is from the North American Aerospace Defense Command (NORAD) [55].
Table A1. Physical parameters of target debris.
Table A1. Physical parameters of target debris.
No.NORAD
Catalog Number
TypeMass
(kg)
ShapeWidth
(m)
Height
(m)
Depth
(m)
Diameter
(m)
Span
(m)
Max. Cross Section (m2)Min. Cross Section (m2)
125733Debris50cone2.91.02.90N/A2.906.611.45
243012Rocket Body1000cylinderN/A1.92N/A2.902.908.635.56
337215Rocket Body1000cylinderN/A7.50N/A2.902.9022.736.60
432959Rocket Body1000cylinderN/A7.50N/A2.902.9022.736.60
525732Rocket Body1000cylinderN/A6.20N/A2.906.2019.156.60
625942Rocket Body1000cylinderN/A7.50N/A2.907.5022.736.60
728059Rocket Body1000cylinderN/A7.50N/A2.907.5022.736.60
827432Rocket Body1000cylinderN/A6.20N/A2.906.2019.156.60
939261Rocket Body1000cylinderN/A7.50N/A2.902.9022.736.60
1032063Rocket Body1000cylinderN/A7.50N/A2.907.5022.736.60
1131114Rocket Body3800cylinderN/A8.40N/AN/A8.4029.488.81
1239203Rocket Body3800cylinderN/A8.40N/AN/A8.429.488.81
1343610Rocket Body3800cylinderN/A9.94N/AN/A9.9434.448.81
1445722Rocket Body3800cylinderN/AN/AN/AN/A7.4926.598.81
1537731Rocket Body3800cylinderN/A8.40N/AN/A8.4029.488.81
1636089Rocket Body3800cylinderN/A8.40N/AN/A8.4029.488.80
1737766Rocket Body3800cylinderN/A8.40N/AN/A8.4029.488.81
1840262Rocket Body3800cylinderN/A8.40N/AN/A8.4029.488.81
1944548Rocket Body4000cylinderN/A8.90N/A3.358.9031.098.81
2041858Rocket Body4006cylinderN/A8.00N/A3.358.0028.218.81
2137932Rocket Body4006cylinderN/A8.00N/A3.358.0028.218.81
2239154Rocket Body4006cylinderN/A8.00N/A3.358.0028.218.81
Table A2. Orbital parameters of target debris.
Table A2. Orbital parameters of target debris.
No.NORAD
Catalog Number
a (km)ei (°) Ω   (°) ω   (°) M (°)
1257337215.450.00143798.91210.47350.699.40
2430127083.690.01271698.73186.11344.7015.04
3372157115.790.00874798.82315.56295.2064.01
4329597130.700.00691698.90131.89282.5887.49
5257327211.300.00364698.90210.33233.36178.63
6259427149.590.00987698.80157.51194.49224.30
7280597095.130.00513098.50100.6186.59274.12
8274327223.480.00501299.08193.56195.95221.49
9392617159.090.00320098.96341.928.18351.99
10320637103.010.00545497.97177.71279.44140.95
11311147205.630.00584698.234.81274.6684.79
12392037087.820.00510198.47359.04345.2514.73
13436107090.600.00909398.57285.13148.48212.19
14457227116.030.00519298.59125.2320.02340.30
15377317024.180.00601197.65180.41152.2720.73
16360897104.230.00814098.20148.55196.06222.06
17377667050.140.0044498.3358.18242.68116.99
18402627035.010.00463398.30330.04137.72222.76
19445487139.020.00060598.13170.97262.8697.19
20418587150.190.00335998.50159.81184.49175.60
21379327196.480.00381198.69111.14108.96262.11
22391547020.850.00206298.4155.4388.48271.88

References

  1. NASA Webpage. Orion Will Go the Distance in Retrograde Orbit During Artemis I. Available online: https://www.nasa.gov/missions/orion-will-go-the-distance-in-retrograde-orbit-during-artemis-i/ (accessed on 2 February 2023).
  2. NASA Webpage. Meet NASA’s Orion Spacecraft. Available online: https://www.nasa.gov/missions/meet-nasas-orion-spacecraft/ (accessed on 2 February 2023).
  3. Broucke, R. NASA Technical Report 32-1168: Periodic Orbits in the Restricted Three-Body Problem with Earth-Moon Masses; Jet Propulsion Laboratory, California Institute of Technology: Pasadena, CA, USA, 1968.
  4. Bezrouk, C.J.; Parker, J. Long duration stability of distant retrograde orbits. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference, San Deiego, CA, USA, 4–7 August 2014. [Google Scholar] [CrossRef]
  5. Strange, N.; Landau, D.; McElrath, T.; Lantoine, G.; Lam, T.; McGuire, M.; Burke, L.; Martini, M.; Dankanich, J. Overview of mission design for NASA asteroid redirect robotic mission concept. In Proceedings of the 33rd International Electric Propulsion Conference, Washington, DC, USA, 6–10 October 2013. [Google Scholar]
  6. Bezrouk, C.J.; Parker, J. Long term evolution of distant retrograde orbits in the Earth-Moon system. Astrophys. Space Sci. 2017, 326, 176. [Google Scholar] [CrossRef]
  7. Welch, C.M.; Parker, J.S.; Buxton, C. Mission Considerations for Transfers to a Distant Retrograde Orbit. J. Astronaut. Sci. 2015, 62, 101–124. [Google Scholar] [CrossRef]
  8. Murakami, N.; Yamanaka, K. Trajectory design for rendezvous in lunar distant retrograde orbit. In Proceedings of the 2015 IEEE Aerospace Conference, Big Sky, MT, USA, 7–14 March 2015; pp. 1–13. [Google Scholar] [CrossRef]
  9. Demeyer, J.; Gurfil, P. Transfer to small distant retrograde orbits. In Proceedings of the American Institute of Physics Conference, Greensboro, NC, USA, 1–2 August 2007; pp. 20–31. [Google Scholar] [CrossRef]
  10. Capdevila, L.; Guzzetti, D.; Howell, K. Various transfer options from Earth into distant retrograde orbits in the vicinity of the Moon. Adv. Astronaut. Sci. 2014, 118, 3659–3678. [Google Scholar]
  11. Ren, J.; Li, M.; Zheng, J. Families of transfers from the Moon to distant retrograde orbits in cislunar space. Astrophys. Space Sci. 2020, 365, 192. [Google Scholar] [CrossRef]
  12. Xu, M.; Xu, S. Exploration of distant retrograde orbits around Moon. Acta Astronaut. 2009, 65, 853–860. [Google Scholar] [CrossRef]
  13. Zhang, R.; Wang, Y.; Zhang, C.; Zhang, H. The transfers from lunar DROs to Earth orbits via optimization in the four body problem. Astrophys. Space Sci. 2021, 366, 49. [Google Scholar] [CrossRef]
  14. Wen, C.; Gao, Y. Reachability study for spacecraft maneuvering from a distant retrograde orbit in the Earth-Moon system. In Proceedings of the 18th Australian International Aerospace Congress, Melbourne, Australia, 24–28 February 2019. [Google Scholar]
  15. Aixue, W.; Chen, Z.; Wang, S.; Zhang, H. Design considerations for access in to Earth-Moon DROs with lunar free-return trajectory. Manned Spacefl. 2022, 28, 81–89. (In Chinese) [Google Scholar] [CrossRef]
  16. Zeng, H.; Li, Z.; Peng, K.; Wang, P.; Huang, Z. Research on application of Earth-Moon NRHO and DRO for Lunar exploration. J. Astronaut. 2020, 41, 910–919. [Google Scholar] [CrossRef]
  17. Capdevila, L.R.; Howell, K.C. A transfer network linking Earth, Moon and the triangular libration point regions in the Earth-Moon system. Adv. Space Res. 2018, 62, 1826–1852. [Google Scholar] [CrossRef]
  18. Morand, V.; Dolado-Perez, J.-C.; Philippe, T.; Handschuh, D.-A. Mitigation rules compliance in low Earth orbit. J. Space Saf. Eng. 2014, 1, 84–92. [Google Scholar] [CrossRef]
  19. PARABOLIC ARC Webpage. NASA ODPO’s Large Constellation Study. Available online: http://34.196.180.244/2018/09/25/nasa-odpos-large-constellation-study/ (accessed on 29 December 2023).
  20. Liou, J.; Johnsonb, N. Instability of the present LEO satellite populations. Adv. Space Res. 2008, 41, 1046–1053. [Google Scholar] [CrossRef]
  21. Wormnes, K.; Le Letty, R.; Summerer, L.; Schonenborg, R.; Dubois-Matra, O.; Luraschi, E.; Cropp, A.; Krag, H.; Delaval, J. ESA technologies for space debris remediation. In Proceedings of the 6th European Conference on Space Debris, ESA Communications ESTEC Noordwijk, Daemstadt, Germany, 22–25 April 2013; pp. 85–92. [Google Scholar]
  22. Botta, E.; Sharf, I.; Misra, A.K. Evaluation of net capture of space debris in multiple mission scenarios. In Proceedings of the 26th AAS/AIAA Space Flight Mechanics Meeting, Napa, CA, USA, 14–18 February 2016. [Google Scholar]
  23. Nikolajsen, J.A.; Kristense, A.S. Self-deployable drag sail folded nine times. Adv. Space Res. 2021, 68, 4242–4251. [Google Scholar] [CrossRef]
  24. Underwood, C.; Viquerat, A.; Schenk, M.; Taylor, B.; Massimiani, C.; Duke, R.; Stewart, B.; Fellowes, S.; Bridges, C.; Aglietti, G.; et al. Inflatesail deorbit flight demonstration results and follow-on dragsail applications. Acta Astronaut. 2019, 162, 344–358. [Google Scholar] [CrossRef]
  25. Bekey, I. Tethers open new space options. Astronaut. Aeronaut. 1983, 21, 32–40. [Google Scholar]
  26. King, L.B.; Parker, G.G.; Deshmukh, S.; Chong, J.-H. Spacecraft Formation-Flying Using Inter-Vehicle Coulomb Forces; Technique Report; NASA Institute for Advanced Concepts: Atlanta, GA, USA, 2002.
  27. Monroe, D.K. Space debris removal using a high-power ground-based laser. In Proceedings of the Space Programs and Technologies Conference and Exhibit, Huntsville, AL, USA, 21–23 September 1993; pp. 276–283. [Google Scholar] [CrossRef]
  28. Kofford, S. System and Method for Creating an Artificial Atmosphere for the Removal of Space Debris. U.S. Patent US2013/0082146A1, 4 April 2013. [Google Scholar]
  29. Andrenucci, M.; Pergola, P.; Ruggiero, A. Active Removal of Space Debris Expanding Foam Application for Active Debris Removal; Technical Report; ESA Advanced Concepts Team: Noordwijk, The Netherlands, 2011. [Google Scholar]
  30. Hughes, J.; Schaub, H. Prospects of using a pulsed electrostatic tractor with nominal geosynchronous conditions. IEEE Trans. Plasma Sci. 2017, 45, 1887–1897. [Google Scholar] [CrossRef]
  31. Hammert, J.; Schaub, H. Debris attitude effects on electrostatic tractor relative motion control performance. In Proceedings of the 2021 AAS/AIAA Astrodynamics Specialist Conference, Big Sky, MT, USA, 8 August 2021. [Google Scholar]
  32. Rubenchik, A.M.; Barty, C.P.J.; Beach, R.J.; Erlandson, A.C.; Caird, J.A. Laser systems for orbital debris removal. In Proceedings of the AIP Conference, Santa Fe, NM, USA, 18–22 April 2010; pp. 347–353. [Google Scholar] [CrossRef]
  33. Phipps, C.R. A laser-optical system to re-enter or lower low Earth orbit space debris. Acta Astronaut. 2014, 93, 418–429. [Google Scholar] [CrossRef]
  34. Gregory, D.A.; Mergen, J.-F. Space Debris Removal Using Upper Atmosphere and Vortex Generator. U.S. Patent US8657235B2, 25 February 2014. [Google Scholar]
  35. Dunn, M.J. Space Debris Removal. U.S. Patent US8800933B2, 12 August 2014. [Google Scholar]
  36. Pergola, P.; Ruggiero, A.; Andrenucci, M.; Olympio, J.; Summerer, L. Expanding foam application for active space debris removal systems. In Proceedings of the 62nd International Astronautical Congress, Cape Town, SA, USA, 27 September–1 October 2010. [Google Scholar]
  37. Castronuovo, M.M. Active space debris removal—A preliminary mission analysis and design. Acta Astronaut. 2011, 69, 848–859. [Google Scholar] [CrossRef]
  38. Tadini, P.; Tancredi, U.; Grassi, M.; Anselmo, L.; Pardini, C.; Francesconi, A.; Branz, F.; Maggi, F.I.; Lavagna, M.I.; DeLuca, L.T.; et al. Active debris multi-removal mission concept based on hybrid propulsion. Acta Astronaut. 2014, 103, 26–35. [Google Scholar] [CrossRef]
  39. Van der Pas, N.; Lousada, J.; Terhes, C.; Bernabeu, M.; Bauer, W. Target selection and comparison of mission design for space debris removal by DL’s advanced study group. Acta Astronaut. 2014, 102, 241–248. [Google Scholar] [CrossRef]
  40. Federici, L.; Zavoli, A.; Colasurdo, G. On the use of A* search for active debris removal mission planning. J. Space Saf. Eng. 2021, 8, 245–255. [Google Scholar] [CrossRef]
  41. De Bonis, A.; Angeletti, F.; Lannelli, P.; Gasbarri, P. Mixed kane-newton multi-body analysis of a dual-arm robotic system for on-orbit servicing missions. Aerotec. Missili Spazio 2021, 100, 375–386. [Google Scholar] [CrossRef]
  42. Schaub, H.; Junkins, J. Analytical Mechanics of Space Systems, 3rd ed.; American Institute of Aeronautics and Astronautics, Inc.: Reston, VA, USA, 2014; Chapter 10; p. 547. [Google Scholar] [CrossRef]
  43. Zhang, R.; Wang, Y.; Zhang, H.; Zhang, C. Transfers from distant retrograde orbits to low lunar orbits. Celest. Mech. Dyn. Astron. 2022, 132, 41. [Google Scholar] [CrossRef]
  44. Zhang, R.; Wang, Y.; Shi, Y.; Zhang, C.; Zhang, H. Performance analysis of impulsive station-keeping strategies for cislunar orbits with the ephemeris model. Acta Astronaut. 2022, 198, 152–160. [Google Scholar] [CrossRef]
  45. Zhang, C.; Zhang, H. Lunar-gravity-assisted low-energy transfer from Earth into DROs. Acta Aeronaut. Astronaut. Sin. 2022, 44, 326507. (In Chinese) [Google Scholar] [CrossRef]
  46. Tselousova, A.; Shirobokov, M.; Trofimov, S. Direct two-impulse transfer from a low-Earth orbit to high circular polar orbits around the Moon. AIP Conf. Proc. 2019, 2171, 130022. [Google Scholar] [CrossRef]
  47. Anthony, W.; Larsen, A.; Butcher, E.A.; Parker, J.S. Impulsive guidance for optimal manifold-based transfers to Earth-Moon L1 halo orbits. In Proceedings of the 23rd AAS/AIAA Spaceflight Mechanics Meeting, Kauai, HI, USA, 10–14 February 2013. [Google Scholar]
  48. Conte, D.; He, G.-W.; Spencer, B.; Melton, G. Trajectory design for LEO to lunar Halo orbits using manifold theory and fireworks optimization. In Proceedings of the 2018 AAS/AIAA Astrodynamics Specialist Conference, Snowbird, UT, USA, 19–23 August 2018. [Google Scholar]
  49. United Nations Office for Outer Space Affairs. Guidelines for the Long-Term Sustainability of Outer Space Activities. Available online: https://www.unoosa.org/res/oosadoc/data/documents/2019/aac_105c_1l/aac_105c_1l_366_0_html/V1805022.pdf (accessed on 2 February 2023).
  50. Shah, V.; Beeson, R.; Coverstone, V. A method for optimizing low-energy transfer in the Earth-Moon system using global transport and genetic algorithms. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference, Long Beach, CA, USA, 13–16 September 2016. [Google Scholar]
  51. Battin, R.H. An Introduction to the Mathematics and Methods of Astrodynamics, revised ed.; American Institute of Aeronautics and Astronautics Press: Reston, VA, USA, 2001; Chapter 11; p. 443. ISBN 978-7-5159-1488-6. [Google Scholar]
  52. Chase, J.P.; Chow, N.; Gralla, E.; Kasdin, N.J. LEO constellation design using the lunar L1 point. In Proceedings of the AAS/AIAA Space Flight Mechanics Meeting, Maui, HI, USA, 8–12 February 2004. [Google Scholar]
  53. Nadoushan, M.J.; Novinzadeh, A.B. Satellite constellation build-up via three-body dynamics. J. Aerosp. Eng. 2013, 228, 155–160. [Google Scholar] [CrossRef]
  54. Carletta, S. A single-launch deployment strategy for lunar constellations. Appl. Sci. 2023, 13, 5104. [Google Scholar] [CrossRef]
  55. CelesTrack Webpage. Satellite Catalog (SATCAT). Available online: https://celestrak.org/satcat/search.php (accessed on 26 February 2024).
Figure 1. Mission scenario. (a) The panorama of flyby mission; (b) the enlarged view of flyby.
Figure 1. Mission scenario. (a) The panorama of flyby mission; (b) the enlarged view of flyby.
Universe 10 00145 g001
Figure 2. Mission scenario. (a) The center of mass rotating frame; (b) the J2000 ECI coordinate.
Figure 2. Mission scenario. (a) The center of mass rotating frame; (b) the J2000 ECI coordinate.
Universe 10 00145 g002
Figure 3. Illustration of flyby mission.
Figure 3. Illustration of flyby mission.
Universe 10 00145 g003
Figure 4. Orbits of target debris, J2000 ECI coordinate.
Figure 4. Orbits of target debris, J2000 ECI coordinate.
Universe 10 00145 g004
Figure 5. (a) Two impulses DRO–Earth transfer orbit in the Earth–Moon Rotating Frame; (b) illustration of Δ V 0 in the Earth–Moon rotating frame.
Figure 5. (a) Two impulses DRO–Earth transfer orbit in the Earth–Moon Rotating Frame; (b) illustration of Δ V 0 in the Earth–Moon rotating frame.
Universe 10 00145 g005
Figure 6. Distribution of all orbits in set S 0 . (a) α 0 versus i * ; (b) β 0 versus i * .
Figure 6. Distribution of all orbits in set S 0 . (a) α 0 versus i * ; (b) β 0 versus i * .
Universe 10 00145 g006
Figure 7. (a) Data of the red box in Figure 6a, the center of mass rotating frame. (b) Data of the red box in Figure 6a, the J2000 ECI frame. (c) Data of the blue box in Figure 6a, the center of mass rotating frame. (d) Data of the blue box in Figure 6a, the J2000 ECI frame. (e) Orbit of the spacecraft when i * = 62.70 , the center of mass rotating frame. (f) Orbit of the spacecraft when i * = 62.70 , J2000 ECI fame. (g) Distribution of | | Δ V 0 | | i 1 . (h) Distribution of | | Δ V 0 | | h p e a r t h . The purple * represents the Earth and the purple + represents the Moon.
Figure 7. (a) Data of the red box in Figure 6a, the center of mass rotating frame. (b) Data of the red box in Figure 6a, the J2000 ECI frame. (c) Data of the blue box in Figure 6a, the center of mass rotating frame. (d) Data of the blue box in Figure 6a, the J2000 ECI frame. (e) Orbit of the spacecraft when i * = 62.70 , the center of mass rotating frame. (f) Orbit of the spacecraft when i * = 62.70 , J2000 ECI fame. (g) Distribution of | | Δ V 0 | | i 1 . (h) Distribution of | | Δ V 0 | | h p e a r t h . The purple * represents the Earth and the purple + represents the Moon.
Universe 10 00145 g007
Figure 8. (a) Distribution of α 1 i 1 . (b) Distribution of β 1 i 1 . (c) Orbit of the data indicated by the red box in (a), the center of mass rotating frame. (d) Orbit of the data indicated by the red box in (a), the J2000 ECI frame. (e) Orbit of the data indicated by the blue box in (b), the center of mass rotating frame. (f) Orbit of the data indicated by the blue box in (b). (g) Distribution of all orbits of set S 1 , ( | | Δ V 0 | | + | | Δ V 1 | | ) versus i 1 . (h) Distribution of all orbits of set S 1 , ( | | Δ V 0 | | + | | Δ V 1 | | ) versus h p e a r t h . The purple * represents the Earth and the purple + represents the Moon.
Figure 8. (a) Distribution of α 1 i 1 . (b) Distribution of β 1 i 1 . (c) Orbit of the data indicated by the red box in (a), the center of mass rotating frame. (d) Orbit of the data indicated by the red box in (a), the J2000 ECI frame. (e) Orbit of the data indicated by the blue box in (b), the center of mass rotating frame. (f) Orbit of the data indicated by the blue box in (b). (g) Distribution of all orbits of set S 1 , ( | | Δ V 0 | | + | | Δ V 1 | | ) versus i 1 . (h) Distribution of all orbits of set S 1 , ( | | Δ V 0 | | + | | Δ V 1 | | ) versus h p e a r t h . The purple * represents the Earth and the purple + represents the Moon.
Universe 10 00145 g008aUniverse 10 00145 g008b
Figure 9. Optimal DRO–Earth transfer orbit. (a) The center of mass rotating frame. (b) The J2000 ECI frame. The purple * represents the Earth and the purple + represents the Moon.
Figure 9. Optimal DRO–Earth transfer orbit. (a) The center of mass rotating frame. (b) The J2000 ECI frame. The purple * represents the Earth and the purple + represents the Moon.
Universe 10 00145 g009
Figure 10. (a) Mission scenario of Phase 2; (b) ranges of the perigee and the apogee of the lunar-resonant orbits.
Figure 10. (a) Mission scenario of Phase 2; (b) ranges of the perigee and the apogee of the lunar-resonant orbits.
Universe 10 00145 g010
Figure 11. Calculation procedure of Phase 2.
Figure 11. Calculation procedure of Phase 2.
Universe 10 00145 g011
Figure 12. Design of multiple transfer orbits.
Figure 12. Design of multiple transfer orbits.
Universe 10 00145 g012
Figure 13. Calculation process when the spacecraft flies without flying by a target.
Figure 13. Calculation process when the spacecraft flies without flying by a target.
Universe 10 00145 g013
Figure 14. HEO transfer orbits: (a) the center of mass rotating frame; (b) J2000 ECI frame.
Figure 14. HEO transfer orbits: (a) the center of mass rotating frame; (b) J2000 ECI frame.
Universe 10 00145 g014
Figure 15. Analysis of HEO orbits: (a) histories of the distance with regard to the relative velocity; (b) impulses; (c) flight times.
Figure 15. Analysis of HEO orbits: (a) histories of the distance with regard to the relative velocity; (b) impulses; (c) flight times.
Universe 10 00145 g015
Figure 16. (a) Schematic of optimal single-impulse transfer; (b) the 1st revolution orbit in the CRTBP system and the two-body system.
Figure 16. (a) Schematic of optimal single-impulse transfer; (b) the 1st revolution orbit in the CRTBP system and the two-body system.
Universe 10 00145 g016
Table 1. Velocity increment budgets for transferring from LEO to different target orbits.
Table 1. Velocity increment budgets for transferring from LEO to different target orbits.
Target Orbits||ΔVtotal||Types of Transfer
DRO3.1934 km/s [45]Lunar gravity assist and weak stability boundary
GEO4.2457 km/s [45]Direct two impulses
PLO3.8 km/s [46]Direct two impulses
NRHO4 km/s [46]Direct two impulses
L1 halo orbit3.61 km/s [47]Two impulses and stable manifold
L2 halo orbit3.36 km/s [48]Two impulses and stable manifold
Table 2. Flyby mission procedure.
Table 2. Flyby mission procedure.
StepsDescription
Step 1At τ ¯ (chosen as 2025/10/1/00:00:00), the spacecraft is located in DRO. The initial phase is 0°. τ ˜ = 2025/11/1/00:00:00.
Step 2At τ 0 ( τ ¯ < τ 0 < τ ˜ ), the spacecraft applies an impulse ΔV0 and leaves DRO, flying towards the Moon.
Step 3The spacecraft enters the Moon’s sphere of influence. At τ 1 , the spacecraft applies an impulse ΔV1 and accomplishes a powered lunar flyby transfer (PLF). The spacecraft enters the Moon–Earth transfer orbit. Then, at τ 2 , the spacecraft near-coplanar flies by a piece of SSO space debris near the perigee.
Step 4The spacecraft flies towards to the apogee. To implement orbit phasing and adjust the orbit, the spacecraft applies an impulse ΔV2 near the apogee. As a result, the spacecraft can near-coplanar fly by the next piece of SSO space debris when the spacecraft is close to the perigee.
Step 5The spacecraft repeats Step 4 until all of the space debris in the database have been flown by at least once.
Table 3. Requirements of the flyby mission.
Table 3. Requirements of the flyby mission.
ParametersMeaningRequirement
Δ τ Days between τ ¯ and τ 0 Δ τ ζ 1
Δ t 0 The transfer time from DRO to the perilune Δ t 0 ζ 1
Δ t 1 The transfer time from the perilune to near the perigee Δ t 1 ζ 1
| | Δ V 0 | | The modulus of Δ V 0 | | Δ V 0 | | ζ 2
| | Δ V 1 | | The modulus of Δ V 1 | | Δ V 1 | | ζ 3
| | Δ V | | The modulus of impulse applied in HEO | | Δ V | | ζ 4
Δ t The period of HEOHEO is resonant with the Moon
h p m o o n The height of perilune h p m o o n h ˜ 1
h p e a r t h The height of perigee h p e a r t h h ˜ 2
| | δ r | | The distance of the flyby ζ 5 | | δ r | | ζ 6
| | δ v | | The magnitude of the relative flyby velocity | | δ v | | ζ 7
Table 4. Parameters of Γ 0 ( X 0 ) and their ergodic ranges.
Table 4. Parameters of Γ 0 ( X 0 ) and their ergodic ranges.
ParametersMeaningRange of Ergodicity
Δ τ Days between τ ¯ and τ 0 The sweep interval is [1, ζ 1 ] days. The sweep stepsize is δ 1 .
| | Δ V 0 | | The modulus of Δ V 0 The sweep interval is [0.005, ζ 2 ] km/s. The sweep stepsize is δ 2 .
( α 0 , β 0 ) Two direction angles of Δ V 0 The sweep interval is [0° 360°]. The sweep stepsize is δ 3 .
Table 5. Orbit parameters of the data indicated by the boxes in Figure 6a,b.
Table 5. Orbit parameters of the data indicated by the boxes in Figure 6a,b.
Colors of Box Δ τ (day) | | Δ V 0 | | (km/s) ( α 0 , β 0 ) (°) a * (km) e * i * (°) Ω * (°) ω * (°) M * (°)
red, Figure 6a230.3650(80,0)224,8180.7838.97348.60319.72240.08
230.3650(100,180)224,8180.7838.97348.60319.72240.08
blue, Figure 6a240.3950(240,180)217,6260.8546.97350.77342.92202.02
240.3950(300,0)217,6260.8546.97350.77342.92202.02
red, Figure 6b110.1850(190,40)221,9270.8257.17351.40161.45226.10
110.1850(350,220)221,9270.8257.17351.40161.45226.10
blue, Figure 6b90.1100(190,130)197,0360.8652.04356.46162.54218.22
90.1100(350.310)197,0360.8652.04356.46162.54218.22
Table 6. Orbit elements when the spacecraft enters and exits the Moon’s SOI.
Table 6. Orbit elements when the spacecraft enters and exits the Moon’s SOI.
Epoch a (km) e i (°) Ω (°) ω (°) M (°)
When spacecraft enters the Moon’s SOI373,8020.6920.987.043.20296.40
When spacecraft exits the Moon’s SOI209,9440.8262.7143.329.93248.97
Table 7. Parameters and their ergodic ranges.
Table 7. Parameters and their ergodic ranges.
ParametersMeaningRange of Ergodicity
Γ 0 ( X 0 ) The first maneuver Γ 0 ( X 0 ) of each G 0 ( Γ 0 ) in set S 0
| | Δ V 1 | | The modulus of Δ V 1 The sweep interval is [0.005, ζ 3 ] km/s. The sweep stepsize is δ 2 .
( α 1 , β 1 ) Two direction angles of Δ V 1 The sweep interval is [ 0 , 360 ] . The sweep stepsize is δ 3 .
Table 8. Orbit parameters of the data indicated by the red and blue boxes in Figure 8a,b.
Table 8. Orbit parameters of the data indicated by the red and blue boxes in Figure 8a,b.
Colors of Box | | Δ V 1 | | (km/s) ( α 1 , β 1 ) (°) a 1 (km) e 1 i 1 (°) Ω 1 (°) ω 1 (°) M 1 (°)
red, Figure 8a0.2000(120, 90)118,5880.9413.69114.27235.9440.00
0.2000(60, 270)118,5880.9413.69114.27235.9440.00
blue, Figure 8a0.2000(240, 90)118,5870.9453.82355.67348.26359.99
0.2000(300, 270)118,5870.9453.82355.67348.26359.99
red, Figure 8b0.1550(230, 170)119,7010.9443.35359.40347.980.00
0.1550(310, 350)119,7010.9443.35359.40347.980.00
blue, Figure 8b0.2000(300, 260)125,2310.9451.05356.39347.44359.99
0.2000(240, 80)125,2310.9451.05356.39347.44359.99
Table 9. Comparison of optimal transfer and flyby orbits for different initial DRO.
Table 9. Comparison of optimal transfer and flyby orbits for different initial DRO.
Parameters | | Δ V 0 | | / m s 1 Δ t 0 /days | | Δ V 1 | | / m s 1 Δ t 1 /days | | Δ V ˜ | | / m s 1 Δ t ˜ /daysFlyby Debris
DRO
3:2 resonance285.027.46185.013.89470.041.3525733
2:1 resonance371.228.02185.73.85556.931.8737766
3:1 resonance365.026.67165.013.92530.040.5932959
Table 10. Parameters of Γ 0 ( X 0 ) and Γ 1 ( X 1 ) , and the results of the flyby.
Table 10. Parameters of Γ 0 ( X 0 ) and Γ 1 ( X 1 ) , and the results of the flyby.
Departure Time φ 0 | | Δ V 0 | | / m s 1 Δ t 0 /days | | Δ V 1 | | / m s 1 Δ t 1 /daysFlyby Debris | | δ r | | /km | | δ v | | / k m s 1
2025/10/31/23:58:43236.70285.027.46185.013.892573369.753.99
Table 11. Parameters of NSGA2 algorithm.
Table 11. Parameters of NSGA2 algorithm.
Population SizeMaximum IterationsSimulates Binary
Crossover Parameters  η c
Polynomial Mutation Parameter  η m Probability of Mutation  p m Probability of Crossover  p c
1500150201000.21
Table 12. Optimization parameters of the NSGA2 algorithm.
Table 12. Optimization parameters of the NSGA2 algorithm.
ParametersMeaning
| | Δ V | | The modulus of the impulse Δ V .
( α , β ) α is the angle between Δ V and the projection of Δ V in the x ^ y ^ plane, and β is the angle between the projection of Δ V in the x ^ y ^ plane and x ^ .
d t 1 Time from the flyby point of the previous revolution to the moment of maneuver in the current revolution.
d t 2 Time from the moment of maneuver to fly by the target debris.
Table 13. Results of transfers and flybys.
Table 13. Results of transfers and flybys.
Revolution of OrbitsTarget Debris | | Δ V | | / m s 1 α /rad β /rad d t 1 /day d t 2 /day | | δ r | | /km | | δ v | | / k m s 1
1257328.746.260.594.036.2799.044.00
2~3/////23.19//
4/6.614.160.114.3710.70//
5/6.924.853.283.036.62//
6457220.164.800.520.009.68100.003.90
7/0.651.076.270.0010.06//
8/4.602.082.234.994.93//
9~11/////31.34//
12/47.164.403.850.1611.84//
13377667.004.781.926.614.4899.914.00
14320630.221.825.900.0011.06100.003.80
153926128.000.933.945.715.68100.004.00
163920329.946.276.289.821.76100.004.00
17/17.164.235.140.009.36//
18~19/////16.87//
204026220.291.042.003.994.3893.043.44
21/11.694.451.708.130.00//
2243610////7.12100.003.23
23/////7.09//
24/27.420.071.583.923.41//
25~26/////13.96//
274301238.551.714.670.0011.59100.003.96
28~31/////46.18//
322743343.861.600.722.958.47100.003.40
332743249.030.122.110.108.42100.003.86
342594212.472.861.744.743.8099.994.00
35/37.874.075.514.024.52//
36~40/////41.25//
41391543.545.382.251.786.61100.003.49
422805916.524.990.414.325.12100.003.59
43/1.804.382.504.654.60//
44/////9.19//
45/51.415.914.264.695.83//
46/////10.60//
47379329.954.444.884.106.31100.004.00
48360895.555.673.254.775.4299.884.00
49377313.712.845.294.046.28100.003.70
50/60.704.220.120.004.90//
51/53.770.870.000.0011.84//
52/////12.05//
53/56.270.814.090.409.65//
54/4.390.995.225.295.23//
55~56/////20.97//
57/6.202.164.5111.930.00//
58/////12.04//
59/4.164.222.135.215.14//
60/2.145.232.975.105.11//
61~62/////21.01//
63/23.315.323.560.007.62//
643295925.773.490.942.235.30100.004.00
65~66/////15.19//
67/48.234.224.804.573.39//
68311149.803.894.224.453.0599.983.88
69~73/////38.36//
74/42.515.212.233.913.92//
75/1.444.243.493.733.73//
76/0.975.184.913.783.74//
77/////7.46//
783920344.645.834.313.194.76100.004.00
79~80/////15.77//
81/31.865.192.493.864.00//
82/10.625.985.973.614.25//
8337215////7.5699.993.79
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, S.; Wang, S. Multiple SSO Space Debris Flyby Trajectory Design Based on Cislunar Orbit. Universe 2024, 10, 145. https://doi.org/10.3390/universe10030145

AMA Style

Zhang S, Wang S. Multiple SSO Space Debris Flyby Trajectory Design Based on Cislunar Orbit. Universe. 2024; 10(3):145. https://doi.org/10.3390/universe10030145

Chicago/Turabian Style

Zhang, Siyang, and Shuquan Wang. 2024. "Multiple SSO Space Debris Flyby Trajectory Design Based on Cislunar Orbit" Universe 10, no. 3: 145. https://doi.org/10.3390/universe10030145

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop