Regional coverage satellite constellation design method based on optimal coverage satisfaction

By calculating the average coverage satisfaction and minimum instantaneous satisfaction of satellite constellations, establishing an optimization model and using the BILP method to solve the problem that traditional design methods are difficult to find global optimal solutions and multi-objective optimization subjectivity, achieving more universal Sexual coverage satisfaction optimization.

CN119210553APending Publication Date: 2024-12-27SOUTHEAST UNIV
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202411149182.9
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2024-08-21
Publication Date
2024-12-27

AI Technical Summary

Technical Problem

Traditional satellite constellation design methods are difficult to find the global optimal solution, and the results of multi-objective optimization problems are subjective, making it difficult to effectively meet the coverage needs of ground stations.

Method used

By calculating the average coverage satisfaction of the constellation, an optimization model based on the minimum instantaneous satisfaction was established, and it was converted into a standard 0-1 integer optimization model. The solution was used to obtain the optimal constellation distribution vector.

Benefits of technology

A more universal satellite constellation design is achieved, single-objective optimization with coverage satisfaction is improved, subjectivity of multi-objective optimization is reduced, and design solutions with stronger physical significance are obtained.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure SMS_1
    Figure SMS_1
  • Figure SMS_2
    Figure SMS_2
  • Figure SMS_5
    Figure SMS_5
Patent Text Reader

Abstract

The invention discloses a regional coverage satellite constellation design method based on optimal coverage satisfaction, and the method comprises the steps: building an average coverage layer number and average coverage satisfaction index model of a system based on an access function, a coverage demand function and a coverage function of a ground station and a constellation distribution vector of a satellite; aiming at the simultaneous optimization of the two indexes, establishing a region coverage constellation design model which takes the maximization of the average coverage satisfaction as a target and takes the condition that the average coverage layer number reaches a set value as a constraint; in the optimization solution stage, an original model is decoupled by using a BIG-M method, an original constellation design model is converted into a standard 0-1 integer optimization model, and then solution is performed by using a BILP method, and the optimal solution is the optimal satellite constellation. Compared with a traditional design scheme, the method has the advantages that indexes with higher physical significance are adopted for single-target optimization modeling, a multi-target optimization algorithm does not need to be adopted, and the model has higher universality compared with the traditional scheme.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The present invention belongs to the technical field of satellite communications, and in particular relates to a method for designing a regional coverage satellite constellation based on optimal coverage satisfaction. Background Art

[0002] At present, traditional satellite constellation design methods are divided into analytical methods and optimization methods. Analytical methods often combine the geometric structure of the constellation with the coverage simulated by STK software for modeling; however, since the goal of constellation design is not unique and is mostly nonlinear and discontinuous, traditional analytical methods are not conducive to finding the global optimal solution. For single-objective optimization problems, traditional optimization methods often use genetic algorithms to design constellations for specific system goals, but their convergence speed is slow and it is difficult to converge to the global optimal solution; for multi-objective optimization problems, previous work used non-dominated sorting genetic algorithms to solve multi-objective optimization problems formulated under actual constraints, but the optimization results are only Pareto optimal sets, and the optimization results still need to be manually selected as design solutions, which is subjective. Summary of the invention

[0003] Purpose of the invention: In order to solve the problems existing in the above-mentioned prior art, the present invention discloses a method for designing a regional coverage satellite constellation based on optimal coverage satisfaction.

[0004] Technical solution: The present invention provides a method for designing a regional coverage satellite constellation based on optimal coverage satisfaction, comprising the following steps:

[0005] Step 1: Calculate the average coverage satisfaction of the common ground track constellation χ(x), where x is the constellation distribution vector, x = (x[0], x[1], …, x[l], …, x[L-1]) T , L represents the number of all discrete moments, l represents the lth discrete moment, x[l] represents the time-shifted position of the satellite relative to the seed satellite at an interval of l discrete moments, and the expression of x[l] is as follows:

[0006]

[0007] l k represents the delay index of the kth satellite;

[0008] Step 2: Based on the average coverage satisfaction of the constellation and the minimum instantaneous satisfaction of the area of ​​interest as a constraint, establish an optimization model M1 based on the maximum average coverage satisfaction;

[0009] Step 3: Convert the model M1 into the decoupled model M2, and then convert the M2 model into the standard 0-1 integer optimization model M3;

[0010] Step 4: Solve M3 to obtain the optimal constellation distribution vector.

[0011] Furthermore, the common ground track constellation in step 1 satisfies the following conditions:

[0012] 1. The orbital period of each satellite is a rational multiple of the Greenwich node period;

[0013] 2. The semi-major axis a, eccentricity e, inclination i and perigee argument ω0 of all satellite orbits at the initial moment are the same;

[0014] 3. Right ascension of the ascending node of the kth satellite Ω k Peace anomaly M k satisfy:

[0015] N P Ω k +N D M k = constant mod(2π),

[0016] Among them, N D is the number of satellite Greenwich node cycles, N P is the number of orbital periods, constant represents a constant, and mod represents the remainder.

[0017] Furthermore, the expression of the constellation average coverage satisfaction χ(x) in step 1 is as follows:

[0018]

[0019] Where L′ represents the total number of discrete moments when all ground stations have communication needs, represents the set of ground stations in the region of interest, represents the total number of ground stations in the target area, U(·) represents the unit step function, ξ j [l′] represents the number of satellite coverage layers of the jth ground station at the discrete moment when the l′th ground station has a communication demand, v 0,j [l′] represents the access function between the seed satellite and the jth ground station at the discrete moment when the l′th ground station has a communication demand. The expression of the access function between the kth satellite and the jth ground station at the lth discrete moment is shown below:

[0020]

[0021] Among them, when k = 0, it means the seed satellite, θ k,j [l] is the elevation angle between the kth satellite and the jth ground station at the lth discrete moment, θ k,j,min [l] is the minimum elevation angle between satellite k and ground station j at discrete time l; Represents the circular convolution operator.

[0022] Furthermore, the expression of the optimization model M1 in step 2 is as follows:

[0023]

[0024] Among them, K represents the intermediate variable, K max represents the maximum number of satellites, σ0[l′] represents the minimum instantaneous satisfaction of the area of ​​interest, σ[l′] represents the instantaneous satisfaction of the area of ​​interest, ζ(x) represents the average number of coverage layers of the ground station, ζ0 represents the minimum average number of coverage layers of the satellite constellation required by the user, and the expressions of ζ(x) and σ[l′] are as follows:

[0025]

[0026]

[0027] Among them, d j [l'] is the instantaneous satisfaction d when l = l' j The value of [l], d j [l]Satisfy the following relationship:

[0028]

[0029] Among them, ξ j [l] is the satellite coverage required by the jth ground station at the lth discrete moment, b j [l] is the coverage function imposed by the constellation on the jth ground station at the lth discrete moment, b j The expression of [l] is:

[0030]

[0031] Wherein, K represents the total number of satellites in the satellite constellation distribution.

[0032] Furthermore, the covering functions of L discrete moments form a covering function vector b j , b j The expression is: b j =(b j [0],b j [1],…,b j [l],…,b j [L-1]) T , b j And the constellation distribution vector x satisfies the following relationship:

[0033] V 0,j x=b j

[0034] Among them, V 0,j represents the access circulant matrix, V 0,j The expression is:

[0035]

[0036] Then b j [l] and x[l] satisfy the following step coupling relationship:

[0037]

[0038] Furthermore, the expression of the model M3 is as follows:

[0039]

[0040] Among them, Q is the augmented vector, Q T =[D T ,x T ], D is the instantaneous coverage satisfaction vector, 0 1×L represents a 1×L-dimensional all-zero matrix, express dimensional all-1 matrix; the expressions of A1, A2, A3, A4, A5, B1, B2, B3, B4, and B5 are as follows:

[0041]

[0042]

[0043]

[0044]

[0045]

[0046] Where I represents the identity matrix, the subscript of I represents the dimension of the identity matrix, and σ0 is the minimum coverage instantaneous satisfaction vector of the region of interest, satisfying σ0=(σ0[0],…,σ0[l′],…,σ0[L′-1]) T , is the new cyclic access matrix after removing the idle time access function from the cyclic access matrix, V j(L′-1)+l′+1 is the access function matrix of the jth ground station at the l′th time, V j(L′-1)+l′+1 is the cyclic access matrix V 0,j The j(L′-1)+l′th row of ξ j =(ξ j[0],…,ξ j [l′],…,ξ j [L′-1]) T , M and ε are both constants.

[0047] Furthermore, the expression of the model M2 is as follows:

[0048]

[0049] Among them, M and ε are constants, Q is the augmented vector, Q T =[D T ,x T ], D is the instantaneous coverage satisfaction vector,

[0050] Furthermore, the step 4 is solved using the BILP method.

[0051] Beneficial effects: The present invention establishes the average coverage layer number and average coverage satisfaction index model of the system based on the access function, coverage demand function, coverage function and satellite constellation distribution vector of the ground station; for the simultaneous optimization of these two indicators, a regional coverage constellation design model is established with the goal of maximizing the average coverage satisfaction and the constraint of the average coverage layer number reaching the set value; in the optimization and solution stage, the BIG-M method is used to decouple the original model, and the original constellation design model is converted into a standard 0-1 integer optimization model, and then the BILP method is used for solution, and the optimal solution is the optimal satellite constellation. Compared with the traditional design scheme, the method of the present invention uses indicators with stronger physical significance for single-objective optimization modeling, without the need to use a multi-objective optimization algorithm, and the model is more universal than the traditional scheme. DETAILED DESCRIPTION

[0052] The exemplary embodiments of the present invention and their description are used to explain the present invention and do not constitute improper limitations of the present invention.

[0053] Consider a constellation mode in which the temporal and spatial positions of all satellites are generated by system calculation so that their sub-satellite point tracks overlap, that is, a single common ground track is used; this type of constellation is called a common ground track constellation. The satellites in the common ground track constellation meet the following three conditions:

[0054] 1) The orbital period of each satellite is a rational multiple of the Greenwich node period;

[0055] 2) The semi-major axis a, eccentricity e, inclination i and perigee argument ω0 of all satellite orbits at t = t0 are the same;

[0056] 3) Right ascension of the ascending node of the kth satellite Ω k Peace anomaly Mk satisfy:

[0057] N P Ω k +N D M k = constant mod(2π),

[0058] Among them, N D is the number of satellite Greenwich node cycles, N P is the number of orbital periods, constant represents a constant, and mod represents the remainder.

[0059] When designing a satellite constellation for regional coverage, the constellation must have spatial and temporal references relative to ground stations and epochs. The virtual satellites in the constellation that convey reference orbital information are defined as seed satellites, and the actual satellites inherit the common orbital characteristics in the orbital elements of the seed satellites. It should be noted that there is no need to have an actual satellite at the seed satellite position; the orbital elements of the seed satellite serve as a reference for the orbital elements of the actual satellites in the system.

[0060] Step 1: Build a system model

[0061] Consider the access function vector between the kth satellite and the jth ground station It stores Boolean information about the satellite access (or visibility) state at each discrete time l∈{0,…,L-1}, where k=0 represents a seed satellite. Therefore, the access function vector v k,j v k,j =(v k,j [0],v k,j [1],…,v k,j [L - 1]) T ,in,(·) T is the transpose of the vector, v k,j [l] is the access function between the kth satellite and the jth ground station at the lth discrete time, satisfying:

[0062]

[0063] Among them, θ k,j [l] is the elevation angle between the kth satellite and the jth ground station at the lth discrete moment; θ k,j,min [l] is the minimum elevation angle threshold between satellite k and ground station j at discrete time l.

[0064] Since there are multiple satellites in a constellation system, the access function vectors must be combined to create a coverage function vector for the ground station; therefore, the coverage function vector is the cumulative sum of the access functions of multiple satellites to the jth ground station; it stores information about the number of satellites seen at each discrete time l. The coverage function vector b j for:

[0065] b j =(b j [0],b j [1],…,b j [L - 1]) T

[0066] Among them, b j [l] is the coverage function of the constellation between the jth ground station at the lth discrete time, which obeys the relationship:

[0067]

[0068] Wherein, K represents the total number of satellites in the satellite constellation distribution.

[0069] The coverage requirement function vector of the jth ground station is recorded as is a non-negative integer vector created by the user according to task requirements. The coverage requirement function vector ξ j satisfy:

[0070] ξ j =(ξ j [0],ξ j [1],…,ξ j [L-1]) T ,

[0071] Among them, ξ j [l] is the number of satellite coverage layers required by the jth ground station at the lth discrete moment.

[0072] Constellation distribution vector represents the time shift of the satellite in orbit relative to the seed satellite. It is a (L×1)-dimensional column vector that satisfies:

[0073] x=(x[0],x[1],…,x[L-1]) T ,

[0074] Where x[l] is the time-shifted position point at l discrete moments relative to the seed satellite, and each element is defined as follows:

[0075]

[0076] Among them, the time positioning index l k is called the delay index of the kth satellite; this is because the time difference lags behind the seed satellite by Δt on the common ground track. k =tstep l k The kth satellite of the constellation distribution vector x is equivalently expressed as l = l k The unit pulse at time t, that is, the kth satellite is running in the system Δt k After time, it will occupy the current position of the seed satellite; whether the jth ground station is visible to the kth satellite at the lth moment is related to its k At each moment, there is the following corresponding relationship between whether the seed satellite is visible:

[0077] v k,j [l] = v 0,j [ll k ].

[0078] The access function vector and the coverage function vector have a cyclic convolution property: in this property, all the access function vectors of the member satellites in a given constellation are periodically the same, that is, the function value moves periodically in a cyclic manner; therefore, any access function v between the kth satellite and the jth ground station k,j Represented as a circular shift seed satellite access function vector v 0,j :

[0079]

[0080] Among them, the matrix P is a (L×L)-dimensional permutation matrix, which has the following expression:

[0081]

[0082] Considering the convolution of the access function and the constellation distribution vector, there is the following relationship:

[0083]

[0084] in, represents a circular convolution operator; or equivalently, this equation is written as:

[0085] V 0,j x=b j

[0086] Among them, the cyclic access matrix is the seed satellite access function vector v corresponding to the jth ground station 0,j The resulting circulant matrix is ​​a special form of the Toeplitz matrix, V 0,j =(v 0,j ,v 1,j ,…,v L-1,j ) is expanded to:

[0087]

[0088] Circulant Matrix V 0,j Each column element value is related to the circular shift seed satellite access function v 0,j Same; expressed in matrix form:

[0089]

[0090] The ground stations are concentrated in the same area of ​​interest, and the set is denoted as The jth ground communication demand at the lth time is ξ j [l] layers; Due to the limitation of the total number of satellites, the constellation cannot meet the communication needs at all times. Therefore, the ground station relaxes its communication demand restrictions and further limits the average number of coverage layers in the region. The average number of coverage layers of the ground station ζ(x) is the coverage function in time (demand time l′=0,1,…,L′-1) and space (ground station ) on average satisfaction:

[0091]

[0092] in, represents the total number of ground stations in the target area; L′ represents the total number of discrete times when all ground stations have communication needs. As a known quantity, it is determined by the communication needs of the ground stations. When calculating the average coverage satisfaction, only discrete moments with communication needs are considered; l′=0,1,…,L′-1 represents the sampling time; when the communication demand of the ground station is continuous coverage, L′=L.

[0093] In order to examine the communication coverage of ground station j at the lth moment, the instantaneous satisfaction d is defined as j [l]:

[0094]

[0095] Written as an expression about the constellation distribution vector x, it satisfies:

[0096]

[0097] Where U(·) represents the unit step function; since when b j [l] = ξ j [l] time j [l] = 1, and covers the demand function ( is a natural number set), so inside the step function U(·) we need to replace ξ j [l] is adjusted to (ξ j [l]-0.5).

[0098] Due to the limitation of the total number of satellites in the satellite constellation, the communication needs of all ground stations may not be fully met. Therefore, it is necessary to introduce two indicators, namely, the minimum instantaneous satisfaction of the region and the average coverage satisfaction of the constellation, to evaluate the performance of the designed satellite system. Consider the following scenario: the ground stations are concentrated in the same area of ​​interest, and the set is denoted as At the lth moment, all ground communication needs are met Due to the limitation of the total number of satellites, the constellation cannot meet the communication needs at all times. Therefore, the ground station relaxes its communication demand restrictions and further limits the instantaneous satisfaction of the region. The instantaneous satisfaction σ[l′] is expressed with respect to the constellation distribution vector x, satisfying:

[0099]

[0100] If at the l′th moment the ground station requires the satellite constellation to reach the regional minimum instantaneous satisfaction σ0[l′], then:

[0101] σ[l′]≥σ0[l′], l′=0,1,…,L′-1.

[0102] The average coverage satisfaction χ(x) is the average of the instantaneous satisfaction σ[l′] of the ground station at the demand time l′=0,1,…,L′-1, which is an expression of the constellation distribution vector x and is defined as:

[0103]

[0104] Step 2: Establish an optimization model based on the maximum average coverage satisfaction

[0105] Consider the following scenario: Due to cost constraints, the maximum number of satellites in the satellite constellation is K max , the minimum instantaneous coverage function required by the user for the ground station is It is required that within a communication demand cycle, the average number of coverage layers of the constellation on the ground station is at least ζ0. At this time, the user hopes that the average coverage satisfaction of the constellation is the best. Then the optimization model M1 is written as:

[0106]

[0107] The first constraint C1 means that the maximum number of satellites in the constellation system is K. max ; The second constraint C2 imposes a binary constraint on the constellation distribution vector x, requiring it to obey an L-dimensional 0-1 distribution to reflect the physical quantization of the satellite; the third constraint C3 requires the ground target point set to meet the minimum instantaneous satisfaction requirement of the region; the fourth constraint C4 gives a limit on the minimum average number of coverage layers.

[0108] The binary integer linear programming (BILP) method is used to solve the optimization model with the maximum average coverage satisfaction. However, it should be noted that the optimization model contains nonlinear terms such as step functions. Therefore, before using the BILP method to solve it, the nonlinear terms in the optimization model need to be processed.

[0109] In the optimization model, the instantaneous coverage satisfaction d at the l′th moment j [l′] has a step coupling relationship with the constellation distribution vector x. In order to facilitate the solution, the coupling relationship needs to be decoupled. The BIG-M method is used to decouple it, and the following constraints are used to perform equivalent replacement of the coupling relationship:

[0110]

[0111] Among them, M is a constant of a large order of magnitude, and M ≥ 10 6 ; ε is a very small number, 0<ε≤10 -6 .

[0112] Since the step function is decoupled using the BIG-M method, a new optimization variable D is added to the optimization model to satisfy where d j =(d j [0],d j [1],…,d j [L′-1]) T Therefore, the optimization variable is transformed from a single x to an augmented vector Q T =[D T ,x T ], the dimension changes from (L′×1) to dimension.

[0113] At this point, the optimized model M1 is transformed into model M2:

[0114]

[0115] Among them, M and ε are both constants.

[0116] First consider the objective function: The objective function is transformed by the BIG-M method into:

[0117]

[0118] in, is a row vector whose element values ​​are all 1, and the subscript is the dimension of the row vector; the matrix f satisfies:

[0119]

[0120] Among them, 0 1×LRepresents a 1×L-dimensional matrix of all zeros.

[0121] Step 3: Convert constraints to standard BILP constraint form

[0122] Consider the constraints: Let the constraints be AQ≤B, satisfying:

[0123]

[0124] Among them, A n ,B n (n=1,…,5) represent the coefficient matrices corresponding to the transformed constraints C1, C3, C4, C5, and C6 respectively; constraints C2 and C7 are binary constraints and do not need to be transformed.

[0125] Convert the constraint C1 into matrix form: The constraint C1 is written as:

[0126] 1 (1×L) x≤K max ,

[0127] The above formula is equivalent to:

[0128]

[0129] therefore,

[0130]

[0131] Convert the constraint C3 into matrix form: The constraint C3 is written as:

[0132]

[0133] The above formula is equivalent to:

[0134]

[0135] Among them, e i =(0,0,…,1,…,0) represents a unit row vector whose i-th element value is 1 and the others are all 0; I (L′×L′) represents the (L′×L′)-dimensional identity matrix; σ0 is the regional minimum coverage satisfaction vector, satisfying σ0=(σ0[0],…,σ0[l′],…,σ0[L′-1]) T ;

[0136] therefore,

[0137]

[0138] Convert the constraint C4 into matrix form: The constraint C4 is written as:

[0139]

[0140] The above formula is equivalent to:

[0141]

[0142] therefore,

[0143]

[0144] in is the new cyclic access matrix after removing the idle time access function from the cyclic access matrix, and its dimension is satisfy: Among them, the matrix V j(L′-1)+l′+1 is the access coefficient of the jth ground station at the l′th time, and its value is equal to the cyclic access matrix V 0,j The value of the j(L′-1)+l′th row of , whose dimension is (1×L).

[0145] Convert constraint C5 into matrix form: constraint C5 is written as:

[0146]

[0147] The above formula is equivalent to:

[0148]

[0149] Among them, the matrix is the demand function vector at the time when there is communication demand, where ξ j =(ξ j [0],…,ξ j [l′],…,ξ j [L′-1]) T ;

[0150] therefore,

[0151]

[0152] Convert constraint C6 into matrix form: Constraint C6 is written as:

[0153]

[0154] The above formula is equivalent to:

[0155]

[0156] therefore,

[0157]

[0158] In summary, the optimization model M2 is transformed into the following model M3:

[0159]

[0160] Among them, the dimension of A is The dimension of B is

[0161] Step 4: Model solution

[0162] The converted standard 0-1 integer optimization model form is solved by using the BILP method, including the implicit method or using a solver, including mosek and gurobi, to solve the model. The model is solved, and the optimal solution obtained is the M1 constellation distribution vector x of the model, thereby obtaining the optimal satellite constellation.

[0163] It should also be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner without contradiction. In order to avoid unnecessary repetition, the present invention will not further describe various possible combinations.

Claims

1. A method for designing a regional coverage satellite constellation based on optimal coverage satisfaction, characterized in that: The steps include: Step 1: Calculate the average coverage satisfaction of the common ground track constellation χ(x), where x is the constellation distribution vector, x = (x[0], x[1], …, x[l], …, x[L-1]) T , L represents the number of all discrete moments, l represents the lth discrete moment, x[l] represents the time-shifted position of the satellite relative to the seed satellite at an interval of l discrete moments, and the expression of x[l] is as follows: l k represents the delay index of the kth satellite; Step 2: Based on the average coverage satisfaction of the constellation and the minimum instantaneous satisfaction of the area of ​​interest as a constraint, establish an optimization model M1 based on the maximum average coverage satisfaction; Step 3: Convert the model M1 into the decoupled model M2, and then convert the M2 model into the standard 0-1 integer optimization model M3; Step 4: Solve M3 to obtain the optimal constellation distribution vector.

2. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 1, characterized in that: The common ground track constellation in step 1 satisfies the following conditions:

1. The orbital period of each satellite is a rational multiple of the Greenwich node period; 2. The semi-major axis a, eccentricity e, inclination i and perigee argument ω0 of all satellite orbits at the initial moment are the same; 3. Right ascension Ω of the ascending node of the kth satellite k Peace anomaly M k satisfy: N P Oh k +N D M k =constant mod(2π), Among them, N D is the number of satellite Greenwich node cycles, N P is the number of orbital periods, constant represents a constant, and mod represents the remainder.

3. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 1, characterized in that: The expression of the constellation average coverage satisfaction χ(x) in step 1 is as follows: Where L′ represents the total number of discrete moments when all ground stations have communication needs, represents the set of ground stations in the region of interest, represents the total number of ground stations in the target area, U(·) represents the unit step function, ξ j [l′] represents the number of satellite coverage layers of the jth ground station at the discrete moment when the l′th ground station has a communication demand, v 0,j [l′] represents the access function between the seed satellite and the jth ground station at the discrete moment when the l′th ground station has a communication demand. The expression of the access function between the kth satellite and the jth ground station at the lth discrete moment is shown below: Among them, when k = 0, it means the seed satellite, θ k,j [l] is the elevation angle between the kth satellite and the jth ground station at the lth discrete moment, θ k,j,min [l] is the minimum elevation angle between satellite k and ground station j at discrete time l; Represents the circular convolution operator.

4. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 3, characterized in that: The expression of the optimization model M1 in step 2 is as follows: C2:x[l]∈{0,1} C4:ζ(x)≥ζ0 Among them, K represents the intermediate variable, K max represents the maximum number of satellites, σ0[l′] represents the minimum instantaneous satisfaction of the area of ​​interest, σ[l′] represents the instantaneous satisfaction of the area of ​​interest, ζ(x) represents the average number of coverage layers of the ground station, ζ0 represents the minimum average number of coverage layers of the satellite constellation required by the user, and the expressions of ζ(x) and σ[l′] are as follows: Among them, d j [l'] is the instantaneous satisfaction d when l = l' j The value of [l], d j [l]Satisfy the following relationship: Among them, ξ j [l] is the satellite coverage required by the jth ground station at the lth discrete moment, b j [l] is the coverage function imposed by the constellation on the jth ground station at the lth discrete moment, b j The expression of [l] is: Wherein, K represents the total number of satellites in the satellite constellation distribution.

5. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 4, characterized in that: The covering functions of L discrete moments form a covering function vector b j , b j The expression is: b j =(b j [0],b j [1],…,b j [l],…,b j [L-1]) T , b j And the constellation distribution vector x satisfies the following relationship: V 0,j x=b j Among them, V 0,j represents the access circulant matrix, V 0,j The expression is: Then b j [l] and x[l] satisfy the following step coupling relationship:

6. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 5, characterized in that: The expression of the model M3 is as follows: stC1:A1Q≤B1 C2:x[l]∈{0,1} C3:A2Q≤B2 C4:A3Q≤B3 C5:A4Q≤B4 C6:A5Q≤B5 Among them, Q is the augmented vector, Q T =[D T ,x T ], D is the instantaneous coverage satisfaction vector, 0 1×L represents a 1×L-dimensional all-zero matrix, express dimensional all-1 matrix; the expressions of A1, A2, A3, A4, A5, B1, B2, B3, B4, and B5 are as follows: Where I represents the identity matrix, the subscript of I represents the dimension of the identity matrix, and σ0 is the minimum coverage instantaneous satisfaction vector of the region of interest, satisfying σ0=(σ0[0],…,σ0[l′],…,σ0[L′-1]) T , is the new cyclic access matrix after removing the idle time access function from the cyclic access matrix, V j(L′-1)+l′+1 is the access function matrix of the jth ground station at the l′th time, V j(L′-1)+l′+1 is the cyclic access matrix V 0,j The j(L′-1)+l′th row of ξ j =(ξ j [0],…,ξ j [l′],…,ξ j [L′-1]) T , M and ε are both constants.

7. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 4, characterized in that: The expression of the model M2 is as follows: C2:x[l]∈{0,1} C4:ζ(x)≥ζ0 Among them, M and ε are constants, Q is the augmented vector, Q T =[D T ,x T ], D is the instantaneous coverage satisfaction vector, 8. The method for designing a regional coverage satellite constellation based on optimal coverage satisfaction according to claim 1, characterized in that: The step 4 is solved by using the BILP method.