Partition satellite mission planning method for earliest completion of regional target coverage

Abstract

The present invention discloses a partition satellite mission planning method for earliest completion of regional target coverage, which includes the following steps: 1. partitioning a to-be-observed rectangular region; 2. allocating observation resources to different regions and selecting certain coverage opportunities and coverage modes thereof to minimize completion time corresponding to a formed coverage solution. The present invention can obtain a satisfactory solution capable of finishing complete coverage of regional targets as early as possible with sufficient observation resources through suitable calculation resources and achieve balance between the consumption of calculation resources and optimality of solutions, so that satellite resources can be fully used to conduct the quick and effective coverage search for important regional targets in a real environment.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of priority from Chinese Patent Application No. 201910957601.4, filed on Oct. 10, 2019. The content of the aforementioned applications, including any intervening amendments thereto, is incorporated herein by reference in its entirety.

TECHNICAL FIELD

The present invention belongs to the field of regional target mission planning for earth observation satellites, and particularly relates to a partition satellite mission planning method for earliest completion of regional target coverage.

BACKGROUND OF THE PRESENT INVENTION

Earth observation satellites (EOS) or imaging satellites (abbreviated as satellites hereinafter) have a main function of observing land, oceans and atmosphere through satellite borne sensors (such as visible light cameras and multi-spectral cameras). Observation requirements are proposed by users from various fields and departments and collected by a ground control center of the satellites. The control center comprehensively formulates an imaging coverage plan of each satellite according to the observation requirements in combination with the application conditions of satellite resources, and generates measurement and control instructions to be transmitted to the satellites via a ground measurement and control station. The satellites make corresponding actions after receiving the instructions to image a specified region. The image data is temporarily saved in a satellite borne hard drive. When the satellites can communicate with the ground station, the image data is transmitted to the ground station. In this process, the link that the ground control center formulates the satellite imaging plan is called satellite mission planning, which is one of key links in the entire satellite application management process.

In the traditional satellite application modes, the satellites formulate the plans separately and perform the imaging missions separately. No collaborative observation is carried out among different satellites. As the number of satellites increases, it is possible to use multiple satellites to collaboratively perform the imaging observation for earth. The collaborative regional imaging of multiple satellites is actually needed.

The satellite can only photograph a strip region of limited length and width at one pass. If the to-be-observed region is large, it is difficult for the satellite to completely observe the entire region at one pass. If the user urgently needs the image data of the region, multiple satellites can be used to perform the collaborative imaging at multiple passes. Since tracks of sub-satellite points of each satellite are not necessarily parallel, the imaging strips may be overlapped. If the satellites are not properly arranged, a large number of strips may be overlapped, causing repeated observation of some regions, while some regions are not observed, which is not conducive to the rational utilization of resources. In order to better use the existing satellite coverage resources, a rational plan should be formulated to solve the above contradiction. It is a typical optimization problem of operations research, which needs to be urgently solved. When the to-be-observed region with this problem is large, a conventional regional target planning method may consume a great amount of calculation resources due to the large scale and high complexity of the problem. This is unfavorable for the application of the method in actual engineering environments.

SUMMARY OF THE PRESENT INVENTION

In order to solve the defects in the prior art, the present invention provides a partition satellite mission planning method for earliest completion of regional target coverage so as to obtain a satisfactory solution capable of completing the complete coverage of regional targets as early as possible with sufficient observation resources through suitable calculation resources, and to achieve balance between the consumption of calculation resources and the optimality of solutions, so that the satellite resources can be fully used to conduct quick and effective coverage search for important regional targets in a real environment.

To solve the technical problems, the present invention adopts the following technical solutions:

Compared with the prior art, the present invention has the following beneficial effects:

1. According to a “divide-and-rule” concept, the present invention proposes a solving strategy based on partition, that is, a large to-be-observed region is divided into a plurality of small subregions; coverage resources are allocated to each subregion; each subregion is solved separately; and the solutions of all subregions are combined to form overall coverage solution, thereby reducing the problem scale and alleviating the optimization difficulty. Different allocation solutions of coverage resources may lead to different overall coverage solutions, so that a simulated annealing method is proposed to search an optimal resource allocation solution.

2. A grid discretization method is used to process the to-be-observed subregions, and the original problem is converted to a problem of covering each grid, so that the coupling with the calculation geometry can be reduced. At the same time, based on grids, a longest basic coverage mode set of each observation opportunity is built as an optional coverage mode set, and decision variables are changed from values in a continuous space to values in a discrete space, thereby reducing the calculation amount.

3. The present invention utilizes a two-phase heuristic algorithm to find a desirable coverage solution for each discretized subregion. At a first phase, a dynamic greed rule is used to find a desirable feasible coverage solution capable of being quickly obtained under the current environment. At a second phase, after a coverage opportunity with maximum end time in the selected coverage opportunities is removed on the basis of the existing solutions, the coverage solution is re-planned to improve the quality of the solutions. The two-phase heuristic algorithm is a polynomial time algorithm. When the total number of the coverage modes is in a rational range, the algorithm is high in solving speed, easy to implement, high in interpretability, stable and reliable, so that in the simulated annealing process of the method, a desirable coverage solution can be rapidly and effectively found to increase the search speed for a solution space of allocation solutions.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart of a partition satellite mission planning method for earliest completion of regional target coverage of the present invention;

FIG. 2a is a schematic diagram of left lower corner oblique coverage of a region target of the present invention;

FIG. 2b is a schematic diagram of right lower corner oblique coverage of the region target of the present invention;

FIG. 2c is a schematic diagram of multi-satellite collaborative coverage of the present invention;

FIG. 3a is a schematic diagram of formation of coverage modes of the present invention;

FIG. 3b is a schematic diagram of a longest coverage mode length of the present invention;

FIG. 3c is a schematic diagram of side-sway of satellites of the present invention;

FIG. 3d is a schematic diagram of a maximum roll angle of the present invention;

FIG. 4a is a schematic diagram of a longest basic coverage mode of a right lower corner oblique type of the present invention;

FIG. 4b is a schematic diagram of the longest basic coverage mode of a left lower corner oblique type of the present invention;

FIG. 5a is a schematic diagram of calculation of a roll angle under one condition of the present invention;

FIG. 5b is a schematic diagram of calculation of a roll angle under another condition of the present invention;

FIG. 6a is a schematic diagram of a feasible upper grid condition 1 of the present invention;

FIG. 6b is a schematic diagram of a feasible upper grid condition 2 of the present invention;

FIG. 6c is a schematic diagram of a feasible upper grid condition 3 of the present invention;

FIG. 6d is a schematic diagram of a feasible upper grid condition 4 of the present invention.

DETAILED DESCRIPTION OF THE PRESENT INVENTION

As shown in FIG. 2a, FIG. 2b and FIG. 2c, in the present embodiment, a partition satellite mission planning method for earliest completion of regional target coverage is used to select a subset from sufficient coverage opportunity set S and to plan corresponding coverage modes for completely covering a to-be-observed region R. The optimization purpose of this method is to minimize the completion time. The coverage opportunity refers to an opportunity for a satellite to perform the coverage observation when passing over the to-be-observed region. As shown in FIG. 2c, an area of the to-be-observed region R is large. The single coverage opportunity can only cover a small part, so that a plurality of coverage opportunities are needed for collaborative observation.

As shown in FIG. 3a, FIG. 3b, FIG. 3c and FIG. 3d, each coverage opportunity s∈s has own attribute and includes sub-satellite point track o_s, height h_s of the satellite to the ground and maximum power-on time m_s, corresponding to a maximum length of coverage mode r_s, maximum roll angle v_s and camera field angle w_s. The satellite orbits around the Earth, and a track formed by sub-satellite points on the ground should be a curve actually, but it can be approximated to a straight line in a small range inside the to-be-observed region. In order to simplify the problem, the sub-satellite point track o_s of the coverage opportunity is considered as the straight line in the present embodiment of the present invention. As shown in FIG. 2a and FIG. 2b, each coverage opportunity has a start time α(s) and an end time β(s) for the to-be-observed region R.

The so-called coverage mode c is a rectangular strip photographed by the satellites and determined by a certain roll angle and on/off time selected by one coverage opportunity.

A coordinate system o-xy is established by adopting any vertex of the to-be-observed rectangular region R as an origin o and two sides adjacent to the origin as axis x and axis y. The coordinates of four vertexes of the to-be-observed rectangular region R are marked as left upper corner vertex LU(R), left lower corner vertex LD(R), right upper corner vertex RU(R) and right lower corner vertex RD(R) respectively. The partition satellite mission planning method includes the following steps:

Step 1: a to-be-observed rectangular region R is partitioned.

Step 1.1: length L_p and width W_p of each subregion are set.

Step 1.2: the length L_R and width W_R of the to-be-observed rectangular region R are calculated according to the coordinates of four vertexes of the to-be-observed rectangular region R.

Step 1.3: L_p is divided by L_R to obtain length section number D_L and remainder R_L of the subregion. If the remainder

$R_L\ge \frac{D_L}{2D_L+1}{L}_p,D_L+1$ is assigned to D_L. W_p is divided by W_R to obtain width section number D_W and remainder R_W of the subregion, and if the remainder

$R_W\ge \frac{D_W}{2D_W+1}{W}_p,D_W+1$ is assigned to D_W.

Step 1.4: L_R is divided by D_L to obtain length L_f of the final subregion, and W_R is divided by D_W to obtain width W_f of the final subregion.

Step 1.5: according to the length L_f and width W_f of the final subregion, the to-be-observed rectangular region R is equally partitioned, thereby obtaining the coordinates of four vertexes of each subregion.

Step 1.5.1: let i=1, the left upper corner coordinate LU(R) of the to-be-observed rectangular region R is used as the left upper corner coordinate LU(n_i) of the i^th subregion n_i.

Step 1.5.2: the abscissa value of the left upper corner coordinate LU(n_i) of the i^th subregion n_i is added with L_f to obtain the right upper corner coordinate RU(n_i) thereof. The ordinate value of the left upper corner coordinate LU(n_i) is added with W_f to obtain the left lower corner coordinate LD(n_i) thereof. The ordinate value of the right upper corner coordinate RU(n_i) is added with f to obtain the right lower coordinate RD(n_i).

Step 1.5.3: if the right lower corner coordinate RD(n_i) of the i^th subregion n_i is equal to the right lower corner coordinate RD(R) of the to-be-observed rectangular region R, step 1.5.6 is carried out; otherwise, step 1.5.4 is carried out.

Step 1.5.4: if the abscissa value of the right upper corner coordinate RU(n_i) of the i^th subregion n_i is equal to the abscissa value of the right upper corner coordinate RU(R) of the to-be-observed rectangular region R, LD(n_(i−DL+1)) is used as the left upper corner coordinate LU(n_i+1) of the (i+1)^th subregion n_i+1; otherwise, the right upper corner coordinate RU(n_i) of the i^th subregion n_i is used as the left upper corner coordinate LU(n_i+1) of the (i+1)^th subregion n_i+1.

Step 1.5.5: after i+1 is assigned to i, step 1.5.2 is carried out to construct the coordinates of four vertexes of a next subregion.

Step 1.5.6: the partitioning is ended, and the coordinates of four vertexes of each subregion are outputted.

Step 1.6: a set N of the subregions is obtained, and the total number of the subregions is |N|=D_L×D_W.

Step 2: observation resources are allocated to different regions, and a certain subset of coverage opportunities and corresponding coverage modes are selected to minimize the completion time of a formed coverage solution. That is, different coverage opportunities are responsible for covering different subregions, while different allocation solutions may finally lead to different coverage solutions which correspond to different completion time. Therefore, a simulated annealing method is used to search an optimum allocation solution here.

Step 2.1: the coverage opportunities are randomly allocated to each subregion to form an initial current allocation solution Y={y^s|∀s∈S} wherein the value of y_s represents the subregion where the coverage opportunity s is allocated.

Step 2.2: the coverage opportunity set S*(Y) selected under the current allocation solution Y, a desirable coverage mode set C*(Y) and corresponding completion time F(Y) are calculated.

Step 2.2.1: different subregions and the set of coverage opportunities allocated to different subregions are used as inputs to construct grids, and a two-phase heuristic algorithm is used to calculate the coverage opportunity set selected for different subregions as well as the desirable coverage solution and corresponding completion time.

Step 2.2.1.1: the i^th subregion is divided into a plurality of consistent square grids to obtain a grid set J_i, wherein a side length of each grid should be minimized. The existing grids are ensured to be completely covered by the coverage modes.

If the to-be-observed region was not partitioned into subregions, the whole large region may form a large number of grids A large number of coverage modes may be generated based on these grids, which consumes vast calculation resources. Some large-scale calculation examples cannot even obtain a solution within a rational time, so that embodiments of the present invention are proposed on the basis of solving strategies based on the partition.

Step 2.2.1.2: according to the divided grid set J_i, a set of longest basic coverage modes is generated for each coverage opportunity in set S_i allocated to the i^th subregion n_i, and jointly form a main set C_i of the optional coverage modes.

As shown in FIG. 4a and FIG. 4b, the longest basic coverage mode refers to that the length is equal to the maximum length r_s, and there are two grids that satisfy the following two conditions:

{circle around (1)} Both grids are completely covered by the coverage mode, that is, the four vertexes are within the area covered by the coverage mode;

{circle around (2)} one grid has a vertex located on a left border of the coverage mode and is named left grid, and the other grid has a vertex located on an upper border of the coverage mode and is named upper grid. The left grid and the upper grid can be the same.

Each coverage opportunity can select infinitely many on/off time and roll angles, so each coverage opportunity may actually correspond to an infinite number of coverage modes. However, all coverage modes cannot be enumerated. Therefore, the solving method of the problem is difficult to design on the basis of the universal set of all coverage modes. The longest basic coverage mode is a special coverage mode and can substitute a majority of non-longest basic coverage modes and is limited in quantity. A majority of non-longest basic coverage modes may be slightly adjusted to be the corresponding longest basic coverage modes without reducing the number of completely covered grids.

Step 2.2.1.2.1: if S_i=Ø, step 2.2.1.2.4 is carried out; otherwise, step 2.2.1.2.2 is carried out.

Step 2.2.1.2.2: one of coverage opportunities s∈S_i is selected, and according to the type of the coverage opportunity s, methods in step 2.2.1.2.2.1 to step 2.2.1.2.2.10 are used to construct a longest basic coverage mode set C_s based on the grid set J_i for the coverage opportunity s.

Step 2.2.1.2.2.1: the type of the coverage opportunity s is judged according to the sub-satellite point track o_s of the coverage opportunity s, and the types include a right lower corner oblique type and a left lower corner oblique type.

Step 2.2.1.2.2.2: let set C_s=Ø.

Step 2.2.1.2.2.3: subsets J_s are selected from the grid set J_i, and each grid in the subsets J_s may be used as the left grid of a longest basic coverage mode of the coverage opportunity s.

Step 2.2.1.2.2.3.1: let J_s=Ø.

Step 2.2.1.2.2.3.2: the set J_i is traversed, and the camera roll angle

${\sigma }_s\left(u\right)=\arctan \frac{d_s^u}{h_s}-\frac{w_s}{2}$ of the coverage opportunity s is calculated when each grid u is used as the left grid of the coverage mode. If the coverage opportunity is the right lower corner oblique type, as shown in FIG. 4a, d_s^u is a distance between the left lower corner vertex LD(u) of the grid u and the sub-satellite point track o_s of the coverage opportunity s. If the coverage opportunity is the left lower corner oblique type, as shown in FIG. 4b, d_s^u is a distance between the left upper corner vertex LU(u) of the grid u and the sub-satellite point track o_s of the coverage opportunity s. If |σ_s(u)|≤v_s, the grid u is added into the subset J_s of the left grid.

Step 2.2.1.2.2.4: if J_s=Ø, step 2.2.1.2.2.10 is carried out; otherwise, after one grid p in J_s is selected, step 2.2.1.2.2.5 is carried out.

Step 2.2.1.2.2.5: according to the grid p, the grid subset J_s^p is selected from the grid set J_i, and when p is used as the left grid, each grid in J_s^p may be used as the upper grid to form a longest basic coverage mode of the coverage opportunity s together with p.

Step 2.2.1.2.2.5.1: let J_s^p=Ø.

Step 2.2.1.2.2.5.2: a strip width η_s(p) of the coverage mode when the grid p is used as the left grid is calculated by using the formula (1):

$\begin{array}{cc}{\eta }_s\left(p\right)=\{\begin{array}{c}{d}_s^p-h_s·\tan \left(\arctan \frac{d_s^p}{h_s}-w_s\right),\arctan \frac{d_s^p}{h_s}>w_s\\ d_s^p+h_s·\tan \left(w_s-\arctan \frac{d_s^p}{h_s}\right),\arctan \frac{d_s^p}{h_s}\le w_s\end{array}& \left(1\right)\end{array}$

As shown in FIG. 5a and FIG. 5b, when the grid u is used as the left grid of the coverage mode, a calculation formula of the camera roll angle σ_s(u) of the coverage opportunity s and the strip width η_s(u) of the coverage mode can both be obtained by elementary geometry. OB is a vertical line, and the length of OB is the height of the satellite to the ground, i.e. h_s. The length of AB is a distance between the straight lines L_u and o_s, i.e. d_s^u. Apparently, ∠ABO is a right angle, so it is easy to get:

$\angle \mathrm{AOB}=\arctan \frac{d_s^u}{h_s},$ and the camera roll angle

${\sigma }_s\left(u\right)=\arctan \frac{d_s^u}{h_s}-\frac{w_s}{2}.$

In FIG. 5a, ∠AOC is a field angle, i.e. w_s, while ∠AOB>w_s. The width of the strip is the length of AC at this moment.

$AC=d_s^u-BC=d_s^u-h_s·\tan \angle \mathrm{BOC}=d_s^u-h_s·\tan \left(\angle \mathrm{AOB}-w_s\right)=d_s^u-h_s·\tan \left(\arctan \frac{d_s^u}{h_s}-w_s\right)$

In FIG. 5b, ∠AOD is a field angle, i.e. w_s. OC is an angular bisector of ∠AOD, i.e.

$\angle AOC=\frac{w_s}{2}.$ The width of the strip is length of AD at this moment.

$\angle \mathrm{AOB}=\arctan \frac{d_s^u}{h_s},\angle \mathrm{DOB}=w_s-\angle \mathrm{AOB},$ so it can obtain: BD=h_s·tan ∠DOB, and it can further obtain:

$\mathrm{AD}=d_s^u+\mathrm{BD}=d_s^u+h_s·\tan \angle \mathrm{BOD}=d_s^u+h_s·\tan \left(w_s-\angle \mathrm{AOB}\right)=d_s^u+h_s·\tan \left(w_s-\arctan \frac{d_s^u}{h_s}\right)$

Step 2.2.1.2.2.5.3: the set J_i is traversed. Whether one grid t satisfies the following four conditions is checked. If all conditions are satisfied, the grid t is added into J_s^p; otherwise, the next grid is checked.

If the coverage opportunity is the right lower corner oblique type, the four conditions are:

{circle around (1)} as shown in FIG. 6a, the straight line L_p∥o_s is made through the left lower corner vertex of the grid p, the straight line L_t′∥o_s is made through the right upper corner vertex of t, and the distance between the straight lines L_p and L_t′ cannot exceed the width η_s(p) of the strip.

{circle around (2)} as shown in FIG. 6b, the straight lines L_p′⊥o_s and L_t⊥o_s are made respectively through the left upper corner vertexes of the grids p and t, and L_t is above L_p′.

{circle around (3)} as shown in FIG. 6c, the straight line L_t∥o_s is made through the left lower corner vertex of the grid t, and the straight line L_t is at the right side of L_p.

{circle around (4)} as shown in FIG. 6d, the straight line L_p⊥o_s is made through the right lower corner vertex of the grid p, and the distance between the straight lines L_p and L_t cannot exceed the length l_s of the strip. “∥” indicates a parallel relation, and “⊥” indicates a vertical relation.

If the coverage opportunity is the left lower corner oblique type, the four conditions are:

{circle around (1)} The straight line L_p∥o_s is made through the left upper corner vertex of the grid p, the straight line L_t′∥o_s is made through the right upper corner vertex of t, and the distance between the straight lines L_p and L_t′ cannot exceed the width η_s(p) of the strip.

{circle around (2)} The straight lines L_p′⊥o_s and L_t⊥o_s are made respectively through the right upper corner vertexes of the grids p and t, and L_t is above L_p′.

{circle around (3)} The straight line L_t∥o_s is made through the left upper corner vertex of the grid t, and the straight line L_t is at the right side of L_p.

{circle around (4)} The straight line L_p⊥o_s is made through the left lower corner vertex of the grid p, and the distance between the straight lines L_p and L_t cannot exceed the length l_s of the strip. “∥” indicates a parallel relation, and “⊥” indicates a vertical relation.

Apparently, the grid p also satisfies the four conditions.

Step 2.2.1.2.2.6: if J_s^p=Ø, step 2.2.1.2.2.9 is carried out; otherwise, after one grid q in J_s^p is selected, step 2.2.1.2.2.7 is carried out.

Step 2.2.1.2.2.7: based on the grids P and q, the longest basic coverage c_s(p,q) belonging to the coverage opportunity s is constructed.

Step 2.2.1.2.2.7.1: if the coverage opportunity is the right lower corner oblique type, the straight line L_p∥o_s is made through the left lower corner vertex LD(p) of the grid p, and the straight line L_q⊥o_s is made through the left upper corner vertex LU(q) of the grid q. If the coverage opportunity is the left lower corner oblique type, the straight line L_p∥o_s is made through the left upper corner vertex LU(p) of the grid p, and the straight line L_q ⊥o_s is made through the right upper corner vertex RU(q) of the grid q.

Step 2.2.1.2.2.7.2: the strip width η_s(p) of the coverage mode when the grid p is used as the left grid is calculated.

Step 2.2.1.2.2.7.3: the straight line L_p is made so as to satisfy:

{circle around (1)} L_p∥L_p;

{circle around (2)} (the distance between L_p and L_p is equal to η_s(p);

{circle around (3)} L_p is at the right side of L_p.

Step 2.2.1.2.2.7.4: the straight line L_q is made so as to satisfy:

{circle around (1)} L_q∥L_q;

{circle around (2)} (the distance between L_q and L_q is equal to l_s;

{circle around (3)} L_q is at the lower side of L_q.

Step 2.2.1.2.2.7.5: four straight lines of L_p, L_p, L_q, L_q are used as four sides of the rectangle, and intersections of the four straight lines are used as four vertexes of the rectangle to construct the coverage mode c_s(p,q).

Step 2.2.1.2.2.8: c_s (p,q) is added into C_s, the grid q is removed from J_s^p, and then step 2.2.1.2.2.6 is carried out.

Step 2.2.1.2.2.9: the grid p is removed from J_s, and then step 2.2.1.2.2.4 is carried out.

Step 2.2.1.2.2.10: the coverage mode set C_s is outputted.

Step 2.2.1.2.3: the coverage mode opportunity s is removed from S_i, and then step 2.2.1.2.1 is carried out.

Step 2.2.1.2.4: let, the coverage mode set C_i is outputted.

Step 2.2.1.3: by utilizing the grid set J_i separated from i_th subregion n_i as well as the coverage opportunity set S_i={s|y_s=n_i, s∈S} allocated to i^th subregion n_i and optional coverage mode set C_i as inputs, the coverage opportunity subset S_i* selected by the i^th subregion n_i is determined by the two-phase heuristic algorithm, and the coverage mode of each coverage opportunity is selected, so that the desirable coverage mode set C_i* is obtained, and the completion time F_i thereof can be obtained at the same time.

Step 2.2.1.3.1: by utilizing the grid set J_i, the coverage opportunity set S_i and the corresponding coverage mode set C_i as inputs, the selected coverage opportunity subsets S_i*, the selected coverage mode sets C_i* and the number A_i of corresponding coverable grids can be calculated with the heuristic algorithm.

Step 2.2.1.3.1.1: all coverage opportunities in S_i are sorted in a non-descending order according to the end time to obtain

$\stackrel{\rightharpoonup }{S_i}=\{s_1,s_2,\dots ,s_{|\stackrel{\rightharpoonup }{s_i|}}\}.$

Step 2.2.1.3.1.2: the state of all grids in J_i is set as “uncovered”, let g=0, and C_i*=Ø.

Step 2.2.1.3.1.3: g+1 is assigned to g, the coverage opportunity s_g is extracted from {right arrow over (S_i)}, and the optional coverage mode set C_sg is traversed. For each coverage mode C_sg, the grid capable of being completely covered by C_sg and with the state of “uncovered” is selected from the grid set J_i and marked as valid grid sets V(c_sg), while the number is marked as |V(c_sg)|.

Step 2.2.1.3.1.4: the coverage mode

$c_{s_g}^{*}=\arg \underset{\forall c_{s_g}\in C_{s_g}}{\max }\left|V\left(c_{s_g}\right)\right|$ covering the most valid grids is selected from C_sg and used as the coverage mode selected by the coverage opportunity s_g. After c_sg* is added into C_i*, the state of all grids in the valid grid sets V(c_sg) of the coverage mode c_sg* is changed as “covered”.

Step 2.2.1.3.1.5: the number of all grids with the state of “covered” in J_i is calculated and marked as A_i.

Step 2.2.1.3.1.6: if A_i=|J_i| or g=|{right arrow over (S_i)}|, step 2.2.1.3.1.7 is carried out; otherwise, step 2.2.1.3.1.3 is carried out.

Step 2.2.1.3.1.7: the coverage opportunity subset S_i*={s₁, s₂, . . . , s_g} as well as the coverage mode set C_i* thereof and the number A_i of covered grids are outputted.

Step 2.2.1.3.2: if A_i<|J_i|, step 2.2.1.3.4 is carried out; otherwise, step 2.2.1.3.3 is carried out.

On the premise of sufficient coverage opportunities, the coverage opportunities that are randomly allocated to each subregion generally can completely cover the subregion, that is, A_i=|J_i|. If some subregions cannot be completely covered, they may be improved in neighbor allocation solutions constructed by an insert operator and an exchange operator at the simulated annealing phase.

Step 2.2.1.3.3: S_i* and C_i* are optimized by using a promotion algorithm.

Step 2.2.1.3.3.1: one coverage opportunity

$s_m=\arg \underset{\forall s\in S_i^{*}}{\max }\beta \left(s\right)$ with maximum end time in S_i* is removed from S_i* to obtain the coverage opportunity set S_i.

Step 2.2.1.3.3.2: by utilizing the grid set J_i, the coverage opportunity set S_i and the optional coverage mode set C_i={c|∀c∀C_s, ∀s∈S_i} corresponding to S_i as inputs, the selected coverage opportunity subsets S_i*, the selected coverage mode set C_i*, and the number A_i of the corresponding coverable grids can be calculated with the heuristic algorithm.

Step 2.2.1.3.3.2.1: the state of all grids in J_i is set as “uncovered”, and C_i*=Ø.

Step 2.2.1.3.3.2.2: let g=1.

Step 2.2.1.3.3.2.3: if g>|C_i|, step a.3.3.2.5 is carried out; otherwise, step a.3.3.2.4 is carried out.

Step 2.2.1.3.3.2.4: the g^th coverage mode in C_i is selected, and the grids capable of being completely covered by the g^th coverage mode and with the state of “uncovered” are selected from the grid set J_i and marked as valid grid set V(g), while its number is marked as |V(g)|. g+1 is assigned to g. Thereafter step a.3.3.2.3 is carried out.

Step 2.2.1.3.3.2.5: if C_i=Ø, step a.3.3.2.8 is carried out; otherwise, step a.3.3.2.6 is carried out.

Step 2.2.1.3.3.2.6: the coverage mode

${\stackrel{\_}{g}}^{*}=\arg \underset{\forall \stackrel{\_}{g}\in \stackrel{\_}{C_i}}{\max }\left|V\left(\stackrel{\_}{g}\right)\right|$ with most valid grids in C_i is selected and used as the coverage mode selected by the corresponding coverage opportunity s*. g* is added into C_i*, and the state of all grids in the valid grid set V(g*) of the coverage mode g* is changed to be “covered”.

Step 2.2.1.3.3.2.7: the coverage mode set C_s* corresponding to the coverage opportunity s* is removed from the optional coverage mode set C_i, and then step a.3.3.2.2 is carried out.

Step 2.2.1.3.3.2.8: the coverage mode set C_i*, selected by each coverage opportunity is outputted, i.e. the final coverage solution. The number of all “covered” grids in J_i is calculated and marked as A_i.

Step 2.2.1.3.3.3: if A_i=|J_i|, S_i* is assigned to S_i* , and C_i* is assigned to C_i*. Thereafter, step 2.2.1.3.3.1 is carried out; otherwise, step 2.2.1.3.4 is carried out.

Step 2.2.1.3.4: the selected coverage opportunity set S_i* and the coverage mode set C_i* are outputted. The maximum value of the end time of the coverage opportunity in S_i* is used as the completion time F_i.

The first phase of the two-phase heuristic algorithm in step 2.2.1.3 utilizes the rule based on dynamic greed. The coverage mode covering most valid grids is selected at each time. Once the coverage opportunity selects one coverage mode, the coverage opportunity may exit the subsequent selection activities. Once the grid is covered in the previous selection, the grid exits the subsequent selection activities. The method is simple, fast and effective and very suitable for rapidly obtaining the desirable feasible coverage solution under the current environment. At the second phase, on the basis of the solution obtained at the first phase, after the selected coverage opportunity with maximum end time is deleted, the complete coverage solution is re-planned to search a better solution space of the solution.

Step 2.2.2: the coverage opportunity sets {S_i*|i=1, 2, . . . , |N|} selected by all subregions are combined to obtain an overall selected coverage opportunity set S*. The coverage solutions {C_i*|i=1, 2, . . . , |N|} of all subregions are combined to obtain an overall coverage solution C* of the to-be-observed rectangular region R. The maximum value of the completion time {F_i|i=1, 2, . . . , |N|} corresponding to the coverage solutions of all subregions is used as the overall completion time F of the to-be-observed rectangular region R.

The solution space of the allocation solution is searched below by utilizing the simulated annealing method to obtain a better allocation solution, so that the completion time corresponding to the formed coverage solution is minimized, which specifically includes step 2.3-step 2.9:

Step 2.3: maximum iteration times k and initial annealing temperature T₁ are set, and the iteration times of the current outer cycle is set as k=1.

Step 2.4: if k>k, step 2.9 is carried out; otherwise, step 2.5 is carried out.

Step 2.5: the insert operator and the exchange operator are used to construct a plurality of neighbor allocation solutions of the current allocation solution Y.

Step 2.5.1: the number of the to-be-constructed neighbor allocation solutions is set as m, and the number of the constructed neighbor allocation solutions is set as m, and m=0.

Step 2.5.2: the insert operator is used to construct the neighbor allocation solution.

Step 2.5.2.1: an idle coverage opportunity ŝ is found from the unselected coverage opportunity set S′(Y)={S/S*(Y)} under the current allocation solution Y. The end time β(ŝ) of the idle coverage opportunity ŝ is smaller than the overall completion time F(Y) but greater than the completion time F_yŝ(Y) that the idle coverage opportunity ŝ is allocated to the subregion y_ŝ. If the idle coverage opportunity cannot be found, step 2.5.2.4 is carried out; otherwise, step 2.5.2.2 is carried out.

Step 2.5.2.2: m+1 is assigned to m, and the idle coverage opportunity ŝ in the current allocation solutions Y is randomly allocated to one subregion, and the end time β(ŝ) of the idle coverage opportunity ŝ is smaller than the completion time of the allocated subregion, thereby forming a neighbor allocation solution Y_m.

Step 2.5.2.3: the idle coverage opportunity ŝ is removed from the coverage opportunity set S′(Y), and then step 2.5.2.1 is carried out.

Step 2.5.2.4: all neighbor allocation solutions obtained by utilizing the insert operator are outputted.

Step 2.5.3: the exchange operator is used to construct the neighbor allocation solution.

Step 2.5.3.1: if m≥m, step 2.5.3.3 is carried out; otherwise, step 2.5.3.2 is carried out.

Step 2.5.3.2: m+1 is assigned to m. Two subregions n_i and n_j are randomly selected from the subregion set N. One coverage opportunity s_i is randomly selected from the coverage opportunity set {s|y_s=n_i} allocated to the subregion n_i. One coverage opportunity s_j is randomly selected from the coverage opportunity set {s|y_s=n_j} allocated to the subregion n_j. s_i and s_j are exchanged in the current allocation solution Y, that is, n_j is assigned to y_si, and n_i is assigned to y_sj, thereby forming a neighbor coverage solution Y_m, and then step 2.5.3.1 is carried out.

Step 2.5.3.3: all neighbor allocation solutions obtained by using the exchange operator are outputted.

Step 2.6: the coverage opportunity set selected under each neighbor allocation solution, desirable coverage mode set and corresponding completion time are calculated separately, and the neighbor allocation solution with minimum completion time is selected as the optimal neighbor allocation solution Y*. The specific steps refer to step 2.2 and are not repeated here.

Step 2.7: The completion time of the optimal neighbor allocation solution Y* is compared with the completion time of the current allocation solution Y. If F(Y*)≤F(Y) the optimal neighbor allocation solution Y* is assigned to the current allocation solution Y, and then step 2.8 is directly carried out; otherwise, step 2.8 is carried out after the step 2.7.1-step 2.7.2.

Step 2.7.1: an inferior solution accept rate

$P_k=\exp \left\{\frac{F\left(Y\right)-F\left({\stackrel{\_}{Y}}^{*}\right)}{T_k}\right\}$ at k^th iteration is calculated.

Step 2.7.2: a random number r_k of 0-1 is generated at k^th times. If r_k≤p_k, the optimal neighbor allocation solution Y* is assigned to the current allocation solution Y.

Step 2.8: k+1 is assigned to k, and the annealing temperature T_k=T_k-1×λ is updated. λ is a cooling coefficient and is a fixed value less than 1 and greater than 0, and then step 2.4 is carried out.

The annealing temperature represents randomness for searching the solution space. If the requirement for the final solution quality is high, λ may be set to be large, and the iteration times may be increased. If the rapid convergence is required, λ is set to be small.

Step 2.9: the current allocation solution Y as well as the overall selected coverage opportunity set S*(Y) under the allocation solution, better coverage modes C*(Y) of each coverage opportunity and corresponding completion time F(Y) are outputted.

Claims

1. A partition satellite mission planning method for earliest completion of regional target coverage, wherein the partition satellite mission planning method comprises: step 1: partitioning a to-be-observed rectangular region R;step 2: allocating observation resources to different regions, and selecting a coverage opportunity and a coverage mode thereof to minimize a completion time of a formed coverage solution, comprising:step 2.1: randomly allocating the coverage opportunity to each subregion to form an initial current allocation solution Y={ys|∇s∈S}, wherein the value of ys represents the subregion where the coverage opportunity s is allocated, and S represents a coverage opportunity set;step 2.2: calculating the coverage opportunity set S*(Y) selected under the current allocation solution Y, a better coverage mode set C*(Y) and corresponding completion time F(Y);step 2.3: setting maximum iteration times k and initial annealing temperature T1, and setting iteration times of a current outer cycle as k=1;step 2.4: if k>k, carrying out step 2.9; otherwise, carrying out step 2.5;step 2.5: using an insert operator and an exchange operator to construct a plurality of neighbor allocation solutions of the current allocation solution Y;step 2.6: calculating the coverage opportunity set selected under each neighbor allocation solution, better coverage mode set and corresponding completion time separately, and selecting the neighbor allocation solution with minimum completion time as the optimal neighbor allocation solution Y*;step 2.7: comparing the completion time of the optimal neighbor allocation solution Y* with the completion time of the current allocation solution Y; if F(Y*)≤F(Y), assigning the optimal neighbor allocation solution Y* to the current allocation solution Y, and then carrying out step 2.8 directly; otherwise, carrying out step 2.8 after the step 2.7.1-step 2.7.2;step 2.7.1: calculating an inferior solution accept rate

2. The partition satellite mission planning method according to claim 1, wherein the step 1 is carried out as follows: step 1.1: setting length Lp and width Wp of each subregion;step 1.2: calculating the length LR and width WR of the to-be-observed rectangular region R according to the coordinates of four vertexes of the to-be-observed rectangular region R;step 1.3: dividing Lp by LR to obtain length section number DL and remainder RL of the subregion; if the remainder

3. The partition satellite mission planning method according to claim 1, wherein the step 2.5 is carried out as follows: step 2.5.1: setting the number of the to-be-constructed neighbor allocation solutions as m, and setting the number of the constructed neighbor allocation solutions as m, and m=0;step 2.5.2: using the insert operator to construct the neighbor allocation solution, comprising step 2.5.2.1-step 2.5.2.4:step 2.5.2.1: finding an idle coverage opportunity ŝ from an unselected coverage opportunity set S′(Y)={S/S*(Y)} under the current allocation solution Y, wherein the end time β(ŝ) of the idle coverage opportunity ŝ is smaller than the overall completion time F(Y) but greater than the completion time Fyŝ(Y) that the idle coverage opportunity ŝ is allocated to the subregion yŝ; if the idle coverage opportunity is not found, carrying out step 2.5.2.4; otherwise, carrying out step 2.5.2.2;step 2.5.2.2: assigning m+1 to m, and randomly allocating the idle coverage opportunity ŝ in the current allocation solutions Y to one subregion, wherein the end time β(ŝ) of the idle coverage opportunity ŝ is smaller than the completion time of the allocated subregion, thereby forming a neighbor allocation solution Ym;step 2.5.2.3: removing the idle coverage opportunity s from the coverage opportunity set S′(Y), and then carrying out step 2.5.2.1;step 2.5.2.4: outputting all neighbor allocation solutions obtained by utilizing the insert operator;step 2.5.3: using the exchange operator to construct the neighbor allocation solution, comprising step 2.5.3.1-step 2.5.3.3:step 2.5.3.1: if m≥m, carrying out step 2.5.3.3; otherwise, carrying out step 2.5.3.2;step 2.5.3.2: assigning m+1 to m; randomly selecting two subregions ni and nj from the subregion set N; randomly selecting one coverage opportunity si from the coverage opportunity set {s|ys=ni} allocated to the subregion ni; randomly selecting one coverage opportunity sj from the coverage opportunity set {s|ys=nj} allocated to the subregion nj; exchanging si and sj in the current allocation solution Y, that is, assigning nj to ysi, and assigning ni to ysj, thereby forming a neighbor coverage solution Ym, and then carrying out step 2.5.3.1;step 2.5.3.3: outputting all neighbor allocation solutions obtained by using the exchange operator.

4. The partition satellite mission planning method according to claim 1, wherein the step 2.2 and the step 2.6 of calculating the selected coverage opportunity set under the given allocation solution, a better coverage mode set selected for each coverage opportunity and corresponding completion time are carried out as follows: step a: using different subregions and the set of coverage opportunities allocated to different subregions as inputs to construct grids, and using a two-phase heuristic algorithm to calculate the coverage opportunity set selected for local part of different subregions as well as the better coverage solution and corresponding completion time;step a.1: dividing the ith subregion ni into a plurality of consistent square grids to obtain a grid set Ji;step a.2: according to the divided grid set Ji, generating a longest basic coverage mode for the coverage opportunity set Si allocated to the ith subregion ni, and forming an optional coverage mode set Ci, whereinthe longest basic coverage mode refers to that the length is equal to the maximum length rs, and the following two types of grids exist: (1) one type of grids is completely covered by the coverage mode, that is, the four vertexes are within the area covered by the coverage mode;(2) one grid has a vertex located on a left border of the coverage mode and is named left grid, and the other grid has a vertex located on an upper border of the coverage mode and is named upper grid, wherein the left grid and the upper grid is overlapped;step a.3: by utilizing the grid set Ji separated from ith subregion ni as well as the coverage opportunity set Si={s|ys=ni, s∈S} allocated to ith subregion ni and the optional coverage mode set Ci as inputs, determining the coverage opportunity subset Si* selected by the ith subregion ni by the two-phase heuristic algorithm, and selecting the coverage mode of each coverage opportunity, in order to obtain the better coverage mode set Ci* and the completion time Fi thereof;step b: combining the coverage opportunity sets {Si*|i=1, 2, . . . , |N|} selected by all subregions to obtain an overall selected coverage opportunity set S*; combining the coverage solutions {Ci*|i=1, 2, . . . , |N|} of all subregions to obtain an overall coverage solution C* of the to-be-observed rectangular region R; and using the maximum value of the completion time {Fi*|i=1, 2, . . . , |N|} corresponding to the coverage solutions of all subregions as the overall completion time F of the to-be-observed rectangular region R.

5. The partition satellite mission planning method according to claim 4, wherein the step a.2 of generating a longest basic coverage mode for the coverage opportunity set according to the grid set and forming the optional coverage mode set is carried out as follows: step a.2.1: if Si=Ø, carrying out step a.2.4; otherwise, carrying out step a.2.2;step a.2.2: selecting one of coverage opportunities s∈Si and according to the type of the coverage opportunity s, using methods in step a.2.2.1 to step a.2.2.10 to construct a longest basic coverage mode set Cs based on the grid set Ji for the coverage opportunity s;step a.2.2.1: judging the type of the coverage opportunity s according to a sub-satellite point track os of the coverage opportunity s, and the types including a right lower corner oblique type and a left lower corner oblique type;step a.2.2.2: making the longest basic coverage mode set as Cs=Ø;step a.2.2.3: selecting subsets Js from the grid set Ji, and using each grid in the subsets Js as the left grid of the longest basic coverage mode of the coverage opportunity s;step a.2.2.3.1: setting Js=Ø;step a.2.2.3.2: traversing the set Ji, and calculating a camera roll angle

6. The partition satellite mission planning method according to claim 5, wherein the step a.2.2.5 is carried out as follows: step a.2.2.5.1: setting Jsp=Ø;step a.2.2.5.2: calculating a strip width ns(p) of the coverage mode when the grid p is used as the left grid by using the formula (1):

7. The partition satellite mission planning method according to claim 5, wherein the step a.2.2.7 is carried out as follows: step a.2.2.7.1: if the coverage opportunity is the right lower corner oblique type, making the straight line Lp∥os through the left lower corner vertex LD(p) of the grid p, and making the straight line Lq⊥os through the left upper corner vertex LU(q) of the grid q; if the coverage opportunity is the left lower corner oblique type, making the straight line Lp∥os is made through the left upper corner vertex LU(p) of the grid p, and making the straight line Lq⊥os through the right upper corner vertex RU(q) of the grid q;step a.2.2.7.2: calculating the strip width ηs(p) of the coverage mode when the grid p is used as the left grid;step a.2.2.7.3: making the straight line Lp in order to satisfy: (1) Lp∥Lp;(2) the distance between Lp and Lp is equal to ηs(p);(3) Lp is at the right side of Lp;step a.2.2.7.4: making the straight line Lq in order to satisfy: (1) Lq∥Lq;(2) the distance between Lq and Lq is equal to ls;(3) Lq is at the lower side of Lq;step a.2.2.7.5: using four straight lines of Lp, Lp, Lq, Lq as four sides of the rectangle, and using intersections of the four straight lines as four vertexes of the rectangle to construct the coverage mode cs(p,q).

8. The partition satellite mission planning method according to claim 4, wherein the step a.3 is carried out as follows: step a.3.1: by utilizing the grid set Ji, the coverage opportunity set Si and the corresponding coverage mode set Ci as inputs, calculating the selected coverage opportunity subsets Si*, the selected coverage mode sets Ci* and the number Ai of corresponding coverable grids with the heuristic algorithm;step a.3.2: if Ai<|Ji|, with |Ji| indicating the number of elements in the set carrying out step a.3.4; otherwise, carrying out step a.3.3;step a.3.3: optimizing Si* and Ci* by using a promotion algorithm;step a.3.3.1: removing one coverage opportunity

9. The partition satellite mission planning method according to claim 8, wherein the step a.3.1 is carried out as follows: step a.3.1.1: sorting all coverage opportunities in Si in a non-descending order according to the end time to obtain

10. The partition satellite mission planning method according to claim 8, wherein the step a.3.3.2 is carried out as follows: step a.3.3.2.1: setting the state of all grids in Ji as “uncovered”, and Ci*=Ø;step a.3.3.2.2: g=1, wherein g indicates the subscript of the gth coverage mode in the coverage mode set Ci*;step a.3.3.2.3: if g>|Ci|, carrying out step a.3.3.2.5; otherwise, carrying out step a.3.3.2.4;step a.3.3.2.4: selecting the gth coverage mode in Ci, and selecting the grids capable of being completely covered by the gth coverage mode and with the state of “uncovered” from the grid set Ji and marking as valid grid set V(g), and marking the number as |V(g)|; assigning g+1 to g; and carrying out step a.3.3.2.3;step a.3.3.2.5: if Ci=Ø, carrying out step a.3.3.2.8; otherwise, carrying out step a.3.3.2.6;step a.3.3.2.6: selecting the coverage mode