Optimization of satellite resource scheduling under regional target coverage conditions

1 Introduction

In recent years, with the continuous development of satellite technology, people have increasingly hoped to be able to personalize the satellite orbit design for a specific task, rather than the traditional general design (Zhang and Feng 2019, Zhang and Wei 2016, Huang et al. 2021b, De Grossi et al. 2021, Quarta et al. 2020). At the same time, due to the influence of satellite height, field of view and other factors, the area covered by a single satellite when performing tasks is generally small, it is difficult to meet the application needs of navigation, meteorological forecasting, space exploration and scientific experiments. With the continuous development of the national economy, to meet the aforementioned application needs, multiple satellites will be formed in accordance with certain laws to form a satellite constellation, compared with the performance of a single satellite, and the satellite constellation has advantages of short revisit interval and good coverage performance (Luders 1961, Draim 1985, 1987, Draim et al. 2002, Mortari and Wilkins 2008, Walker 1982, Adams and Rider 1987, Yan et al. 2016, Guo et al. 2021).

In the constellation design of the area coverage mission, due to the diversity and flexibility of the area to be covered, there is currently no universal design method and conclusion to obtain the desired results (Wu and Wu 2007). Therefore, based on the current mature evolutionary algorithm, it is also an efficient and meaningful constellation design method to search for constellation parameters that can meet the constraint requirements by reasonably constructing an optimization model (Hu et al. 2021, Meziane-Tani et al. 2016, Savitri et al. 2017, Appel et al. 2014, Huang et al. 2021a). In this paper, the characteristics of common ground trajectory constellations that meet certain revisit interval constraints are analyzed (Wu et al. 1999) using the optimized model of the constructed regional coverage constellation design. The orbit parameters of the first datum satellite are first determined, and then other parameters in the common ground trajectory constellation are determined according to the revisit interval constraint and the characteristics of the area to be covered (Wu and Wu 2007). Moreover, in view of the emergence of some contingencies, it is sometimes necessary to shorten the re-visit interval of observations based on the existing constellations of common ground trajectories, so that updated observational data can be obtained in a shorter time period. This requires dispatching some existing satellites to complement the designed common ground trajectory constellations (Zeng 2014). In this paper, a simple and effective two-pulse maneuvering scheme is used to construct an optimized model of satellite mission scheduling to meet the needs of burst missions.

Finally, since the waters of the Northwest Pacific Ocean are the source of frequent tropical cyclone formations, posing a great threat to the lives and property of people in neighbouring countries, it is necessary to conduct surveillance observations of this area (Wu and Wu 2006, Lin et al. 2022, Wei et al. 2016). In such a mission requirement, through the optimization model and method of regional coverage constellation design and existing satellite scheduling proposed in this paper, the regular observation task of this region is successfully completed, which shows the effectiveness of the optimization model and method and has certain practical use value.

2 Optimization model construction for co-ground trajectory constellation design

2.1 Analysis of the characteristics of constellations of common ground trajectories

A common ground trajectory constellation is a special type of constellation in which all moons in a constellation follow an unchanging ground trajectory (Wu et al. 1999), its substellar trajectory repeats regularly and more compactly, making it suitable for coverage of a specified area. Because the substellar coverage zone of the satellite in the near-polar orbit presents a nearly “straight up and straight down” shape characteristics at the middle and low latitudes, it provides a certain foundation for simplifying the system management and control of the common ground trajectory constellation. Therefore, this paper uses near-polar orbit satellites to construct a common ground trajectory constellation, and its orbital inclination is generally between 8 5 ° and 9 0 ° .

In recent years, some domestic and foreign literature have given the analysis of the parameters and characteristics of the common ground trajectory constellation (Pennoni 1994, Wu and Wu 2007), but the corresponding parameter characteristic analysis has not been given for the application of near-polar orbit satellites with revisit interval constraints in the design of the co-ground trajectory constellation. For the sake of narrative, a mathematical description of the target area is given here. A bounded arbitrary shape of the target area, along the tangent of its boundary to make a smallest area of rectangular area, is taken, and then, the rectangular area is evenly dispersed into a target point, expressed by a combination ( λ , ϕ ) of latitude and longitude. The degree of discrete density, that is, the longitude interval and latitude interval between the two target points, is found, and this is further determined by the width of the cover zone at the substellar point of the satellite in the specific question. Therefore, the region can be represented by a combination of latitude and longitude of a series of discrete target points. As shown in Figure 1, the target region in is mathematically represented as ( λ i j , ϕ i j ) = ( λ 0 + i × Δ λ , ϕ 0 + j × Δ ϕ ) , i = { 0 , 1 , … , i max } , j = { 0 , 1 , … , j max } , and the units are expressed in degrees.

Figure 1

Schematic of a mathematical description of the target area.

To make full use of the good substratum trajectory repetition characteristics of the common ground trajectory constellation, and thus achieve continuous coverage of specific areas, regression orbits are usually used in the design of co-ground trajectory constellations. The previous studies (Duan 2006, Yang et al. 2016) presented the mathematical relationships that the regression orbit should satisfy:

(1) N T Ω = D T e .

The intersection day T e is defined as the time interval at which the zero longitude line on the earth rotates 1 week relative to the upcross of the satellite orbit, and the intersection period T Ω is defined as the time interval between the satellites passing through the orbital upbringing point twice in a row. The regression orbit is the orbit of the satellite after passing through the intersection day, the orbit around the earth (after the intersection period), and the orbit of the substellar point repeats. The regression orbit is the orbit of the satellite after passing through D intersection days, rotating N circles around the Earth (after N intersection periods), and the orbit of the substellar point is repeated.

According to the derivation of literature (Wu and Wu 2007), in a common ground trajectory constellation using a regression orbit, the upshoot equatorial warp difference between adjacent orbital planes and the phase difference between adjacent satellites satisfy a simple linear relationship:

(2) Δ M = Δ Ω × N D .

The difference between the right ascension of the ascending node is Δ Ω = Ω 2 − Ω 1 , and the difference in phase is Δ M = M 2 − M 1 , which indicates that the phase of satellite 2 lags behind satellite 1.

In Figure 1, we suppose that there are two near-polar orbiting satellites forming a common ground trajectory constellation. The right ascension of the ascending node Ω 1 and mean anomaly M 1 of the first satellite are known, and at the moment t 1 , its substellar trajectory begins to cover the longitude band λ 1 specified in the graph. To revisit the longitude band after the Δ t time interval, the right ascension of the ascending node Ω 2 and mean anomaly M 2 of the second satellite need to satisfy the following mathematical relationships.

(3) Ω 2 = Ω 1 + Δ t 3,600 ω e ′ − 36 0 ° × n .

In the formula, ω e ′ is the average rotational angular velocity of the earth, but the unit used here is not the international unit, but deg/h. The unit of the right ascension of the ascending node Ω is deg, and the unit of Δ t is s. n is an integer, and its function is to guarantee 0 ° ⩽ Ω < 36 0 ° by appropriately subtracting the integer multiple of 36 0 ° . Substitute Eq. (3) into Eq. (2) to obtain:

(4) M 2 = M 1 − Δ t 3,600 ω e ′ × N D + 36 0 ° × n .

To cover the longitude band λ 2 adjacent to the longitude band λ 1 , only the right ascension of the ascending node of satellite 1 and satellite 2 needs to be adjusted to obtain the orbital parameters of satellite 3 and satellite 4 covering the longitude band, satisfying the following mathematical relationship:

(5) Ω 3 = Ω 1 + ( λ 2 − λ 1 ) Ω 4 = Ω 2 + ( λ 2 − λ 1 ) .

2.2 Analysis optimized design of orbital parameters of reference satellites

In this paper, considering only the gravitational field of the earth center and the role of J2-term perturbation, according to the literature (Meng and Dai 2005, Zeng 2014), to minimize the error caused by the influence of perturbation on orbital parameters on the model, it is best to use a circular orbit, that is, eccentricity e = 0 . For the selection of orbital altitude, the following two aspects are mainly considered, one is the impact of the Van Allen radiation belt on the spaceborne electronic device, and the other is the impact of the earth’s atmospheric resistance on the life of the satellite. The Van Allen radiation inner belt is in a height range of approximately 1,500–5,000 km, and to avoid adverse effects on satellite spaceborne electronics, it is best to design the orbital altitude to be lower than the minimum altitude of the Van Allen radiation inner belt. At the same time, to minimize the impact of atmospheric resistance on the orbital life of satellites, the orbital height should not be too low, generally designed to be greater than 500 km.

Based on the previous analysis, the most important thing is to first determine the orbital parameters of the first reference satellite. As shown in Figure 2, we refer to the process of the datum near-polar orbit satellite covering the longitude band λ 1 for the first time in a day through the target area, which is called the upstream process, and the time elapsed is called the upward time t up . Similarly, the process of covering adjacent longitude band through the target area for the second time in a day by the reference satellite is called the downward process, and the corresponding time is called the downward time t down . The period of time that elapses from the end of the upside to the beginning of the downside is called the uplink and downlink interval t interval .

Figure 2

Schematic diagram of the optimization design process for a near-polar-orbit reference satellite.

In the J2000 inertial coordinate system, the orbital dynamics equation for the orbital motion of a satellite that takes only term perturbations around the Earth is stated as follows:

(6) x ¨ = − μ e r 3 x + A J 2 μ e 15 x z 2 r 7 − 3 x r 5 y ¨ = − μ e r 3 y + A J 2 μ e 15 y z 2 r 7 − 3 y r 5 z ¨ = − μ e r 3 z + A J 2 μ e 15 z 3 r 7 − 9 z r 5 .

In the formula, A J 2 = 1 2 J 2 R e 2 , r = x 2 + y 2 + z 2 , μ e is the Earth’s gravitational constant and R e is the average radius of the Earth. x , y , z are coordinate components of the satellite in the inertial coordinate system. Suppose the initial moment of the epoch is T 0 , the number of orbital roots of the satellite is σ = ( a , e , i , Ω , ω , M ) , within the Δ t time interval from the initial moment, the satellite’s substellar point first ascended to a location near the boundary target point ( λ 1 , ϕ 0 ) . In the geostationary system, the spherical distance between the point under the satellite star and the target point is d 1 . Then, orbital dynamics Eq. (6) is used, it predicts the coordinates of the substellar point at which the satellite reaches the location near the boundary target point ( λ 1 , ϕ 0 + j max Δ ϕ ) after the uplink time t up , at which time the spherical distance from the target point is d 2 . Then the coordinates of the substellar point where the satellite reaches the location near the target point ( λ 2 , ϕ 0 + j max Δ ϕ ) after the up-down interval time t interval is predicted, and the spherical distance from the substellar point to the target point is d 3 . Finally, it predicts the coordinates of the substellar point at which the satellite reaches the location near the boundary target point ( λ 2 , ϕ 0 ) after the downlink time t down , at which time the spherical distance from the target point is d 4 .

In summary, the optimal design of the reference satellite orbit parameters is equivalent to considering the following constrained optimization problems. Considering that the reference satellite adopts a near-polar circular orbit with regression characteristics, the eccentricity and the argument of perigee are both set to zero. Therefore, the design variables can be expressed as the following vectors:

(7) X 1 = [ a , i , Ω , M , t up , t interval , t down ] T .

There are two main constraints to consider, one is to meet the constraints of the maximum range of satellite exploration:

(8) d 1 , d 2 , d 3 , d 4 ⩽ D max ,

where D max represents the maximum distance from the satellite observation target. The second needs to satisfy the constraints of the regression orbit, that is, to satisfy the mathematical relation Eq. (1). The optimized objective function is expressed as follows:

(9) min f ( X 1 ) = d 1 + d 2 + d 3 + d 4 .

3 Optimization model construction for satellite scheduling

An important feature of the satellite dispatch problem is that there is a constraint of the visible time window in the regional coverage task, and only within the visible time window can the scheduling task be completed more economically. Therefore, an important principle followed by satellite scheduling is to try to choose satellites whose time window covering the target area at the sub-satellite point is roughly consistent with the time range required in the scheduling task.

As shown in Figure 3, it may be advisable to select the bar area as the area to be covered when performing the scheduling task, and the shape of the area is a long strip of longitude band λ 1 acrossing latitude ϕ 0 – ϕ 4 .

Figure 3

Schematic diagram of the existing satellite optimization scheduling process.

Based on the aforementioned time window principle, it is assumed that satellite is selected as the satellite to be dispatched. The orbital parameter of the satellite at the epoch T 0 is σ 1 , but its substellar point coverage band does not completely cover the entire strip area a r e a 1 . To simplify the problem, the satellite’s propulsion uses instantaneous velocity pulses. Suppose that the component of the applied velocity pulse in the J2000 geocentric equatorial inertial coordinate system is Δ v x , Δ v y , Δ v z , time before and after the speed pulse is applied are, respectively, t − and t + , and the position and velocity of the satellite in the inertial frame satisfies the following equation:

In the aforementioned formula, x t − , y t − , z t − , x ˙ t − , y ˙ t − , z ˙ t − represent the position and velocity component of the satellite in the inertial coordinate system before the velocity pulse is applied. What’s more, x t + , y t + , z t + , x ˙ t + , y ˙ t + , and z ˙ t + represent the position and velocity component of the satellite in the inertial coordinate system after the velocity pulse is applied.

In this paper, a simple and effective two-pulse motorized orbit change scheme is adopted. The basic process is as follows: Suppose that the two orbital maneuver times from the epoch time T are, respectively, t 1 and t 2 , using the orbital dynamics Eq. (6). First, the position r 1 − and speed v 1 − of the satellite after the elapsed time t 1 are predicted, using the pulsed orbit change Eq. (10), it can obtain the position r 1 + and velocity v 1 + after applying the instantaneous velocity pulse, where r 1 − = r 1 + . Similarly, it is possible to obtain the position r 2 − and velocity v 2 − before the second maneuver change time t 2 , and it is possible to obtain the position r 2 + and velocity v 2 + after t 2 , where r 2 − = r 2 + . After the second maneuver, the satellite’s sub-satellite point will completely cover the bar area a r e a 1 within the required time interval [ t 3 , t 4 ] . To constrain satellite s a t 1 to be able to fully cover the strip area a r e a 1 within the required “transit” time interval [ t 3 , t 4 ] after the maneuver, this article uses the gold separation method to conduct a one-dimensional search for time interval [ t 3 , t 4 ] , the minimum spherical distance d 5 from the target point ( λ 1 , ϕ 0 ) when passing through the bar area a r e a 1 and the minimum spherical distance d 6 from the target point ( λ 1 , ϕ 4 ) when starting to leave the bar area a r e a 1 are obtained. The boundary constraint described earlier can be satisfied as long as both d 5 and d 6 are less than the maximum distance D max of the satellite observation target. The iteration steps are as follows.

• Step 1: Determine the initial parameters. Use the time interval [ t 3 , t 4 ] as the initial search interval. Choose iterative precision δ > 0 and calculate the values for the first two heuristic moments:

  (11) λ = t 4 − 0.618 ( t 4 − t 3 ) μ = t 3 + 0.618 ( t 4 − t 3 ) .

• Step 2: Compare spherical distance sizes. If d ( λ ) < d ( μ ) then proceed to Step 3, otherwise proceed to Step 4.

• Step 3: If ∣ t 4 − t 3 ∣ ⩽ δ , then the iteration is stopped, the output is a spherical distance d t 3 + t 4 2 from the target point ( λ 1 , ϕ 0 ) . Otherwise, we obtain

  (12) t 4 = μ μ = λ d ( μ ) = d ( λ ) λ = t 4 − 0.618 ( t 4 − t 3 ) .

  Calculate d ( λ ) and proceed to Step 2.

• Step 4: If ∣ t 4 − t 3 ∣ ⩽ δ , then the iteration is stopped, and the output is a spherical distance d t 3 + t 4 2 from the target point ( λ 1 , ϕ 0 ) . Otherwise, we obtain

  (13) t 3 = λ λ = μ d ( λ ) = d ( μ ) μ = t 3 − 0.618 ( t 4 − t 3 ) .

  Calculate d ( μ ) and proceed to Step 2.

The spherical distance d t 3 + t 4 2 from the target point ( λ 1 , ϕ 0 ) is the minimum spherical distance d 5 required. Similarly, using the iterative process of the aforementioned golden section method, the minimum spherical distance d 6 from the boundary target point ( λ 1 , ϕ 4 ) can be obtained.

In summary, the satellite dispatch problem can be considered as the following constraint optimization problem. The design variables are as follows:

where Δ v 1 represents the magnitude of the pulse velocity applied at the first maneuver at time t 1 , α 1 represents the angle between the applied pulse velocity Δ v 1 and the x -axis of the inertial coordinate system, and β 1 indicates the angle between the projection of the pulse velocity Δ v 1 in the y o z coordinate plane and the y axis. Then, the three components of the velocity pulse Δ v 1 in the inertial coordinate system can be expressed as follows:

Similarly, the three components Δ v 2 x , Δ v 2 y , and Δ v 2 z of the velocity pulse Δ v 1 applied at the second maneuver moment t can also be expressed.

There are three main aspects of the constraints to be considered. One is that the above boundary constraints need to be met:

Second, according to the analysis in Section 2, the orbital height h needs to meet a certain range of constraints:

Third, the satellite orbit after the mission dispatch also needs to meet the regression characteristics, and the equation constraint of the regression orbit can meet the mathematical relationship Eq. (1). The optimized objective function is expressed as follows:

4 Examples and analysis

Some scholars have conducted statistical analysis of the characteristics of tropical cyclones in the northwest Pacific Ocean in the past 50 years from 1949 to 1996 (Liu and Liu 2012) and concluded that most of the stronger typhoons originated from the ocean surface east of 12 5 ° east longitude and more than 90% of tropical cyclones in the northwest Pacific Ocean originate from latitudes 6– 2 2 ° north, as shown in Figure 4. Therefore, this paper mainly selects this part of the sea area with a range of 125– 13 0 ° east longitude and 6– 2 2 ° north latitude as the area to be covered for task design. According to the target area representation method described in Section 1, the areas of 125– 13 0 ° east longitude and 6– 2 2 ° north latitude are uniformly dispersed into 102 target points according to the condition that it meets D max = 50 km , mathematically expressed as ( λ i j , ϕ i j ) = ( 125 + i , 6 + j ) , i = { 0 , 1 , … , 5 } and j = { 0 , 1 , … , 16 } .

Figure 4

The distribution map of the annual average frequency of tropical cyclones in the Northwest Pacific from 1990 to 2009, design coverage area for the task within the red box (Liu and Liu 2012).

The mission starts at 0:00 on January 1, 2020. Through the optimization of the design, the observation mission of the aforementioned areas can be realized by determining the common ground trajectory constellation scheme composed of nine small satellites. The difference Δ Ω between adjacent track polygons can be determined to be 120.32867 6 ° by Eq. (3) based on the 8 h revisit interval constraint. Here, the first satellite covering the longitude belt of 125 and 12 6 ° east longitude is selected as the reference satellite, and the orbit parameters of the entire common ground trajectory constellation are obtained, as shown in Table 1(a)–(c).

If in the prime of tropical cyclone activity, existing satellites need to be dispatched to supplement the designed common ground trajectory constellations to reduce revisit intervals to capture tropical cyclone activity in a timely manner. The orbital parameters of the four existing satellites to be dispatched are shown in Table 2. Here, take the strip of 6– 9 ° north latitude in the longitude belt of 12 7 ° east as an example, and according to the time window of selecting the satellite to be dispatched, the time window of the target area covered by the sub-satellite point of the Old_1 satellite is more consistent with the time range required in the scheduling task, and the fuel consumption is the least after the simulation comparison. Therefore, the Old_1 satellite is selected as the satellite to be scheduled.

Using the existing satellite dispatch optimization model constructed in Section 2, two-pulse maneuvering orbit changes for Old_1 satellite. It can obtain the orbital parameters of the satellite at 3:46:11.862173 on January 1, 2020, as presented in Table 3. After the dispatch of Old_1 satellite, the sub-satellite point will cover the strip area of 6– 2 2 ° north latitude in the 12 7 ° longitude belt east longitude at about 4 o’clock every day, which reduces the observation interval of the area to a certain extent.

5 Conclusion

In this paper, the constellation optimization design and satellite optimization scheduling problem under the condition of regional target coverage are studied. To cover the area most frequently monitored by tropical cyclone activity in the northwestern Pacific Ocean, the mission was designed for this area, which ranged from 125 to 13 0 ° east longitude and 6 to 2 2 ° north latitude. Simulation analysis results show that the designed low-polar orbit common ground trajectory constellation and the optimized dispatch satellite can provide good coverage performance for the region, and also verify the effectiveness of the regional coverage constellation design optimization model and satellite dispatch optimization model constructed in this paper. Currently, it is also more efficient and meaningful to construct a reasonable optimization model based on evolutionary algorithm for regional coverage constellation design.