Area coverage constellation orbit optimization design method

Abstract

The invention discloses an area coverage constellation orbit optimization design method. The method comprises the steps that a plurality of target points in a target area are selected and numbered; searching regression orbits based on a genetic algorithm to obtain orbital elements of the regression orbits; calculating a time interval of a common sub-satellite point trajectory constellation; and according to the orbital elements and the time intervals of the regression orbits, obtaining the orbital elements of the common sub-satellite point trajectory constellation. According to the method, the constellation orbit design with any precision requirement is realized, and the design efficiency is improved.

Description

Technical Field

[0001] The invention relates to a method for optimizing the design of a regional coverage constellation orbit, and belongs to the technical field of constellation orbit design. Background Art

[0002] According to the different target areas, satellite constellations can be divided into global coverage constellations, latitude band coverage constellations and regional coverage constellations. Global coverage constellations can cover the entire world; latitude band coverage constellations can cover the strip area between two equal latitudes in the northern and southern hemispheres; regional coverage constellations only need to cover a fixed ground area, emphasizing coverage of specific areas. In terms of orbit design of regional coverage constellations, there is currently little research because the problem requires a huge number of variables to be designed and the design method is highly dependent on the target area.

[0003] Reference [1] adopted a design method for constructing a satellite constellation based on multiple common subsatellite point trajectories to meet the observation needs of 225 target points in a specific area on the ground. However, the method adopted in this reference requires solving a type of optimization problem each time when calculating the semi-major axis of the orbit, which greatly increases the complexity of the design process. In addition, the design accuracy of the method adopted in this reference is limited, and designers are required to manually fine-tune the calculation results, which makes it difficult to ensure the optimality of the design results. Reference [2] proposed a constellation optimization design method based on low-density grid coverage analysis, but this method requires optimization of 8 variables and the optimization process is too complicated. Reference [3] proposed a weighted coverage of key areas in a specific area and integrated the genetic ant colony algorithm for optimization and solution, but after integrating the two heuristic algorithms, more manual settings are required. parameters, which puts higher requirements on designers. In addition, compared with a single optimization algorithm, the algorithm design process after fusion is more complicated, and the results obtained are not substantially improved. Reference [4] proposes a method for solving the orbit parameters of a single satellite based on a particle swarm algorithm and solving the orbit parameters of a constellation based on an analytical method. However, this method only contains a common sub-satellite point trajectory constellation and cannot guarantee full coverage of the target area. Patent [5] discloses a method for designing a low-orbit common ground trajectory constellation covering a local area, but this method lacks a global search process for the number of orbital elements and is difficult to guarantee a global optimal solution. Patent [6] discloses a satellite constellation design method for covering multiple target areas with discrete distributions, but the accuracy of the orbital elements solved by this method is limited and still cannot solve the problem in reference [1].

[0004] [1] MA Jian, MENG Yazhe, ZHU Xiaolong, et al. Optimal design method for low-orbit satellite constellation for intensive observation in specific areas[J]. SCIENCE CHINA: Technological Sciences, 2018, 48(02): 170-184.

[0005] [2] Ma Yuanye, Zheng Jianhua, Li Mingtao. Optimal design of revisit constellation based on low-density grid coverage analysis[J]. Chinese Journal of Space Science, 2020, 40(04): 554-561.

[0006] [3] Hua Bing, Wang Peiyuan, Wu Yunhua, et al. Regional coverage satellite constellation optimization algorithm based on genetic ant colony algorithm[J]. Journal of National University of Defense Technology, 2023, 45(02):167-171.

[0007] [4] Liu Ming, Yang Ming, Gao Xing, et al. Design of continuous coverage constellation for irregular regions based on particle swarm[J]. Aerospace Control, 2020, 38(02): 31-37.

[0008] [5] Zhang Xiang, Lei Jizhao, Zhao Hang, et al. A method for designing a low-orbit common ground track constellation covering a local area [P]. Chongqing: CN202110284531.8, 2022-06-10.

[0009] [6]Qiao Penghao, Qian Mingjing, Li Jianqing. A satellite constellation design method for discretely distributed multi-target area coverage[P]. Beijing: CN202311536225.4, 2024-05-03. Summary of the invention

[0010] The technical problem solved by the present invention is: to overcome the deficiencies of the prior art, to provide a method for optimizing the design of a regional coverage constellation orbit, and to achieve constellation orbit design with arbitrary precision requirements.

[0011] The technical solution of the present invention is:

[0012] The present invention discloses a method for optimizing the orbit of a regional coverage constellation, comprising:

[0013] Select several target points in the target area and number them;

[0014] Search the regression orbit based on the genetic algorithm and obtain the orbital elements of the regression orbit;

[0015] Calculate the time interval of the common subsatellite point trajectory constellation;

[0016] According to the orbital elements and time interval of the regression orbit, the orbital elements of the common subsatellite point trajectory constellation are obtained.

[0017] Furthermore, in the above method, the plurality of target points in the target area are selected and numbered, and the specific method is as follows:

[0018] Calculate the semi-geocentric angle of a single satellite's coverage area;

[0019] Select a target point every half geocentric angle in longitude and latitude in the target area;

[0020] Number each target point.

[0021] Furthermore, in the above method, the genetic algorithm is used to search for the regression orbit to obtain the number of orbital elements of the regression orbit. The specific method is:

[0022] S31, setting the search accuracy of orbital inclination i, ascending node right ascension Ω, and mean anomaly M; setting the search range of orbital inclination, ascending node right ascension, and mean anomaly;

[0023] S32. Generate S groups of regression orbits according to the search accuracy and search range, each group including N regression orbits, and obtain the orbital inclination i, the right ascension Ω of the ascending node, the mean anomaly M, the eccentricity e and the argument of perigee ω of each regression orbit;

[0024] S33, solving the semi-major axis a according to the orbital inclination i and the eccentricity e;

[0025] S34, calculating the number of target points within the coverage of the constellation according to the semi-major axis a, orbit inclination i, right ascension of ascending node Ω, mean anomaly M, eccentricity e and argument of perigee ω;

[0026] S35, converting the orbital inclinations, right ascensions of ascending nodes, and mean anomalies of all regression orbits in each group of regression orbits from decimal to binary and sequentially connecting them to form a binary code for each group of regression orbits, thereby obtaining S groups of binary codes for the population;

[0027] S36, randomly changing or exchanging the binary codes of any two groups of regression tracks to obtain new binary codes;

[0028] S37, converting the new binary code into decimal to obtain S new groups of regression orbits, and calculating the number of target points that each group of regression orbits can cover;

[0029] S38, judging whether there is a regression track that can cover all target points according to the number of target points, if so, proceeding to step S40; otherwise, proceeding to step S39;

[0030] S39, determine whether the population generation number reaches the set upper limit, if yes, proceed to step S40; if no, repeat steps S35 to S38;

[0031] S40, determine whether the search accuracy reaches the required accuracy; if so, proceed to step S41; otherwise, update the search accuracy and search range; repeat steps S32 to S39;

[0032] S41, let N = N-1, repeat steps S31 to S40 until the minimum N that can cover all target points is obtained; record the number of orbital elements of the minimum N regression orbits (ak ,e k ,i k ,Ω k ,ω k ,M k ),k=1,2,…,N.

[0033] Furthermore, in the above method, the updating of the search precision and search range is specifically as follows: setting the search precision of the orbital inclination i, the ascending node right ascension Ω, and the mean anomaly M to 0.1 times of the original precision; setting the search range of the orbital inclination to [i op -Δi,i op +Δi], the search range of the ascending node right ascension is set to [Ω op -ΔΩ,Ω op +ΔΩ]; the search range of the mean anomaly angle is set to [M op -ΔM,M op +ΔM]; where i op is the orbital inclination angle searched in step S38 or S39, Δi is the search accuracy of the original orbital inclination angle; Ω op is the right ascension of the ascending node obtained by searching in step S38 or S39, ΔΩ is the search accuracy of the right ascension of the original ascending node, M op is the mean anomaly angle obtained by searching in step S38 or S39, and ΔM is the search accuracy of the original mean anomaly angle.

[0034] Furthermore, in the above method, in step S34, the number of target points within the coverage of the constellation is calculated by:

[0035] S341, calculating the orbital elements at the current moment;

[0036] S342, converting the orbital elements to positions in the geocentric coordinate system;

[0037] S343. Calculate the sub-satellite position of the j-th satellite at the current moment according to the position in the geocentric coordinate system: Where R is the position in the geocentric coordinate system, a is the semi-major axis, and R e is the average radius of the Earth;

[0038] S344. Calculate the semi-geocentric angle Λ of each satellite coverage area j , j = 1, 2, ..., N;

[0039] S345, calculate the sub-satellite position r of any target point k and any satellite j Dj The geocentric angular distance between kj ; Among them, r Tk is the position vector of the kth target point in the geocentric coordinate system; k = 1, 2, ..., N total, N total is the total number of target points.

[0040] S346. If for any j, there is Λ kj >Λ j , then the target point k is not within the coverage of the constellation, otherwise the target point k is within the coverage of the constellation;

[0041] S347, repeat S341 to S346 until the regression period is reached, and count the number of target points N within the coverage range of the constellation. fg .

[0042] Furthermore, in the above method, the time interval for calculating the common subsatellite point trajectory constellation is specifically:

[0043]

[0044] Among them, Δt is the time interval between two adjacent satellite sub-satellite points passing the same target point. is the round-up symbol, T hg is the regression period of the regression orbit, T cf is the maximum revisit time, N s is the number of satellites contained in a common subsatellite point trajectory constellation.

[0045] Furthermore, in the above method, the orbital elements of the common subsatellite point trajectory constellation are obtained according to the orbital elements and time interval of the regression orbit, and the specific method is:

[0046] S61, the orbital root number (a k ,e k ,i k ,Ω k ,ω k ,M k ), k = 1, 2, ..., N as the orbital element number of the first satellite in each common subsatellite point trajectory constellation; N is the regression orbit number;

[0047] S62, satellite S k,n The orbital elements are converted to satellite S k,n Position vector and velocity vector in the geocentric inertial coordinate system; where S k,n is the nth satellite in the kth common subsatellite point trajectory constellation;

[0048] S63, according to satellite S k,n The position vector and velocity vector in the geocentric inertial coordinate system are obtained as follows: k,n+1 Position vector and velocity vector in the geocentric inertial coordinate system; where t0 is the initial time, Δt is the time interval of the subsatellite point trajectory constellation;

[0049] S64, the satellite S at time t0+Δt k,n+1 The position vector and velocity vector in the geocentric inertial coordinate system are transformed into the satellite S at time t0+Δt k,n+1 The number of orbital elements of

[0050] S65, for satellite S at time t0+Δt k,n+1 The orbital elements of the satellite S are obtained by integrating forward for time Δt. k,n+1 The number of orbital elements at time t0;

[0051] S66, repeat steps S61 to S65 to obtain N s The number of orbital elements of satellites at time t0 forms a system containing N·N s The regional coverage constellation of satellites; where N is the number of regression orbits, N s is the number of satellites contained in a common subsatellite point trajectory constellation.

[0052] Furthermore, in the above method, in step S63, at time t0+Δt, satellite S k,n+1 The position vector and velocity vector in the geocentric inertial coordinate system are:

[0053]

[0054] Among them, M EI (t0+Δt) is the coordinate transformation matrix from the geocentric coordinate system to the geocentric inertial coordinate system at time t0+Δt, ω E is the Earth's rotation angular velocity vector, r Ik,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The position in the geocentric inertial coordinate system, v Ik,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The velocity in the Earth-centered inertial coordinate system, r Ek,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The position in the geocentric coordinate system, v Ek,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The velocity in the geocentric coordinate system, r Ek,n (t0) is the satellite S at time t0 k,n The position in the geocentric coordinate system, v Ek,n (t0) is the satellite S at time t0 k,n The velocity in the geocentric coordinate system.

[0055] Further, in the above method, in step S33, the semi-major axis a is solved as follows:

[0056]

[0057] h min +R e ≤a≤h max +R e

[0058] Among them, ω E is the Earth's rotation angular rate, μ is the Earth's gravitational constant, J2 is the second-order band harmonic coefficient, P and Q are mutually prime natural numbers, h min is the minimum required orbital height, h max is the required maximum orbital height; e is the eccentricity, and i is the orbital inclination.

[0059] Furthermore, in the above method, in step S341, the number of orbital elements is calculated, specifically:

[0060]

[0061] Where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, M is the mean anomaly, μ is the Earth's gravitational constant, and R e is the average radius of the earth, and J2 is the second-order harmonic coefficient.

[0062] Furthermore, in the above method, the semi-geocentric angle of the coverage range of a single satellite is specifically:

[0063]

[0064] Among them, R e is the average radius of the Earth, h min is the minimum value required for the constellation orbit height, γ d is the distance from the satellite to the edge of the ground coverage area, η is the satellite half-viewing angle, and γ is the auxiliary angle.

[0065] The beneficial effects of the present invention and the prior art are:

[0066] (1) The method disclosed in the present invention can realize constellation orbit design with arbitrary precision requirements by continuously improving the search precision of orbit elements and reducing the search range of orbit elements, thus solving the problem of limited search precision and the need for manual adjustment in the prior art;

[0067] (2) The method disclosed in the present invention provides the necessary conditions that the semi-major axis needs to meet under the earth's J2 perturbation, solves the problem of optimizing the semi-major axis in the relevant literature, and improves the design efficiency;

[0068] (3) The method disclosed in the present invention solves the problem that the existing technology may not be able to achieve full coverage of the target area by calculating the semi-geocentric angle of the satellite coverage range to select the discrete spacing of the target point.

[0069] (4) The present invention solves the problems in the prior art such as limited design accuracy, a large number of optimization variables or design parameters, and difficulty in covering the entire target area. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 It is a flow chart of regression track search based on genetic algorithm of the present invention. DETAILED DESCRIPTION

[0071] The present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0072] The present invention discloses a method for optimizing the orbit of a regional coverage constellation, comprising:

[0073] Select several target points in the target area and number them;

[0074] Search the regression orbit based on the genetic algorithm and obtain the orbital elements of the regression orbit;

[0075] Calculate the time interval of the common subsatellite point trajectory constellation;

[0076] According to the orbital elements and time interval of the regression orbit, the orbital elements of the common subsatellite point trajectory constellation are obtained.

[0077] Preferably, several target points in the target area are selected and numbered, and the specific method is as follows:

[0078] Calculate the semi-geocentric angle of a single satellite's coverage area;

[0079] Select a target point every half geocentric angle in longitude and latitude in the target area;

[0080] Number each target point.

[0081] Preferably, the regression orbit is searched based on a genetic algorithm to obtain the number of orbital elements of the regression orbit. The specific method is:

[0082] S31, setting the search accuracy of orbital inclination i, ascending node right ascension Ω, and mean anomaly M; setting the search range of orbital inclination, ascending node right ascension, and mean anomaly;

[0083] S32. Generate S groups of regression orbits according to the search accuracy and the search range, each group containing N regression orbits, and obtain the orbital inclination i, the right ascension Ω of the ascending node, the mean anomaly M, the eccentricity e and the argument of perigee ω of each regression orbit;

[0084] S33, solving the semi-major axis a according to the orbital inclination i and the eccentricity e;

[0085] S34. Calculate the number of target points within the coverage of the constellation according to the semi-major axis a, the orbit inclination i, the right ascension of the ascending node Ω, the mean anomaly M, the eccentricity e and the argument of perigee ω;

[0086] S35, converting the orbital inclinations, right ascensions of ascending nodes, and mean anomalies of all regression orbits in each group of regression orbits from decimal to binary and sequentially connecting them to form a binary code for each group of regression orbits, thereby obtaining S groups of binary codes for the population;

[0087] S36, randomly changing or exchanging the binary codes of any two groups of regression tracks to obtain new binary codes;

[0088] S37, converting the new binary code into decimal to obtain S new groups of regression orbits, and calculating the number of target points that each group of regression orbits can cover;

[0089] S38, judging whether there is a regression track that can cover all target points according to the number of target points, if so, proceeding to step S40; otherwise, proceeding to step S39;

[0090] S39, determine whether the population generation number reaches the set upper limit, if yes, proceed to step S40; if no, repeat steps S35 to S38;

[0091] S40, determine whether the search accuracy reaches the required accuracy; if so, proceed to step S41; otherwise, update the search accuracy and search range; repeat steps S32 to S39;

[0092] S41, let N = N-1, repeat steps S31 to S40 until the minimum N that can cover all target points is obtained; record the number of orbital elements of the minimum N regression orbits (a k ,e k ,i k ,Ω k ,ω k ,M k ),k=1,2,…,N.

[0093] Preferably, the search accuracy and search range are updated, specifically: the search accuracy of the orbital inclination i, the ascending node right ascension Ω, and the mean anomaly M are set to 0.1 times the original accuracy; the search range of the orbital inclination is set to [i op -Δi,i op +Δi], the search range of the ascending node right ascension is set to [Ω op -ΔΩ,Ω op +ΔΩ]; the search range of the mean anomaly angle is set to [M op -ΔM,M op +ΔM]; where i opis the orbital inclination angle searched in step S38 or S39, Δi is the search accuracy of the original orbital inclination angle; Ω op is the right ascension of the ascending node obtained by searching in step S38 or S39, ΔΩ is the search accuracy of the right ascension of the original ascending node, M op is the mean anomaly angle obtained by searching in step S38 or S39, and ΔM is the search accuracy of the original mean anomaly angle.

[0094] Preferably, in step S34, the number of target points within the coverage of the constellation is calculated, and the specific method is:

[0095] S341, calculating the orbital elements at the current moment;

[0096] S342, converting the orbital elements to positions in the geocentric coordinate system;

[0097] S343. Calculate the sub-satellite position of the j-th satellite at the current moment according to the position in the geocentric coordinate system: Where R is the position in the geocentric coordinate system, a is the semi-major axis, and R e is the average radius of the Earth;

[0098] S344. Calculate the semi-geocentric angle Λ of each satellite coverage area j , j = 1, 2, ..., N;

[0099] S345, calculate the sub-satellite position r of any target point k and any satellite j Dj The geocentric angular distance between kj ; Among them, r Tk is the position vector of the kth target point in the geocentric coordinate system; k = 1, 2, ..., N total , N total is the total number of target points.

[0100] S346. If for any j, there is Λ kj >Λ j , then the target point k is not within the coverage of the constellation, otherwise the target point k is within the coverage of the constellation;

[0101] S347, repeat S341 to S346 until the regression period is reached, and count the number of target points N within the coverage range of the constellation. fg .

[0102] Preferably, the time interval of the common subsatellite point trajectory constellation is calculated as follows:

[0103]

[0104] Among them, Δt is the time interval between two adjacent satellite sub-satellite points passing the same target point. is the round-up symbol, T hg is the regression period of the regression orbit, T cf is the maximum revisit time, N s is the number of satellites contained in a common subsatellite point trajectory constellation.

[0105] Preferably, the orbital elements of the common subsatellite point trajectory constellation are obtained according to the orbital elements and time interval of the regression orbit, and the specific method is:

[0106] S61, the orbital root number (a k ,e k ,i k ,Ω k ,ω k ,M k ), k = 1, 2, ..., N as the orbital element number of the first satellite in each common subsatellite point trajectory constellation; N is the regression orbit number;

[0107] S62, satellite S kn The orbital elements are converted to satellite S kn Position vector and velocity vector in the geocentric inertial coordinate system; where S k,n is the nth satellite in the kth common subsatellite point trajectory constellation;

[0108] S63, according to satellite S k,n The position vector and velocity vector in the geocentric inertial coordinate system are obtained as follows: k,n+1 Position vector and velocity vector in the geocentric inertial coordinate system; where t0 is the initial time, Δt is the time interval of the subsatellite point trajectory constellation;

[0109] S64, the satellite S at time t0+Δt k,n+1 The position vector and velocity vector in the geocentric inertial coordinate system are transformed into the satellite S at time t0+Δt k,n+1 The number of orbital elements of

[0110] S65, for satellite S at time t0+Δt k,n+1 The orbital elements of the satellite S are obtained by integrating forward for time Δt. k,n+1 The number of orbital elements at time t0;

[0111] S66, repeat steps S61 to S65 to obtain N s The number of orbital elements of satellites at time t0 forms a system containing N·N s The regional coverage constellation of satellites; where N is the number of regression orbits, N s is the number of satellites contained in a common subsatellite point trajectory constellation.

[0112] Preferably, in step S63, at time t0+Δt, satellite Sk,n+1 The position vector and velocity vector in the geocentric inertial coordinate system are:

[0113]

[0114] Among them, M EI (t0+Δt) is the coordinate transformation matrix from the geocentric coordinate system to the geocentric inertial coordinate system at time t0+Δt, ω E is the Earth's rotation angular velocity vector, r Ik,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The position in the geocentric inertial coordinate system, v Ik,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The velocity in the Earth-centered inertial coordinate system, r Ek,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The position in the geocentric coordinate system, v Ek,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The velocity in the geocentric coordinate system, r Ek,n (t0) is the satellite S at time t0 k,n The position in the geocentric coordinate system, v Ek,n (t0) is the satellite S at time t0 k,n The velocity in the geocentric coordinate system.

[0115] Preferably, in step S33, the semi-major axis a is solved as follows:

[0116]

[0117] h min +R e ≤a≤h max +R e

[0118] Among them, ω E is the Earth's rotation angular rate, μ is the Earth's gravitational constant, J2 is the second-order band harmonic coefficient, P and Q are mutually prime natural numbers, h min is the minimum required orbital height, h max is the required maximum orbital height; e is the eccentricity, and i is the orbital inclination.

[0119] Preferably, in step S341, the number of orbital elements is calculated as follows:

[0120]

[0121] Where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, M is the mean anomaly, μ is the Earth's gravitational constant, and Re is the average radius of the earth, and J2 is the second-order harmonic coefficient.

[0122] Preferably, the semi-geocentric angle of the coverage range of a single satellite is specifically:

[0123]

[0124] r d =R e cosγ+(h min +R e )cosη

[0125]

[0126] Among them, R e is the average radius of the Earth, h min is the minimum value required for the constellation orbit height, γ d is the distance from the satellite to the edge of the ground coverage area, η is the satellite half-viewing angle, and γ is the auxiliary angle.

[0127] Example

[0128] like Figure 1 As shown, the present invention provides a method for optimizing the orbit of a regional coverage constellation, comprising the following steps:

[0129] Step 1: discretize the target area.

[0130] (1) According to the minimum value h required by the constellation orbit height min And the auxiliary angle γ is calculated from the half field of view angle η of a single satellite:

[0131]

[0132] Where R e =6371004m is the average radius of the earth

[0133] (2) Calculate the distance r between the edge of a single satellite and the satellite d :

[0134] r d =R e cosγ+(h min +R e )cosη

[0135] (3) Calculate the semi-geocentric angle Λ of the coverage area of ​​a single satellite

[0136]

[0137] (4) Select a target point every Λ in the longitude and latitude of the target area, take all the target points as the objects to be covered, and number each target point, record as target point 1, target point 2, ..., target point N total .

[0138] Step 2: Search the regression trajectory based on genetic algorithm.

[0139] (1) Set the number of regression orbits N contained in a constellation;

[0140] (2) The search precision of orbit inclination i, ascending node right ascension Ω, and mean anomaly M is set to 1°, the search range of orbit inclination is set to [-90°, 90°], and the search ranges of ascending node right ascension and mean anomaly are set to [0, 360°];

[0141] (3) Generate S groups of regression orbits, each group contains N regression orbits. The orbital inclination i, the right ascension Ω of the ascending node, and the mean anomaly M of each regression orbit are random values ​​generated according to the set search accuracy and search range. The eccentricity e and the argument of perigee ω are both 0. The semi-major axis a is solved according to the following formula:

[0142]

[0143] in

[0144]

[0145] In the above formula, ω E =7.292115×10 -5 rad / s is the Earth's rotation angular velocity, μ = 3.986004418 × 10 14 m 3 / s 2 is the Earth’s gravitational constant, J2=0.1082629832258×10 -2 is the second-order band harmonic coefficient, P and Q are mutually prime natural numbers, and by changing the values ​​of P and Q, the semi-major axis a is within the required range, that is, h min +R e ≤a≤h max +R e ,h min is the minimum required orbital height, h max is the maximum required orbital height.

[0146] (4) Each group of regression tracks is regarded as an individual, and S groups of regression tracks together constitute a population. In a regression cycle, the target point numbers that each individual can cover are counted, and the number of covered target points is regarded as the individual's fitness. The method for calculating the number of target points covered by a certain physical fitness is:

[0147] I) Calculate the orbital elements (a, e, i, Ω, ω, M) according to the following formula:

[0148]

[0149] II) Convert the orbital elements (a, e, i, Ω, ω, M) to the position and velocity (R, V) in the geocentric coordinate system. Then the sub-satellite position of the jth satellite at the current moment is

[0150]

[0151] III) Calculate the semi-geocentric angle Λ of each satellite coverage area according to step 1 j (j=1,2,…,N);

[0152] IV) Calculate any target point k (k = 1, 2, ..., N) total ) and the sub-satellite position r of any satellite j Dj The geocentric angular distance between kj :

[0153]

[0154] where r Tk is the position vector of the kth target point in the geocentric coordinate system;

[0155] V) Determine whether the target point k is within the coverage of the constellation. If for any j,

[0156] Λ kj >Λ j

[0157] Then the target point k is not within the coverage of the constellation, otherwise the target point k is within the coverage of the constellation;

[0158] VI) repeating steps I)-V) until the regression cycle is reached;

[0159] VII) Count the number of target points N within the coverage of the constellation fg .

[0160] (5) converting the orbital inclination, right ascension of the ascending node, and mean anomaly of all the regression orbits in each regression orbit group from decimal to binary and then concatenating them in sequence to form the binary code of the individual;

[0161] (6) After obtaining the S groups of binary codes of the population, randomly change or exchange the binary codes of the individuals and perform crossover and mutation genetic algorithm operations on the population;

[0162] (7) Convert the binary code of each individual after step (6) back to decimal to obtain new S individuals, and calculate the number of target points that each individual can cover according to step (4);

[0163] (8) If the optimal individual in the population can cover all target points, or the population generation number reaches the upper limit, then stop the iteration; otherwise, continue to generate the next generation of population;

[0164] (9) Repeat steps (3)-(8) up to a maximum number of repetitions k max ;

[0165] (10) Determine whether the search accuracy reaches the required accuracy. If it reaches the required accuracy, proceed to step (11). If not, set the search accuracy of the orbital inclination i, the right ascension of the ascending node Ω, and the mean anomaly M to 0.1 times the original accuracy, and set the search range to the original search accuracy superimposed on the last calculation result of step (9), that is, set the search range of the orbital inclination to [i op -Δi,i op +Δi], where i op is the orbit inclination calculated for the last time in step (9), Δi is the original search accuracy, and the search range of the ascending node right ascension is set to [Ω op -ΔΩ,Ω op +ΔΩ], where Ω op is the right ascension of the ascending node calculated for the last time in step (9), ΔΩ is the original search accuracy, and the search range of the mean anomaly angle is set to [M op -ΔM,M op +ΔM], where M op is the orbital inclination calculated for the last time in step (9), and ΔM is the original search accuracy. Repeat steps (3)-(9);

[0166] (11) Determine whether N has reached the minimum value. If the optimal individual can cover all target points, set N = N-1 and repeat steps (2)-(10) until the minimum N that can cover all target points is obtained. When N is the minimum, the corresponding optimal individual is the optimal regression orbit. The number of orbital elements of these N regression orbits is recorded as (a k ,e k ,i k ,Ω k ,ω k ,M k ),k=1,2,…,N.

[0167] Step 3: Calculate the distribution interval of the common subsatellite point trajectory constellation.

[0168] The maximum revisit time T cf , as the time interval T between the subsatellite points of two adjacent satellites in the common subsatellite point trajectory constellation passing through the same target point g , thus obtaining the number of satellites N that should be deployed on the corresponding common subsatellite point trajectory constellation s :

[0169]

[0170] In the formula is the round-up symbol, T hg is the regression period of the regression orbit. Therefore, the time interval Δt between two adjacent satellite sub-satellite points passing the same target point is:

[0171]

[0172] Step 4: Calculate the number of orbital elements of the common subsatellite point trajectory constellation.

[0173] (1) Substitute the orbital elements (a) obtained in step 2 k ,e k ,i k ,Ω k ,ω k ,M k ), k = 1, 2, ..., N as the orbital element number of the first satellite in each subsatellite point trajectory constellation;

[0174] (2) Name the nth satellite in the kth common subsatellite point trajectory constellation as satellite S k,n ;

[0175] (3) S k,1 Orbital element number (a k ,e k ,i k ,Ω k ,ω k ,M k ) is converted to the position vector r in the geocentric inertial coordinate system Ik,1 (t0) and velocity v Ik,1 (t0), thus obtaining its position vector r in the geocentric coordinate system Ek,1 (t0) and velocity vector v Ek,1 (t0):

[0176]

[0177] Among them, M IE (t0) is the coordinate transformation matrix from the geocentric inertial coordinate system to the geocentric coordinate system at time t0, ω E is the Earth's rotation angular velocity vector.

[0178] (4) Satellite S k,n The position vector r in the geocentric coordinate system at time t0+(n-1)Δt Ek,n (t0+(n-1)Δt) and satellite S k,1 The position vector r in the geocentric coordinate system at time t0 Ek,1 (t0) coincides, and the satellite Sk,n The velocity vector v in the geocentric coordinate system at time t0+(n-1)Δt Ek,n (t0+(n-1)Δt) and satellite S k,n The velocity vector v in the geocentric coordinate system at time t0 Ek,1 (t0) is the same, that is,

[0179]

[0180] (5) Satellite S k,n The position vector and velocity vector in the geocentric inertial coordinate system are:

[0181]

[0182] Among them, M EI (t0+Δt) is the coordinate transformation matrix from the geocentric coordinate system to the geocentric inertial coordinate system at time t0+Δt;

[0183] (6) Satellite S k,n The position vector and velocity vector in the geocentric inertial coordinate system are converted into orbital elements, and the orbital perturbation equation containing the J2 term is integrated forward for (n-1)Δt to obtain the satellite S k,n The six elements of the orbit at time t0 (a k,n ,e k,n ,i k,n ,Ω k,n ,ω k,n ,M k,n ). This is to numerically integrate the differential equations in step 2 (4) I), and the total duration of the integration is -(n-1)△t.

[0184] (7) Similarly, we can calculate the number of orbital elements (a k ,e k ,i k ,Ω k ,ω k ,M k ), k=1,2,…,N expands to N s satellites, and we get a total of N·N s Regional coverage constellation of satellites.

[0185] Although the content of the present invention has been described in detail through the above preferred embodiments, it should be appreciated that the above description should not be considered as a limitation of the present invention. After reading the above content, it will be apparent to those skilled in the art that various modifications and substitutions of the present invention will occur. Therefore, the protection scope of the present invention should be limited by the appended claims.

[0186] The contents not described in detail in the specification of the present invention belong to the common knowledge of the professionals in this field.

Claims

1. A method for optimizing the orbit of a regional coverage constellation, characterized in that: include: Select several target points in the target area and number them; Search the regression orbit based on the genetic algorithm and obtain the orbital elements of the regression orbit; Calculate the time interval of the common subsatellite point trajectory constellation; According to the orbital elements and time interval of the regression orbit, the orbital elements of the common subsatellite point trajectory constellation are obtained.

2. The method for optimizing the regional coverage constellation orbit design according to claim 1, characterized in that: The method of selecting and numbering several target points in the target area is as follows: Calculate the semi-geocentric angle of a single satellite's coverage area; Select a target point every half geocentric angle in longitude and latitude in the target area; Number each target point.

3. The method for optimizing the regional coverage constellation orbit design according to claim 1, characterized in that: The genetic algorithm is used to search for the regression orbit and obtain the number of orbital elements of the regression orbit. The specific method is as follows: S31, setting the search accuracy of orbital inclination i, ascending node right ascension Ω, and mean anomaly M; setting the search range of orbital inclination, ascending node right ascension, and mean anomaly; S32. Generate S groups of regression orbits according to the search accuracy and search range, each group including N regression orbits, and obtain the orbital inclination i, the right ascension Ω of the ascending node, the mean anomaly M, the eccentricity e and the argument of perigee ω of each regression orbit; S33, solving the semi-major axis a according to the orbital inclination i and the eccentricity e; S34, calculating the number of target points within the coverage of the constellation according to the semi-major axis a, orbit inclination i, ascending node right ascension Ω, mean anomaly M, eccentricity e and argument of perigee ω; S35, converting the orbital inclinations, right ascensions of ascending nodes, and mean anomalies of all regression orbits in each group of regression orbits from decimal to binary and sequentially connecting them to form a binary code for each group of regression orbits, thereby obtaining S groups of binary codes for the population; S36, randomly changing or exchanging the binary codes of any two groups of regression tracks to obtain new binary codes; S37, converting the new binary code into decimal to obtain S new groups of regression orbits, and calculating the number of target points that each group of regression orbits can cover; S38, judging whether there is a regression track that can cover all target points according to the number of target points, if so, proceeding to step S40; Otherwise, go to step S39; S39, determine whether the population generation number reaches the set upper limit, if so, proceed to step S40; If not, repeat steps S35 to S38; S40, determining whether the search accuracy reaches the required accuracy; if so, proceeding to step S41; Otherwise, update the search accuracy and search range; repeat steps S32 to S39; S41, let N = N-1, repeat steps S31 to S40, until the minimum N that can cover all target points is obtained; record the number of orbital elements of the minimum N regression orbits (a k ,e k ,i k ,Ω k ,ω k ,M k ),k=1,2,…,N.

4. The method for optimizing the regional coverage constellation orbit design according to claim 3, characterized in that: The updating of the search precision and search range is specifically as follows: the search precision of the orbital inclination i, the ascending node right ascension Ω, and the mean anomaly angle M is set to 0.1 times the original precision; the search range of the orbital inclination is set to [i op -Δi,i op +Δi], the search range of the ascending node right ascension is set to [Ω op -ΔΩ,Ω op +ΔΩ]; the search range of the mean anomaly angle is set to [M op -ΔM,M op +ΔM]; where i op is the orbital inclination angle searched in step S38 or S39, Δi is the search accuracy of the original orbital inclination angle; Ω op is the right ascension of the ascending node obtained by searching in step S38 or S39, ΔΩ is the search accuracy of the right ascension of the original ascending node, M op is the mean anomaly angle obtained by searching in step S38 or S39, and ΔM is the search accuracy of the original mean anomaly angle.

5. The method for optimizing the regional coverage constellation orbit design according to claim 3, characterized in that: In step S34, the number of target points within the coverage of the constellation is calculated, and the specific method is: S341, calculating the orbital elements at the current moment; S342, converting the orbital elements to positions in the geocentric coordinate system; S343. Calculate the sub-satellite position of the j-th satellite at the current moment according to the position in the geocentric coordinate system: Where R is the position in the geocentric coordinate system, a is the semi-major axis, and R e is the average radius of the Earth; S344. Calculate the semi-geocentric angle Λ of each satellite coverage area j , j = 1, 2, ..., N; S345, calculate the sub-satellite position r of any target point k and any satellite j Dj The geocentric angular distance between kj ; Among them, r Tk is the position vector of the kth target point in the geocentric coordinate system; k = 1, 2, ..., N total , N total is the total number of target points. S346. If for any j, there is Λ kj >Λ j , then the target point k is not within the coverage of the constellation, otherwise the target point k is within the coverage of the constellation; S347, repeat S341 to S346 until the regression period is reached, and count the number of target points N within the coverage range of the constellation. fg .

6. The method for optimizing the regional coverage constellation orbit design according to claim 1, characterized in that: The time interval for calculating the common subsatellite point trajectory constellation is specifically: Among them, Δt is the time interval between two adjacent satellite sub-satellite points passing the same target point. is the round-up symbol, T hg is the regression period of the regression orbit, T cf is the maximum revisit time, N s is the number of satellites contained in a common subsatellite point trajectory constellation.

7. The method for optimizing the regional coverage constellation orbit design according to claim 6, characterized in that: The orbital elements and time interval of the regression orbit are used to obtain the orbital elements of the common subsatellite point trajectory constellation. The specific method is as follows: S61, the orbital root number (a k ,e k ,i k ,Ω k ,ω k ,M k ), k = 1, 2, ..., N as the orbital element number of the first satellite in each common subsatellite point trajectory constellation; N is the regression orbit number; S62, satellite S k,n The orbital elements are converted to satellite S k,n Position vector and velocity vector in the geocentric inertial coordinate system; where S k,n is the nth satellite in the kth common subsatellite point trajectory constellation; S63, according to satellite S k,n The position vector and velocity vector in the geocentric inertial coordinate system are obtained as follows: k,n+1 Position vector and velocity vector in the geocentric inertial coordinate system; where t0 is the initial time, Δt is the time interval of the subsatellite point trajectory constellation; S64, the satellite S at time t0+Δt k,n+1 The position vector and velocity vector in the geocentric inertial coordinate system are transformed into the satellite S at time t0+Δt k,n+1 The number of orbital elements of S65, for satellite S at time t0+Δt k,n+1 The orbital elements of the satellite S are obtained by integrating forward for time Δt. k,n+1 The number of orbital elements at time t0; S66, repeat steps S61 to S65 to obtain N s The number of orbital elements of satellites at time t0 forms a system containing N·N s The regional coverage constellation of satellites; where N is the number of regression orbits, N s is the number of satellites contained in a common subsatellite point trajectory constellation.

8. The method for optimizing the regional coverage constellation orbit design according to claim 7, characterized in that: In step S63, at time t0+Δt, satellite S k,n+1 The position vector and velocity vector in the geocentric inertial coordinate system are: Among them, M EI (t0+Δt) is the coordinate transformation matrix from the geocentric coordinate system to the geocentric inertial coordinate system at time t0+Δt, ω E is the Earth's rotation angular velocity vector, r Ik,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The position in the geocentric inertial coordinate system, v Ik,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The velocity in the Earth-centered inertial coordinate system, r Ek,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The position in the geocentric coordinate system, v Ek,n+1 (t0+Δt) is the satellite S at time t0+Δt k,n+1 The velocity in the geocentric coordinate system, r Ek,n (t0) is the satellite S at time t0 k,n The position in the geocentric coordinate system, v Ek,n (t0) is the satellite S at time t0 k,n The velocity in the geocentric coordinate system.

9. The method for optimizing the regional coverage constellation orbit design according to claim 3, characterized in that: In step S33, the semi-major axis a is solved as follows: h min +R e ≤a≤h max +R e Among them, ω E is the Earth's rotation angular rate, μ is the Earth's gravitational constant, J2 is the second-order band harmonic coefficient, P and Q are mutually prime natural numbers, h min is the minimum required orbital height, h max is the required maximum orbital height; e is the eccentricity, and i is the orbital inclination.

10. The method for optimizing the regional coverage constellation orbit design according to claim 4, characterized in that: In step S341, the number of orbital elements is calculated, specifically: Where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the right ascension of the ascending node, ω is the argument of perigee, M is the mean anomaly, μ is the Earth's gravitational constant, and R e is the average radius of the earth, and J2 is the second-order harmonic coefficient.

11. The method for optimizing the regional coverage constellation orbit design according to claim 2, characterized in that: The semi-geocentric angle of a single satellite’s coverage is: r d =R e cosγ+(h min +R e )cosη Among them, R e is the average radius of the Earth, h min is the minimum value required for the constellation orbit height, γ d is the distance from the satellite to the edge of the ground coverage area, η is the satellite half-viewing angle, and γ is the auxiliary angle.

Images