Orbital Maneuver Optimization of Earth Observation Satellites Using an Adaptive Differential Evolution Algorithm

: Earth observation satellite (EOS) systems often encounter emergency observation tasks oriented to sudden disasters (e


Introduction
Earth observation satellite (EOS) systems can acquire images of the Earth's surface via their remote sensing instruments. Due to the advantages such as large-scale observation coverage and high observation frequency, EOSs have been widely implemented to monitor and observe disasters such as earthquakes, floods, landslides, and debris flow [1][2][3]. Although the number of EOSs is continuously increasing, there are still several limitations to satisfy all kinds of user requirements. For example, when an earthquake occurs, EOSs are required to take ground images urgently to provide timely support for rescue operations. However, EOSs in their regular orbits may not be able to observe the earthquake area timely or clearly. Thus, an appropriately selected satellite needs to be transferred to a new orbit to provide better coverage properties, which is termed the orbit maneuver optimization problem.
Generally, the orbit maneuver optimization problem can be treated as a kind of orbit design problem [4][5][6]. Numerous studies have been carried out to investigate the orbit design problem. For example, Graham et al. [7] studied a minimum-time Earthorbit transfers optimization problem using low-thrust propulsion with eclipsing. They developed an initial guess generation method to construct a useful guess and analyzed the approximate place where the spacecraft enters and exits the Earth's shadow. A similar problem was addressed by Wang et al. [8], who adopted a convex optimization method. Zhang et al. [9] investigated a minimum-fuel optimization problem using lowthrust in the circular restricted three-body scenario. By considering actuation uncertainties, Mohammadi et al. [10] proposed a robust optimization approach for the impulsive orbit transfers optimization problem. In their study, the genetic algorithm, Monte-Carlo sampling, and surrogate model are combined to balance the optimization accuracy and time. Cheng et al. [11] developed a real-time optimal control approach based on multiscale deep neural networks for the orbit transfer problem of the solar sail spacecraft. In a recent study, Morante et al. [12] proposed a multi-objective optimization approach for an orbit-raising optimization problem, in which chemical, electrical, and hybrid trajectories are considered.
However, most of the orbit design problems aim to find an optimal orbit for improving orbit performance (e.g., coverage time and fuel consumption) [4,5,13]. Those studies assume that the satellite flies along a fixed orbit without orbit maneuvers and consider orbit elements as decision variables. For the cases in which orbit maneuvers are considered, there are few existing studies that mainly focus on reconfiguration problems of satellite constellations [14][15][16]. For example, McGrath et al. [17] presented a satellite constellation reconfiguration problem, in which a restricted low-thrust Lambert rendezvous scenario was included. Soleymani et al. [18] investigated an optimal mission planning problem of the reconfiguration process of satellite constellations. They applied a combination of particle swarm optimization and genetic algorithm to find the optimal departure and arrival positions of each satellite. He et al. [19] developed a physical programming method together with a genetic algorithm, to solve a multi-objective satellite constellation reconfiguration problem for disaster monitoring purposes. Wang et al. [20] proposed a hybrid-resampling particle swarm optimization method for an agile satellite constellation design problem, in which different types of sensors, the attitude maneuver of sensors, and different coverage performance indices are considered. To satisfy the requirements of emergency observation, a recent study proposed by Hu et al. [21] carried out a multiobjective optimization framework for the satellite constellation optimization problem.
It can be concluded that although many relevant studies have been published, the orbit maneuver optimization problem that optimizes maneuvers of a satellite is still a minor branch of orbit design problems and is rarely investigated. Hence, in this study, we make effort to address the orbit maneuver optimization problem from a scheduling perspective. Specifically, since a satellite can transfer its orbit by conducting an impulsive maneuver at a specific time instance and the maneuver result would affect the orbit performance, it would be crucial to determine the reasonable magnitude and direction of the impulse, as well as the maneuver moment. Different from most of the previous studies that aim to determine the promising position (i.e., orbit elements) of a satellite, our study optimizes the orbit maneuver in terms of velocity increments for an impulsive maneuver and the maneuver moment. Meanwhile, our study considers multiple satellites and the most suitable satellite would be selected to execute the task according to scheduling results.
On the other hand, to improve the service quality, diverse user requirements are being considered in the orbit maneuver optimization in recent years. For instance, since the fuel capacity is limited and the remaining fuel affects the lifetime of a satellite, some users may require a low fuel consumption solution. In case of some emergency tasks that need to be accomplished at all costs, the users may want the satellite to respond to observation requests as quickly as possible. Further, in some rescue operations, the orbit altitude is the optimization objective since an appropriate orbit altitude that can provide higher ground resolution is crucial. Therefore, we build three models that, respectively, optimize three objectives, including response time, ground resolution, and fuel consumption to satisfy diverse user requirements. Meanwhile, since we focus on EOS, specific constraints such as the resolution constraint are included in models.
Since the studied problem considers orbit maneuvers at every second as decision variables, the search space would be very large. Meanwhile, specific constraints of EOS would increase the difficulties of solving the problem. All of the above reasons propose challenges for solving the problem. In this regard, evolutionary algorithms would be a promising solution method owing to their powerful and effective search capabilities. Previously, evolutionary algorithms have been widely employed to address the orbit design problem. The algorithms used mainly include particle swarm optimization [20,22,23], genetic algorithms [16,21,24], and hybrid algorithms [5,25]. For example, Shirazi [25] applied a hybridization of the genetic algorithm and simulated annealing to a multiobjective orbit maneuver optimization problem. Based on the particle swarm optimization (PSO) algorithm, Pontani et al. [22] solved four kinds of impulsive orbital transfer problems, focusing on the optimization of impulsive transfers between two coplanar and non-coplanar, circular and elliptic orbits, respectively. Yao et al. [26] investigated the application of an improved DE algorithm on an orbit design problem by adding self-adaption and stochastic mechanisms. To optimize coverage-related metrics, as well as the number and semi-major axes of satellites in multiple constellations, Hitomi et al. [27] proposed a variable-length chromosome-based evolutionary algorithm.
Particularly, our studied problem can be treated as a continuous optimization problem. As a simple and efficient evolutionary algorithm, especially for continuous optimization, differential evolution (DE) which has rarely been implemented by previous related studies would be a promising candidate for addressing our problem. However, due to the well-known no-free-lunch theorem [28], the same optimization algorithm with the same configurations may have different performances on different problems. We have three models with different constraints and different optimization objectives, which propose challenges for optimizers. Moreover, DE highly depends on the configuration of genetic strategies and control parameters [29]. It would be time-consuming to find effective combinations of configurations to obtain high-quality solutions on different optimization models by using the same algorithm. Previously, many techniques have been developed to relieve this issue, such as ensemble and adaption techniques [28,30,31]. In this study, we implement the adaption technique to DE. Specifically, we form the genetic strategies and parameters of DE into a directed acyclic graph, in which each path indicates a combination of the genetic strategies and parameters. As the pheromone trails and property always enable the ant colony to find a reasonable path from the graph, an ant colony optimization (ACO) is adopted to search for effective combinations during the evolution. The hybridization of ACO and DE exhibits the effective search capability of DE that has been proved in previous studies [4,26,32]. Furthermore, it can dynamically optimize the algorithm configurations to improve the adaptive capability of DE, such that higher-quality solutions can be obtained for all three optimization models.
In summary, this paper has the following contributions.
(i) We investigate the orbit maneuver optimization problem considering diverse user requirements. In the problem, a satellite is selected from a set of satellites and transferred to a new orbit based on appropriate maneuvers (i.e., the velocity increment and maneuver moment) to respond to an emergency observation request. By analyzing orbit coverage and dynamics, we build three optimization models that optimize response time, ground resolution, and fuel consumption, respectively, to satisfy different user requirements.
(ii) To solve the proposed optimization models, we implement an adaptive differential evolution based on graph search. In the algorithm, key algorithm components (i.e., crossover strategies, mutation strategies, and control parameters) are formed into a directed acyclic graph and an ACO is adopted to find reasonable combinations of configurations during the evolution. The implemented algorithm is a hybrid of ACO and DE, therefore it is named ACODE.
(iii) We conduct simulation experiments to verify the efficiency of ACODE. The ACODE is compared with three representative algorithms including EPSDE [33], CSO [34], and SLPSO [35] in simulation scenarios where multiple EOSs are requested to observe a ground target. The simulation results show the superiority of ACODE.
This paper is organized as follows. Section 2 details the orbit coverage and dynamics analysis, as well as three optimization models. Sections 3 and 4 introduce the solution method and simulation experiments, respectively. Finally, the conclusions are remarked by Section 5.

Problem Description
In this section, we elaborate on the orbit maneuver optimization problem based on orbit coverage and dynamics calculations, followed by three optimization models with different optimization objectives (i.e., response time, fuel consumption, and ground resolution). As a part of the satellite system design, orbit maneuver optimization is associated with many complicated environmental factors. Therefore, some reasonable assumptions are adopted to simplify the problem.
(i) There are some perturbations (e.g., atmospheric drag, solar radiation pressure, and third body effects) that have negative impacts on the operation of the satellite. We only consider J 2 perturbation of Earth oblateness in the model, which is a common assumption in existing studies on orbit design problems [6,32,36].
(ii) Assume that the sensor equipped on each satellite is visible to the ground target when the satellite flies in the sunshine and the sunshine time is from 6:00 to 18:00 local time. Further, the other factors that may affect the imaging such as clouds and weather conditions, as well as the altitude of ground targets are assumed to be negligible.
(iii) Each satellite is assumed to be independent. Therefore, the orbit maneuver of a satellite does not affect the flying of another satellite.
(iv) The time required by the satellite to process task information and start the rocket engine is assumed to be negligible.
(v) The ground target is assumed to be a point target. Hence, the ground target can be imaged by the satellite once the satellite passes over it.
Main notations used in this section are displayed in Table 1.

Orbit Coverage Analysis
The visibility between a satellite and a ground target depends on many factors, such as the location of the ground target (i.e., longitude and latitude), the orbit elements, and the field of view (FOV) of the satellite. To conduct the orbit coverage analysis, we assume that the Earth is a round body, the orbit is approximately circular, and the FOV on the ground is rectangular as in [32,37]. The ground target is visible to the satellite when it lies in the FOV, which can be determined by calculating the longitudes and latitudes of four vertices. A typical satellite coverage on the Earth is shown in Figure 1. As Figure 1 shows, the Earth's angular radius λ h defines the half-ground range that may visible to the satellite, which can be expressed by where R E is the Earth's radius, and r sat is the distance between the Earth's center and the satellite. The slant range to the horizon, ρ h , can be written as However, in practical applications, due to some limitations such as the imaging angle of the sensor and sunshine conditions, the actual half ground range would be smaller than λ h . Hence, a general expression for the slant range to any point, ρ, can be expressed by [38] where γ is the intermediate angle and η is the boresight angle of the satellite (i.e., half of the sensor angle). Afterward, the half-ground range from the subsatellite point can be calculated by For the satellite equipped with a scanning sensor, the geometry of the FOV is no longer symmetrical about the subsatellite point, requiring more processing to obtain the ground range angle. Given the center boresight angle η c of the satellite, the maximum and minimum ground-range angles from the subsatellite point can be obtained. Specifically, the maximum and minimum boresight angles can be written as [38] Then, the maximum and minimum Earth's angular radiuses (i.e., λ max and λ min ) can be obtained according to Equations (5)- (7). Since the sensor of the satellite can be rotated on multiple axes, in this study we assume that the sensor half-angle equals η max and the Earth's angular radius equals λ max for convenience. Define the latitude and longitude of the subsatellite point as [lat s , lon s ], the latitudes and longitudes of the four vertices of the FOV can be calculated by According to the latitude and longitude information of the FOV, the latitude and longitude information of the ground target, the positions of the satellite at each moment, as well as the right ascension of Greenwich at the initial moment, we can obtain the key orbit performance indices of a satellite [39], such as the response time [32,40]. The response time is defined as the time required from when a request is received to observe a ground target until the satellite can observe it. The method that assesses whether the target lies in the FOV at moment t, as well as the response time t imag can be found in [32]. Moreover, the calculation method of the latitude and longitude of the subsatellite point at each moment is introduced in the next section.

Orbit Dynamics Model
The position of a satellite in its orbit can be obtained by using six orbit elements, including semimajor axis a, eccentricity e, inclination i, longitude of ascending node Ω, argument of perigee ω, and true anomaly θ, as Figure 2 shows. By using the orbit elements, we can calculate the position of the satellite, as well as the latitude and longitude of each subsatellite point at each moment. In this section, we briefly introduce the calculation methods, and more detailed derivation steps can be found in [41]. Given a satellite flying around the Earth, it is well-known that the period P for an orbit of the satellite is calculated by where µ is the Earth's gravitational parameter. Then, the time t 0 since perigee at the initial epoch can be calculated as where M is the initial mean anomaly. According to Kepler's equation, M can be calculated by where E is the eccentric anomaly, which yields the relation with true anomaly θ as Given a time change ∆t, the longitude of ascending node Ω, argument of perigee ω at the moment t = t 0 + ∆t can be expressed by whereΩ andω are time variations of Ω and ω, which are determined by J 2 perturbation of Earth oblateness. The expressions ofΩ andω are written aṡ where J 2 = 1.083 × 10 −3 . The orbit elements are updated by repeating Equations (9)-(13) at each moment t. Meanwhile, the newly found true anomaly θ at the moment t can be used to calculate the state vector of the satellite in the perifocal coordinate coordinate system (PQW). The satellite position vector r O and velocity vector v O in the PQW frame can be expressed by where h is the angular momentum of the satellite, yielding a relation with the semimajor axis a as below Particularly, the position vectors r O can be transformed to the Earth-centered inertial (ECI) frame through the transformation matrix R I/O (C ≡ cos and S ≡ sin) written as by where r I is the satellite position in the ECI frame. Meanwhile, r I can be expressed in the Rotating Earth-fixed frame by [4] which can be written in vector notation Define a notation r = √ X 2 + Y 2 + Z 2 , the latitude and longitude of the subsatellite point can be calculated by Based on the above equations, the position of the satellite, as well as the latitude and longitude of each subsatellite point at each moment can be obtained.
Furthermore, in this study, the orbit maneuver focuses on how to move a satellite in the same plane, which can be treated as a co-orbital rendezvous problem [42]. In the co-orbital rendezvous problem, two satellites are assumed to be located in the same orbit and one satellite maneuvers its orbit by two-impulse Hohmann transfer to catch up with the other one. Therefore, a velocity increment ∆v at the moment t m is considered to calculate the orbit elements before and after maneuvering based on orbit equations.

Formulation of the Optimization Problem
As mentioned above, the optimization problem aims to find appropriate velocity increment ∆v and maneuver moment t m to obtain a reasonable scheduling scheme for transferring the satellite. Hence, ∆v and t m can be considered as decision variables of the optimization problem, and they are constrained by Since ∆v is associated with the capacity of fuel, which is the key parameter of the satellite remained lifetime, constraint (25) ensures that the velocity increment is limited to a reasonable range. Here ∆v max represents the maximum allowed velocity increment and the negative value indicates the velocity in the reverse direction. Constraint (26) defines the range of a feasible maneuver moment (when the satellite starts its rocket engine). Furthermore, there are other constraints introduced in the following.
To obtain sufficient information from a single observation result, constraint (27) guarantees that the ground resolution is smaller than the required resolution, in which the ground resolution is associated with the satellite altitude over the ground target divided by the horizontal number of pixels. When the satellite is transferring its orbit, the change of orbit altitude should be limited to a reasonable range to ensure the stable operation of the satellite. As the satellite altitude is typically between 250 and 1300 km [43], constraint (28) is carried out to limit the satellite altitude after maneuvering. Constraint (29) is used to ensure that the observation moment lies in the sunshine time window. Furthermore, as timeliness is crucial for emergency observation tasks, we use constraint (30) to limit the maximum response time of the satellite.
In practical applications, users require different solutions depending upon the purpose. To satisfy diverse user requirements, we build three models considering response time, ground resolution, and fuel consumption as objectives, respectively, which are written as (29).

s.t. Constraints
(32) The calculation processes of objectives are as follows. During the optimization process, appropriate decision variables (i.e., velocity vector increment ∆v and maneuver moment t m ) will be searched by the algorithm under the constraints mentioned in Equations (31)- (33). Since the satellite conducts an impulsive maneuver whose direction and magnitude are determined by ∆v at moment t m to transfer its orbit, the new state velocity vector at moment t m can be determined by decision variables. Then, the state velocity vector of the satellite can be transformed into the position vector represented by orbit elements by using Equations (19) and (20). As the satellite flies around the Earth, the position vector of the satellite changes with time. The changes in position vector can be tracked by Kepler's equation coupled with orbit equations mentioned in Equations (8)- (18). Meanwhile, according to the position vector, the subsatellite point at the same moment can be obtained by Equations (21)- (24). The FOV of the satellite is determined by the subsatellite point at the same moment according to Equations (1)-(7) introduced in the orbit coverage analysis. When a FOV covers the target point, it indicates that the satellite can observe this target point. Thus, the first objective f 1 is calculated by the difference between the maneuver time t m and the time instance when the satellite can observe the target. In the second objective f 2 , the satellite altitude is the distance between the satellite and the subsatellite point when the satellite can observe the target. As to the third objective f 3 , it is determined by the decision variable ∆v directly.

Adaptive Differential Evolution Algorithm Based on Graph Search
Differential evolution (DE) is an efficient population-based stochastic optimization approach for solving optimization problems over continuous space, and many variants of DE have been implemented in engineering fields [30,44,45]. In this study, we conduct problem-specific modifications on the framework of an adaptive DE, named ACODE, first proposed in [31] that concerned data clustering problems, to solve the orbit maneuver optimization problem. The ACODE can be treated as a hybridization of DE and ACO, which will be detailed in this section after a brief introduction to the classical DE.

Classical DE Algorithm
Typically, the DE includes four basic steps [46]: Initialization, mutation, crossover, and selection.
(i) Initialization. This step randomly creates an initial population consisting of N individuals. When the iteration number G = 0, the i-th individual is initialized in the search space constrained by the minimum bound X min = {x 1 min , x 2 min , . . . , x D min } and the maximum bound X max = {x 1 max , x 2 max , . . . , x D max }, according to the following method where rand(0, 1) is a uniformly distributed number within [0, 1] and D is the number of dimensions.
(ii) Mutation. The mutation operation perturbs a target vector X i,G from the current generation to obtain a donor vector V i,G , which can be written as where F is the scaling factor, which is a positive control parameter for scaling the difference vectors. The indices r i 1 , r i 2 , and r i 3 are mutually exclusive integers randomly chosen from the range [1, N] and they are different from the base vector index i.
(iii) Crossover. The crossover operation can improve the diversity of the population by exchanging the components of the donor vector V i,G with the target vector X i,G to form the There are two kinds of commonly used crossover strategies, including exponential (i.e., two-point modulo) and binomial (i.e., uniform). The exponential crossover makes the trial vector contains a sequence of consecutive components taken from the parent vector. The structure of the trial vector can be expressed by where j n is a modulo function with modules D, k and L are two integers randomly chosen from [1, D]. On the other hand, the binomial strategy can be outlined as where CR is the crossover rate and j rand is a randomly chosen index lying in the interval [1, D].
(iv) Selection. The selection operation determines whether the target or the trial vector survives to the next generation according to the objective function, which is described as Once an initial population is created, the mutation, crossover, and selection strategies are repeated until a stopping criterion is satisfied to obtain promising solutions. It should be noted that different mutation strategies demarcate a DE scheme from other schemes. Except for the mutation strategy introduced above, there are some other well-known mutation strategies, such as "DE/best/1", "DE/best/2", "DE/rand/2", "DE/rand-tobest/1", "DE/current-to-pbest/1", "DE/current-to-rand/1", etc. [46].

ACODE Algorithm
The performance of DE highly depends on four key components, i.e., mutation strategy, crossover strategy, scaling factor F, and crossover rate CR [31]. We transform these components into a directed acyclic graph and implement an ant colony optimization-based adaptive DE algorithm to conduct the optimization process.

Directed Acyclic Graph Formed by Configurations
An example of the directed acyclic graph formed by configurations is shown in Figure 3. The graph includes five levels, of which a virtual point lies in the first level to gather all ants and the remaining four levels represent four key components (i.e., mutation strategy, crossover strategy, scaling factor F, and crossover rate CR), respectively. Every node in each level indicates a candidate configuration, and the nodes of two adjacent levels are fully connected. A path that starts from the first level and terminates at the fifth level can be treated as a combination of four candidate configurations, as the blue path in Figure 3 shows. Mathematically, the directed acrylic graph can be described by Φ = {V, E}, where V is the set of nodes and E ⊆ V × V indicates directed arcs. The pheromone trail on the arcs connecting v ∈ V to adjacent nodes in the next level is recorded by pheromone vector B v . Hence, the length of B v depends on the number of nodes in the next level. According to empirical considerations [31,33,46], the candidate configurations we used in this paper are summarized in Table 2

Framework of ACODE
The framework of ACODE is displayed in Algorithm 1. The algorithm starts with the initialization of a population P with size N, configuration matrix M, iteration counter g, and pheromone matrix B g (line 1). Particularly, the number of individuals in P equals the number of ants. Then, the algorithm runs until the stopping criterion is satisfied (lines 2-11). Every ant at each iteration is utilized to find a reasonable combination of configurations from the graph, and the combination is recorded in M (line 4). M i indicates the configuration combination of the i-th individual and it is implemented to evolve individual x i,g (line 5). The offspring u i,g is compared with x i,g to determine which solution is preserved into the next generation (lines 6-9). Finally, the pheromone matrix B g+1 that would be used in the next generation is updated according to the fitness information in line 10.

Solution Representation and Initialization
As above-mentioned, we consider the velocity increment ∆v and maneuver moment t m as decision variables. During the optimization process, the decision variables are used to calculate the positions of the satellite expressed by orbit elements, while the objectives are determined by position vectors, as introduced in Section 2. The velocity increment ∆v is the magnitude of the change in the velocity vector, which can be represented by three velocities (i.e., ∆v x , ∆v y , and ∆v z ) on three axes of the Cartesian coordinate system. Here the X axis is directed to the eccentricity vector, Z axis is in the direction of the satellite's angular momentum which lies perpendicular to the orbital plane, and the Y axis completes the right-hand set of co-ordinate axis. Therefore, a chromosome should be composed of velocity increments in three dimensions and the moment when the maneuver occurs. Figure 4 shows the representation of a chromosome, where x i is the expression vector while t m , ∆v x , ∆v y , and ∆v z are decision variables. Since the satellite conducts a coplanar maneuver, the velocity increment ∆v z always equals 0. Nevertheless, we still include ∆v z in the chromosome for computation convenience. In addition, we generate the initial population randomly and the boundaries of decision variables are determined by constraints (25) and (26). Since the search space of each variable is large and the performance of DE algorithm is seriously influenced by the diversity of the initial population, the initial population should be uniformly distributed in the search space. We use the Latin hypercube sampling (LHS) to generate the initial population. The LHS is a statistical method that can generate a quasi-random sampling distribution, which has been widely applied in other studies to obtain a high-quality initial population [47].

Parameters Adaption Based on ACO
Based on the directed acyclic graph, we conduct the parameter adaption by utilizing an ACO method. Specifically, each ant searches for a reasonable combination (i.e., path) for configuring an individual according to the pheromone trail on each arc in the graph. Given a node v from which an ant departs, there would be h candidate arcs that can be chosen. The probability p j for picking arc j is written as where B g v,j is the pheromone trail on arc j with the starting node v at the g-th iteration. The process of parameter adaption is shown in Algorithm 2. In the algorithm, an ant departs from the virtual node v 1 and travels through four nodes in the remaining four levels. At the l-th level, all probabilities for choosing arcs connecting starting node v l with all nodes in the next level are calculated based on Equation (39) and recorded in P l (line 3). Then, roulette wheel selection (i.e., RouletteWheel()) is adopted to choose an arc j that determines the node v l+1 (i.e., end node of arc j) at the (l + 1)-th level (lines 4-5). The roulette wheel selection is a well-known stochastic selection method, in which the probability for the selection of an arc is proportional to the pheromone trails on it. The above steps are repeated until a path consisting of four arcs is obtained. This algorithm is embedded into Algorithm 1 by executing once for each individual.

Pheromone Update
In Algorithm 1, the pheromone trails of the whole graph at the generation g is recorded by a pheromone matrix B g . The pheromone trail on each arc is updated at the end of each iteration by the following method where ∆τ g v,j is the pheromone increment on arc j with starting node v at the g-th iteration, P g j is the set of individuals who use the configurations corresponding to arc j at the gthe iteration, and ρ is the evaporation rate. Equation (40) indicates that the pheromone increment on an arc is determined by accumulated fitness improvements of individuals who passed this arc divided by that of all individuals. The pheromone trail B g+1 v,j on arc j with starting node v at the (g + 1)-th iteration is updated by pheromone increment ∆τ g v,j and pheromone trail B g v,j at the g-the iteration, as well as evaporation rate ρ in Equation (41). Further, to avoid premature convergence, the Max-Min ant system [48] is implemented in this study to limit the pheromone level on each arc within a range [0.1, 0.9].

Computational Experiments
To demonstrate the efficiency of ACODE on the proposed problem, simulation experiments are conducted in this section. All algorithms are coded in Matlab and run on a 64-bit Windows OS with Intel Core(TM) i5-8265U, 1.6 GHz, and 8 GB RAM.

Scenario Settings
We assume a set of scenarios in which three satellites are requested to observe four ground targets within 12 h (from 1 December 2020 14:00:00 to 2 December 2020 02:00:00, Beijing time). At the initial moment (i.e., 1 December 2020 14:00:00), each ground target is invisible to each satellite. Once an observation task is received, an appropriate satellite would be selected from these satellites to undertake orbital maneuvers to accomplish the task according to users' requirements. The initial orbital elements are displayed in Table 3, where the first column is the satellite ID and the other columns indicate semimajor axis a, inclination i, right ascension of the ascending node Ω, eccentricity e, argument of perigee ω, and mean anomaly M. The four ground targets are randomly located in low-latitude, mid-latitude, high-latitude, and higher-latitude areas, and their geographical information is summarized in Table 4. The maximum scanning angle of the satellite, maximum response time, minimum ground resolution, and maximum velocity increment predefined by users are set to 45 • , 12 h, 2 m, and 300 m/s, respectively. Moreover, the maximum number of fitness evaluations (FEs) of the algorithm is set to 50,000 and the evaporation rate ρ is set to 0.8 according to pre-experiments.   Table 4. Geographical information of ground targets.

Simulation Results
The ACODE is implemented to solve the three optimization models with different optimization objectives based on the generated scenarios. The experiment results are summarized in Table 5, in which the columns indicate scenarios, satellites selected to accomplish observation tasks, maneuver moment t m , velocity increments (∆v x and ∆v y ), and objective values ( f 1 , f 2 , and f 3 ) of the three optimization models. Particularly, the scenario index is composed of the ground target ID and optimization objective. For instance, T1O1 means in this scenario the satellites are requested to observe ground target 1 and the optimization objective is f 1 involved by the first optimization model. Here f 1 , f 2 , and f 3 are response time, ground resolution, and fuel consumption, respectively. Note that the minimum fuel consumption is represented by minimum velocity increment, as discussed in Section 2. Although only one objective is considered in each scenario, we provide the values of the other two objectives corresponding to the optimal solution of the scenario. The value of the optimized objective considered in each scenario is in boldface.
The results in Table 5 indicate that all scenarios can be well-addressed by ACODE. Furthermore, it can be observed that huge differences in objective values can be obtained if we execute the same observation task based on different optimization models. For example, T1O1, T1O2, and T1O3 are three scenarios in which the satellites are requested to observe ground target 1 with three different optimization objectives, respectively. The solution of T1O1 selects satellite 3 to execute the task and earns the minimum response time while yielding the poorest ground resolution compared with solutions of T1O2 and T1O3. More specifically, the solution of T1O1 decreases the response time by up to 84.44% and increases the ground resolution by up to 200.12% compared with the solutions of T1O2 and T1O3. Meanwhile, its fuel consumption is a little less than the solution of T1O2 that optimizes the ground resolution and much more than the solution of T1O3 that aims to find the minimum fuel consumption.  1 The values in boldface are optimized objectives in each scenario.

Algorithm Comparisons
To further demonstrate the superiority of ACODE, we compare it with three wellknown evolutionary algorithms in existing studies, i.e., EPSDE [33], CSO [34], and SLPSO [35]. Particularly, EPSDE is an ensemble-based DE algorithm, in which a pool of mutation strategies along with a pool of corresponding control parameters compete to produce offspring individuals. CSO is a competitive swarm optimizer inspired by particle swarm optimization. In CSO, a pairwise competition mechanism is used to update the position of the particle that loses the competition by learning from the winner. Similarly, SLPSO adopts social learning mechanisms for particle swarm optimization. Meanwhile, a dimensiondependent parameter control method is embedded into the SLPSO to ease the burden of parameter settings.
The comparison results are summarized in Table 6, in which the last four columns are best, worst, mean, and standard deviation values of three optimization objectives over 10 runs obtained by all algorithms. Note that the fuel consumption is represented by the value of velocity increment. For each scenario, the best results are in boldface. Wilcoxon rank-sum tests with a significance level of 0.05 are used for the significance tests. It can be found that ACODE significantly outperforms EPSDE, CSO, and SLPSO in almost all scenarios, in terms of response time, ground resolution, and fuel consumption. Especially, the superiority of ACODE is more significant when it optimizes orbital maneuvers for observing ground target 1 in terms of response time (scenario T1O1) and ground target 2 in terms of fuel consumption (scenario T2O3). It should be noted that EPSDE is similar to ACODE, as EPSDE ensembles a set of mutation strategies and corresponding control parameters in DE. Hence, EPSDE shows similar performance compared with ACODE for observing ground target 1 in terms of ground resolution (scenario T1O2). Nevertheless, ACODE is superior to EPSDE in other scenarios. The reasons can be twofold. First, EPSDE only ensembles mutation strategies and corresponding control parameters, while crossover strategies and corresponding control parameters that also can affect algorithm performance are not involved. On the contrary, ACODE considers both mutation strategies, crossover strategies, and their control parameters. Second, each component of EPSDE conducts the adaption independently while ACODE configures all components in a holistic manner, which is also the difference between ACODE and ensemble-based algorithms.

Experiments with Insufficient Satellite Resources
In the above sections, we calculate the solution of every satellite and select the most appropriate satellite out of three satellites to observe the ground target. The simulation results are obtained by the algorithm with sufficient satellite resources. However, since some satellites may be occupied by other tasks that cannot be interrupted when emergencies occur (i.e., some satellites may be infeasible for executing the observation task), it is interesting to investigate the impact of insufficient satellite resources on the orbital maneuver scheme and algorithm performance. Hence, this section analyzes the experiment results with different numbers of satellites based on the scenarios generated by removing satellites from the scenarios introduced in Section 4.1. Specifically, the one-satellite scenarios in this section preserve satellite 1, the two-satellite scenarios preserve satellite 1 and satellite 2, and the three-satellite scenarios are the same as before.
The simulation results are presented in Figure 5. Each figure indicates the simulation results for observing the same ground target, and the same color means the simulation results in the scenarios that consider the same optimization model. Moreover, to understand the trade-off among three objectives, we normalize all results into [0, 1], and a smaller value indicates a better solution in a direction. Since we generate the scenarios by removing satellites from the scenarios that already have solutions in Section 4.2, some scenarios would have the same solution as before. For example, the solution schemes of three scenarios that observe target 1 while optimizing the ground resolution with different numbers of satellites select satellite 1 to execute the task. Hence, the results of scenarios T1O2-onesatellite, T1O2-two-satellite, and T1O2-three-satellite are the same, as Figure 5a shows. On the other hand, other solutions indicate that the number of satellites significantly affects the algorithm results. For example, scenarios T1O1-one-satellite, T1O1-two-satellite, and T1O1-three-satellite select three different satellites to execute the task, respectively. To observe ground target 1 with the aim of optimizing response time, satellite 2 is selected in the two-satellite scenario and the response time is increased by 249.93% compared with the solution of the three-satellite scenario, as Figure 5a shows. Furthermore, it can be found that with the increase in the number of satellites, the value of the optimization objective that corresponds to each optimization model can be significantly improved. However, the trade-off results among the three objectives show that the improvement on one objective may not always promote the improvement of other objectives. For example, the ground resolution for observing ground target 2 is significantly improved as the number of satellites increases from 1 to 3, while the fuel consumption is still very high and the response time is even increased, as Figure 5b shows.

Conclusions
In this paper, we investigate the orbital maneuver optimization problem of Earth observation satellites oriented to emergency tasks. Based on the analysis of orbit coverage and dynamics, we propose three kinds of optimization models that aim to, respectively, optimize response time, ground resolution, and fuel consumption, to satisfy diverse user requirements. Meanwhile, we implement an adaptive differential evolution algorithm based on graph search to solve the proposed optimization problems, which is named ACODE. The main feature of ACODE is to form the key components of DE into a directed acyclic graph and adopt an ACO method to search for combinations of these components from the graph, thereby adaptively configuring reasonable components for DE. The key components considered in this paper include mutation strategies, crossover strategies, as well as their corresponding control parameters, both of which can affect the performance of DE.
Finally, computational experiments are conducted to verify the proposed three optimization models and ACODE. The simulation results show that all simulation scenarios that consider different optimization objectives can be well-addressed by ACODE. Comparison experiments are also carried out to demonstrate the superiority of ACODE on the proposed problem. The comparison results indicate that ACODE is superior to three well-known algorithms (i.e., EPSDE, CSO, and SLPSO). Further, we find that insufficient satellite resources would affect the efficiency of the orbital maneuver scheme and algorithm.
In future studies, we would like to investigate the multi-objective optimization algorithm that can optimize the three optimization objectives simultaneously for better decision making operations.