# Large region targets observation scheduling by multiple satellites using
resampling particle swarm optimization

Yi Gu, Chao Han, Yuhan Chen, Shenggang Liu, and Xinwei Wang Yi Gu, Chao Han,
and Shenggang Liu are with the School of Astronautics, Beihang University,
Beijing, 100191, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>liushg@buaa.edu.cn).Yuhan Chen is with the China Satellite Network Innovation
Co., Ltd, Beijing, 100029, China (e-mail: chenyuhan@buaa.edu.cn).Xinwei Wang
is with the Department of Transport and Planning, Delft University of
Technology, Delft 2628 CD, The Netherlands (e-mail:
xinwei.wang.china@gmail.com). (Corresponding authors: Shenggang Liu, Xinwei
Wang.)

###### Abstract

The last decades have witnessed a rapid increase of Earth observation
satellites (EOSs), leading to the increasing complexity of EOSs scheduling. On
account of the widespread applications of large region observation, this paper
aims to address the EOSs observation scheduling problem for large region
targets. A rapid coverage calculation method employing a projection reference
plane and a polygon clipping technique is first developed. We then formulate a
nonlinear integer programming model for the scheduling problem, where the
objective function is calculated based on the developed coverage calculation
method. A greedy initialization-based resampling particle swarm optimization
(GI-RPSO) algorithm is proposed to solve the model. The adopted greedy
initialization strategy and particle resampling method contribute to
generating efficient and effective solutions during the evolution process. In
the end, extensive experiments are conducted to illustrate the effectiveness
and reliability of the proposed method. Compared to the traditional particle
swarm optimization and the widely used greedy algorithm, the proposed GI-RPSO
can improve the scheduling result by 5.42% and 15.86%, respectively.

###### Index Terms:

large region targets, multiple satellites, observation scheduling, resampling
particle swarm optimization

##  I Introduction

Earth observation satellites (EOSs) play a crucial role in acquiring images of
specific areas of the Earth [1]. Taking advantages of space-based location and
rapid revisit ability, EOSs can collect information on the Earth’s surface
conveniently by the observation payload. The last decades have seen a growing
trend towards observation requests for the large region target, whose boundary
size is much larger than the imaging width of EOSs. Thus cooperative
observation by multiple EOSs is necessary for large region target observation,
which is essential to environmental monitoring, crop survey, sea
reconnaissance, and other practical areas [2]. The large region targets
observation scheduling problem has attracted increasing attention due to its
significant economic and social value. Meanwhile, the growing number of EOSs
[3, 4] makes it possible for the cooperation observation by multiple EOSs.

According to the relative size between the target and the field of view of the
satellite, the observation targets on the Earth can be divided into point
targets and area targets. To further distinguish from the area target which
can be fully observed during only one pass of satellite [5], we define the
large region target which can only be completely observed by multiple EOSs
cooperative observations. Figure 1 shows the satellite observation for point
and area targets on the Earth. An EOS could generate a long observation strip
whose width depends on the altitude and the field of view of the satellite.
Generally, the size of point targets and the small area target is less than
the width of the observation strip and they can be successfully observed
during one imaging. On the contrary, the large region target enclosed by the
red boundary can not be covered by the observation strip depicted in Figure 1.
This study will focus on the large region targets observation scheduling
problem by multiple EOSs.

Figure 1: Schematic of the imaging for different targets on the Earth.

###  I-A Related work

TABLE I: Overview of the EOS scheduling with respect to region target. Satellite | The classification of observation targets
---|---
Single region | Multiple regions
Single EOS |  Exact [6, 7], Heuristic [8, 9] | $-$
Multiple EOSs |  MOEA [10, 11], Heuristic [12, 13, 14, 15] |  MOEA [16, 17], Heuristic [18, 19, 20, 21, 22, 23]

A considerable amount of literature has been published on the point targets
observation scheduling problem over the past decades. Existing algorithms to
solve the Earth observation scheduling problem can be divided into three
types: the exact algorithm [24, 25, 26], heuristic [27, 28, 29, 30], and
multi-objective evolutionary algorithm (MOEA) [31]. Although region targets
can be decomposed into point targets, the scale of the scheduling problem
would increase dramatically with respect to large region targets. Meanwhile,
the dramatic increase in computational complexity prevents the application of
the observation scheduling method for point targets in the scheduling problem
of observing large region targets.

Existing studies on the region target observation scheduling problem by EOSs
are comparatively limited. In order to make a clear description of the
relevant literature, we summarize previous studies in Table I in line with
applied scheduling algorithms. The single satellite scheduling problem for a
region target was firstly addressed in [8], where a set covering model was
established. Lemaître and Verfaillie [6] introduced an integer linear
programming for a single agile satellite scheduling problem and solved the
constraint programming model with a commercial solver. Mancel and Lopez [7]
adopted a column generation algorithm to solve the integer programming model
for scheduling the PLEIADES satellite of France. The observation coverage time
was optimized by the particle swarm optimization (PSO) algorithm hybridized
with a differential evolutionary algorithm. To accomplish the observation for
an area target by a single satellite, Du $et\ al$. [9] transformed the
original problem into a path planning problem for visiting nodes and developed
an improved ant colony algorithm to solve it. Vazquez $et\ al$. [32] proposed
an automated satellite antenna assignment algorithm for observation data
transmission.

Due to the great advantage of multiple EOSs cooperation, most research has
been conducted on the multiple EOSs scheduling problem. For observing a single
region target, Wang $et\ al$. [10] constructed a multi-objective optimization
model and adopted a strength Pareto evolutionary algorithm to obtain the
scheduling result. For disaster emergency response, Niu $et\ al$. [11] solved
the satellite scheduling problem of a large region task by using a multi-
objective genetic algorithm. Four optimization objectives including the
imaging finish time, the coverage, the mean spatial resolution, and the
average slewing angle have been taken into consideration. The 2008 Wenchuan
earthquake was taken as a real disaster scenario to evaluate the performance
of the algorithm. For the rapid observation of a target area in an emergency,
Zhu $et\ al$. [12] have taken the satellite scheduling and orbital transfer
problems into consideration, simultaneously. In order to maximize the
observation coverage, Chen $et\ al$. [13] proposed a hybrid optimization
method based on the PSO and genetic algorithm. Perea $et\ al$. [14] modeled
the scheduling problem as a set covering problem and presented an adaptive
search algorithm to solve it. He $et\ al$. [15] constructed an integer
programming model for the area target scheduling problem, which aims at
maximizing the coverage ratio while minimizing the response time by
considering the weighted sum. A balanced heuristic has been presented to solve
the model and the constant-factor performance guarantee can be provided.

As for the multi-region observation, Li $et\ al$. [16] introduced the
preference incorporation to the region targets observation scheduling problem
and established a multi-objective optimization model considering the profit,
quality, and timeliness simultaneously. In order to maximize the target area
coverage and minimize the observation resource utilization, Chen $et\ al$.
[17] adopted a non-dominated sorting genetic algorithm to solve the model.
Wang $et\ al$. [18] established a nonlinear model and presented a heuristic
for solving the disaster monitoring observation scheduling problem by four
satellites. For the polygon region observation request, Zhu $et\ al$. [19]
proposed a three-phase solution method that integrates with a dynamic greedy
algorithm and a tabu search algorithm. Moreover, the scheduling problem for
the synthetic aperture radar satellite, which owns a specific field of view,
has been researched for area observing [20, 21]. Focusing on large region
target observation, Xu $et\ al$. [22] and E $et\ al$. [23] developed improved
genetic algorithms to solve the multi-satellite scheduling problem, aiming at
improving the summed observation profit.

The calculation of observation coverage is a significant component in the
observation scheduling of large region targets. Two calculation methods are
mainly utilized to determine the coverage in previous studies. The grid
discretization-based area calculation (GDAC) method is widely adopted owing to
its good performance in precision [22, 23]. While concerning large region
targets, we have to strike a balance between the computational efficiency and
the calculation accuracy of the GDAC. Meanwhile, the polygon clipping
technique for plane polygon intersection operation has also been utilized to
conduct the coverage analysis [33]. It is worth pointing out that a general
projection from a spherical to a plane, serving as the basis of clipping, will
result in relatively large deformation.

###  I-B Objective and contributions

Although the above research has been carried out on area target observation by
EOSs, few studies pay sufficient attention to the large region target except
for [11, 17, 22, 23]. A possible explanation might be that there exist
extensive candidate tasks derived from a large region target, which makes the
coverage calculation and scheduling process more difficult. However, there
still exist two challenges for large region targets observation scheduling.
First, existing studies generally calculate the plane intersection between the
large region and observation strips instead of the spherical intersection. The
equidistant cylindrical projection in [23] is reasonable when the region
target is small. With respect to the large region target or high latitude
region, the area boundary would be greatly deformed during the projection
between plane and sphere. This results in inaccurate calculation of the
intersection, and even affects the successful execution of the observation.
Second, the approximate coverage calculation by counting the number of grid
points is relatively time-consuming. The calculation complexity would be
sharply increased when larger region targets are required to observe.

The large region targets observation scheduling is additionally complicated by
extensive observation strips and the considerable solution space. This study,
therefore, dedicates to investigating a rapid and accurate coverage
calculation method and then establishing an efficient scheduling algorithm for
the large region target observation by multiple EOSs. Unlike existing studies,
a coverage calculation method based on a polygon clipping technique [34] is
proposed to substitute the GDAC method. To eliminate the calculation bias
between a plane and a spherical surface, a projection reference plane is
defined for each large region target and the center projection method is
adopted. The precise intersection point between the region target and
satellite strips can then be obtained through the proposed projection
transformation and inverse transformation. The presented projection reference
plane integrated polygon clipping technique can be denoted as the amended
polygon clipping technique (APCT) in this study. On the basis of the proposed
coverage calculation method, a nonlinear integer programming model is
constructed for maximizing the total observation area and satisfying
engineering constraints. Given the multiple EOSs scheduling has been proved
NP-hard [23, 35], we propose an improved resampling particle swarm
optimization (RPSO) algorithm to solve the established model. The algorithm
based on particle searching has been demonstrated as well-performance in the
large-scale optimization problem [36]. The greedy initialization strategy is
designed to enhance the convergence speed while the presented resampling
method has the potential to obtain better solutions. Moreover, a particle
reconstruction method has been adopted to handle several engineering
constraints for practical applications.

The main contributions of this study can be summarized as twofold: (1) We
propose an accurate coverage calculation method based on the specially
designed projection reference plane and polygon clipping technique, which
ensures efficient objective function calculation of the scheduling problem.
(2) An improved RPSO scheduling algorithm, integrated with the greedy
initialization strategy and the particle reconstruction method, has been
developed to enhance the large region targets scheduling result.

The remaining part of this paper proceeds as follows. Section II describes the
region coverage calculation method. In Section III, the nonlinear integer
mathematical model is constructed in line with necessary assumptions and
notations. Section IV presents the overall framework of RPSO, followed by a
detailed design of each algorithm component. Experiment setting and simulation
results are provided in Section V. Finally, the paper ends with conclusions in
Section VI.

##  II Coverage calculation method based on polygon clipping

With respect to the large region targets observation problem, we aim to
maximize the total coverage for the regions to be observed [23]. Efficient and
accurate coverage calculation of each area is thus essential. In this section,
we first describe the observation scheduling problem for large region targets,
and then define the projection reference plane. Based on the introduced
polygon clipping technique, the coverage calculation method can be finally
provided.

###  II-A Problem description

Figure 2: Schematic of satellite imaging and observation strips.

As seen in Figure 2, an EOS is flying over a large region target, and the
visible time window for the target can be determined according to the area
boundaries and the satellite maneuverability. Generally, the field of view
angle of the payload is limited, but the satellite can adjust the orientation
of the payload through attitude maneuvering so as to increase the observable
area. Taking the satellite rolling ability into consideration, the visible
area along with the satellite ground track can be uniquely defined. The
boundary of the visible area along the ground track corresponds to the visible
start time and end time, respectively. Meanwhile, the boundaries of the
visible area perpendicular to the ground track are tangent to the boundary of
the region target. Available observation strips are discretely distributed
within the visible area and the roll angle interval between two adjacent
observation strips is denoted as $\delta\phi$. Note that only one observation
strip can be collected during one passing of a satellite. As shown in Figure
2, each observation strip owns a specific observation duration which
corresponds to varying degrees of coverage. Thus the selection of an
observation strip plays a key role in the scheduling problem.

Note that although Figure 2 is displayed on a plane surface for the sake of
simplification, the division of observation strips and coverage calculation
actually are performed on the spherical surface. By dividing the observation
strip according to a certain time duration, each independent observation strip
can be approximated as a spherical polygon on the surface of the Earth.
Meanwhile, large region targets in this study are defined as spherical
polygons of the Earth. For regions with other different definitions, they can
generally be approximated as spherical polygons through discretization. Then
the coverage calculation can be transformed into the problem of intersection
and merging of spherical polygons. The intersection strip represents the
intersection area of the observation strip and the large region target. Based
on the polygon clipping technique, boundary points after merging different
intersection strips can be obtained to support the coverage calculation.

###  II-B Projection reference plane

Calculating the coverage of satellite observation strips plays a critical role
in the observation scheduling of large region targets. There exist two main
methods for computing intersection strips: plane polygons and spherical
polygons. Owing to the fact that the intersection computing based on spherical
polygons is time-consuming and inefficient [22], the intersection on plane
polygons is utilized in this study.

To reduce the distortion during the projection transformation from sphere to
plane, a projection reference plane is defined for each large region target.
This means that observation strips and region targets will be projected on the
reference plane for intersection and merging computing. Meanwhile, a central
projection rule which takes the Earth center as the projection reference point
is adopted. Thus the spherical polygon boundary can be accurately projected as
a straight line on the projection plane.

Figure 3: Schematic of projection from Earth surface to reference plane.

As shown in Figure 3, two spherical polygons in different colors are projected
onto the reference plane. The blue polygon indicates a large region target
while the red one could denote an observation strip. Points $B$ and $D$
represent an area boundary point and an intersection of two polygons,
respectively. Note that all intersection points could be determined accurately
in this way. The projection reference plane should be determined as close to
the large region target as possible to reduce deformation. A general framework
to select the projection reference plane for each large region target is
proposed as follows.

  1. 1.

The set of boundary points of a large region target is denoted as $P$, and
$|P|$ indicates the set size. Coordinates of each boundary point in the Earth-
Centered Earth-Fixed (ECEF) coordinate system can be denoted as
$(x_{p},y_{p},z_{p}),p\in P$. The arithmetic average center of these boundary
points can be referred to as $C$. Denote the vector from Earth center to point
$C$ as $\bm{R}_{c}=(x_{c},y_{c},z_{c})$, which can be formulated as

$\displaystyle\left\\{{\begin{array}[]{*{20}{c}}x_{c}=\sum\limits_{p\in{P}}x_{p}/|P|\\\
y_{c}=\sum\limits_{p\in{P}}y_{p}/|P|\\\
z_{c}=\sum\limits_{p\in{P}}z_{p}/|P|\end{array}}\right.$ (4)

  2. 2.

To describe the reference plane, a cartesian coordinate system $S_{c}$ is
established with point $C$ as the origin. Three unit basis vectors of $S_{c}$
can be denoted as

$\displaystyle\left\\{{\begin{array}[]{*{20}{c}}\bm{i}_{c}=\dfrac{\bm{z}_{f}\times\bm{R}_{c}}{\|\bm{z}_{f}\times\bm{R}_{c}\|}\\\
\bm{j}_{c}=\bm{k}_{c}\times\bm{i}_{c}\\\
\bm{k}_{c}=\bm{R}_{c}/\|\bm{R}_{c}\|\end{array}}\right.$ (8)

where $\bm{z}_{f}$ indicates the unit vector of the $z$-axis in ECEF.

  3. 3.

The projection reference plane can be set as $z=\|\bm{R}_{c}\|$ in $S_{c}$.
For arbitrary spherical polygon boundary point $B$, denote
$\bm{R}_{B}=(x_{B},y_{B},z_{B})$ as the vector from Earth center to point $B$
in ECEF coordinate system. Then, point $B$ can be projected onto the reference
plane as
$(x_{B},y_{B},z_{B})\rightarrow(\widetilde{x}_{B},\widetilde{y}_{B})$, and the
projection coordinate in $S_{c}$ can be calculated as

$\displaystyle\left\\{{\begin{array}[]{*{20}{c}}\widetilde{x}_{B}=\dfrac{\bm{i}_{c}\cdot\bm{R}_{B}}{\|\bm{i}_{c}\cdot\bm{R}_{B}\|}\cdot\|\bm{R}_{c}\|\\\
\widetilde{y}_{B}=\dfrac{\bm{j}_{c}\cdot\bm{R}_{B}}{\|\bm{j}_{c}\cdot\bm{R}_{B}\|}\cdot\|\bm{R}_{c}\|\end{array}}\right.$
(11)

Then the intersection computing based on the projection reference plane can be
conducted utilizing the polygon clipping technique.

###  II-C Polygon clipping technique

Recent decades have seen an extensive application of polygon clipping in
industrial manufacturing and computer-aided design [37, 38]. In this study,
both the observation strip of EOSs and the large region target are not
restricted as the convex polygon. Therefore, the general efficient polygon
clipping technique [39, 34] which supports convex and concave polygons
simultaneously is utilized to conduct the intersection and merging computing.
The C++ version of the polygon clipping code can be obtained from the website:
http://www.angusj.com.

The utilization of the polygon clipping technique mainly consists of two
stages. As depicted in Figure 4, all intersection points between observation
strips and large region targets on the projection reference plane can be
obtained in the first stage. Intersection points among all observation strips
could be calculated in the second stage. Then the union of observation strips
can be generated, which is the area with the red color boundary in Figure 4.
Meanwhile, coordinates of all intersection points in the reference coordinate
system $S_{c}$ can be determined.

Figure 4: Schematic of the intersection and merging of observation strips.

###  II-D Strip coverage calculation

In order to transform the result of intersection and merging operation on the
plane to spherical, it is required to perform an inverse transformation. After
generating boundary points of the union observation strip on the reference
plane, an inverse transformation should be conducted for arbitrary
intersection point $\widetilde{D}$ in Figure 3. It is not difficult to know
that point $D$ is the result of the inverse transformation, which is the
intersection between the Earth’s surface and the line crossing both the
Earth’s center and point $\widetilde{D}$. Denote coordinates of point
$\widetilde{D}$ on the reference plane as
$(\widetilde{x}_{D},\widetilde{y}_{D})$, and the vector from the Earth center
to point $\widetilde{D}$ can be described as

$\bm{R}_{\widetilde{D}}=\widetilde{x}_{D}\bm{i}_{c}+\widetilde{y}_{D}\bm{j}_{c}+\|\bm{R}_{c}\|\bm{k}_{c}$
(12)

Projected coordinates of $\bm{R}_{\widetilde{D}}$ in ECEF can be denoted as
$[x_{Df},y_{Df},z_{Df}]^{T}$, and corresponding coordinates of point $D$ on
the Earth can be formulated as

$\displaystyle\left\\{{\begin{array}[]{*{20}{c}}x_{D}=x_{Df}/Nor\\\
y_{D}=y_{Df}/Nor\\\ z_{D}=z_{Df}/Nor\end{array}}\right.$ (16)
$Nor=\sqrt{{{\left({\frac{{{x_{Df}}}}{{{R_{E}}}}}\right)}^{2}}+{{\left({\frac{{{y_{Df}}}}{{{R_{E}}}}}\right)}^{2}}+{{\left({\frac{{{z_{Df}}}}{{{R_{E}}}}}\right)}^{2}}}$
(17)

where $R_{E}$ indicates the equatorial radius of the Earth. Subsequently, the
calculation method of spherical polygon area [40] can be employed to determine
the strip coverage.

##  III Large Region Targets Observation Scheduling by EOSs

###  III-A Assumptions

This study focuses on solving the observation scheduling problem of large
region targets by multiple EOSs. The following necessary assumptions and
simplifications are adopted for practical applications.

  1. 1.

This study takes large region targets into account which can not be completely
observed in one pass of a satellite.

  2. 2.

The EOS mentioned in this paper refers to the non-agile satellite without
pitching maneuverability.

  3. 3.

Only one observation can be conducted by each satellite at the same time and
the observation task cannot be interrupted.

  4. 4.

The satellite attitude remains unchanged during the observation to ensure
imaging quality.

  5. 5.

Supposing that there exist enough ground stations in the scenario, the data
download, and instruction upload fall outside of the scope of our study.

###  III-B Notations

TABLE II: Notations. Targets
---
$T$ |  Set of large region targets, $T=\\{1,...,|T|\\}$ and $|T|$ is the target set size
$i$ |  Target index, $i\in T$
$i_{s}$ |  Target index of denoting the successor observation of $i$, $i_{s}\in T\cup\\{|T|+1\\}$, in which $|T|+1$ is a dummy target
$i_{1},i_{2}$ |  Target index, $i_{1}\in T,i_{2}\in T$
$A_{i}$ |  The area of large region target $i$, $i\in T$
Satellies & Orbits |
$S$ |  Set of satellites, $S=\\{1,...,|S|\\}$ and $|S|$ is the satellite set size
$j$ |  Satellite index, $j\in S$
$O_{j}$ |  Set of orbits of satellite $j$, $O_{j}=\\{1,...,|O_{j}|\\}$ and $|O_{j}|$ is the orbit set size, $j\in S$
$k,k_{1},k_{2}$ |  Orbit index, $k,k_{1},k_{2}\in O_{j}$
$M_{jk},E_{jk}$ |  Memory and energy capacity of satellite $j$ on orbit $k$, $j\in S,k\in O_{j}$
$et_{jk}$ |  Energy consumption for each unit angle of attitude transition of satellite $j$ on orbit $k$, $j\in S,k\in O_{j}$
$mo_{jk},eo_{jk}$ |  Memory and energy consumption for unit time observation of satellite $j$ on orbit $k$, $j\in S,k\in O_{j}$
$[VS_{ijk},VE_{ijk}]$ |  Visible time window of large region target $i$ by satellite $j$ on orbit $k$, $i\in T,j\in S,k\in O_{j}$
$\delta\phi$ |  Real parameter, indicting the discrete granularity interval of roll angle
Observation strips |
$L_{ijk}$ |  Set of feasible observation strips, $L_{ijk}=\\{1,...,|L_{ijk}|\\}$ and $|L_{ijk}|$ is the number of observation strips of target $i$ by satellite $j$ on orbit $k$, $i\in T,j\in S,k\in O_{j}$
$l,l_{1},l_{2}$ |  Strip index, $l\in L_{ijk}$, $l_{1}\in L_{i_{1}jk_{1}}$, $l_{2}\in L_{i_{2}jk_{2}}$, $i,i_{1},i_{2}\in T$, $j\in S$, $k,k_{1},k_{2}\in O_{j}$
$l_{s}$ |  Strip index of denoting the successor observation of $l$, $l_{s}\in L_{i_{s}jk}$, $i_{s}\in T\cup\\{|T|+1\\},j\in S,k\in O_{j}$
$[OS_{ijk}^{l},OE_{ijk}^{l}]$ |  Observation time window of strip $l$ within target $i$ by satellite $j$ on orbit $k$, $l\in L_{ijk},i\in T,j\in S,k\in O_{j}$
$od_{ijk}^{l}$ |  The observation duration of strip $l$ within target $i$ by satellite $j$ on orbit $k$, $l\in L_{ijk},i\in T,j\in S,k\in O_{j}$
$st_{j}^{l_{1}l_{2}}$ |  Attitude transition time between observation strips $l_{1}$ and $l_{2}$ of satellite $j$, $j\in S$, $l_{1}\in L_{i_{1}jk_{1}}$, $l_{2}\in L_{i_{2}jk_{2}}$
$se_{ijk}^{ll_{s}}$ |  Energy consumption in the transition from strip $l$ to $l_{s}$ of satellite $j$ on orbit $k$, $i\in T,j\in S,k\in O_{j}$, $l\in L_{ijk},l_{s}\in L_{i_{s}jk}$
Decision variable
$x_{ijk}^{l}$ |  Binary decision variable. $x_{ijk}^{l}=1$ if strip $l$ within target $i$ is scheduled to be observed by satellite $j$ on orbit $k$, otherwise $x_{ijk}^{l}=0$, $l\in L_{ijk},i\in T,j\in S,k\in O_{j}$

For ease of reference, some notations in this study are summarized in Table
II. Let $T$ and $S$ be the set of large region targets and satellites,
respectively. For each satellite $j\in S$, $O_{j}$ indicates the set of orbits
of satellite $j$ during the given scenario. Each orbit $k\in O_{j}$ is
associated with five parameters: $M_{jk}$ indicting the maximum memory
capacity, $E_{jk}$ denoting the maximum energy capability, $et_{jk}$ denoting
the energy consumption for a unit angle of attitude transition, $mo_{jk}$ and
$eo_{jk}$ representing the memory and energy consumption for unit time
observation, respectively. Moreover, the visible time window
$[VS_{ijk},VE_{ijk}]$ denotes visible start and end time of large region
target $i$ by satellite $j$ on orbit $k$, corresponding to a specific visible
area. A visible area can be obtained according to the intersection of large
region target $i$ and the satellite’s field of view within
$[VS_{ijk},VE_{ijk}]$. Then the visible area can be divided into several
observation strips relying on $\delta\phi$, namely the discrete interval of
the roll angle. Let $L_{ijk}$ denote the set of available observation strips
formulated by the division of the visible area. The observation time window
$[OS_{ijk}^{l},OE_{ijk}^{l}]$ can be calculated for each strip $l\in L_{ijk}$.
The observation duration of strip $l$ by satellite $j$ on orbit $k$ is
described as $od_{ijk}^{l}=OE_{ijk}^{l}-OS_{ijk}^{l}$.

When a satellite accomplishes an observation for strip $l_{1}\in
L_{i_{1}jk_{1}}$, an attitude transition process is needed before an
observation for strip $l_{2}\in L_{i_{2}jk_{2}}$. Notice the required
transition time $st_{j}(i_{1},k_{1},l_{1},i_{2},k_{2},l_{2})$ from $l_{1}$ to
$l_{2}$ is also related to target index $i_{1},i_{2}$ and orbit index
$k_{1},k_{2}$. We denote the attitude transition time between adjacent
observation strips as $st_{j}^{l_{1}l_{2}}$ for simplicity afterward. Taking
the attitude stabilization process into account [41, 42], the attitude
transition time can be formulated as

$st_{j}^{l_{1}l_{2}}=\displaystyle{\frac{|\phi_{i_{2}}^{l_{2}}-\phi_{i_{1}}^{l_{1}}|}{v_{j}^{roll}}}{\rm{+}}\left\\{{\begin{array}[]{*{20}{c}}5&{}\hfil&{|\phi_{i_{2}}^{l_{2}}-\phi_{i_{1}}^{l_{1}}|\leq
15^{\circ}}\\\
{10}&{}\hfil&{15^{\circ}<|\phi_{i_{2}}^{l_{2}}-\phi_{i_{1}}^{l_{1}}|\leq
40^{\circ}}\\\
{15}&{}\hfil&{40^{\circ}<|\phi_{i_{2}}^{l_{2}}-\phi_{i_{1}}^{l_{1}}|}\end{array}}\right.$
(18)

where $v_{j}^{roll}$ represents the attitude maneuvering angular velocity of
the roll axis for satellite $j$, $\phi_{i_{2}}^{l_{2}}$ and
$\phi_{i_{1}}^{l_{1}}$ indicate roll angles for observing strips $l_{1}$ and
$l_{2}$, respectively. The attitude maneuvering process from successive
observation strips $l$ to $l_{s}$ by satellite $j$ on orbit $k$ consumes
energy $se_{jk}(i,i_{s},l,l_{s})$, which can be simplified as
$se_{ijk}^{ll_{s}}$. In order to maintain satellite energy security, we assume
that the roll angle of EOSs at the beginning of each orbital period is zero
when calculating the total energy consumption on each orbit. If target $i$ is
the initial or ending observation on orbit $k$, therefore, corresponding
energy consumption during attitude transition can be defined as
${et_{jk}\cdot|\phi_{i}^{l}|}$. The energy consumption during each transition
process can be calculated as

$se_{ijk}^{ll_{s}}=\left\\{{\begin{array}[]{*{20}{c}}{et_{jk}\cdot|\phi_{i}^{l}|}&{}\hfil&{i\text{
is the initial observation}}\\\ {et_{jk}\cdot|\phi_{i}^{l}|}&{}\hfil&{i\text{
is the ending observation}}\\\
{et_{jk}\cdot|\phi_{i_{s}}^{l_{s}}-\phi_{i}^{l}|}&{}\hfil&{\text{otherwise}}\end{array}}\right.$
(19)

where $\phi_{i}^{l}$ and $\phi_{i_{s}}^{l_{s}}$ represent observation roll
angles for strips $l\in L_{ijk}$ and $l_{s}\in L_{i_{s}jk}$, respectively.

###  III-C Mathematical model

Based on the above notations and equations, the mathematical model for the
large region targets observation scheduling problem by multiple EOSs can be
constructed as follows.

max $\displaystyle\sum\limits_{i\in{T}}\frac{\xi_{i}}{A_{i}}\cdot
f\Big{(}{\color[rgb]{0,0,0}\bigcup_{{j\in{S},k\in{O_{j}},l\in{L_{ijk}}}}x_{ijk}^{l}}\Big{)}$
(20) subject to $\displaystyle\sum\limits_{l\in{L_{ijk}}}x_{ijk}^{l}\leq 1$
$\displaystyle\forall i\in T,j\in S,k\in O_{j}$ (21)
$\displaystyle\sum\limits_{i\in{T}}\sum\limits_{l\in{L_{ijk}}}(x_{ijk}^{l}\cdot
od_{ijk}^{l}\cdot mo_{jk})\leq M_{jk}$ $\displaystyle\forall j\in S,k\in
O_{j}$ (22)
$\displaystyle\sum\limits_{i\in{T}}\sum\limits_{l\in{L_{ijk}}}x_{ijk}^{l}\cdot(od_{ijk}^{l}\cdot
eo_{jk}+se_{ijk}^{ll_{s}})\leq E_{jk}$ $\displaystyle\forall j\in S,k\in
O_{j}$ (23) $\displaystyle x_{ijk}^{l}\in\\{0,1\\}$ $\displaystyle\forall i\in
T,j\in S,k\in O_{j}$ (24) $\displaystyle\big{\\{}x_{i_{1}jk_{1}}^{l_{1}}\cdot
x_{i_{2}jk_{2}}^{l_{2}}\cdot(OS_{i_{2}jk_{2}}^{l_{2}}-$ $\displaystyle
OE_{i_{1}jk_{1}}^{l_{1}}-st_{j}^{l_{1}l_{2}})\geq 0|\text{if
}OS_{i_{2}jk_{2}}^{l_{2}}\geq OE_{i_{1}jk_{1}}^{l_{1}}\big{\\}}$ (25)
$\displaystyle\forall i_{1},i_{2}\in T,j\in S,k_{1},k_{2}\in O_{j},l_{1}\in
L_{i_{1}jk_{1}},l_{2}\in L_{i_{2}jk_{2}}$

  * •

The objective function (20) is to maximize the total area of the observed
region. The collection of scheduled observation strips within target $i$ is
denoted as
$\Big{(}{\color[rgb]{0,0,0}\bigcup_{{j\in{S},k\in{O_{j}},l\in{L_{ijk}}}}x_{ijk}^{l}}\Big{)}$,
and the total area can be calculated by an implicit function $f$, which
corresponds to the coverage calculation method illustrated in Section II. A
coefficient factor $\xi_{i}$ is defined to balance the area difference of
different regions, which can be described as

$\xi_{i}=\frac{A_{i}}{\sum\limits_{i\in{T}}A_{i}}\times 100$ (26)

Thus the maximum value of the objective function (20) can be normalized as 100
regardless of the number of regional targets in the scenario. Note we can also
prioritize large region targets by introducing a weight coefficient $w_{i}$
for each large region as $\xi_{i}^{{}^{\prime}}=\xi_{i}\cdot w_{i}$. We may
also vary the objective function, and set a hard constraint for observing an
entire large region target, which is out of the scope of this work.

  * •

Constraints (21) regulate that for each region target, at most one available
observation strip can be observed on every satellite orbit, which is
guaranteed by assumptions 3 and 4.

  * •

Constraints (22) restrict that the memory consumption for observations is
always within the memory capacity on each orbit.

  * •

Constraints (23) indicate that the total energy consumption from target
observing and attitude maneuvering cannot exceed the energy security range. As
defined in Equation (19), the initial and ending attitude maneuvering process
on each orbit is considered relative to the case where the roll angle is 0.
Energy constraints (23), namely a relaxation of actual operation constraints,
contribute to modeling the energy consumption on each orbit and enhancing the
energy safety margin.

  * •

Constraints (24) point out that the value of the decision variable should be
restricted as 0 or 1.

  * •

Constraints (25) reveal that the attitude transition time
$st_{j}^{l_{1}l_{2}}$ of satellite $j$ should be less than the observation
time interval between strips $l_{1}$ and $l_{2}$ on arbitrary orbits
$k_{1},k_{2}\in O_{j}$. The attitude transition constraint will be verified
for the scheduling sequence of each satellite.

##  IV Greedy Initialization based RPSO Scheduling Algorithm

Taking considerable observation strips and the highly combinatorial
characteristic of the multi-satellite scheduling problem into account, we have
developed an improved RPSO algorithm to solve the large-scale combinatorial
optimization problem. The PSO algorithm is adopted in this study owing to its
advantage in searching the strip which contributes most to the overall profit
within all available observation strips. In order to improve the solution
efficiently, several modifications in the adaptability of PSO have been
proposed. First, a greedy initialization strategy is integrated to enhance the
fast convergence of the solution with almost no increase in the calculation.
Then an individual reconstruction method is designed to handle complex
operation constraints for generating feasible scheduling solutions.
Furthermore, a resampling mechanism is introduced to the PSO algorithm which
will effectively improve the quality of the solution.

###  IV-A Framework

The primary framework of the improved RPSO algorithm is summarized in Figure
5. Firstly, each scheduling solution can be described as a particle by
introducing the integer coding of particles. Meanwhile, most particles are
randomly generated except for one particle created based on a greedy strategy.
Subsequently, a resampling procedure will be conducted if the resampling
condition is satisfied. The current position and velocity of particle $n$,
namely ${\bm{x}}_{n}^{t}$ and ${\bm{v}}_{n}^{t}$, will be resampled to
overcome the defect of particle hysteresis [43, 44]. According to the PSO
algorithm based on the constriction coefficient [45], the new velocity
$\bm{v}_{n}^{t+1}$ and position $\bm{x}_{n}^{t+1}$ for each particle $n$ will
be updated as follows.

${\bm{v}}_{n}^{t+1}=\chi\left({{\bm{v}}_{n}^{t}+{c_{1}}{r_{1}}\left({{\bm{p}}_{n}^{t}-{\bm{x}}_{n}^{t}}\right)+{c_{2}}{r_{2}}\left({{\bm{p}}_{g}^{t}-{\bm{x}}_{n}^{t}}\right)}\right)$
(27)
$\bm{x}_{n}^{t+1}=\bm{x}_{n}^{t}+\left\lfloor{\bm{v}_{n}^{t+1}}\right\rfloor$
(28)

where $r_{1}$ and $r_{2}$ are random real numbers within [0,1],
${\bm{p}}_{n}^{t}$ indicates the current best position of particle $n$ and
${\bm{p}}_{g}^{t}$ represents the current global best solution. Rounding the
velocity $\left\lfloor{\bm{v}_{n}^{t+1}}\right\rfloor$ could ensure that the
particle position remains as an integer. The constriction coefficient $\chi$
can be determined as

$\chi=\frac{2}{{\left|{2-\varphi-\sqrt{{\varphi^{2}}-4\varphi}}\right|}}$ (29)
$\varphi={c_{1}}+{c_{2}},\varphi>4$ (30)

Generally, $\varphi$ is set as 4.1 while $c_{1}$ and $c_{2}$ are fixed as
2.05.

Figure 5: The framework of the scheduling algorithm for large region targets.

After the velocity and position update process, some particles may not satisfy
the constraints presented in Section III-C. In order to ensure that all
particles are feasible solutions, a particle reconstruction method is adopted
in this study by referring to the individual reconfiguration in [23].
Successively, the objective function (20) will be calculated for each particle
based on the coverage calculation method described in Section II. At the same
time, the local and global optimum solutions will be updated to
${\bm{p}}_{n}^{t+1}$ and ${\bm{p}}_{g}^{t+1}$. Afterward, the termination
condition should be verified and the best solution can be output when the
stopping criterion is met. The detailed algorithm framework is described in
the following sections.

###  IV-B Coding and initialization

Note that the satellite attitude remains unchanged within the observation
duration, which means that only one observation strip can be selected during
one pass of each EOS. An integer coding strategy is adopted to describe each
particle, where each gene position corresponds to a visible time window
$[VS_{ijk},VE_{ijk}]$. As shown in Figure 6, the index value of each
observation strip is assigned to the corresponding gene. By arranging
according to satellite observations, each particle can represent the
scheduling scheme of all satellites. The number of genes within a particle
equals the total number of available observation opportunities. The value of
each gene can take an integer varying from 0 to $|L_{ijk}|$, where the zero
value indicates that no observation strip is scheduled in this visible time
window.

Figure 6: The coding structure of each particle.

Compared to a completely random initialization method, a greedy initialization
strategy for one particle is introduced in our method. A particle is generated
by the greedy strategy while other particles are randomly initialized. In each
satellite crossing, the observation strip with the largest coverage area,
corresponding to the maximum observation duration, will be selected to observe
in the greedy algorithm. This strategy can guarantee a relatively good initial
solution, as well as contribute to the convergence speed of the solution.

###  IV-C Particle reconstruction method

Because constraints (21) and (24) are already satisfied in the previous
heuristic process, there exists no need to handle these constraints in the
particle reconstruction procedure. For each generated particle, constraints
(22)-(23) and (25) still should be verified. As shown in Figure 7, boxes along
the $Time$ axis indicate various observation tasks with different observation
durations, where the same color represents that these tasks belong to the same
satellite or orbit. Attitude transition constraints will be verified for each
satellite in subplot (a) while memory and energy constraints should be
confirmed on each orbit as depicted in subplot (b). Note that attitude
transition constraints of $Sat_{2}$ are all satisfied owing to the sufficient
time intervals between adjacent tasks. With respect to $Sat_{1}$, the time
interval $\Delta t$ between observation tasks S1A and S1B is too nervous to
accomplish the attitude transition. In subplot (b), $Orb_{1}$ and $Orb_{2}$
are within the same orbit set $O_{j}$. Considering the memory and energy
consumption, tasks on $Orb_{1}$ can be successfully observed without violating
any constraint while $Orb_{2}$ can not. Therefore, some observation tasks on
$Orb_{2}$ should be removed to guarantee all resource constraints.

Figure 7: Schematic of particle reconstruction for constraints handling.
Algorithm 1 Particle reconstruction algorithm

0: Current scheduling solution $Particle_{n}$

0: The updated scheduling solution

1: for each satellite $j\in S$ do

2: while (1) do

3: Obtain current scheduling scheme of satellite $j\in S$ according to
$Particle_{n}$;

4: Verify the attitude transition constraint for arbitrary two consecutive
observation tasks;

5: if all constraints are satisfied then

6: break;

7: else

8: for observation strips $\\{l_{1},l_{2}\\}$ not satisfying constraint (25)
do

9: Adjust the observation strip $l_{2}$ to $l_{2}^{{}^{\prime}}$ and update
$Particle_{n}$;

10: Verify the attitude transition constraint for the updated
$\\{l_{1},l_{2}^{{}^{\prime}}\\}$;

11: if the constraint is not satisfied then

12: Remove the later observation task and update $Particle_{n}$;

13: end if

14: end for

15: end if

16: end while

17: for each orbit $k\in O_{j}$ do

18: Obtain observation tasks of satellite $j$ on orbit $k$ according to
$Particle_{n}$;

19: Verify memory constraints (22) and energy constraints (23) on orbit $k$;

20: if all constraints are satisfied then

21: break;

22: else

23: while (1) do

24: Remove the task with the minimum observation duration and update
$Particle_{n}$;

25: Execute Steps 19-21;

26: end while

27: end if

28: end for

29: end for

The detailed process of the particle reconstruction method is presented as
Algorithm 1. After generating a particle, namely a scheduling solution, the
proposed reconstruction algorithm should be conducted for each satellite $j\in
S$. In order to decouple constraints (25) and constraints (22)-(23), the
attitude transition constraint will be handled at first and followed by memory
and energy constraints. The transition constraint will be verified for
arbitrary two consecutive observation tasks. Denote $\\{l_{1},l_{2}\\}$ as two
observation strips of two consecutive tasks. When constraint (25) is not
satisfied for $\\{l_{1},l_{2}\\}$, a new strip $l_{2}^{{}^{\prime}}$ would be
selected. As can be seen from Figure 7(a), $\phi_{A}$ remains unchanged while
$\phi_{B}$ will be adjusted to make $\Delta t$ as large as possible. Then the
attitude transition constraint will be verified for
$\\{l_{1},l_{2}^{{}^{\prime}}\\}$. The observation task S1B will be removed if
constraint (25) is still not satisfied. It is worth noting that whole attitude
transition constraints for satellite $j$ will be verified again because the
adjustment of each task may affect the transition time of adjacent tasks. The
loop will be stopped if all transition constraints are satisfied.

Subsequently, the verification for constraints (22) and (23) on each orbit
$k\in O$ will be conducted. Both memory and energy consumption are all related
to the number of observation tasks and the length of observation duration. By
referring to [23], we will remove the task with the minimum observation
duration at first, namely the less coverage. The observation task on the
current orbit will be removed in this order until memory and energy
constraints are satisfied simultaneously. Finally, the particle satisfying all
constraints can be generated by the reconstruction algorithm.

###  IV-D Resampling strategy

As a core of the proposed improved RPSO algorithm, the resampling strategy
contributes to generating a better scheduling solution within given iteration
steps. In the standard PSO algorithm, the particle velocity will be updated as
Equation (27) and then the new particle position can be generated according to
Equation (28). Generally, the gene value of each newly generated particle
position may be out of the valid range $[0,|L_{ijk}|]$. When the gene value is
larger than $|L_{ijk}|$, it will be fixed as the upper boundary. Meanwhile,
the gene value is set as 0 if a negative number has been assigned to it.
Similarly, each gene value of particle velocity will be limited within the
valid range in the same way.

However, there exists a defect in the standard PSO. Although most particles
will converge towards a certain local optimum, some particles may move to
several non-optimal regions owing to the random number introduced in the PSO
algorithm. This part of particles contributes little to the optimization
process as well as wastes computing resources during the process when they
approach better particles. Therefore, it is necessary to solve the moving
hysteresis defect by the resampling strategy.

Figure 8: Illustration of the resampling strategy. Algorithm 2 Particle
resampling algorithm

0: The position $\bm{x}_{n}^{t}$, velocity $\bm{v}_{n}^{t}$, the objective
function value $F(n)$ for each particle $n\in\\{1,2,...,N\\}$, the current
global best solution ${\bm{p}}_{g}^{t}$;

0: The resampled particle position and particle velocity;

1: if the reampling criteria is met then

2: Initialization: $N^{m}\leftarrow 0$;

3: for each particle $n$ do

4:
${q_{n}}=\dfrac{1}{{\sqrt{2\sigma\pi}}}\exp\left({-\dfrac{{({F(n)-{\bm{p}}_{g}^{t}})^{2}}}{{2\sigma}}}\right)$

5: ${Q_{n}}=q_{n}/\sum\limits_{n=1}^{N}{q_{n}}$

6: end for

7: $\bar{Q}_{n}=\dfrac{1}{N}\sum\limits_{n=1}^{N}{q_{n}}$;

8: for $n\in\\{1,2,...,N\\}$ do

9: if ${Q_{n}}>\bar{Q}_{n}$ then

10: $N_{c}=\left\lfloor{Q_{n}/\bar{Q}_{n}}\right\rfloor$

11: Copy the current particle $N_{c}$ times and set the weight as
$\bar{Q}_{n}$;

12: $N^{m}=N^{m}+N_{c}$

13: end if

14: end for

15: Inherit $N-N^{m}$ particles from $\\{1,2,...,N\\}$ and set the weight as
$\bar{Q}_{n}$;

16: end if

As shown in Figure 8, each dot represents a particle while the arrow indicates
the velocity of each particle. All particles are located in an elliptical area
which denotes the feasible region. The weight of each particle is calculated
and assigned during the process from subplot (a) to (b), where the larger dot
indicates the particle with a higher weight and the hollow circle represents
the particle with a low weight. The resampling process is illustrated in the
variation from subplot (b) to (c), where the particle with a larger weight is
copied twice and the particle with a small weight is discarded to keep the
total number of particles unchanged. After resampling, all particle weights
are restored to the same level. Although the newly generated particles are
completely the same in Figure 8(c), the velocity in the next time step would
be different owing to the random parameters in Equation (27). These generated
particles will search in different directions from the same location, and then
more particles can search for solutions in better areas, which would make up
the defect of particle hysteresis to a certain extent.

The main steps for the resampling algorithm are shown in Algorithm 2. The
resampling criteria here is that the resampling process will be performed
twice after every ten iterations by referring to [36]. Each particle
$n\in\\{1,2,...,N\\}$ would be given a weight value $q_{n}$ according to its
objective function value $F(n)$, which can be determined by Equation (20). The
parameter $\sigma$ indicates the sample variance of ${F(n)-{\bm{p}}_{g}^{t}}$
and $Q_{n}$ represents the unitary weight value of particle $n$. For each
particle $n$, if ${Q_{n}}>\bar{Q}_{n}$ is satisfied, particle $n$ will be
copied $N_{c}$ times as depicted in Figure 8(c) and the weight for each
particle is set as the average value $\bar{Q}_{n}$. The total number of
particles obtained by copying is denoted as $N^{m}$. There are still $N-N^{m}$
particles missing compared to the initial situation. Finally, these particles
will be inherited from original unreplicated particles in descending order of
weight. Particles with smaller weights will be discarded in this procedure.
Through the resampling process, more particles will be concentrated in a
better area, which will facilitate the optimization of the solution.

##  V Computational Experiments

Extensive simulations are conducted to evaluate the proposed coverage
calculation method and the improved RPSO algorithm in this section. Three
objectives are assigned to the numerical experiments: verify the accuracy of
the coverage calculation, evaluate the effectiveness of the improved RPSO
algorithm, and conduct the sensitivity analysis.

###  V-A Data generation

The scheduling algorithm for large region targets is coded in C++, and
experiments are conducted on a desktop with a CPU of 3.20 GHz and a memory of
16.0 GB. Because there lacks the benchmark for the large region targets
observation scheduling problem, we have designed the simulation instances by
referring to [22, 23]. The detailed description of the scenario is as follows.

Preliminary simulation shows that the proposed method can achieve the optimal
solution in small instances. In order to verify the performance of the
proposed method on a practical scale, seven large region targets distributed
in different latitudes of the world are considered in the experiment. The
longitude and latitude information of the seven region targets are provided in
Table III, where the area is provided by the System Tool Kit (STK) software
according to the geographical location of each region. These regional targets
are assumed to be of equal priority in the observation requirement from
customers.

TABLE III: Basic information of seven large region targets. Targets | Latitude & Longitude |  Area($km^{2}$)
---|---|---
$T_{1}$ | (5.33,30.26),(5.33,9.32),(-5.17,9.32),(-5.17,30.26) | 2732840.4
$T_{2}$ | (-0.53,-50.5),(-0.53,-75.26),(-9.06,-75.26),(-9.06,-50.5) | 2628413.6
$T_{3}$ | (32.58,114.91),(32.58,105.09),(22.09,105.09),(22.09,114.91) | 1129008.3
$T_{4}$ | (20.557,111.658),(10.52,111.658),(10.52,119.35),(20.557,119.355) | 916658.2
$T_{5}$ | (34.67,178.48),(24.44,178.48),(24.44,168.15),(34.67,168.15) | 1134530.3
$T_{6}$ | (35.5,85.5),(25.8,85.5),(25.8,95.5),(35.5,95.5) | 1030102.7
$T_{7}$ | (70.5,21.5),(60.8,21.5),(60.8,33.8),(70.5,33.8) | 608088.1

The simulation experiment is conducted from the same start time of 2021-04-07
00:00:00. Twenty satellites in the Chinese satellite platform are selected,
and initial orbital parameters and the maximal roll angle $\phi_{max}$ of
involved EOSs are shown in Table IV. Meanwhile, other performance parameters
of satellites and orbits are set as the same and have been shown in Table V.

TABLE IV: Specific orbital elements and maneuvering ability of the selected satellites. $ID$ |  $a$(km) | $e$ |  $i$(∘) |  $\Omega$ (∘) |  $\omega$ (∘) |  $M$ (∘) |  $\phi_{max}$(∘)
---|---|---|---|---|---|---|---
Sat1 | 7015.341 | 0.001108 | 97.8860 | 339.1370 | 155.8560 | 20.9230 | 30
Sat2 | 7002.399 | 0.001291 | 97.8560 | 54.1270 | 157.1140 | 306.0990 | 35
Sat3 | 7126.453 | 0.001121 | 98.5290 | 91.1400 | 90.2280 | 268.9880 | 30
Sat4 | 7015.610 | 0.000861 | 97.9960 | 140.7450 | 163.0120 | 219.4300 | 30
Sat5 | 7149.376 | 0.001479 | 98.7240 | 117.3840 | 140.6560 | 242.8550 | 32
Sat6 | 6876.495 | 0.001410 | 97.3820 | 96.3900 | 155.0710 | 56.6480 | 30
Sat7 | 6997.446 | 0.002413 | 97.8490 | 215.1140 | 109.3060 | 186.4790 | 30
Sat8 | 7001.883 | 0.001245 | 97.9200 | 235.5390 | 90.2690 | 114.1530 | 30
Sat9 | 7001.856 | 0.001260 | 97.8600 | 229.4870 | 90.0140 | 354.3360 | 30
Sat10 | 6869.771 | 0.001486 | 97.3560 | 350.2770 | 77.9030 | 271.3230 | 30
Sat11 | 6887.457 | 0.001387 | 97.5150 | 161.0020 | 90.5810 | 271.2620 | 30
Sat12 | 6886.176 | 0.001263 | 97.3050 | 300.6480 | 87.0550 | 9.4970 | 30
Sat13 | 6887.253 | 0.001410 | 97.4430 | 170.5300 | 87.5400 | 79.7450 | 30
Sat14 | 7618.359 | 0.002772 | 99.9940 | 240.6510 | 47.7570 | 255.9060 | 30
Sat15 | 7578.863 | 0.000758 | 100.4130 | 7.0110 | 60.8300 | 337.8330 | 30
Sat16 | 7579.543 | 0.000695 | 100.2430 | 301.4140 | 79.4160 | 138.3540 | 30
Sat17 | 7578.169 | 0.000451 | 100.3980 | 345.9620 | 75.0290 | 250.0000 | 30
Sat18 | 7575.721 | 0.001826 | 100.3020 | 293.3490 | 77.6880 | 332.1050 | 30
Sat19 | 7142.139 | 0.000764 | 98.3110 | 301.8740 | 186.2800 | 54.1290 | 30
Sat20 | 7142.139 | 0.000780 | 101.3580 | 223.8170 | 82.2040 | 23.1210 | 30
TABLE V: Constraint parameters for each satellite and orbit. Parameter | Value
---|---
$M_{jk}$ | 30000 MB
$mo_{jk}$ | 150 MB/s
$E_{jk}$ | 40000 J
$eo_{jk}$ | 100 J/s
$et_{jk}$ |  150 J/∘
$v_{j}^{roll}$ |  3 ∘/s

The particle population size generally varies from 20 to 50 in practical
applications [45, 36], and we set the particle size as 50. According to
Section IV-B, the position range of each gene value is [0, $L_{ijk}$], and the
corresponding velocity range is set as
[-$\left\lfloor|L_{ijk}|/10\right\rfloor$,
$\left\lfloor|L_{ijk}|/10\right\rfloor$]. The parameter $\delta\phi=0.001rad$
is initialized to obtain discrete observation strips.

###  V-B Results and discussions

####  V-B1 Coverage calculation

First, the proposed grid-less coverage calculation method in Section II is
verified by comparing the GDAC method and the STK software. As shown in Figure
9, large region targets and observation strips located at different latitudes
are selected in STK for the comparison experiment. Areas enclosed by the blue
border are three large region targets from low latitude to high latitude,
which correspond to $T_{1}$, $T_{6}$, and $T_{7}$ in Table III, respectively.
Observation strips are represented as the area enclosed by the red border and
can be generated arbitrarily to verify the coverage calculation result. All
intersection points between observation strips and region targets can be
obtained accurately by the APCT in this paper. As depicted in Figure 9, we
have added the intersection point to STK according to the latitude and
longitude information computed by APCT. It is not difficult to find that these
intersection points are exactly located on the region boundary, which
indicates that intersection calculation results are consistent with the two-
dimensional display results of STK. The intersection area consists of black
boundaries connecting intersection points, and then the coverage calculation
can be compared.

Figure 9: Three large regions at different latitudes for coverage calculation
comparison.

The coverage can be directly obtained based on the APCT in this paper.
According to the geographical location information of the large region target
and the intersection area, STK can provide the corresponding acreage value,
and then the coverage can be calculated. The GDAC method would count the
number of discrete grid points within the intersection area. The different
dispersion of grid points can be adjusted by STK and the coordinates of each
point are output for the counting by GDAC. Comparison results are shown in
Table VI, where column Dis indicates the different dispersion of the grid and
column Cov indicates the coverage calculated by different methods. What stands
out in the table is that the absolute error between APCT and STK is
exclusively small in all three instances, which are completely negligible in
practical engineering. The program running time for almost all cases by GDAC
is higher than that of APCT in addition to one result highlighted by the
underscore. Meanwhile, the coverage calculation result by GDAC is throughout
not as accurate as APCT, although its running time has exceeded 1000 times
longer than APCT when dispersion equals 0.05∘. This might be explained by the
inherent defect of GDAC that grid points close to the boundary are not
completely inside or outside the area. Then the coverage calculation result is
accurate only when the dispersion is small enough, which results in an
exhaustively increase in running time. Closer inspection of Table VI shows the
relative errors between GDAC and STK are 1.91%, 6.62%, and 9.90% for the three
instances when the dispersion is 1.0∘. This indicates that the coverage
calculation error is expected to decline as the region area increases under
the same degree of dispersion. In conclusion, the presented APCT in this paper
can calculate the coverage accurately within a shorter time, which is
definitely an improvement to the calculation method by grid points.

TABLE VI: Comparison of coverage calculation results of different methods. Method | Dis/∘ | Low latitude region |  | Middle latitude region |  | High latitude region
---|---|---|---|---|---|---
Cov/$\%$ |  Time/$ms$ |  |  Cov/$\%$ |  Time/$ms$ |  |  Cov/$\%$ |  Time/$ms$
GDAC | 1.000 | 16.4021 | 4173.00 |  | 26.3889 | 493.00 |  | 28.8889 | 1210.00
| 0.500 | 16.4903 | 25189.00 |  | 23.8683 | 3390.00 |  | 27.0370 | 3780.00
| 0.100 | 16.7409 | 344398.00 |  | 24.7928 | 61455.00 |  | 26.2234 | 35950.00
| 0.050 | 16.7490 | 1962054.00 |  | 24.6708 | 528109.00 |  | 26.2438 | 127513.00
STK | $-$ | 16.7219 | $-$ |  | 24.7493 | $-$ |  | 26.2875 | $-$
APCT | $-$ | 16.7217 | 1186.00 |  | 24.7494 | 600.00 |  | 26.2908 | 631.00

####  V-B2 Large region targets observation scheduling

The performance of the improved greedy initialization-based RPSO (GI-RPSO)
algorithm is evaluated by comparing the greedy initialization-based PSO (GI-
PSO) algorithm, the RPSO algorithm with no greedy strategy, and the standard
PSO algorithm. The scenario duration is set as 1 day and the maximal iteration
number is initialized as 1000 to validate the algorithm performance. Without
loss of generality, 10 runs for each algorithm have been conducted and the
average objective function value of Equation (20) within the iteration can be
calculated and depicted in Figure 10. Owing to the introduced greedy
initialization strategy, the initial solution of GI-RPSO and GI-PSO is the
same as the greedy result. It is apparent that there is a sharp rise in the
scheduling result for all four algorithms with the increase in iteration
numbers. Though the random initial solutions of PSO and RPSO are much worse
than the greedy result, an opposite situation would arise when the iteration
number exceeds 100. What is interesting in this figure is the rapid
convergence of the scheduling result of GI-RPSO and GI-PSO, which could be
attributed to the introduction of the greedy strategy. From the data point in
Figure 10, it is apparent that the resampling method could contribute to the
improvement of the scheduling solution. Average objective function values of
GI-RPSO and RPSO are consistently higher than that of GI-PSO and PSO,
respectively. Closer inspection of the figure shows that the result of GI-PSO
would converge quickly, which would be trapped in the local optima. A higher
objective function value can be obtained by GI-RPSO because the resampling
strategy could contribute to searching for a better solution. Considering that
the greedy algorithm is widely utilized in practical engineering [22], the
comparison experiment with the greedy method can validate the efficiency of
the proposed algorithm in this study. Meanwhile, the GI-RPSO performs the best
among all five algorithms, which emphasizes the ability of the improved RPSO
algorithm in the large region targets observation scheduling problem. It can
be concluded that the proposed resampling method and greedy initialization
strategy perform well in the improved RPSO algorithm.

Figure 10: Performance comparison of several algorithms at different
iterations.

The distribution of scheduling results of 10 runs at different iteration
numbers is shown in Figure 11, where the green triangle symbol and the orange
horizontal line in the box represent the average and the median of 10 results,
respectively. Each box indicates the distributions, and different surface
colors represent various iteration numbers in line with the legend. The
program execution time of the scheduling algorithm for each run is denoted as
the purple square dot. Overall, there exists no big difference in scheduling
results among different runs. Besides, the introduction of the greedy strategy
might enhance the stability of the algorithm, which could be found by the
comparison of scheduling results between PSO and GI-PSO. It is not difficult
to find that the running times of different algorithms at the same iteration
number are analogous, which indicates that proposed strategies are not time-
consuming. The program running time increases with the growth of iteration
numbers, and the maximal running time does not exceed 50 seconds, which shows
the excellent efficiency of the proposed algorithm.

Figure 11: The comparison result of different algorithms in different
iterations.

To further evaluate the performance of various algorithms, fifteen unique
target sets are generated randomly to conduct the comparison experiment. Each
target set contains fifteen regions, including the previous seven regions and
eight newly generated targets. In line with [15], these eight targets are
randomly generated as a square region with a width of 4∘, located within the
latitude range [60∘S, 60∘N] and the longitude range [180∘W, 180∘E]. For each
target set, 10 runs are conducted and the average results in different
iterations are recorded. The statistical results of fifteen target sets are
presented in Figure 12. We can observe a steady increase in the objective
function value as the iteration number grows with respect to all four
algorithms. The performance of the greedy algorithm is provided by the results
of the first iteration of GI-PSO and GI-RPSO. Note that the PSO and RPSO will
exceed the greedy algorithm after about 50 iterations. The results of GI-PSO
and GI-RPSO at 200 iterations are better than the ultimate performance of PSO
and RPSO, demonstrating the promotion effect of the greedy-initialization
strategy. Meanwhile, it can be seen that the average results of GI-RPSO
outperform other algorithms, which validates the effectiveness of the proposed
heuristic. Considering the performance of GI-RPSO improves little when the
iteration number varies from 600 to 1000, the maximal iteration will be set as
600 in the following experiments.

Figure 12: The statistical result of various target sets in different
algorithms.

The verification for the proposed particle reconstruction method is conducted
by experiments with different constraint parameters. As shown in Table VII,
the value of $mo_{jk}$ can be adjusted and corresponding simulation results by
the GI-RPSO algorithm in different iterations are recorded. Parameters
$mo_{jk}$ vary from 100 to 200, and column IFN indicates the initial feasible
number of particles. It can be found that for the instance with
$mo_{jk}=100$MB/s, all initially generated particles are feasible and the
reconstruction is not performed. When the value of $mo_{jk}$ increases to 125,
there exist 7 initial particles that do not meet constraints. Meanwhile, the
particle generated by the greedy strategy is still feasible in this instance
according to the simulation result at the first iteration. If $mo_{jk}$
exceeds 150 MB/s, however, there exists no feasible particle at the
initialization procedure. Then all particles will be changed to the solution
satisfying all constraints by utilizing the particle reconstruction method.
What can be clearly seen in Table VII is the steady decline of the objective
function value at each iteration with the increase of $mo_{jk}$. There is a
dramatic decline in the average value of the objective function when $mo_{jk}$
varies from 150 to 200, which might be explained by the fact that the energy
constraints (22) would be violated seriously in these instances. It can be
concluded that the constraint is evaluated in this experiment and the proposed
particle reconstruction method performs efficiently.

TABLE VII: Comparison of simulation results under different constraint parameters. $mo_{jk}$ | IFN |  Simulation results in different iteration numbers
---|---|---
1 | 100 | 200 | 300 | 400 | 500 | 600
100 | 50 | 71.574 | 82.704 | 83.555 | 83.719 | 83.880 | 83.973 | 84.007
125 | 43 | 71.574 | 82.145 | 82.935 | 83.186 | 83.325 | 83.368 | 83.392
150 | 0 | 64.022 | 73.449 | 74.128 | 74.326 | 74.466 | 74.635 | 74.712
175 | 0 | 39.343 | 47.620 | 49.119 | 50.076 | 50.620 | 51.072 | 51.390
200 | 0 | 18.485 | 25.608 | 26.757 | 27.458 | 27.887 | 28.223 | 28.446
Figure 13: Sensitivity analysis of coefficients $c_{1}$ and $c_{2}$.

The coefficients $c_{1}$ and $c_{2}$ in Equation (27) could impact the
optimization result of the proposed heuristic. To quantify the impact of
$c_{1}$ and $c_{2}$, a sensitivity analysis has been conducted with respect to
the proposed GI-RPSO algorithm. The seven large regions are scheduled to
observe by the selected 20 satellites, and the number of iterations is set as
600. Following [36], $c_{1}$ and $c_{2}$ are initially both set as 2.05 to
satisfy the constraint that the sum of $c_{1}$ and $c_{2}$ exceeds 4. We then
vary two coefficients from 2.05 to 3.15, with a step of 0.1 for the
experiment. For each combination of $c_{1}$ and $c_{2}$, ten runs are
conducted. Therefore we entirely have $12\times 12\times 10=1440$ runs. The
simulation results are shown in Figure 13, where the objective function value
is reported by a surface plot with the horizontal axes $c_{1}$ and $c_{2}$.
Clearly the objective function value has a decreasing trend with the growth of
$c_{1}$ and $c_{2}$, although results fluctuate in a few cases. It can be
observed that the highest objective function value is obtained when $c_{1}$
and $c_{2}$ are both set as 2.05, which is consistent with the recommended
values in [45].

In order to evaluate the influence of discrete granularity of the roll angle,
a sensitivity analysis of $\delta\phi$ has been conducted. Simulation results
under two different scenario durations are shown in Figure 14, where the
iteration number is taken as 600. Different values of $\delta\phi$ varying
from 0.1 to $0.0005rad$, as depicted in the legend, are chosen for the
comparative experiment. What is interesting in this figure is the steady
growth of the objective function value with the decrease of $\delta\phi$. This
result may be explained by the fact that more available observation strips can
be provided when the value of $\delta\phi$ is set as a smaller value. The
increased solution space may potentially contain better scheduling results,
and then the improved RPSO could search for a better solution compared to the
instance with a higher value of $\delta\phi$. Figure 14 shows that the
difference of results under $\delta\phi$ equals 0.001 and 0.0005 is not
obvious, which indicates that $\delta\phi=0.001rad$ is fine enough to obtain
an excellent scheduling result. Therefore the value of $\delta\phi$ is set as
$0.001rad$, namely about 0.057∘, for the corresponding experiments in this
study. Meanwhile, it is not difficult to find that the program running time
changes insignificantly with the variation of $\delta\phi$, which indicates
the excellent efficiency of the proposed scheduling algorithm. According to
the data provided on the website of Natural Resources Satellite Remote Sensing
Cloud Service Platform [46], the attitude pointing accuracy of the high-
resolution multi-mode satellite is better than 0.01∘. Owing to the fact that a
considerable degree of pointing accuracy can be achieved by a practical
satellite, the parameter value set in this article is reasonable.

Figure 14: Sensitivity analysis of $\delta\phi$ under different scenario
durations.

##  VI Conclusions

We have addressed the large region targets observation scheduling problem by
multiple EOSs in this work. An efficient coverage calculation method has been
developed based on the polygon clipping technique. Meanwhile, an improved RPSO
scheduling algorithm has been proposed to solve the multiple large region
targets scheduling problem, where the greedy initialization strategy and the
resampling method perform significant roles. Comparative experiments show that
the proposed coverage calculation method with respect to the large region
exhibits better performance than the GDAC method in terms of calculation
accuracy and efficiency. This approach will be useful in improving the
calculation ability of large region coverage problems in practical
engineering. In the scheduling process, the greedy initialization strategy has
been proven effective in enhancing the convergence speed and solution quality
of the optimization algorithm. The resampling method contributes to searching
for a better solution, and the proposed GI-RPSO outperforms other heuristics,
especially the greedy algorithm widely used in practice. Extensive experiments
have been conducted to verify the effectiveness and stability of the proposed
algorithm.

Further studies regarding the role of agile EOSs in the large region targets
observation would be worthwhile. Compared to the non-agile EOSs, agile
satellites possess stronger attitude maneuverability, leading to more
observation opportunities. In addition, the cloud occlusion to optical sensors
could be taken into consideration for future research.

## References

  * [1] X. Wang, G. Wu, L. Xing, and W. Pedrycz, “Agile Earth observation satellite scheduling over 20 years: formulations, methods and future directions,” _IEEE Systems Journal_ , vol. 15, no. 3, pp. 3881–3892, 2021.
  * [2] X. Wang, Y. Gu, G. Wu, and J. R. Woodward, “Robust scheduling for multiple agile Earth observation satellites under cloud coverage uncertainty,” _Computers & Industrial Engineering_, vol. 156, p. 107292, 2021.
  * [3] Y. Gu, C. Han, Y. Chen, and W. W. Xing, “Mission replanning for multiple agile Earth observation satellites based on cloud coverage forecasting,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 15, pp. 594–608, 2021.
  * [4] G. Wu, Q. Luo, X. Du, Y. Chen, P. N. Suganthan, and X. Wang, “Ensemble of meta-heuristic and exact algorithm based on the divide and conquer framework for multi-satellite observation scheduling,” _IEEE Transactions on Aerospace and Electronic Systems_ , 2022.
  * [5] Z. Chang, Z. Zhou, L. Xing, and F. Yao, “Integrated scheduling problem for Earth observation satellites based on three modeling frameworks: an adaptive bi-objective memetic algorithm,” _Memetic Computing_ , vol. 13, no. 2, pp. 203–226, 2021.
  * [6] M. Lemaîitre and G. Verfaillie, “Daily management of an Earth observation satellite: comparison of ILOG solver with dedicated algorithms for valued constraint satisfaction problems,” 1997.
  * [7] C. Mancel and P. Lopez, “Complex optimization problems in space systems,” in _13th International Conference on Automated Planning & Scheduling (ICAPS’03)_, 2003.
  * [8] J. T. Walton, “Models for the management of satellite-based sensors,” Ph.D. dissertation, Massachusetts Institute of Technology, 1993.
  * [9] B. Du, S. Li, Y. She, W. Li, H. Liao, and H. Wang, “Area targets observation mission planning of agile satellite considering the drift angle constraint,” _Journal of Astronomical Telescopes, Instruments, and Systems_ , vol. 4, no. 4, p. 047002, 2018.
  * [10] J. Wang, J. Li, J. Chen, Y. GUO, and N. Jing, “An approach for multiobjective uniting imaging scheduling of Earth observing satellites,” _Journal of Astronautics_ , vol. 28, pp. 354–193, 2007.
  * [11] X. Niu, H. Tang, and L. Wu, “Satellite scheduling of large areal tasks for rapid response to natural disaster using a multi-objective genetic algorithm,” _International journal of disaster risk reduction_ , vol. 28, pp. 813–825, 2018.
  * [12] K.-J. Zhu, J.-F. Li, and H.-X. Baoyin, “Satellite scheduling considering maximum observation coverage time and minimum orbital transfer fuel cost,” _Acta Astronautica_ , vol. 66, no. 1-2, pp. 220–229, 2010.
  * [13] Y. Chen, D. Zhang, M. Zhou, and H. Zou, “Multi-satellite observation scheduling algorithm based on hybrid genetic particle swarm optimization,” in _Advances in Information Technology and Industry Applications_. Springer, 2012, pp. 441–448.
  * [14] F. Perea, R. Vazquez, and J. Galan-Viogue, “Swath-acquisition planning in multiple-satellite missions: An exact and heuristic approach,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 51, no. 3, pp. 1717–1725, 2015.
  * [15] L. He, B. Liang, J. Li, and M. Sheng, “Balancing coverage and response time in area target scheduling for satellite networks,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 6, pp. 6848–6853, 2020.
  * [16] L. Li, Y. Wang, H. Trautmann, N. Jing, and M. Emmerich, “Multiobjective evolutionary algorithms based on target region preferences,” _Swarm and Evolutionary Computation_ , vol. 40, pp. 196–215, 2018.
  * [17] Y. Chen, M. Xu, X. Shen, G. Zhang, Z. Lu, and J. Xu, “A multi-objective modeling method of multi-satellite imaging task planning for large regional mapping,” _Remote Sensing_ , vol. 12, no. 3, p. 344, 2020.
  * [18] P. Wang, G. Reinelt, P. Gao, and Y. Tan, “A model, a heuristic and a decision support system to solve the scheduling problem of an Earth observing satellite constellation,” _Computers & Industrial Engineering_, vol. 61, no. 2, pp. 322–335, 2011.
  * [19] W. Zhu, X. Hu, W. Xia, and H. Sun, “A three-phase solution method for the scheduling problem of using Earth observation satellites to observe polygon requests,” _Computers & Industrial Engineering_, vol. 130, pp. 97–107, 2019.
  * [20] H. Kim and Y. K. Chang, “Mission scheduling optimization of SAR satellite constellation for minimizing system response time,” _Aerospace Science and Technology_ , vol. 40, pp. 17–32, 2015.
  * [21] ——, “Optimal mission scheduling for hybrid synthetic aperture radar satellite constellation based on weighting factors,” _Aerospace Science and Technology_ , vol. 107, p. 106287, 2020.
  * [22] Y. Xu, X. Liu, R. He, and Y. Chen, “Multi-satellite scheduling framework and algorithm for very large area observation,” _Acta Astronautica_ , vol. 167, pp. 93–107, 2020.
  * [23] E. Zhibo, R. Shi, L. Gan, H. Baoyin, and J. Li, “Multi-satellites imaging scheduling using individual reconfiguration based integer coding genetic algorithm,” _Acta Astronautica_ , vol. 178, pp. 645–657, 2021.
  * [24] M. Lemaître, G. Verfaillie, F. Jouhaud, J.-M. Lachiver, and N. Bataille, “Selecting and scheduling observations of agile satellites,” _Aerospace Science and Technology_ , vol. 6, no. 5, pp. 367–381, 2002.
  * [25] J. Wang, E. Demeulemeester, and D. Qiu, “A pure proactive scheduling algorithm for multiple Earth observation satellites under uncertainties of clouds,” _Computers & Operations Research_, vol. 74, pp. 1–13, 2016.
  * [26] X. Wang, G. Song, R. Leus, and C. Han, “Robust Earth observation satellite scheduling with uncertainty of cloud coverage,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 56, no. 3, pp. 2450–2461, 2020.
  * [27] X.-W. Wang, Z. Chen, and C. Han, “Scheduling for single agile satellite, redundant targets problem using complex networks theory,” _Chaos, Solitons & Fractals_, vol. 83, pp. 125–132, 2016.
  * [28] X. Liu, G. Laporte, Y. Chen, and R. He, “An adaptive large neighborhood search metaheuristic for agile satellite scheduling with time-dependent transition time,” _Computers & Operations Research_, vol. 86, pp. 41–53, 2017.
  * [29] G. Peng, R. Dewil, C. Verbeeck, A. Gunawan, L. Xing, and P. Vansteenwegen, “Agile Earth observation satellite scheduling: An orienteering problem with time-dependent profits and travel times,” _Computers & Operations Research_, vol. 111, pp. 84–98, 2019.
  * [30] X. Wang, C. Han, R. Zhang, and Y. Gu, “Scheduling multiple agile Earth observation satellites for oversubscribed targets using complex networks theory,” _IEEE Access_ , vol. 7, pp. 110 605–110 615, 2019.
  * [31] Y. Du, L. Xing, J. Zhang, Y. Chen, and Y. He, “MOEA based memetic algorithms for multi-objective satellite range scheduling problem,” _Swarm and Evolutionary Computation_ , vol. 50, p. 100576, 2019.
  * [32] R. Vazquez, F. Perea, and J. G. Vioque, “Resolution of an antenna–satellite assignment problem by means of integer linear programming,” _Aerospace Science and Technology_ , vol. 39, pp. 567–574, 2014.
  * [33] H. Wang and S. Bai, “A versatile method for target area coverage analysis with arbitrary satellite attitude maneuver paths,” _Acta Astronautica_ , 2022\.
  * [34] K. Max, “Computer graphics and geometric modelling: Implementation and modelling,” 2005.
  * [35] F. Xhafa and A. W. Ip, “Optimisation problems and resolution methods in satellite scheduling and space-craft operation: a survey,” _Enterprise Information Systems_ , vol. 15, no. 8, pp. 1022–1045, 2021.
  * [36] X. Wang, H. Zhang, S. Fan, and H. Gu, “Coverage control of sensor networks in IoT based on RPSO,” _IEEE Internet of Things Journal_ , vol. 5, no. 5, pp. 3521–3532, 2018.
  * [37] F. Martinez, A. J. Rueda, and F. R. Feito, “A new algorithm for computing boolean operations on polygons,” _Computers & Geosciences_, vol. 35, no. 6, pp. 1177–1185, 2009.
  * [38] L. J. Simonson, “Industrial strength polygon clipping: A novel algorithm with applications in VLSI CAD,” _Computer-Aided Design_ , vol. 42, no. 12, pp. 1189–1196, 2010.
  * [39] B. R. Vatti, “A generic solution to polygon clipping,” _Communications of the ACM_ , vol. 35, no. 7, pp. 56–63, 1992.
  * [40] R. G. Chamberlain and W. H. Duquette, “Some algorithms for polygons on a sphere,” Pasadena, CA: Jet Propulsion Laboratory, Tech. Rep., 2007.
  * [41] C. Han, Y. Gu, G. Wu, and X. Wang, “Simulated annealing based heuristic for multiple agile satellites scheduling under cloud coverage uncertainty,” _arXiv preprint arXiv:2003.08363_ , 2020.
  * [42] X. Wang, C. Han, Y. Yu, and Z. Yang, “A modified adaptive backstepping method for shaft deflection tracking control of magnetically suspended momentum wheel with nonlinear magnetic torque,” _Journal of the Franklin Institute_ , vol. 355, no. 7, pp. 3274–3287, 2018.
  * [43] E. R. Beadle and P. M. Djuric, “A fast-weighted bayesian bootstrap filter for nonlinear model state estimation,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 33, no. 1, pp. 338–343, 1997.
  * [44] J. S. Liu, R. Chen, and T. Logvinenko, “A theoretical framework for sequential importance sampling with resampling,” in _Sequential Monte Carlo methods in practice_. Springer, 2001, pp. 225–246.
  * [45] M. Clerc and J. Kennedy, “The particle swarm-explosion, stability, and convergence in a multidimensional complex space,” _IEEE transactions on Evolutionary Computation_ , vol. 6, no. 1, pp. 58–73, 2002.
  * [46] Natural Resources Satellite Remote Sensing Cloud Service Platform, “High-resolution multi-mode satellite,” http://www.sasclouds.com/english/satellite/chinese/multimode, accessed December 15, 2021.