# Reconfigurable Satellite Constellation Design for Disaster Monitoring Using Physical Programming

Authors: Xiangyue He; Haiyang Li; Luyi Yang; Jian Zhao
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

Data collection by satellites during and after a natural disaster is of great significance. In this work, a reconfigurable satellite constellation is designed for disaster monitoring, and satellites in the constellation are made to fly directly overhead of the disaster site through orbital transfer. By analyzing the space geometry relations between satellite orbit and an arbitrary disaster site, a mathematical model for orbital transfer and overhead monitoring is established. Due to the unpredictability of disasters, target sites evenly spaced on the Earth are considered as all possible disaster scenarios, and the optimal reconfigurable constellation is designed with the intention to minimize total velocity increment, maximum and mean reconfiguration time, and standard deviation of reconfiguration times for all target sites. To deal with this multiobjective optimization, a physical programming method together with a genetic algorithm is employed. Numerical results are obtained through the optimization, and different observation modes of the reconfigurable constellation are analyzed by a specific case. Superiority of our design is demonstrated by comparing with the existing literature, and excellent observation performance of the reconfigurable constellation is demonstrated.


---

## Body

## 1. Introduction

A global or regional observation satellite constellation is used for some emergency or significant missions such as disaster monitoring, and usually, the observation pattern of a specific constellation is restricted by both on-board instruments and geometric configuration of satellites [1, 2]. However, natural disasters such as earthquakes, forest fires, or extreme weathers occur randomly, and usually, data collection is achieved only when satellites happen to fly over the disaster location [3, 4]. To get timely and adequate observation for the disaster area, orbital maneuvers of satellites are required. Therefore, configuration of a constellation needs to be adjusted due to the change of observation mission, and reconfiguration of the constellation is demanded.The satellite constellation reconfiguration is defined as “a deliberate change of the relative arrangements of satellites in a constellation by addition or subtraction of satellites and orbital maneuvering in order to achieve desired changes in coverage or capacity” [5]. Usually, optimization of constellation reconfiguration includes final constellation configuration design and optimal orbital transfer, and the purposes are to minimize the total fuel consumption as well as the reconfiguration time [6, 7]. Based on the satellite constellation reconfiguration, the concept of the reconfigurable constellation (ReCon) is proposed [8, 9]. Operation of a ReCon comprises the following two modes: global observation mode (GOM) for normal operations and regional observation mode (ROM) for contingent responses, and due to the flexibility of the ReCon, timely and adequate observation for uncertain areas of interest can be achieved [10].In recent years, plenty of research about constellation reconfiguration has been carried out. De Weck et al. [5] used the auction algorithm to optimize the reconfiguration of constellation, and the optimal orbital transfer strategy for all satellites was designed with the purpose of minimizing the total fuel consumption. Reconfiguration time is not considered as a design metric in their work. Hong et al. [11] introduced a concept and a baseline design of a dual-mode disaster monitoring constellation, and a normal mode as well as a disaster mode is considered. However, the constellation needs extra on-orbit propellant depots to provide plenty of fuel for the mode transformation. Chen et al. [12] investigated the observation for some particular targets on the ground when an emergent requirement arises in a short period and proposed a variable-size multiobjective differential evolution method to deal with the reconfiguration of constellation, where the main idea of the reconfiguration is to adjust the mean anomaly of satellites. The optimization is implemented with known initial constellations and target sites in their work. Fakoor et al. [13, 14] used the hybrid Invasive Weed Optimization/Particle Swarm Optimization algorithm to deal with the constellation reconfiguration. Optimal assignment of the satellites and optimal orbital transfer were combined in one step to simplify the reconfiguration process. Yet, reconfiguration time is not considered as an optimization objective in their work. Hoskins et al. [3] aimed at forest fire monitoring and introduced the L-shaped method and stochastic programing approach to design the initial constellation and maneuver sequences for reconfiguration. The direct overhead monitoring for the disaster site was achieved after the reconfiguration, and results proved that their design performs better than the existing observation satellites.As for regional or specific target observation, repeating ground track orbit is a usual choice, which has been proved to have better partial coverage properties than the unrepeated one [15, 16]. A ground track is considered as repeating or periodic if it repeats at a fixed time interval, and the interval is usually an integral multiple of one nodal day [17–20]. Hence, sites which are located on the ground track can be regularly observed. In addition, since the on-board sensors benefit from minimizing the distance between the satellite and the site to be monitored, an optimal way to achieve data collection for the disaster area is to make satellites fly directly over the target site [3]. Therefore, to achieve the timely and adequate observation for a natural disaster site, the repeating ground track orbit, of which the ground track passes the site directly, can be chosen as the target orbit for each satellite in the initial constellation. After the reconfiguration, each satellite in the final constellation is able to fly directly over the target site for data collection in a fixed period, and the period can be designed as one or more nodal days according to the requirement of a specific mission.Previous investigations mentioned above considered plenty of approaches to optimize the reconfiguration process and also the use of repeating orbit for partial coverage. However, most of them were carried out with a known initial constellation, a known final constellation, or a set of known observation locations, and few of them considered the reconfiguration for monitoring the unpredictable disaster site as well as the design of a reconfigurable initial constellation. Besides, most of the previous researches took fuel consumption as the optimization target, and few of them considered the optimization of reconfiguration time or tradeoff between them. In addition, none of the current constellations with the task of observing the Earth can achieve active maneuvers to collect data for the natural disaster or other significant events [3]. Therefore, the design of a ReCon according to the specific constraint of fuel consumption, reconfiguration time, and unpredicted disaster sites is of great significance.In this paper, we aim to design a ReCon, which can achieve timely and adequate observation for an arbitrary and unpredictable disaster site through orbital maneuvers. Satellites in the constellation flying directly over the disaster area for data collection are the purpose of the reconfiguration. Due to the unpredictability of disaster, target sites evenly spaced on the Earth are considered as all possible disaster scenarios. For the optimal ReCon, several optimization objectives are considered, which are (1) minimizing fuel consumption, which is equivalent to minimizing the total velocity increment during reconfiguration; (2) minimizing maximum and mean reconfiguration time to ensure timely observation for the target site; and (3) minimizing standard deviation of reconfiguration times to ensure similar reconfiguration ability for all target sites. To deal with this multiobjective optimization, a physical programming method together with a genetic algorithm is employed. The preferred design metrics in physical programming are considered according to the specific requirement or preference of the designer, which makes the design much more feasible, practical, and reasonable. Furthermore, GOM and ROM of the ReCon are analyzed to demonstrate the superior observation performance for normal global scan as well as some emergency disaster monitoring. As a whole, the main contributions of our work can be summarized as follows: firstly, a mathematical model for orbital transfer and overhead monitoring is established by analyzing the space geometry relations between satellite orbit and an arbitrary disaster site. Secondly, a ReCon for unpredictable disaster monitoring is designed, and several optimization objectives are considered to optimize the reconfiguration process as well as to ensure that the ReCon has similar reconfiguration ability to observe any possible disaster site. Thirdly, two observation modes are presented in the ReCon, and normal global scan as well as emergency disaster monitoring with excellent observation performance can be achieved, respectively.The rest of this paper is organized as follows: The detailed problem formulation process is described in Section2. An overview of the physical programming method is proposed in Section 3. Numerical results of a design example and comparison with existing literature are given in Section 4. Finally, conclusions are presented in Section 5.


## 2. Problem Formulation

Problem formulation is presented in this section. To simplify the formulation process, as well as to make the model easier to establish and understand, some reasonable and necessary assumptions ought to be made. Firstly, maneuvers for orbital transfer are impulsive, and this is reasonably accurate during the initial mission planning [19]. Secondly, the specific sensors are not modeled, since all sensors (e.g., visual spectrum cameras, infrared cameras, and radar sounders) benefit from minimizing the range between the satellite and the area being monitored [3]. Mathematic models for orbital transfer and overhead monitoring, as well as the description of the optimization problem for the optimal ReCon design, are presented as follows.


### 2.1. Model Overview

#### 2.1.1. Problem Description

The aim of our work is to design a ReCon for disaster monitoring. There are two observation modes for the ReCon, which are GOM for normal operations and ROM for contingent responses [10]. When an unpredictable disaster occurs, the observation model of the constellation can be converted to ROM through reconfiguration, and satellites in the constellation are able to fly overhead the disaster site for monitoring through orbital maneuvers. In this work, constellation with GOM is defined as initial constellation, while constellation with ROM is the target constellation, and initial constellation is designed to achieve optimal reconfiguration for monitoring any possible disaster. Features of the initial and target constellation, as well as the orbital transfer, can be described as follows.(1) Initial Constellation. For the initial constellation with GOM, the number of satellites is K, which is also the number of orbit planes. Eccentricity and argument of perigee are both set to be 0 for each orbit. Radius r1 and inclination i are set to be the same for each orbit. Right ascension of the ascending node (RAAN) and argument of latitude (AOL) of orbit kk=1,2,⋯,K are Ωk and uk, respectively. It should be mentioned that for all satellites, Ωk are chosen to be evenly distributed along the equator to ensure the uniformity of the observation, and for the first satellite, Ω1=0; thus, we have Ωk=2πk−1/K. uk is the angle between the ascending node and the satellite, and 0≤uk<2π. As a result, the design variables of the initial constellation are K, i, r1, and uk, where k=1,2,⋯,K.J2 perturbation of the Earth is considered in this work. For the initial orbit, assume that change rates of the RAAN, the argument of perigee, and mean anomaly are Ω̇1, ω̇1, and Ṁ1, respectively. n1 denotes the mean angular velocity of the satellite. Then, the nodal period Tp1 of the initial orbit can be defined as
(1)Tp1=2πn1+Ṁ1+ω̇1.Change ratesΩ̇1, ω̇1, and Ṁ1 can be expressed as
(2)ω̇1=ξn14−5sin2i,Ω̇1=−2ξn1cosi,Ṁ1=−ξn13sin2i−2,where ξ=3Re2J2/4r12, Re denotes the average radius of the Earth, and the value of J2 is 1.082626×10−3 [21].Also, for the nodal period of initial orbit, the corresponding angular velocityn1p, which is different from n1, is defined as
(3)np1=2πTp1.(2) Target Constellation. As for the target constellation with ROM, the number of orbit planes, the RAAN, and the inclination stay the same as the initial constellation, for only coplanar orbital transfer is considered. The type of each target orbit is considered to be the repeating ground track circular orbit, and radius r2 of the orbit is determined by the number of revolutions and number of nodal days in the repetition period. Also, circular orbit means constant distance between the satellite and the Earth, which can ensure the uniformity of observation.A ground track is considered as repeating or periodic if it repeats at a fixed time interval allowing a cyclic observation of the Earth [1]. Let Tp2 be the nodal period of the target orbit and Td be the nodal day. Np and Nd are positive integers, which are prime to one another. For the repeating ground track orbit, it should be guaranteed that the period Tp2 precisely matches the period Td, and this can be described as a mathematical condition, which is the following expression [20]:
(4)Tr=NpTp2=NdTd,where Tr is a fixed time interval, which is also the period of repetition.For the target repeating ground track orbit, assume that change rates of the RAAN, the argument of perigee, and mean anomaly areΩ̇2, ω̇2, and Ṁ2, respectively. n2 denotes the mean angular velocity of the satellite, ωe denotes the angular velocity of the Earth, and μ denotes the geocentric gravitational constant. Then, the period Td and Tp2 can be defined as
(5)Td=2πωe−Ω̇2,Tp2=2πn2+Ṁ2+ω̇2.Also, radiusr2 can be calculated by [1]
(6)r23.5+b1r22+b2=0,where
(7)b1=−μqωe,b2=3J2Re2μ2qωe1−4cos2i+qcosi,q=NpNd.For the nodal period of orbit, the corresponding angular velocitynp2, which is different from n2, is defined as
(8)np2=2πTp2.In this work, the repetition period of one nodal day is considered to ensure the high frequency of observation, which also means thatNd=1. Also, it is convenient to adjust our model for other kinds of repetition periods.It should be mentioned thatnp1 and np2 will be used as the angular velocity of the initial and target orbit in the following analysis.(3) Orbital Transfer. In this work, Hohmann transfer, as the optimal two-impulse coplanar circular orbital transfer, is introduced. The total velocity increment Δv of Hohmann transfer is determined by the radius of the initial circular orbit r1 and radius of the target circular orbit r2. Let nT=r2/r1, and μ denotes the geocentric gravitational constant; then, Δv can be expressed by
(9)Δv=μr12nTnT+11−1nT+1nT−1.The transfer timeΔTh is a half of the period of the elliptical transfer orbit, which is defined as
(10)ΔTh=πaT3μ,where aT is the semimajor axis of the transfer orbit and aT=r1+r2/2.


#### 2.1.2. Mathematical Modeling

Due to the unpredictability of disaster, target sites around the Earth’s surface in step of 1 degree are set as possible scenarios, and all satellites in the constellation are considered to monitor these target sites through orbital maneuvers. To illustrate the model mathematically, one arbitrary target site on the ground and an orbit with one satellite are considered as the basic scenario. The main task of this subsection is to calculate the specific transition time of each satellite monitoring each target site in the basic scenario and then to form the whole constellation reconfiguration process. It should be noted that transition time refers to the time from the mission beginning to the satellite finishing the transfer, while reconfiguration time represents the maximum transition time of all satellites in the constellation.The ascending pass of a satellite is first considered for the monitoring mission. The basic scenario is presented in Figure1. In the figure, OEXYZ is the J2000 inertial system. The intersection of the X-axis and the equator is O, and the intersection of the Z-axis and the Earth is P. Ground track of the orbit shown in the figure is defined on the nonrotational Earth in the inertial system, which is also the projection of the orbit. It is a circle and would rotate along the Z-axis at a regular angular velocity due to the perturbation. T is the target site.Figure 1
Monitoring during ascending pass.Suppose that the longitude and latitude of the target siteT are λ and φ, respectively, and the current time is assumed to be t0, which is also the mission start time. The corresponding right ascension and declination of the target site at t0 moment are λ0 and φ0, respectively, where 0≤λ0<2π and −π/2≤φ0≤π/2. As for the satellite, inclination of the orbit is i, RAAN at t0 moment is Ωt0, and eccentricity as well as argument of perigee are both set to be 0.In the inertial system, target siteT is rotating along the Z-axis at a regular angular velocity, which is equal to the Earth rotating angular velocity. Let i0=mini,π/2−i be the latitude limit of the ground track; then, only target sites which satisfy −i0≤φ0≤i0 can intersect the ground track and get the chance to be directly monitored by the satellite flying overhead, and hence, only this kind of target sites are considered in the following analysis.After a specific time, target siteT will first intersect the ground track at point T1 due to the rotation of the Earth. In this paper, this specific time is defined as the “minimum distance 1” between the target site and the ground track, which is marked as ΔTl1min, and is determined by TT1⌢, where the number “1” represents the ascending pass (the same below). Only after at least ΔTl1min then the satellite could possibly get the chance to fly directly over the target site T during the ascending pass. Obviously, the target site will intersect the ground track every single nodal day after the first intersection. Therefore, if the satellite cannot finish the orbital transfer before the target site intersects the ground track, then it will get another one nodal day time to wait and to perform the orbital transfer, and so on. Here, the number of nodal days is marked as D1, and it is defined as “waiting days” in this paper.For pointT1 on the ground track, the declination is φ0 and the right ascension is assumed to be λ1 at t0 moment. Make an orthodrome crossing point P and T and intersecting with the equator at point M. Also, make an orthodrome crossing point P and T1 and intersecting with the equator at point M1. Let u be the geocentric angle of NT1⌢. Then, in the spherical right triangle NT1M1, u and λ1 can be expressed as
(11)u=arcsinsinφ0sini,(12)λ1=arcsintanφ0tani+Ωt0,where i≠π/2. For the case that i=π/2, it is obvious that u=φ0 and λ1=Ωt0.Therefore, the defined “minimum distance 1”ΔTl1min can be calculated by
(13)ΔTl1min=λ1−λ02πωe−Ω̇,where the function X2π ensures that X falls into the interval of 0,2π. ωe is the angular velocity of the Earth, and Ω̇ is the change rate of the RAAN.Here, we takeΔTl1 as the time limit for the satellite to perform the orbital transfer. According to the analysis above, ΔTl1 can be expressed as
(14)ΔTl1=ΔTl1min+D1Td,where D1=0,1,2,⋯, is the number of waiting days. Td is the time interval of a nodal day.As for the descending pass, similar analysis can be carried out. As shown in Figure2, the “minimum distance 2” ΔTl2min can be determined by TT2⌢, where the number “2” represent the descending pass.Figure 2
Monitoring during descending pass.LetP2 be the midpoint of T1T2⌢. Make an orthodrome crossing point P and P2 and intersecting with the ground track and equator at point P1 and P3, respectively. Also, make an orthodrome crossing point P and T2 and intersecting with the equator at point M2. Declination φ0 and right ascension λ1 of point T1 are already obtained. For point T2 on the ground track, the declination is φ0 and the right ascension is assumed to be λ2 at t0 moment. Thus, for point P1 on the ground track, the declination is i0 and the right ascension is λ1+λ2/2. For point P2, the declination is φ0 and the right ascension is also λ1+λ2/2.Then, in the spherical right triangleNP3P1, it is easy to obtain that
(15)λ2=2arccosseci+Ωt0−λ1,where i≠π/2, and for the case that i=π/2, it is obvious that λ2=λ1+π.Therefore, the defined “minimum distance 2”ΔTl2min can be calculated by
(16)ΔTl2min=λ2−λ02πωe−Ω̇.Similarly, we takeΔTl2 as the time limit for the satellite to perform the orbital transfer. According to the analysis above, ΔTl2 can be expressed as
(17)ΔTl2=ΔTl2min+D2Td,where D2=0,1,2,⋯, is the number of waiting days.OnceΔTl1 and ΔTl2 are determined, the next task is to obtain the transition time of the satellite. In this paper, we define the required position of the satellite in the target orbit and the position of the satellite in the initial orbit as “target position” PT and “initial position” PI, respectively. Obviously, these two positions are changing over time, and the main idea of the transition is to wait for a suitable moment to achieve the Hohmann transfer.As illustrated in Figure3, r1 denotes the radius of the initial circular orbit and r2 denotes the radius of the target circular orbit. Let np1 and np2 be the angular velocity of the initial and target orbits, respectively, which are calculated by equations (3) and (8). To calculate the transition time, the direct overhead observation during the ascending pass of a satellite is analyzed below, and for the descending pass, analysis can be carried out in the same way.Figure 3
Orbital transfer.Att0 moment, AOL of the initial position is assumed to be u1, which is a given parameter. After ΔTl1, the target position is right above the ground point T1 presented in Figure 1, and the corresponding AOL is u. Therefore, at t0 moment, the AOL u2 of the target position can be expressed by
(18)u2=arcsinsinφ0sini−np2ΔTl12π.Then, the phase differenceΔu0 between the initial and target position at t0 moment can be obtained, which is
(19)Δu0=u1−u22π.As shown in Figure3, if Hohmann transfer is to be performed, the phase difference ΔuT of the initial and target position need to be restricted as
(20)ΔuT=np2ΔTh−π2π,where ΔTh is the transfer time defined by equation (10). Also, the target position is required to be PH, as illustrated in Figure 3.The difference betweenΔu0 and ΔuT is gradually changed due to the difference between np1 and np2. Define Δuc as the phase difference that needs to be reduced, then
(21)Δuc=Δu0−ΔuT2π.For the case thatr1>r2, the target position is “chasing” the initial position, until the target position is ΔuT behind the initial position, and then, the first maneuver can be implemented. Define the “chasing” time as ΔTc, which is also the time for waiting before the orbital transfer, then ΔTc can be expressed by
(22)ΔTc=Δucnp2−np1.Therefore, the total transition timeΔT1 can be expressed as
(23)ΔT1=ΔTc+ΔTh.For the case thatr1<r2, the analysis can be carried out in a similar way. Also, for the descending pass observation, similar analysis can be done, and the total transition time ΔT2 can be obtained.So far, time limitsΔTl1 and ΔTl2 and total transition times ΔT1 and ΔT2 are obtained, and the problem turns to be judging whether the satellite can achieve the orbital transfer within the time limit of ΔTl1 or ΔTl2. It should be mentioned that in order to ensure the quality of observation, the satellite needs a time interval to prepare for the monitoring mission after the last impulse, and the interval in this work is set to be a half of the target orbit period 2π/np2 [3].In conclusion, for the basic scenario that one satellite overhead monitors one arbitrary target site during ascending pass, the calculation process can be summarized as the following seven steps.Step 1.
Input target siteλ0,φ0, parameters t0, i, Ωt0, u1, r1, and r2 of initial and target orbits, and set D1=0.Step 2.
Calculate minimum distanceΔTl1min by equations (11)–(13).Step 3.
Calculate time limitΔTl1 by equation (14).Step 4.
Calculate total transition timeΔT1 by equations (18)–(23).Step 5.
Judge whetherΔT1+π/np2<ΔTl1. If the inequality above holds, then go to Step 6. Otherwise, D1=D1+1 and return back to Step 3.Step 6.
Calculate total velocity incrementΔv by equation (9).Step 7.
OutputΔv, ΔT1, and D1.Similarly, analysis for descending pass can be carried out, and the mathematical model for one satellite monitoring one target can be concluded as Figure4.Figure 4
Mathematical model for one satellite monitoring one target.As a result, when consideringK satellites in the constellation maneuvering for monitoring S target sites around the Earth’s surface, the same analysis given in Figure 4 can be carried out individually, and S×K sets of Δv, ΔT, and D can be obtained. Finally, the total reconfiguration model is established.


### 2.2. Optimization Problem

Constellation design is generally a complicated optimization problem with plenty of design variables; especially when there are several design objectives, the computational process can be very time-consuming [22] [23]. In this paper, the optimal design is simplified by restricting some design variables. Reasons and details for the simplification are analyzed as follows.Firstly, radiusr1 or altitude h1 and inclination i of each initial circular orbit are the same. The RAAN Ω of each orbit is evenly distributed along the equator, and the inclination i is determined by the highest latitude of the region that needs to be monitored.Secondly, radiusr2 or altitude h2 of the target orbit is determined by the number of revolutions in the period of repetition, and it is a given reference for the design of the altitude h1 of the initial orbit.Thirdly, the repetition period of the target orbit is one nodal day, and there is only one satellite in each orbit plane. The number of satellitesK, which is also the number of orbit planes, is determined by the specific observation frequency requirements and the economic conditions of the designer.Therefore, once the number of satellitesK, the inclination i, and the altitude h2 of the target orbit are determined by the designer, design variables would be reduced to the altitude h1 of the initial orbit and the AOL uk each initial orbit, where k=1,2,⋯,K. For the first satellite, it is reasonable to set u1=0, since only phase differences between each satellite matter.Based on the analysis above, our design can be concluded as designing the optimal altitude of the initial orbit and phase differences between each orbit plane, and the purpose is to minimize the fuel consumption, maximum reconfiguration time, mean transition time of all satellites for all targets, the standard deviation of the reconfiguration time for all targets, and maximum waiting days for all target sites.Before formulating the optimization problem, symbols of design variables and design metrics need to be defined. Assume that there areS (S≥1) target sites on the Earth’s surface and K (K≥2) satellites in the constellation. Design variables are the altitude h1 of the initial orbit and the phase differences Δuk between satellite k and the first satellite, where 0≤Δuk≤2π for k=2,⋯,K. For an arbitrary target site s (s=1,2,⋯,S) and an arbitrary satellite k (k=1,2,⋯,K), the impulse Δv, transition time ΔTs,k, and waiting days Ds,k can be obtained based on the mathematical model described before. It should be mentioned that for a constellation, the reconfiguration time is defined as the maximum transition time of all satellites. Therefore, reconfiguration time ΔTs of the constellation for target site s can be defined as
(24)ΔTs=maxΔTs,1,ΔTs,2,⋯,ΔTs,K.Also, the maximum waiting days of the constellation for the target sites can be expressed as
(25)Ds=maxDs,1,Ds,2,⋯,Ds,K.In conclusion, the optimization problem for ReCon design can be formulated as(26)find:minfix,i=1,2,3,4,5subjectto:xLB≤x≤xUB,where design metrics fi are defined as
(27)f1=Δv,f2=maxΔT1,ΔT2,⋯,ΔTS,f3=1S1K∑s=1S∑k=1KΔTs,k,f4=∑s=1SΔTs−ΔT¯2S−1,f5=maxD1,D2,⋯,DS,where ΔT¯ denotes the mean value of reconfiguration times for all target sites, which can be expressed by
(28)ΔT¯=1S∑s=1SΔTs.x denotes the vector of design variables, and it is expressed as
(29)x=h1,Δu2,Δu3,⋯,ΔuK,xLB and xUB are the lower and upper bound, and they are determined by the designer. For the altitude h1, hmin and hmax are defined as the lower and upper bound, while for the phase differences Δuk (k=2,⋯,K), 0 and 2π are the lower and upper bound. Thus, xLB and xUB can be expressed as
(30)xLB=hmin,0,0,⋯,0,xUB=hmax,2π,2π,⋯,2π,where numbers of “0” and “2π” are both K−1.


## 2.1. Model Overview

### 2.1.1. Problem Description

The aim of our work is to design a ReCon for disaster monitoring. There are two observation modes for the ReCon, which are GOM for normal operations and ROM for contingent responses [10]. When an unpredictable disaster occurs, the observation model of the constellation can be converted to ROM through reconfiguration, and satellites in the constellation are able to fly overhead the disaster site for monitoring through orbital maneuvers. In this work, constellation with GOM is defined as initial constellation, while constellation with ROM is the target constellation, and initial constellation is designed to achieve optimal reconfiguration for monitoring any possible disaster. Features of the initial and target constellation, as well as the orbital transfer, can be described as follows.(1) Initial Constellation. For the initial constellation with GOM, the number of satellites is K, which is also the number of orbit planes. Eccentricity and argument of perigee are both set to be 0 for each orbit. Radius r1 and inclination i are set to be the same for each orbit. Right ascension of the ascending node (RAAN) and argument of latitude (AOL) of orbit kk=1,2,⋯,K are Ωk and uk, respectively. It should be mentioned that for all satellites, Ωk are chosen to be evenly distributed along the equator to ensure the uniformity of the observation, and for the first satellite, Ω1=0; thus, we have Ωk=2πk−1/K. uk is the angle between the ascending node and the satellite, and 0≤uk<2π. As a result, the design variables of the initial constellation are K, i, r1, and uk, where k=1,2,⋯,K.J2 perturbation of the Earth is considered in this work. For the initial orbit, assume that change rates of the RAAN, the argument of perigee, and mean anomaly are Ω̇1, ω̇1, and Ṁ1, respectively. n1 denotes the mean angular velocity of the satellite. Then, the nodal period Tp1 of the initial orbit can be defined as
(1)Tp1=2πn1+Ṁ1+ω̇1.Change ratesΩ̇1, ω̇1, and Ṁ1 can be expressed as
(2)ω̇1=ξn14−5sin2i,Ω̇1=−2ξn1cosi,Ṁ1=−ξn13sin2i−2,where ξ=3Re2J2/4r12, Re denotes the average radius of the Earth, and the value of J2 is 1.082626×10−3 [21].Also, for the nodal period of initial orbit, the corresponding angular velocityn1p, which is different from n1, is defined as
(3)np1=2πTp1.(2) Target Constellation. As for the target constellation with ROM, the number of orbit planes, the RAAN, and the inclination stay the same as the initial constellation, for only coplanar orbital transfer is considered. The type of each target orbit is considered to be the repeating ground track circular orbit, and radius r2 of the orbit is determined by the number of revolutions and number of nodal days in the repetition period. Also, circular orbit means constant distance between the satellite and the Earth, which can ensure the uniformity of observation.A ground track is considered as repeating or periodic if it repeats at a fixed time interval allowing a cyclic observation of the Earth [1]. Let Tp2 be the nodal period of the target orbit and Td be the nodal day. Np and Nd are positive integers, which are prime to one another. For the repeating ground track orbit, it should be guaranteed that the period Tp2 precisely matches the period Td, and this can be described as a mathematical condition, which is the following expression [20]:
(4)Tr=NpTp2=NdTd,where Tr is a fixed time interval, which is also the period of repetition.For the target repeating ground track orbit, assume that change rates of the RAAN, the argument of perigee, and mean anomaly areΩ̇2, ω̇2, and Ṁ2, respectively. n2 denotes the mean angular velocity of the satellite, ωe denotes the angular velocity of the Earth, and μ denotes the geocentric gravitational constant. Then, the period Td and Tp2 can be defined as
(5)Td=2πωe−Ω̇2,Tp2=2πn2+Ṁ2+ω̇2.Also, radiusr2 can be calculated by [1]
(6)r23.5+b1r22+b2=0,where
(7)b1=−μqωe,b2=3J2Re2μ2qωe1−4cos2i+qcosi,q=NpNd.For the nodal period of orbit, the corresponding angular velocitynp2, which is different from n2, is defined as
(8)np2=2πTp2.In this work, the repetition period of one nodal day is considered to ensure the high frequency of observation, which also means thatNd=1. Also, it is convenient to adjust our model for other kinds of repetition periods.It should be mentioned thatnp1 and np2 will be used as the angular velocity of the initial and target orbit in the following analysis.(3) Orbital Transfer. In this work, Hohmann transfer, as the optimal two-impulse coplanar circular orbital transfer, is introduced. The total velocity increment Δv of Hohmann transfer is determined by the radius of the initial circular orbit r1 and radius of the target circular orbit r2. Let nT=r2/r1, and μ denotes the geocentric gravitational constant; then, Δv can be expressed by
(9)Δv=μr12nTnT+11−1nT+1nT−1.The transfer timeΔTh is a half of the period of the elliptical transfer orbit, which is defined as
(10)ΔTh=πaT3μ,where aT is the semimajor axis of the transfer orbit and aT=r1+r2/2.


### 2.1.2. Mathematical Modeling

Due to the unpredictability of disaster, target sites around the Earth’s surface in step of 1 degree are set as possible scenarios, and all satellites in the constellation are considered to monitor these target sites through orbital maneuvers. To illustrate the model mathematically, one arbitrary target site on the ground and an orbit with one satellite are considered as the basic scenario. The main task of this subsection is to calculate the specific transition time of each satellite monitoring each target site in the basic scenario and then to form the whole constellation reconfiguration process. It should be noted that transition time refers to the time from the mission beginning to the satellite finishing the transfer, while reconfiguration time represents the maximum transition time of all satellites in the constellation.The ascending pass of a satellite is first considered for the monitoring mission. The basic scenario is presented in Figure1. In the figure, OEXYZ is the J2000 inertial system. The intersection of the X-axis and the equator is O, and the intersection of the Z-axis and the Earth is P. Ground track of the orbit shown in the figure is defined on the nonrotational Earth in the inertial system, which is also the projection of the orbit. It is a circle and would rotate along the Z-axis at a regular angular velocity due to the perturbation. T is the target site.Figure 1
Monitoring during ascending pass.Suppose that the longitude and latitude of the target siteT are λ and φ, respectively, and the current time is assumed to be t0, which is also the mission start time. The corresponding right ascension and declination of the target site at t0 moment are λ0 and φ0, respectively, where 0≤λ0<2π and −π/2≤φ0≤π/2. As for the satellite, inclination of the orbit is i, RAAN at t0 moment is Ωt0, and eccentricity as well as argument of perigee are both set to be 0.In the inertial system, target siteT is rotating along the Z-axis at a regular angular velocity, which is equal to the Earth rotating angular velocity. Let i0=mini,π/2−i be the latitude limit of the ground track; then, only target sites which satisfy −i0≤φ0≤i0 can intersect the ground track and get the chance to be directly monitored by the satellite flying overhead, and hence, only this kind of target sites are considered in the following analysis.After a specific time, target siteT will first intersect the ground track at point T1 due to the rotation of the Earth. In this paper, this specific time is defined as the “minimum distance 1” between the target site and the ground track, which is marked as ΔTl1min, and is determined by TT1⌢, where the number “1” represents the ascending pass (the same below). Only after at least ΔTl1min then the satellite could possibly get the chance to fly directly over the target site T during the ascending pass. Obviously, the target site will intersect the ground track every single nodal day after the first intersection. Therefore, if the satellite cannot finish the orbital transfer before the target site intersects the ground track, then it will get another one nodal day time to wait and to perform the orbital transfer, and so on. Here, the number of nodal days is marked as D1, and it is defined as “waiting days” in this paper.For pointT1 on the ground track, the declination is φ0 and the right ascension is assumed to be λ1 at t0 moment. Make an orthodrome crossing point P and T and intersecting with the equator at point M. Also, make an orthodrome crossing point P and T1 and intersecting with the equator at point M1. Let u be the geocentric angle of NT1⌢. Then, in the spherical right triangle NT1M1, u and λ1 can be expressed as
(11)u=arcsinsinφ0sini,(12)λ1=arcsintanφ0tani+Ωt0,where i≠π/2. For the case that i=π/2, it is obvious that u=φ0 and λ1=Ωt0.Therefore, the defined “minimum distance 1”ΔTl1min can be calculated by
(13)ΔTl1min=λ1−λ02πωe−Ω̇,where the function X2π ensures that X falls into the interval of 0,2π. ωe is the angular velocity of the Earth, and Ω̇ is the change rate of the RAAN.Here, we takeΔTl1 as the time limit for the satellite to perform the orbital transfer. According to the analysis above, ΔTl1 can be expressed as
(14)ΔTl1=ΔTl1min+D1Td,where D1=0,1,2,⋯, is the number of waiting days. Td is the time interval of a nodal day.As for the descending pass, similar analysis can be carried out. As shown in Figure2, the “minimum distance 2” ΔTl2min can be determined by TT2⌢, where the number “2” represent the descending pass.Figure 2
Monitoring during descending pass.LetP2 be the midpoint of T1T2⌢. Make an orthodrome crossing point P and P2 and intersecting with the ground track and equator at point P1 and P3, respectively. Also, make an orthodrome crossing point P and T2 and intersecting with the equator at point M2. Declination φ0 and right ascension λ1 of point T1 are already obtained. For point T2 on the ground track, the declination is φ0 and the right ascension is assumed to be λ2 at t0 moment. Thus, for point P1 on the ground track, the declination is i0 and the right ascension is λ1+λ2/2. For point P2, the declination is φ0 and the right ascension is also λ1+λ2/2.Then, in the spherical right triangleNP3P1, it is easy to obtain that
(15)λ2=2arccosseci+Ωt0−λ1,where i≠π/2, and for the case that i=π/2, it is obvious that λ2=λ1+π.Therefore, the defined “minimum distance 2”ΔTl2min can be calculated by
(16)ΔTl2min=λ2−λ02πωe−Ω̇.Similarly, we takeΔTl2 as the time limit for the satellite to perform the orbital transfer. According to the analysis above, ΔTl2 can be expressed as
(17)ΔTl2=ΔTl2min+D2Td,where D2=0,1,2,⋯, is the number of waiting days.OnceΔTl1 and ΔTl2 are determined, the next task is to obtain the transition time of the satellite. In this paper, we define the required position of the satellite in the target orbit and the position of the satellite in the initial orbit as “target position” PT and “initial position” PI, respectively. Obviously, these two positions are changing over time, and the main idea of the transition is to wait for a suitable moment to achieve the Hohmann transfer.As illustrated in Figure3, r1 denotes the radius of the initial circular orbit and r2 denotes the radius of the target circular orbit. Let np1 and np2 be the angular velocity of the initial and target orbits, respectively, which are calculated by equations (3) and (8). To calculate the transition time, the direct overhead observation during the ascending pass of a satellite is analyzed below, and for the descending pass, analysis can be carried out in the same way.Figure 3
Orbital transfer.Att0 moment, AOL of the initial position is assumed to be u1, which is a given parameter. After ΔTl1, the target position is right above the ground point T1 presented in Figure 1, and the corresponding AOL is u. Therefore, at t0 moment, the AOL u2 of the target position can be expressed by
(18)u2=arcsinsinφ0sini−np2ΔTl12π.Then, the phase differenceΔu0 between the initial and target position at t0 moment can be obtained, which is
(19)Δu0=u1−u22π.As shown in Figure3, if Hohmann transfer is to be performed, the phase difference ΔuT of the initial and target position need to be restricted as
(20)ΔuT=np2ΔTh−π2π,where ΔTh is the transfer time defined by equation (10). Also, the target position is required to be PH, as illustrated in Figure 3.The difference betweenΔu0 and ΔuT is gradually changed due to the difference between np1 and np2. Define Δuc as the phase difference that needs to be reduced, then
(21)Δuc=Δu0−ΔuT2π.For the case thatr1>r2, the target position is “chasing” the initial position, until the target position is ΔuT behind the initial position, and then, the first maneuver can be implemented. Define the “chasing” time as ΔTc, which is also the time for waiting before the orbital transfer, then ΔTc can be expressed by
(22)ΔTc=Δucnp2−np1.Therefore, the total transition timeΔT1 can be expressed as
(23)ΔT1=ΔTc+ΔTh.For the case thatr1<r2, the analysis can be carried out in a similar way. Also, for the descending pass observation, similar analysis can be done, and the total transition time ΔT2 can be obtained.So far, time limitsΔTl1 and ΔTl2 and total transition times ΔT1 and ΔT2 are obtained, and the problem turns to be judging whether the satellite can achieve the orbital transfer within the time limit of ΔTl1 or ΔTl2. It should be mentioned that in order to ensure the quality of observation, the satellite needs a time interval to prepare for the monitoring mission after the last impulse, and the interval in this work is set to be a half of the target orbit period 2π/np2 [3].In conclusion, for the basic scenario that one satellite overhead monitors one arbitrary target site during ascending pass, the calculation process can be summarized as the following seven steps.Step 1.
Input target siteλ0,φ0, parameters t0, i, Ωt0, u1, r1, and r2 of initial and target orbits, and set D1=0.Step 2.
Calculate minimum distanceΔTl1min by equations (11)–(13).Step 3.
Calculate time limitΔTl1 by equation (14).Step 4.
Calculate total transition timeΔT1 by equations (18)–(23).Step 5.
Judge whetherΔT1+π/np2<ΔTl1. If the inequality above holds, then go to Step 6. Otherwise, D1=D1+1 and return back to Step 3.Step 6.
Calculate total velocity incrementΔv by equation (9).Step 7.
OutputΔv, ΔT1, and D1.Similarly, analysis for descending pass can be carried out, and the mathematical model for one satellite monitoring one target can be concluded as Figure4.Figure 4
Mathematical model for one satellite monitoring one target.As a result, when consideringK satellites in the constellation maneuvering for monitoring S target sites around the Earth’s surface, the same analysis given in Figure 4 can be carried out individually, and S×K sets of Δv, ΔT, and D can be obtained. Finally, the total reconfiguration model is established.


## 2.1.1. Problem Description

The aim of our work is to design a ReCon for disaster monitoring. There are two observation modes for the ReCon, which are GOM for normal operations and ROM for contingent responses [10]. When an unpredictable disaster occurs, the observation model of the constellation can be converted to ROM through reconfiguration, and satellites in the constellation are able to fly overhead the disaster site for monitoring through orbital maneuvers. In this work, constellation with GOM is defined as initial constellation, while constellation with ROM is the target constellation, and initial constellation is designed to achieve optimal reconfiguration for monitoring any possible disaster. Features of the initial and target constellation, as well as the orbital transfer, can be described as follows.(1) Initial Constellation. For the initial constellation with GOM, the number of satellites is K, which is also the number of orbit planes. Eccentricity and argument of perigee are both set to be 0 for each orbit. Radius r1 and inclination i are set to be the same for each orbit. Right ascension of the ascending node (RAAN) and argument of latitude (AOL) of orbit kk=1,2,⋯,K are Ωk and uk, respectively. It should be mentioned that for all satellites, Ωk are chosen to be evenly distributed along the equator to ensure the uniformity of the observation, and for the first satellite, Ω1=0; thus, we have Ωk=2πk−1/K. uk is the angle between the ascending node and the satellite, and 0≤uk<2π. As a result, the design variables of the initial constellation are K, i, r1, and uk, where k=1,2,⋯,K.J2 perturbation of the Earth is considered in this work. For the initial orbit, assume that change rates of the RAAN, the argument of perigee, and mean anomaly are Ω̇1, ω̇1, and Ṁ1, respectively. n1 denotes the mean angular velocity of the satellite. Then, the nodal period Tp1 of the initial orbit can be defined as
(1)Tp1=2πn1+Ṁ1+ω̇1.Change ratesΩ̇1, ω̇1, and Ṁ1 can be expressed as
(2)ω̇1=ξn14−5sin2i,Ω̇1=−2ξn1cosi,Ṁ1=−ξn13sin2i−2,where ξ=3Re2J2/4r12, Re denotes the average radius of the Earth, and the value of J2 is 1.082626×10−3 [21].Also, for the nodal period of initial orbit, the corresponding angular velocityn1p, which is different from n1, is defined as
(3)np1=2πTp1.(2) Target Constellation. As for the target constellation with ROM, the number of orbit planes, the RAAN, and the inclination stay the same as the initial constellation, for only coplanar orbital transfer is considered. The type of each target orbit is considered to be the repeating ground track circular orbit, and radius r2 of the orbit is determined by the number of revolutions and number of nodal days in the repetition period. Also, circular orbit means constant distance between the satellite and the Earth, which can ensure the uniformity of observation.A ground track is considered as repeating or periodic if it repeats at a fixed time interval allowing a cyclic observation of the Earth [1]. Let Tp2 be the nodal period of the target orbit and Td be the nodal day. Np and Nd are positive integers, which are prime to one another. For the repeating ground track orbit, it should be guaranteed that the period Tp2 precisely matches the period Td, and this can be described as a mathematical condition, which is the following expression [20]:
(4)Tr=NpTp2=NdTd,where Tr is a fixed time interval, which is also the period of repetition.For the target repeating ground track orbit, assume that change rates of the RAAN, the argument of perigee, and mean anomaly areΩ̇2, ω̇2, and Ṁ2, respectively. n2 denotes the mean angular velocity of the satellite, ωe denotes the angular velocity of the Earth, and μ denotes the geocentric gravitational constant. Then, the period Td and Tp2 can be defined as
(5)Td=2πωe−Ω̇2,Tp2=2πn2+Ṁ2+ω̇2.Also, radiusr2 can be calculated by [1]
(6)r23.5+b1r22+b2=0,where
(7)b1=−μqωe,b2=3J2Re2μ2qωe1−4cos2i+qcosi,q=NpNd.For the nodal period of orbit, the corresponding angular velocitynp2, which is different from n2, is defined as
(8)np2=2πTp2.In this work, the repetition period of one nodal day is considered to ensure the high frequency of observation, which also means thatNd=1. Also, it is convenient to adjust our model for other kinds of repetition periods.It should be mentioned thatnp1 and np2 will be used as the angular velocity of the initial and target orbit in the following analysis.(3) Orbital Transfer. In this work, Hohmann transfer, as the optimal two-impulse coplanar circular orbital transfer, is introduced. The total velocity increment Δv of Hohmann transfer is determined by the radius of the initial circular orbit r1 and radius of the target circular orbit r2. Let nT=r2/r1, and μ denotes the geocentric gravitational constant; then, Δv can be expressed by
(9)Δv=μr12nTnT+11−1nT+1nT−1.The transfer timeΔTh is a half of the period of the elliptical transfer orbit, which is defined as
(10)ΔTh=πaT3μ,where aT is the semimajor axis of the transfer orbit and aT=r1+r2/2.


## 2.1.2. Mathematical Modeling

Due to the unpredictability of disaster, target sites around the Earth’s surface in step of 1 degree are set as possible scenarios, and all satellites in the constellation are considered to monitor these target sites through orbital maneuvers. To illustrate the model mathematically, one arbitrary target site on the ground and an orbit with one satellite are considered as the basic scenario. The main task of this subsection is to calculate the specific transition time of each satellite monitoring each target site in the basic scenario and then to form the whole constellation reconfiguration process. It should be noted that transition time refers to the time from the mission beginning to the satellite finishing the transfer, while reconfiguration time represents the maximum transition time of all satellites in the constellation.The ascending pass of a satellite is first considered for the monitoring mission. The basic scenario is presented in Figure1. In the figure, OEXYZ is the J2000 inertial system. The intersection of the X-axis and the equator is O, and the intersection of the Z-axis and the Earth is P. Ground track of the orbit shown in the figure is defined on the nonrotational Earth in the inertial system, which is also the projection of the orbit. It is a circle and would rotate along the Z-axis at a regular angular velocity due to the perturbation. T is the target site.Figure 1
Monitoring during ascending pass.Suppose that the longitude and latitude of the target siteT are λ and φ, respectively, and the current time is assumed to be t0, which is also the mission start time. The corresponding right ascension and declination of the target site at t0 moment are λ0 and φ0, respectively, where 0≤λ0<2π and −π/2≤φ0≤π/2. As for the satellite, inclination of the orbit is i, RAAN at t0 moment is Ωt0, and eccentricity as well as argument of perigee are both set to be 0.In the inertial system, target siteT is rotating along the Z-axis at a regular angular velocity, which is equal to the Earth rotating angular velocity. Let i0=mini,π/2−i be the latitude limit of the ground track; then, only target sites which satisfy −i0≤φ0≤i0 can intersect the ground track and get the chance to be directly monitored by the satellite flying overhead, and hence, only this kind of target sites are considered in the following analysis.After a specific time, target siteT will first intersect the ground track at point T1 due to the rotation of the Earth. In this paper, this specific time is defined as the “minimum distance 1” between the target site and the ground track, which is marked as ΔTl1min, and is determined by TT1⌢, where the number “1” represents the ascending pass (the same below). Only after at least ΔTl1min then the satellite could possibly get the chance to fly directly over the target site T during the ascending pass. Obviously, the target site will intersect the ground track every single nodal day after the first intersection. Therefore, if the satellite cannot finish the orbital transfer before the target site intersects the ground track, then it will get another one nodal day time to wait and to perform the orbital transfer, and so on. Here, the number of nodal days is marked as D1, and it is defined as “waiting days” in this paper.For pointT1 on the ground track, the declination is φ0 and the right ascension is assumed to be λ1 at t0 moment. Make an orthodrome crossing point P and T and intersecting with the equator at point M. Also, make an orthodrome crossing point P and T1 and intersecting with the equator at point M1. Let u be the geocentric angle of NT1⌢. Then, in the spherical right triangle NT1M1, u and λ1 can be expressed as
(11)u=arcsinsinφ0sini,(12)λ1=arcsintanφ0tani+Ωt0,where i≠π/2. For the case that i=π/2, it is obvious that u=φ0 and λ1=Ωt0.Therefore, the defined “minimum distance 1”ΔTl1min can be calculated by
(13)ΔTl1min=λ1−λ02πωe−Ω̇,where the function X2π ensures that X falls into the interval of 0,2π. ωe is the angular velocity of the Earth, and Ω̇ is the change rate of the RAAN.Here, we takeΔTl1 as the time limit for the satellite to perform the orbital transfer. According to the analysis above, ΔTl1 can be expressed as
(14)ΔTl1=ΔTl1min+D1Td,where D1=0,1,2,⋯, is the number of waiting days. Td is the time interval of a nodal day.As for the descending pass, similar analysis can be carried out. As shown in Figure2, the “minimum distance 2” ΔTl2min can be determined by TT2⌢, where the number “2” represent the descending pass.Figure 2
Monitoring during descending pass.LetP2 be the midpoint of T1T2⌢. Make an orthodrome crossing point P and P2 and intersecting with the ground track and equator at point P1 and P3, respectively. Also, make an orthodrome crossing point P and T2 and intersecting with the equator at point M2. Declination φ0 and right ascension λ1 of point T1 are already obtained. For point T2 on the ground track, the declination is φ0 and the right ascension is assumed to be λ2 at t0 moment. Thus, for point P1 on the ground track, the declination is i0 and the right ascension is λ1+λ2/2. For point P2, the declination is φ0 and the right ascension is also λ1+λ2/2.Then, in the spherical right triangleNP3P1, it is easy to obtain that
(15)λ2=2arccosseci+Ωt0−λ1,where i≠π/2, and for the case that i=π/2, it is obvious that λ2=λ1+π.Therefore, the defined “minimum distance 2”ΔTl2min can be calculated by
(16)ΔTl2min=λ2−λ02πωe−Ω̇.Similarly, we takeΔTl2 as the time limit for the satellite to perform the orbital transfer. According to the analysis above, ΔTl2 can be expressed as
(17)ΔTl2=ΔTl2min+D2Td,where D2=0,1,2,⋯, is the number of waiting days.OnceΔTl1 and ΔTl2 are determined, the next task is to obtain the transition time of the satellite. In this paper, we define the required position of the satellite in the target orbit and the position of the satellite in the initial orbit as “target position” PT and “initial position” PI, respectively. Obviously, these two positions are changing over time, and the main idea of the transition is to wait for a suitable moment to achieve the Hohmann transfer.As illustrated in Figure3, r1 denotes the radius of the initial circular orbit and r2 denotes the radius of the target circular orbit. Let np1 and np2 be the angular velocity of the initial and target orbits, respectively, which are calculated by equations (3) and (8). To calculate the transition time, the direct overhead observation during the ascending pass of a satellite is analyzed below, and for the descending pass, analysis can be carried out in the same way.Figure 3
Orbital transfer.Att0 moment, AOL of the initial position is assumed to be u1, which is a given parameter. After ΔTl1, the target position is right above the ground point T1 presented in Figure 1, and the corresponding AOL is u. Therefore, at t0 moment, the AOL u2 of the target position can be expressed by
(18)u2=arcsinsinφ0sini−np2ΔTl12π.Then, the phase differenceΔu0 between the initial and target position at t0 moment can be obtained, which is
(19)Δu0=u1−u22π.As shown in Figure3, if Hohmann transfer is to be performed, the phase difference ΔuT of the initial and target position need to be restricted as
(20)ΔuT=np2ΔTh−π2π,where ΔTh is the transfer time defined by equation (10). Also, the target position is required to be PH, as illustrated in Figure 3.The difference betweenΔu0 and ΔuT is gradually changed due to the difference between np1 and np2. Define Δuc as the phase difference that needs to be reduced, then
(21)Δuc=Δu0−ΔuT2π.For the case thatr1>r2, the target position is “chasing” the initial position, until the target position is ΔuT behind the initial position, and then, the first maneuver can be implemented. Define the “chasing” time as ΔTc, which is also the time for waiting before the orbital transfer, then ΔTc can be expressed by
(22)ΔTc=Δucnp2−np1.Therefore, the total transition timeΔT1 can be expressed as
(23)ΔT1=ΔTc+ΔTh.For the case thatr1<r2, the analysis can be carried out in a similar way. Also, for the descending pass observation, similar analysis can be done, and the total transition time ΔT2 can be obtained.So far, time limitsΔTl1 and ΔTl2 and total transition times ΔT1 and ΔT2 are obtained, and the problem turns to be judging whether the satellite can achieve the orbital transfer within the time limit of ΔTl1 or ΔTl2. It should be mentioned that in order to ensure the quality of observation, the satellite needs a time interval to prepare for the monitoring mission after the last impulse, and the interval in this work is set to be a half of the target orbit period 2π/np2 [3].In conclusion, for the basic scenario that one satellite overhead monitors one arbitrary target site during ascending pass, the calculation process can be summarized as the following seven steps.Step 1.
Input target siteλ0,φ0, parameters t0, i, Ωt0, u1, r1, and r2 of initial and target orbits, and set D1=0.Step 2.
Calculate minimum distanceΔTl1min by equations (11)–(13).Step 3.
Calculate time limitΔTl1 by equation (14).Step 4.
Calculate total transition timeΔT1 by equations (18)–(23).Step 5.
Judge whetherΔT1+π/np2<ΔTl1. If the inequality above holds, then go to Step 6. Otherwise, D1=D1+1 and return back to Step 3.Step 6.
Calculate total velocity incrementΔv by equation (9).Step 7.
OutputΔv, ΔT1, and D1.Similarly, analysis for descending pass can be carried out, and the mathematical model for one satellite monitoring one target can be concluded as Figure4.Figure 4
Mathematical model for one satellite monitoring one target.As a result, when consideringK satellites in the constellation maneuvering for monitoring S target sites around the Earth’s surface, the same analysis given in Figure 4 can be carried out individually, and S×K sets of Δv, ΔT, and D can be obtained. Finally, the total reconfiguration model is established.


## 2.2. Optimization Problem

Constellation design is generally a complicated optimization problem with plenty of design variables; especially when there are several design objectives, the computational process can be very time-consuming [22] [23]. In this paper, the optimal design is simplified by restricting some design variables. Reasons and details for the simplification are analyzed as follows.Firstly, radiusr1 or altitude h1 and inclination i of each initial circular orbit are the same. The RAAN Ω of each orbit is evenly distributed along the equator, and the inclination i is determined by the highest latitude of the region that needs to be monitored.Secondly, radiusr2 or altitude h2 of the target orbit is determined by the number of revolutions in the period of repetition, and it is a given reference for the design of the altitude h1 of the initial orbit.Thirdly, the repetition period of the target orbit is one nodal day, and there is only one satellite in each orbit plane. The number of satellitesK, which is also the number of orbit planes, is determined by the specific observation frequency requirements and the economic conditions of the designer.Therefore, once the number of satellitesK, the inclination i, and the altitude h2 of the target orbit are determined by the designer, design variables would be reduced to the altitude h1 of the initial orbit and the AOL uk each initial orbit, where k=1,2,⋯,K. For the first satellite, it is reasonable to set u1=0, since only phase differences between each satellite matter.Based on the analysis above, our design can be concluded as designing the optimal altitude of the initial orbit and phase differences between each orbit plane, and the purpose is to minimize the fuel consumption, maximum reconfiguration time, mean transition time of all satellites for all targets, the standard deviation of the reconfiguration time for all targets, and maximum waiting days for all target sites.Before formulating the optimization problem, symbols of design variables and design metrics need to be defined. Assume that there areS (S≥1) target sites on the Earth’s surface and K (K≥2) satellites in the constellation. Design variables are the altitude h1 of the initial orbit and the phase differences Δuk between satellite k and the first satellite, where 0≤Δuk≤2π for k=2,⋯,K. For an arbitrary target site s (s=1,2,⋯,S) and an arbitrary satellite k (k=1,2,⋯,K), the impulse Δv, transition time ΔTs,k, and waiting days Ds,k can be obtained based on the mathematical model described before. It should be mentioned that for a constellation, the reconfiguration time is defined as the maximum transition time of all satellites. Therefore, reconfiguration time ΔTs of the constellation for target site s can be defined as
(24)ΔTs=maxΔTs,1,ΔTs,2,⋯,ΔTs,K.Also, the maximum waiting days of the constellation for the target sites can be expressed as
(25)Ds=maxDs,1,Ds,2,⋯,Ds,K.In conclusion, the optimization problem for ReCon design can be formulated as(26)find:minfix,i=1,2,3,4,5subjectto:xLB≤x≤xUB,where design metrics fi are defined as
(27)f1=Δv,f2=maxΔT1,ΔT2,⋯,ΔTS,f3=1S1K∑s=1S∑k=1KΔTs,k,f4=∑s=1SΔTs−ΔT¯2S−1,f5=maxD1,D2,⋯,DS,where ΔT¯ denotes the mean value of reconfiguration times for all target sites, which can be expressed by
(28)ΔT¯=1S∑s=1SΔTs.x denotes the vector of design variables, and it is expressed as
(29)x=h1,Δu2,Δu3,⋯,ΔuK,xLB and xUB are the lower and upper bound, and they are determined by the designer. For the altitude h1, hmin and hmax are defined as the lower and upper bound, while for the phase differences Δuk (k=2,⋯,K), 0 and 2π are the lower and upper bound. Thus, xLB and xUB can be expressed as
(30)xLB=hmin,0,0,⋯,0,xUB=hmax,2π,2π,⋯,2π,where numbers of “0” and “2π” are both K−1.


## 3. Physical Programming Method

A physical programming method is an efficient way to deal with multiobjective problems, of which the main idea is to convert the problem into a single-objective problem by using preference functions that capture the designer’s preferences [24, 25]. Tradeoffs between all design metrics are considered during problem formulation, and Pareto design which is in the preferred region of the Pareto surface is generated. Consequently, a compromise solution is obtained.


### 3.1. Preference Function Description

Preference functions are defined with the intention to reflect the preferences of the designer, and each preference function corresponds to a design metric. Assume that there arensc design metrics, and each of them is marked as gi; also, value of the corresponding preference function is marked as Higi, where i=1,2,⋯,nsc. The less the value of preference function Higi, the more satisfied the design metric gi. It should be noticed that gi is directly determined by the vector of design variables x, and hence, gi can also be expressed as function gix.Generally, preferences in physical programming are divided into the following four main classes: Class 1, smaller is better; Class 2, larger is better; Class 3, value is better; and Class 4, range is better. Each class is divided into soft and hard types as well, such as Class 1S and Class 1H, where the suffix letters S and H represent soft type and hard type, respectively. For the soft one, value of the preference function is gradually changed in the feasible region, yet for the hard one, value of the preference function in the feasible region stays the same, which is set to be 0. As a whole, preference functions can be concluded to be eight basic types, as listed in Figure5.Figure 5
Preference function classification for physical programming.Furthermore, for soft types, design metrics can be divided into several regions. Here, we take Class 1S as an example, and the rest of soft types can be defined in the similar pattern. Regions for Class 1S can be listed as follows:Region 1: highly desirable range (gi≤gi1).Region 2: desirable range (gi1≤gi≤gi2).Region 3: tolerable range (gi2≤gi≤gi3).Region 4: undesirable range (gi3≤gi≤gi4).Region 5: highly undesirable range (gi4≤gi≤gi5).Region 6: unacceptable range (gi>gi5),wheregi1 to gi5 are the region end points, which are given by the designer according to the specific preferences. These six regions are given in Figure 6.Figure 6
Preference function regions for Class 1S.


### 3.2. Preference Function Formulation

Preference functions are formulated based on the six interval divisions of the design metrics. By designing the value and the slope of preference function at each region end point, and fitting the value of preference function in each region, the specific preference function can be obtained.Class 1S is illustrated as an example. Letk be the region number, and k=1,2,⋯,5. As for region 6, the preference function is a linear extension of the end point of region 5. Value and slope of the preference function at the k region end point gik are Hik and sik, respectively.To calculate the valueHik, the modification ΔHik of the preference function in the k region is introduced. For the first region, the value and the modification of the preference function are usually set to be 0.1, and the rests can be calculated by the following equations:
(31)Hi1=ΔHi1=0.1,ΔHik=βnscΔHik−1,k=2,3,4,5,Hik=Hik−1+ΔHik,k=2,3,4,5,where β is a constant determined by the designer and β>5.To calculate the slopesik, the average slope s¯ik of the preference function in k region is introduced, which is defined as
(32)s¯ik=ΔHikλik.Then, the slope of the first region can be expressed as(33)si1=αs¯i2,where α is a constant determined by the designer and α>0.Slopes for the adjacent region can be expressed as(34)sik=134s¯ik−sik−1+8α3s¯ik−sik−1,k=2,3,4,5.Therefore, the value of preference function in each region can be obtained, which is illustrated as follows.Whenk=1,
(35)Hi=Hi1expsi1Hi1gi−gi1,gi≤gi1.Whenk=2,3,4,5,
(36)Hi=T0ξikgik−1+T1ξikgik+T2λik,ξiksik−1+T3λik,ξiksik,where
(37)ξik=gi−gik−1gik−gik−1,λik=gik−gik−1,and functions T0 to T3 are expressed as follows:
(38)T0ξ≡12ξ4−12ξ−14−2ξ+32,T1ξ≡−12ξ4+12ξ−14+2ξ−12,T2ξ,λ≡λ18ξ4−38ξ−14−12ξ+38,T3ξ,λ≡λ18ξ4−18ξ−14−12ξ+18.As a result, the preference function of Class 1S is formulated. Other types of the preference functions can be formulated in a similar way.Finally, the aggregate objective function is formed, and the multiobjective function is consequently converted to a single-objective function, which can be expressed as(39)minfx=log101nsc∑i=1nscHigix,where x is the vector of design variables.In this paper,nsc=5, and gi is the design metric fi proposed in equation (27). For the classification of design metrics, velocity increment, maximum reconfiguration time, mean transition time, and the standard deviation of the reconfiguration time for all targets belong to Class 1S, while maximum waiting days belong to Class 1H.


## 3.1. Preference Function Description

Preference functions are defined with the intention to reflect the preferences of the designer, and each preference function corresponds to a design metric. Assume that there arensc design metrics, and each of them is marked as gi; also, value of the corresponding preference function is marked as Higi, where i=1,2,⋯,nsc. The less the value of preference function Higi, the more satisfied the design metric gi. It should be noticed that gi is directly determined by the vector of design variables x, and hence, gi can also be expressed as function gix.Generally, preferences in physical programming are divided into the following four main classes: Class 1, smaller is better; Class 2, larger is better; Class 3, value is better; and Class 4, range is better. Each class is divided into soft and hard types as well, such as Class 1S and Class 1H, where the suffix letters S and H represent soft type and hard type, respectively. For the soft one, value of the preference function is gradually changed in the feasible region, yet for the hard one, value of the preference function in the feasible region stays the same, which is set to be 0. As a whole, preference functions can be concluded to be eight basic types, as listed in Figure5.Figure 5
Preference function classification for physical programming.Furthermore, for soft types, design metrics can be divided into several regions. Here, we take Class 1S as an example, and the rest of soft types can be defined in the similar pattern. Regions for Class 1S can be listed as follows:Region 1: highly desirable range (gi≤gi1).Region 2: desirable range (gi1≤gi≤gi2).Region 3: tolerable range (gi2≤gi≤gi3).Region 4: undesirable range (gi3≤gi≤gi4).Region 5: highly undesirable range (gi4≤gi≤gi5).Region 6: unacceptable range (gi>gi5),wheregi1 to gi5 are the region end points, which are given by the designer according to the specific preferences. These six regions are given in Figure 6.Figure 6
Preference function regions for Class 1S.


## 3.2. Preference Function Formulation

Preference functions are formulated based on the six interval divisions of the design metrics. By designing the value and the slope of preference function at each region end point, and fitting the value of preference function in each region, the specific preference function can be obtained.Class 1S is illustrated as an example. Letk be the region number, and k=1,2,⋯,5. As for region 6, the preference function is a linear extension of the end point of region 5. Value and slope of the preference function at the k region end point gik are Hik and sik, respectively.To calculate the valueHik, the modification ΔHik of the preference function in the k region is introduced. For the first region, the value and the modification of the preference function are usually set to be 0.1, and the rests can be calculated by the following equations:
(31)Hi1=ΔHi1=0.1,ΔHik=βnscΔHik−1,k=2,3,4,5,Hik=Hik−1+ΔHik,k=2,3,4,5,where β is a constant determined by the designer and β>5.To calculate the slopesik, the average slope s¯ik of the preference function in k region is introduced, which is defined as
(32)s¯ik=ΔHikλik.Then, the slope of the first region can be expressed as(33)si1=αs¯i2,where α is a constant determined by the designer and α>0.Slopes for the adjacent region can be expressed as(34)sik=134s¯ik−sik−1+8α3s¯ik−sik−1,k=2,3,4,5.Therefore, the value of preference function in each region can be obtained, which is illustrated as follows.Whenk=1,
(35)Hi=Hi1expsi1Hi1gi−gi1,gi≤gi1.Whenk=2,3,4,5,
(36)Hi=T0ξikgik−1+T1ξikgik+T2λik,ξiksik−1+T3λik,ξiksik,where
(37)ξik=gi−gik−1gik−gik−1,λik=gik−gik−1,and functions T0 to T3 are expressed as follows:
(38)T0ξ≡12ξ4−12ξ−14−2ξ+32,T1ξ≡−12ξ4+12ξ−14+2ξ−12,T2ξ,λ≡λ18ξ4−38ξ−14−12ξ+38,T3ξ,λ≡λ18ξ4−18ξ−14−12ξ+18.As a result, the preference function of Class 1S is formulated. Other types of the preference functions can be formulated in a similar way.Finally, the aggregate objective function is formed, and the multiobjective function is consequently converted to a single-objective function, which can be expressed as(39)minfx=log101nsc∑i=1nscHigix,where x is the vector of design variables.In this paper,nsc=5, and gi is the design metric fi proposed in equation (27). For the classification of design metrics, velocity increment, maximum reconfiguration time, mean transition time, and the standard deviation of the reconfiguration time for all targets belong to Class 1S, while maximum waiting days belong to Class 1H.


## 4. Numerical Results

In this section, simulation based on the proposed model is carried out. Initial conditions and preferences of design metrics in physical programming are set up by referring to the literature [3], and also, results are compared with that obtained in the literature. Furthermore, the GOM and ROM of the ReCon are analyzed and compared through a specific case, and the superior observation properties of the ReCon are demonstrated.


### 4.1. Simulation Setup

Initial conditions for the simulation are presented as follows. The inclination is selected as 90° for the global monitoring, and the number of satellitesK is set to 5, which are the same as the BLASTOFF constellation presented in literature [3]. Target sites around the Earth’s surface in step of 1 degree are considered, ignoring the north and south poles, where satellites will pass by directly without any maneuvering. That is to say, 64440 (179×360) target sites in total are considered as possible disaster sites. For the final repeating ground track orbit, parameters introduced in Section 2 are listed in Table 1.Table 1
Parameters of the final repeating ground track orbit.
ParameterValueUnitTr1436.07MinuteNp15—Tp295.74MinuteNd1—Td1436.07Minuteh2547.883kmAs for design variables and design metrics, specific ranges are determined according to literature [3]. In the referenced literature, minimizing the change in the orbital period is introduced as a means to minimizing the fuel consumption. The minimum period change is 12.125 minutes for 5 satellites, and the corresponding velocity increment for each satellite is about 62.961 m/s. Considering that the constellation in our design needs to go back to GOM after disaster monitoring, two Hohmann transfers are required for each satellite. Thus, if half of 62.961 m/s is used for the transfer, then the initial orbital altitude is about 605.8 km. Therefore, it is reasonable to set boundaries of the initial altitude as hmin=580 km and hmax=620 km, and the corresponding velocity increments are about 17.5 m/s and 39.2 m/s. Also, the RAAN and the AOL of the first satellite, which is also the anchor satellite, are both set to 0, for only the relative positions of satellites matter.Based on the boundaries of the initial altitude, the differences of angular velocity of the initial and target orbits are determined, and the transfer ability can be expressed as equation (22). Therefore, it is reasonable to set the boundaries of the maximum reconfiguration time as 100 and 200 hours, and the maximum waiting days are set to 6 days. Also, preferred value of the average reconfiguration time is set to 48 hours according to the BLASTOFF constellation in the literature [3]. As for the standard deviation of the reconfiguration time, 0 is a preferred value to ensure that the constellation has the similar ability to monitor any target site through orbital maneuver. As a whole, feasible regions of design metrics can be designed, and details of the region end points for Class 1S are provided in Table 2.Table 2
Details of the region end points for Class 1S.
ObjectiveUnitHighly desirableDesirableTolerableUndesirableHighly undesirablegi1gi2gi3gi4gi5g1m/s1720253240g2Hour100125150175200g3Hour2436486072g4Hour012243648A genetic algorithm is employed for the optimization. Parameters used for the genetic algorithm are listed in Table3.Table 3
Parameters of the genetic algorithm.
ParameterValuePopulation size50Generations100Crossover fraction0.8Migration fraction0.2Function tolerance1e−5


### 4.2. Optimization Results

Based on the introduced physical programming method and genetic algorithm, optimization is carried out, and results are given as follows.Change of the best and mean fitness value during the optimization are presented in Figure7, and these two values tend to be convergent within the given tolerance after 83 generations.Figure 7
Best and mean fitness in each generation.Values of design variables are listed in Table4.Table 4
Values of design variables.
Variableh1Δu2Δu3Δu4Δu5Value597.766216.15372.357287.738144.232UnitkmdegdegdegdegFrom the results presented in Table4, we can find that each phase difference is nearly an integral of 72. Actually, due to the uniformly distributed target sites and the uniformly distributed orbits in this problem, it is easy to form an intuitive conclusion that satellites should be distributed uniformly, and the optimization results exactly verified this conclusion.Furthermore, if we consider the walker-δ constellation, of which configuration parameters can be expressed as N/P/F, i, h [22], and set N=P=5, F=3, i=90°, and h=h1, then we can find that results presented in Table 4 are quite close to the proposed walker-δ constellation. Actually, it is easy to verify that the proposed walker-δ constellation, which can also be marked as configuration 5/5/3 [26], is the optimal result for our design by our program. Besides, uniform performance of the walker-δ constellation for the whole global can also consequently be verified in our work.Values of design metrics are given in Table5.Table 5
Values of design metrics.
Variableg1g2g3g4g5Value27.172150.01050.66122.5626Unitm/sHourHourHourDayFrom the results presented in Table5, we can find that the final velocity increment for each satellite is 27.172 m/s, and considering that the constellation is required to turn back to the GOM, the velocity increment for each satellite during the whole process is 54.344 m/s, which is about 86.3% of that of the BLASTOFF. It should be mentioned that the velocity increment is optimized by the designer’s preference, and hence, the value can be smaller at the cost of the increment of the reconfiguration time.Also, the maximum reconfiguration time is about 150 hours, while the mean transition time is only 50 hours. That is to say, for any target site on the ground, the whole constellation can finish the reconfiguration in 6 days, and the average transition time for each satellite is about 2 days. Therefore, satellites in the constellation can achieve the orbital transfer gradually within 6 days for any target on the global, and the target site can be overhead monitored during the whole reconfiguration process. Besides, the standard deviation of reconfiguration times for all target sites is minimized to about 22.6 hours, which is much less than the maximum reconfiguration time, and similar reconfiguration ability for all target sites can be achieved.In addition, a repeating ground track orbit is chosen to be the target orbit of all satellites in our design, which means the disaster site can always be directly overhead monitored by the whole constellation after the reconfiguration. Once the observation mission ends, the constellation can convert back to the GOM in an inverse process.


### 4.3. Case Study

In this subsection, a practical case is studied to demonstrate the superiority and effectiveness of the designed ReCon. A monitoring mission for a forest fire in Liangshan, Sichuan, China, is introduced. The GOM and ROM of the Recon, as well as the BLASTOFF proposed in literature [3], are all considered for the monitoring mission, and specific observation performances of these constellations are compared.The mission starting time is 2020.3.30.10:00:00 UTC, and the geographic latitude and longitude of the target site are 28.531°N and 101.241°E. For the ReCon, RAAN and AOL of the anchor satellite at the starting time are both set to 0, and orbital altitude and phase differences of the rest satellites are provided in Table4. For the BLASTOFF, specific parameters at the starting time can be found in Table 2 of literature [3]. As for the camera on the satellite, specific parameters are set according to the GF-1 [27]. High-resolution modes of the camera are considered for the overhead fire monitoring, and the width of the satellite swath on the ground is 60 km.Observation sequences of different constellations for the target site are presented in Figure8. The X-axis represents the time from the mission beginning, and the unit is hour. Each triangle along the X-axis denotes an observation, and the x-value of the vertex of a triangle represents the monitoring moment, since the observation duration is quite short. The figure of GOM shows the observation sequences of the initial constellation, while the figure of ROM presents the observation sequences of the ReCon during and after the reconfiguration.Figure 8
Observation sequences of different constellations.The total number of observations varies with the time from the mission beginning for each constellation, and the specific change trend can be concluded based on the observation sequences, as illustrated in Figure9.Figure 9
Total observation numbers of different constellations.From Figures8 and 9, we can find that the total number of observations of the ROM and that of the BLASTOFF are much more than that of the GOM, and hence, observation for the disaster area after reconfiguration is more adequate and frequent. Moreover, the ROM can achieve 4 extra observations within 2 days compared to the BLASTOFF, and thus, observation of the ROM is much more timely, although the total number of observations of the ROM is 2 less than that of the BLASTOFF. Besides, observation of the ROM becomes regular (5 times per day as the BLASTOFF) after about 100 hours, and this regular observation would remain until the mission ends without extra orbital transfer of any satellite in the constellation.In the foregoing case, each satellite in the ROM monitors the target site during the ascending pass, and the orbit of each satellite is a repeating orbit with the same repeating period. Therefore, all satellites in the ROM are flying along the same repeating ground track, and this constellation is also the so-called Flower constellation [17], which is marked as ROM1 in this work. Actually, each satellite is able to monitor the target site during the descending pass as illustrated in our mathematical mode. If the initial condition of the constellation changes, the final ROM might be different from the ROM1. Here, we set the AOL of the anchor satellite at the starting time as 180°, and the rest conditions stay the same. Then, we can get a new observation mode, which is marked as ROM2. In the ROM1, all satellites observe the target site during ascending pass, and the AOL of each satellite stays the same with that of the anchor one, while in the ROM2, there is one satellite observing the target site during the descending pass, and its AOL is 56.740° behind the others’. Specific configurations of the GOM, ROM1, and ROM2 are presented in Figure 10.Figure 10
Configurations of different constellations.
(a)
GOM(b)
ROM1(c)
ROM2Based on the presented configurations, ground tracks and average observation numbers per day of the GOM, the ROM1, and the ROM2 can be analyzed, as shown in Figure11, where the red point represents the target site. It should be mentioned that average observation numbers are calculated by the method presented by Razoumny [28, 29]. From the figure, we can find that in the GOM, satellites can evenly scan the entire region within a latitude band without bias, while in the ROM, satellites can scan the region near a fixed ground track frequently. Specifically, in the GOM, satellites can observe the target site about 0.25 times a day, which means that the target site could only get one observation in about 4 days, while all areas in the latitude band can be monitored. In the ROM, satellites can observe the target site 5 times a day, and particularly in the ROM1, only the region around the target site can get 5 observations in one day, and most of the areas cannot be monitored at all. As a whole, the ReCon can scan the Earth uniformly in GOM as a normal operation, and when disaster occurs, it can turn to ROM, and satellites can scan the disaster area adequately and frequently.Figure 11
Ground tracks and average observation numbers per day of different constellations.
(a)
GOM(b)
ROM1(c)
ROM2To further analyze different properties of the GOM and the ROM, the low-resolution mode of the GF-1 camera is considered, of which the swath width is 800 km. Similarly, average observation numbers per day of the GOM, the ROM1, and the ROM2 are calculated, and results are given in Figure12. With a wider swath of each satellite, the GOM can achieve about 3 times of observations per day for most of areas globally, while the ROM can monitor the target site 10 times per day. For the GOM, observation numbers per day for a specific latitude is still almost the same, while for both ROM1 and ROM2, areas around the target are the main observation object. Therefore, different modes of the ReCon show excellent observation performances for both global observation and disaster monitoring.Figure 12
Average observation numbers per day of different constellations.
(a)
GOM(b)
ROM1(c)
ROM2In addition, to verify the accuracy of the proposed mathematical model, distances between the target site and the satellite ground tracks in the ROM1 and the ROM2 are presented in Figure13, which are also the enlarged views of the ground tracks shown in Figures 11(b) and 11(c). It should be noted that the orbital elements and ground tracks of each satellite in the ROM1 and the ROM2 are calculated through numerical simulation, while in the numerical simulation, the transition times and the moments to implement the orbital transfer are given by the proposed mathematical model. From Figure 13, we can find that the minimum arc distance between the target site and the satellite ground tracks is less than 0.01 degrees, which is about a few hundred meters and does not affect the overhead observation. Therefore, the analytical calculations of the transition times and the moments to implement the orbital transfer are accurate enough for the initial observation mission planning.Figure 13
Enlarged views of the ground tracks of the ROM1 and the ROM2.
(a)
ROM1(b)
ROM2


## 4.1. Simulation Setup

Initial conditions for the simulation are presented as follows. The inclination is selected as 90° for the global monitoring, and the number of satellitesK is set to 5, which are the same as the BLASTOFF constellation presented in literature [3]. Target sites around the Earth’s surface in step of 1 degree are considered, ignoring the north and south poles, where satellites will pass by directly without any maneuvering. That is to say, 64440 (179×360) target sites in total are considered as possible disaster sites. For the final repeating ground track orbit, parameters introduced in Section 2 are listed in Table 1.Table 1
Parameters of the final repeating ground track orbit.
ParameterValueUnitTr1436.07MinuteNp15—Tp295.74MinuteNd1—Td1436.07Minuteh2547.883kmAs for design variables and design metrics, specific ranges are determined according to literature [3]. In the referenced literature, minimizing the change in the orbital period is introduced as a means to minimizing the fuel consumption. The minimum period change is 12.125 minutes for 5 satellites, and the corresponding velocity increment for each satellite is about 62.961 m/s. Considering that the constellation in our design needs to go back to GOM after disaster monitoring, two Hohmann transfers are required for each satellite. Thus, if half of 62.961 m/s is used for the transfer, then the initial orbital altitude is about 605.8 km. Therefore, it is reasonable to set boundaries of the initial altitude as hmin=580 km and hmax=620 km, and the corresponding velocity increments are about 17.5 m/s and 39.2 m/s. Also, the RAAN and the AOL of the first satellite, which is also the anchor satellite, are both set to 0, for only the relative positions of satellites matter.Based on the boundaries of the initial altitude, the differences of angular velocity of the initial and target orbits are determined, and the transfer ability can be expressed as equation (22). Therefore, it is reasonable to set the boundaries of the maximum reconfiguration time as 100 and 200 hours, and the maximum waiting days are set to 6 days. Also, preferred value of the average reconfiguration time is set to 48 hours according to the BLASTOFF constellation in the literature [3]. As for the standard deviation of the reconfiguration time, 0 is a preferred value to ensure that the constellation has the similar ability to monitor any target site through orbital maneuver. As a whole, feasible regions of design metrics can be designed, and details of the region end points for Class 1S are provided in Table 2.Table 2
Details of the region end points for Class 1S.
ObjectiveUnitHighly desirableDesirableTolerableUndesirableHighly undesirablegi1gi2gi3gi4gi5g1m/s1720253240g2Hour100125150175200g3Hour2436486072g4Hour012243648A genetic algorithm is employed for the optimization. Parameters used for the genetic algorithm are listed in Table3.Table 3
Parameters of the genetic algorithm.
ParameterValuePopulation size50Generations100Crossover fraction0.8Migration fraction0.2Function tolerance1e−5


## 4.2. Optimization Results

Based on the introduced physical programming method and genetic algorithm, optimization is carried out, and results are given as follows.Change of the best and mean fitness value during the optimization are presented in Figure7, and these two values tend to be convergent within the given tolerance after 83 generations.Figure 7
Best and mean fitness in each generation.Values of design variables are listed in Table4.Table 4
Values of design variables.
Variableh1Δu2Δu3Δu4Δu5Value597.766216.15372.357287.738144.232UnitkmdegdegdegdegFrom the results presented in Table4, we can find that each phase difference is nearly an integral of 72. Actually, due to the uniformly distributed target sites and the uniformly distributed orbits in this problem, it is easy to form an intuitive conclusion that satellites should be distributed uniformly, and the optimization results exactly verified this conclusion.Furthermore, if we consider the walker-δ constellation, of which configuration parameters can be expressed as N/P/F, i, h [22], and set N=P=5, F=3, i=90°, and h=h1, then we can find that results presented in Table 4 are quite close to the proposed walker-δ constellation. Actually, it is easy to verify that the proposed walker-δ constellation, which can also be marked as configuration 5/5/3 [26], is the optimal result for our design by our program. Besides, uniform performance of the walker-δ constellation for the whole global can also consequently be verified in our work.Values of design metrics are given in Table5.Table 5
Values of design metrics.
Variableg1g2g3g4g5Value27.172150.01050.66122.5626Unitm/sHourHourHourDayFrom the results presented in Table5, we can find that the final velocity increment for each satellite is 27.172 m/s, and considering that the constellation is required to turn back to the GOM, the velocity increment for each satellite during the whole process is 54.344 m/s, which is about 86.3% of that of the BLASTOFF. It should be mentioned that the velocity increment is optimized by the designer’s preference, and hence, the value can be smaller at the cost of the increment of the reconfiguration time.Also, the maximum reconfiguration time is about 150 hours, while the mean transition time is only 50 hours. That is to say, for any target site on the ground, the whole constellation can finish the reconfiguration in 6 days, and the average transition time for each satellite is about 2 days. Therefore, satellites in the constellation can achieve the orbital transfer gradually within 6 days for any target on the global, and the target site can be overhead monitored during the whole reconfiguration process. Besides, the standard deviation of reconfiguration times for all target sites is minimized to about 22.6 hours, which is much less than the maximum reconfiguration time, and similar reconfiguration ability for all target sites can be achieved.In addition, a repeating ground track orbit is chosen to be the target orbit of all satellites in our design, which means the disaster site can always be directly overhead monitored by the whole constellation after the reconfiguration. Once the observation mission ends, the constellation can convert back to the GOM in an inverse process.


## 4.3. Case Study

In this subsection, a practical case is studied to demonstrate the superiority and effectiveness of the designed ReCon. A monitoring mission for a forest fire in Liangshan, Sichuan, China, is introduced. The GOM and ROM of the Recon, as well as the BLASTOFF proposed in literature [3], are all considered for the monitoring mission, and specific observation performances of these constellations are compared.The mission starting time is 2020.3.30.10:00:00 UTC, and the geographic latitude and longitude of the target site are 28.531°N and 101.241°E. For the ReCon, RAAN and AOL of the anchor satellite at the starting time are both set to 0, and orbital altitude and phase differences of the rest satellites are provided in Table4. For the BLASTOFF, specific parameters at the starting time can be found in Table 2 of literature [3]. As for the camera on the satellite, specific parameters are set according to the GF-1 [27]. High-resolution modes of the camera are considered for the overhead fire monitoring, and the width of the satellite swath on the ground is 60 km.Observation sequences of different constellations for the target site are presented in Figure8. The X-axis represents the time from the mission beginning, and the unit is hour. Each triangle along the X-axis denotes an observation, and the x-value of the vertex of a triangle represents the monitoring moment, since the observation duration is quite short. The figure of GOM shows the observation sequences of the initial constellation, while the figure of ROM presents the observation sequences of the ReCon during and after the reconfiguration.Figure 8
Observation sequences of different constellations.The total number of observations varies with the time from the mission beginning for each constellation, and the specific change trend can be concluded based on the observation sequences, as illustrated in Figure9.Figure 9
Total observation numbers of different constellations.From Figures8 and 9, we can find that the total number of observations of the ROM and that of the BLASTOFF are much more than that of the GOM, and hence, observation for the disaster area after reconfiguration is more adequate and frequent. Moreover, the ROM can achieve 4 extra observations within 2 days compared to the BLASTOFF, and thus, observation of the ROM is much more timely, although the total number of observations of the ROM is 2 less than that of the BLASTOFF. Besides, observation of the ROM becomes regular (5 times per day as the BLASTOFF) after about 100 hours, and this regular observation would remain until the mission ends without extra orbital transfer of any satellite in the constellation.In the foregoing case, each satellite in the ROM monitors the target site during the ascending pass, and the orbit of each satellite is a repeating orbit with the same repeating period. Therefore, all satellites in the ROM are flying along the same repeating ground track, and this constellation is also the so-called Flower constellation [17], which is marked as ROM1 in this work. Actually, each satellite is able to monitor the target site during the descending pass as illustrated in our mathematical mode. If the initial condition of the constellation changes, the final ROM might be different from the ROM1. Here, we set the AOL of the anchor satellite at the starting time as 180°, and the rest conditions stay the same. Then, we can get a new observation mode, which is marked as ROM2. In the ROM1, all satellites observe the target site during ascending pass, and the AOL of each satellite stays the same with that of the anchor one, while in the ROM2, there is one satellite observing the target site during the descending pass, and its AOL is 56.740° behind the others’. Specific configurations of the GOM, ROM1, and ROM2 are presented in Figure 10.Figure 10
Configurations of different constellations.
(a)
GOM(b)
ROM1(c)
ROM2Based on the presented configurations, ground tracks and average observation numbers per day of the GOM, the ROM1, and the ROM2 can be analyzed, as shown in Figure11, where the red point represents the target site. It should be mentioned that average observation numbers are calculated by the method presented by Razoumny [28, 29]. From the figure, we can find that in the GOM, satellites can evenly scan the entire region within a latitude band without bias, while in the ROM, satellites can scan the region near a fixed ground track frequently. Specifically, in the GOM, satellites can observe the target site about 0.25 times a day, which means that the target site could only get one observation in about 4 days, while all areas in the latitude band can be monitored. In the ROM, satellites can observe the target site 5 times a day, and particularly in the ROM1, only the region around the target site can get 5 observations in one day, and most of the areas cannot be monitored at all. As a whole, the ReCon can scan the Earth uniformly in GOM as a normal operation, and when disaster occurs, it can turn to ROM, and satellites can scan the disaster area adequately and frequently.Figure 11
Ground tracks and average observation numbers per day of different constellations.
(a)
GOM(b)
ROM1(c)
ROM2To further analyze different properties of the GOM and the ROM, the low-resolution mode of the GF-1 camera is considered, of which the swath width is 800 km. Similarly, average observation numbers per day of the GOM, the ROM1, and the ROM2 are calculated, and results are given in Figure12. With a wider swath of each satellite, the GOM can achieve about 3 times of observations per day for most of areas globally, while the ROM can monitor the target site 10 times per day. For the GOM, observation numbers per day for a specific latitude is still almost the same, while for both ROM1 and ROM2, areas around the target are the main observation object. Therefore, different modes of the ReCon show excellent observation performances for both global observation and disaster monitoring.Figure 12
Average observation numbers per day of different constellations.
(a)
GOM(b)
ROM1(c)
ROM2In addition, to verify the accuracy of the proposed mathematical model, distances between the target site and the satellite ground tracks in the ROM1 and the ROM2 are presented in Figure13, which are also the enlarged views of the ground tracks shown in Figures 11(b) and 11(c). It should be noted that the orbital elements and ground tracks of each satellite in the ROM1 and the ROM2 are calculated through numerical simulation, while in the numerical simulation, the transition times and the moments to implement the orbital transfer are given by the proposed mathematical model. From Figure 13, we can find that the minimum arc distance between the target site and the satellite ground tracks is less than 0.01 degrees, which is about a few hundred meters and does not affect the overhead observation. Therefore, the analytical calculations of the transition times and the moments to implement the orbital transfer are accurate enough for the initial observation mission planning.Figure 13
Enlarged views of the ground tracks of the ROM1 and the ROM2.
(a)
ROM1(b)
ROM2


## 5. Conclusion

In this paper, a reconfigurable constellation, which can achieve timely and adequate observation for an unpredictable disaster site through orbital maneuvers, is successfully designed. By analyzing the space geometry relations between an arbitrary satellite orbit and disaster site, a mathematical model for orbital transfer and overhead monitoring is established. Aiming at minimizing both total velocity increment and reconfiguration time, optimization of the initial constellation and the reconfiguration process are combined and carried out, and a physical programming method along with a genetic algorithm is employed for the multiobjective optimization problem. Preferences in physical programming are set by referring to literature [3], and results are compared to demonstrate the superiority of our design. Also, the GOM and ROM of the ReCon as well as the BLASTOFF are analyzed and compared, and excellent observation performances of the ReCon are demonstrated.As for the results of the optimization, the optimal reconfigurable constellation is extremely close to the walker-δ constellation. The reconfiguration velocity increment for each satellite is 54.344 m/s, which is about 86.3% of that in literature [3]. The maximum reconfiguration time is about 150 hours, while the mean transition time is only 50 hours, which means satellites can finish transition gradually, and the target site can be overhead monitored during the whole reconfiguration process. Also, the standard deviation of reconfiguration times for all target sites is minimized to less than one day to achieve similar reconfiguration ability for all target sites. A specific disaster case is analyzed, and the observation performance of the GOM, the ROM, and the BLASTOFF is compared. The GOM can achieve even observation for the global, while the ROM can monitor the target site frequently and adequately and much more timely than the BLASTOFF. Besides, different properties of the GOM and the ROM are further analyzed and compared, and results indicate that the ReCon has excellent observation performances for both global observation and disaster monitoring.Furthermore, it should be mentioned that the preferences in our model are designed according to the designer, and if less velocity increment or less reconfiguration time is required, it is easy to adjust our model and obtain the preferred optimal results. Future work would be to monitor 2 or more disaster sites, as well as to monitor the moving disaster site through constellation reconfiguration. This research can provide basic idea and method for the preliminary design of reconfigurable constellation for disaster monitoring.



---

*Source: 8813685-2020-09-01.xml*