A heterogeneous satellite constellation configuration correction method based on stable earth coverage performance

Abstract

The application discloses a heterogeneous satellite constellation configuration correction method based on stable earth coverage performance. The method comprises the following steps: constructing a first query table and a second query table; for each orbit plane, using a working orbit parameter optimization method, combining the first query table and the second query table, analyzing and calculating dynamic visibility conditions and a visibility list, and then correcting the working semi-major axis and the latitude argument AoL of each satellite on the orbit plane; for each orbit plane, according to the working semi-major axis correction value of the satellite on the orbit plane, combining the first query table, calculating the right ascension of the ascending node RAAN correction value of the satellite; and the working semi-major axis correction value, the latitude argument AoL correction value and the right ascension of the ascending node RAAN correction value of each satellite on all orbit planes constitute the configuration parameters of the heterogeneous satellite constellation. The application restrains the drift of the earth coverage index of the heterogeneous satellite constellation, reduces the constellation control frequency, and avoids the additional increase of the number of satellites due to the drift of the coverage index and the multiple artificial iterative design.

Description

Technical Field

[0001] This invention relates to a method for correcting the configuration of heterogeneous constellations, and more particularly to a method for correcting the configuration of heterogeneous constellations with Earth overlay. Background Technology

[0002] Heterogeneous constellations refer to constellations composed of satellites or formations with different functions, appearances, physical structures, and payloads. Such constellations can integrate different types of satellites, such as those for communication, navigation, and remote sensing. Through information interconnection and resource integration, they can improve the constellation's operational efficiency, functional density, cost-effectiveness, and inter-satellite collaboration capabilities, thereby enhancing the system's resilience and conserving orbital slot resources. This represents an inevitable trend in current and future constellation development.

[0003] Compared to traditional isomorphic constellations, the differences in the windward surface area ratio and sensor coverage of individual satellites in heterogeneous constellations are the main factors affecting constellation configuration design and stability. The varying coverage of each satellite means that the analysis of constellation coverage capabilities and the design of its configuration cannot rely on the homogenization, uniformity, and linear superposition approaches of traditional constellations. Furthermore, the differences in drift characteristics among heterogeneous satellites cause drift in inter-satellite configurations, coverage characteristics, and coverage indicators, increasing the number of satellites and control costs compared to traditional isomorphic constellation design and control methods. Among these, coverage drift caused by latitudinal argument (AoL) drift in coplanar following configurations is the most severe, while semi-major axis drift and right ascension of the ascending node (RAAN) drift also adversely affect the stability of constellation configuration and coverage indicators.

[0004] Constellation design methods can be broadly categorized based on their design principles into three types: those based on classical configurations, analytical methods based on geometric principles, and those based on modern optimization methods and intelligent algorithms. Classical configuration-based methods commonly include the Walker constellation and its derivatives, the Iridium constellation, etc. Analytical methods based on geometric principles include the coverage zone method and the ground trajectory method, which have rigorous mathematical proofs and can be flexibly customized according to requirements, typically for special configurations such as return orbits. Design methods based on modern optimization methods and intelligent algorithms focus on establishing and solving optimization models, rather than directly arranging orbits using sub-satellite point projections. However, current extensive research has only addressed specific coverage problems of homogeneous constellations, failing to solve the analytical analysis and evaluation of coverage performance for heterogeneous constellations under varying satellite coverage capabilities, the drift problem of configuration and coverage indicators caused by differences in satellite drift trajectories, and the multi-faceted solution to heterogeneous constellation design problems related to coverage effectiveness, configuration stability, and coverage indicator stability. Furthermore, the performance analysis and design of heterogeneous constellations rely on simulation iterations, significantly increasing the complexity of manual design.

[0005] In the field of heterogeneous constellation design, some general methods have been provided for the design of heterogeneous constellation configurations. However, the heterogeneous meaning involved is singular and applicable to a variety of heterogeneous satellite payloads. However, the stability drift of coverage indicators caused by the difference in the surface mass ratio of heterogeneous satellites has not been taken into account, nor has the constellation design problem based on the stability of ground coverage performance when heterogeneous payloads and heterogeneous surface mass ratios act simultaneously been considered.

[0006] In summary, traditional constellation configuration design methods typically aim to minimize costs or maximize coverage metrics. However, because they do not adequately consider the configuration stability and coverage stability of heterogeneous constellations, directly applying these methods to heterogeneous constellations will lead to coverage metrics drifting over time. This necessitates the addition of extra satellites or constellation configuration maintenance, thereby increasing the manufacturing and control costs of the constellation. Summary of the Invention

[0007] To address the problem of effectively optimizing the configuration of heterogeneous constellations for Earth coverage, this invention proposes a method for correcting the configuration of heterogeneous satellite constellations based on the stability of Earth coverage performance. This method not only achieves analytical calculation of coverage indicators for heterogeneous constellations through geometric methods, but also solves the problem of Earth coverage indicator drift in medium and low Earth orbit (LEO) heterogeneous constellations by correcting satellite orbital parameters. It also addresses the issues of high constellation control frequency or additional satellite numbers caused by Earth coverage indicator drift. This invention provides an effective analytical method for optimizing the configuration design of LEO circular heterogeneous Earth coverage constellations, avoiding multiple manual iterative designs.

[0008] The technical solution adopted in this invention is:

[0009] The method of the present invention specifically includes the following steps:

[0010] Step S1: Construct a first lookup table and a second lookup table. The first lookup table contains the semi-major axis attenuation parameters and RAAN change rate of satellites with different areal-to-mass ratios under each characteristic semi-major axis; the second lookup table contains the coverage geocentric angle of satellites with different coverage semi-subtraction angles under each characteristic semi-major axis, as well as the initial RAAN-AoL visibility conditions and initial RAAN visibility feature range for each ground target at the initial time.

[0011] Step S1 includes:

[0012] First, coverage data and orbital drift data under different characteristic semi-major axes are obtained; the coverage data includes the geocentric angle of coverage of satellites with different coverage semi-angles under each characteristic semi-major axis; the orbital drift data includes the semi-major axis attenuation parameters and RAAN change rate of satellites with different surface-to-mass ratios under each characteristic semi-major axis, and the semi-major axis attenuation parameters include quadratic semi-major axis attenuation parameters and linear semi-major axis attenuation parameters.

[0013] Finally, combining the satellite-ground visibility conditions, the initial RAAN-AoL visibility condition set and the initial RAAN visibility characteristic range set of each ground target of the satellite are obtained.

[0014] Among them, the initial RAAN-AoL visibility conditions and the initial RAAN visibility characteristic range of the satellite for each ground target are calculated through the following process: First, according to the characteristic semi-major axis and the covering semi-aperture angle of the satellite, one or two sets of initial RAAN visibility conditions of the satellite for the ground target are calculated. Each set of initial RAAN visibility conditions mainly consists of boundary conditions (two boundary points), a center point and a visibility range. Subsequently, for each set of initial RAAN visibility conditions, according to the boundary conditions and the center point, several RAAN characteristic points including the boundary condition Ω L 、Ω R and the center point Ω M are generated through discretization, and the two AoL boundary points corresponding to each RAAN characteristic point are calculated. Polynomial fitting is performed on all RAAN characteristic points and the corresponding AoL boundary points to obtain the initial RAAN-AoL visibility conditions.

[0015] The process of one or two sets of initial RAAN visibility conditions of the satellite for the ground target specifically refers to: When the satellite satisfies φ n <i - θ for the ground target n, the satellite has two sets of initial RAAN visibility conditions for the ground target n. The RAAN difference between the two sets of boundary points in the two sets of initial RAAN visibility conditions is the same, but the values of the two sets of boundary points are different; when the satellite satisfies i - θ < φ n <θ for the ground target n, the satellite has one set of initial RAAN visibility conditions for the ground target n. Among them, φ n represents the latitude of the ground target n, θ represents the covering geocentric angle of the satellite, and i represents the orbital inclination of the orbital plane to which the satellite belongs.

[0016] Step S2: For each orbital plane, adopt the working orbit parameter optimization method, combine the first query table and the second query table, analyze and calculate the dynamic visibility conditions and the visibility list, and then correct the working semi-major axis and the latitude argument of latitude AoL of each satellite on the orbital plane. Finally, the corrected values of the working semi-major axis and the latitude argument of latitude AoL of each satellite on the orbital plane are obtained.

[0017] In the said step S2, the working orbit parameter optimization method includes:

[0018] Step S2.1: For each satellite on the orbital plane, through the analytical method, combine the first query table and the second query table to obtain the initial working semi-major axis, the initial AoL, the covering geocentric angle and the maximum AoL offset of the satellite.

[0019] The said analytical method includes:

[0020] For each satellite on the orbital plane, the satellite's semi-major axis attenuation parameter is obtained from the first lookup table based on the satellite's area-to-mass ratio and nominal semi-major axis; the satellite's coverage geocentric angle is obtained from the second lookup table based on the satellite's coverage semi-angle and nominal semi-major axis.

[0021] The initial working semi-major axis is calculated based on the satellite's nominal semi-major axis and semi-major axis attenuation parameters;

[0022] The nominal AoL of each satellite is obtained based on the total number of satellites in the orbital plane;

[0023] The maximum AoL offset is obtained based on the satellite's initial working semi-major axis, nominal semi-major axis, and semi-major axis attenuation parameters.

[0024] The initial AoL of the satellite is obtained based on the satellite's nominal AoL and maximum AoL offset.

[0025] Step S2.2: If the maximum AoL offset of each satellite on the orbital plane is less than or equal to twice the coverage geocentric angle, proceed to step S2.3; otherwise, proceed to step S2.5.

[0026] Step S2.3: Based on the initial working semi-major axis and initial AoL of all satellites on the orbital plane, obtain the coverage index of the orbital plane within the mission cycle by analytically calculating the dynamic visibility conditions and visibility list.

[0027] Step S2.3 includes:

[0028] Step S2.3.1: For each satellite on the orbital plane, based on the satellite's surface-to-mass ratio, coverage half-angle, and initial working half-major axis, obtain the satellite's RAAN change rate, coverage geocentric angle, and initial RAAN-AoL visibility conditions and initial RAAN visibility feature range for each ground target from the first lookup table and the second lookup table, and then obtain the maximum concentrated visibility period of the satellite for each ground target.

[0029] The maximum concentrated visible period of the satellite for each ground target is expressed as follows:

[0030]

[0031] In the formula, Γ Ω0,p,s A1 represents the maximum concentrated visible period of satellite s on orbital plane p relative to ground target n; A2 represents the initial offset coefficient; t represents the long-term drift ratio coefficient. Γ1 or t Γ2 t represents the duration of the single maximum concentrated visible period. B +t Γ1 or t B +tΓ2 The periodicity of the maximum concentrated visible time is represented by τ1, τ2, or τ3, which represents the long-term drift time. φ n Let θ represent the latitude of the ground target n, θ represent the geocentric angle of the coverage of satellite s, and i represent the orbital inclination of orbital plane p.

[0032] Step S2.3.2: For each satellite on the orbital plane, obtain a list of visible times for each ground target based on the maximum concentrated visible time period of the satellite for each ground target and the initial AoL.

[0033] Step S2.3.2 includes:

[0034] Based on the satellite's surface-to-mass ratio, coverage half-angle, and initial working half-major axis, the orbital attenuation coefficient, RAAN rate of change, and the initial RAAN-AoL visibility conditions of the satellite on the ground target are obtained from the first and second lookup tables.

[0035] Iterate through each moment in the maximum concentrated visibility period, and obtain the RAAN and AoL values ​​of the satellite at the current moment based on the satellite's orbital attenuation coefficient, RAAN change rate, and initial AoL. Map the initial RAAN-AoL visibility conditions of the satellite to the ground target to the current moment to obtain the RAAN-AoL visibility conditions of the satellite to the ground target at the current moment. Determine the visibility based on the RAAN value, AoL value, and RAAN-AoL visibility conditions of the satellite to the ground target at the current moment: if the RAAN value is within the RAAN visibility conditions and the AoL value is within the AoL visibility conditions, then the satellite is visible to the ground target at the current moment; otherwise, the satellite is not visible to the ground target at the current moment.

[0036] By integrating the visibility assessment results of all times within the maximum concentrated visibility period, a list of times when the satellite is visible to ground targets is obtained.

[0037] Step S2.3.3: Merge the list of visible times of all satellites in the orbital plane to the same ground target to obtain a list of visible time periods of all satellites in the orbital plane to the ground target.

[0038] Step S2.3.4: Based on the list of visible time periods for each ground target by all satellites in the orbital plane, calculate the coverage index and coverage index consistency for each ground target by all satellites in the orbital plane.

[0039] The coverage index for each ground target by all satellites in the orbital plane is the average visibility duration or average revisit time of all satellites in the orbital plane for the ground target.

[0040] Step S2.3.5: Obtain the average coverage index and the consistency of the average coverage index of all satellites on the orbital plane over all ground targets, which together form the coverage index of the orbital plane.

[0041] Step S2.4: If the coverage index of the orbital plane meets the preset conditions, the initial working semi-major axis and initial AoL obtained in step S2.1 are used as the correction values ​​for the working semi-major axis and latitude argument AoL of the satellite, respectively; otherwise, proceed to step S2.5. In specific implementation, the preset conditions can be determined according to the coverage performance requirements in the mission design phase.

[0042] Step S2.5: Using an optimization algorithm, the initial working semi-major axis and initial AoL of all satellites on the orbital plane are used as variables to be optimized, and the coverage index of the orbital plane within the mission cycle is minimized as the optimization objective. During the optimization process, the coverage index is obtained according to step S2.3, and finally the working semi-major axis correction value and latitude argument AoL correction value of each satellite on the orbital plane are obtained.

[0043] Step S3: For each orbital plane, based on the working semi-major axis correction value of each satellite on the orbital plane and in conjunction with the first lookup table, calculate the right ascension (RAAN) correction value of the ascending node of each satellite.

[0044] Step S3 includes the following steps:

[0045] For each satellite on the orbital plane, the orbital attenuation parameters are obtained from the first lookup table based on the satellite's working semi-major axis correction value and surface mass ratio, and then the RAAN change of the satellite during the mission cycle is calculated.

[0046] Calculate the average RAAN change of all satellites on the orbital plane during the mission period;

[0047] The difference between the nominal RAAN and the average RAAN change of the orbital plane is obtained to obtain the RAAN correction value of the orbital plane.

[0048] For each satellite on the orbital plane, the difference between the RAAN correction value on the orbital plane and the change in RAAN of the satellite during the mission period is obtained to obtain the RAAN correction value of the right ascension of the ascending node of the satellite.

[0049] Step S4: The working semi-major axis correction value, latitude argument AoL correction value, and ascending node right ascension RAAN correction value of each satellite on all orbital planes constitute the configuration parameters of the heterogeneous satellite constellation.

[0050] The beneficial effects of this invention are:

[0051] 1. This invention, based on traditional uniform constellations such as the Walker constellation, obtains differentiated working orbit parameters for each satellite by compensating and correcting three types of orbital parameters: semi-major axis, AoL, and RAAN. This increases the stability of the Earth coverage index of the low-to-medium orbit circular heterogeneous constellation composed of satellites with different non-surface ratios and different coverage half-angles. Stability is represented by the consistency of coverage indexes, i.e., the coverage indexes remain close to the design values ​​throughout the mission cycle. Commonly used coverage indexes include: revisit time, transit time, number of transits, coverage overlap, and the time required for complete imaging of the covered area.

[0052] 2. This invention achieves the beneficial effects of curbing the drift of ground coverage indicators in heterogeneous ground coverage constellations, reducing constellation control frequency, avoiding additional increases in the number of satellites due to coverage indicator drift, and avoiding multiple manual iterative designs. Attached Figure Description

[0053] Figure 1 This is a schematic diagram showing the trajectory of the satellite and the coverage of ground stations.

[0054] Figure 2 The relationship between the shape of the visible range of the right ascension-latitude argument of the ascending node and the orbital inclination and the latitude of the ground station.

[0055] Figure 3 A schematic diagram showing the geometric relationship between the satellite coverage half-angle and the coverage geocentric angle.

[0056] Figure 4 This refers to the visible range of the right ascension of the ascending node and the period of maximum concentrated visibility.

[0057] Figure 5 This shows the change in the visible range of the right ascension-latitude argument of the ascending node over time.

[0058] Figure 6 The satellite trajectory is obtained by analytically solving the satellite orbital parameters.

[0059] Figure 7 This is a schematic diagram of the overall process of the method of the present invention.

[0060] Figure 8 This is a flowchart illustrating the method of the present invention. Detailed Implementation

[0061] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. The specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of the invention in any way.

[0062] The method of the present invention specifically includes the following steps:

[0063] Step S1: Construct the first query table and the second query table; the purpose of this step is to save subsequent computing resources. Among them, the first query table contains the semi-major axis decay parameters and RAAN change rates of satellites with different area-to-mass ratios at each characteristic semi-major axis. The second query table contains the covered geocentric angles of satellites with different covered half angles at each characteristic semi-major axis, as well as the initial RAAN-AoL visibility conditions and the initial RAAN visibility characteristic ranges of each satellite for each ground target.

[0064] Specifically, step S1 includes:

[0065] First, obtain the coverage data and orbit drift data at different characteristic semi-major axes; the coverage data includes the covered geocentric angles of satellites with different covered half angles at each characteristic semi-major axis; the orbit drift data includes the semi-major axis decay parameters and RAAN change rates of satellites with different area-to-mass ratios at each characteristic semi-major axis, and the semi-major axis decay parameters include quadratic semi-major axis decay parameters and linear semi-major axis decay parameters;

[0066] Finally, combined with the satellite-ground visibility conditions, obtain the set of initial RAAN-AoL visibility conditions and the set of initial RAAN visibility characteristic ranges of each satellite for each ground target. Among them, the RAAN-AoL visibility condition refers to the right ascension of the ascending node RAAN and the argument of latitude AoL conditions that need to be satisfied for the satellite and the ground target to be visible.

[0067] Furthermore, the initial RAAN-AoL visibility conditions and the initial RAAN visibility characteristic ranges of each satellite for each ground target can be calculated through the following process:

[0068] First, according to the characteristic semi-major axis a of the satellite at the initial moment r,x and the covered half angle α y , calculate one or two sets of initial RAAN visibility conditions of the satellite for the ground target. Each set of initial RAAN visibility conditions (Ω L , Ω R , Ω M , △L) is mainly composed of the boundary conditions Ω A , Ω D , the center point Ω M and the visible range △L. Specifically: when the satellite satisfies φ n <i - θ for the ground target n, the satellite has two sets of initial RAAN visibility conditions for the ground target n. The RAAN difference between the two boundary points in the two sets of initial RAAN visibility conditions is the same, but the values of the two boundary points are different; when the satellite satisfies i - θ < φ n <θ for the ground target n, the satellite has one set of initial RAAN visibility conditions for the ground target n. Among them, φ nLet θ represent the latitude of the ground target n, θ represent the geocentric angle of the satellite's coverage, and i represent the orbital inclination of the satellite's orbital plane.

[0069] Finally, for each set of initial time RAAN visibility conditions (Ω) L Ω R Ω M Based on the boundary conditions and the center point, several boundary conditions Ω are generated through discretization. L Ω R and the center point Ω M The RAAN feature points are included, and the two AoL boundary points u1 and u2 corresponding to each RAAN feature point are calculated. Polynomial fitting is performed on all RAAN feature points and their corresponding AoL boundary points to obtain an initial RAAN-AoL visibility condition f. n (Ω t=0 (u1, u2), the initial RAAN-AoL visibility condition is a function f of the visible polygon. n (Ω t=0 (u1, u2).

[0070] Furthermore, satellite semi-major axis attenuation data can be obtained through the following process: using simulation methods, simulate the trajectory of a satellite located on a specific characteristic semi-major axis and having a specific surface-to-mass ratio within a mission duration T, and extract the semi-major axis change data and the RAAN change ΔΩ. x,z A quadratic fitting method was used to fit the semi-major axis variation data to obtain the quadratic and linear semi-major axis attenuation parameters. The RAAN variation ΔΩ was then used. x,z Divide the RAAN rate of change by the simulation task duration T.

[0071] Step S2: For each orbital plane, the working orbital parameter optimization method is adopted. Combining the first lookup table and the second lookup table, the dynamic visibility conditions and visibility list are calculated and analyzed, and then the working semi-major axis and latitudinal argument AoL of each satellite on the orbital plane are corrected.

[0072] In step S2, the optimization method for working trajectory parameters includes:

[0073] Step S2.1: For each satellite on the orbital plane, the initial working semi-major axis, initial AoL, coverage geocentric angle, and maximum AoL offset of the satellite are obtained by analytical method, combined with the first lookup table and the second lookup table.

[0074] The analytical method is as follows:

[0075] For each satellite on the orbital plane, the satellite's semi-major axis attenuation parameter is obtained from the first lookup table based on the satellite's surface-to-mass ratio and nominal semi-major axis, and the satellite's coverage geocentric angle is obtained from the second lookup table based on the satellite's coverage semi-angle and nominal semi-major axis.

[0076] The initial working semi-major axis is calculated based on the satellite's nominal semi-major axis and semi-major axis attenuation parameters;

[0077] The maximum AoL offset is obtained based on the satellite's initial working semi-major axis, nominal semi-major axis, and semi-major axis attenuation parameters.

[0078] Based on the total number of satellites in orbit, the nominal AoL of each satellite is calculated using the following formula:

[0079] u s =[2(s-1)π] / S p , s∈[1,S p ]

[0080] In the formula, u s 'Indicates the nominal AoL, S of satellite s p This represents the total number of satellites in the orbital plane where satellite s is located;

[0081] The initial AoL of each satellite is calculated using the following formula, based on the nominal AoL and maximum AoL offset of each satellite:

[0082] u s =u s '-△u max,s / 2

[0083] In the formula, u s Denotes the initial AoL, u of satellite s s 'Indicates the nominal AoL of satellite s, Δu max,s This represents the maximum AoL offset of satellite s.

[0084] Step S2.2: If the maximum AoL offset of each satellite on the orbital plane is less than or equal to twice its own coverage geocentric angle, that is, the conditions for using the analytical method are met, then proceed to step S2.3; otherwise, proceed to step S2.5.

[0085] Step S2.3: Based on the initial working semi-major axis and initial AoL of all satellites on the orbital plane, obtain the coverage index of the orbital plane within the mission cycle by analytically calculating the dynamic visibility conditions and visibility list. The coverage index of the orbital plane includes the average coverage index of all satellites on the orbital plane over all ground targets and the consistency of the average coverage index.

[0086] Step S2.3 includes:

[0087] Step S2.3.1: For each satellite on the orbital plane, based on the satellite's surface-to-mass ratio, coverage half-angle, and initial working half-major axis, obtain the satellite's RAAN change rate, coverage geocentric angle, and initial RAAN-AoL visibility conditions and initial RAAN visibility feature range for each ground target from the first lookup table and the second lookup table, and then obtain the maximum concentrated visibility period of the satellite for each ground target.

[0088] The maximum period of concentrated visibility of a satellite for each ground target is expressed as:

[0089]

[0090] In the formula, Γ Ω0,p,s A1 represents the maximum concentrated visible period of satellite s on orbital plane p relative to ground target n; A2 represents the initial offset coefficient; t represents the long-term drift ratio coefficient. Γ1 or t Γ2 t represents the duration of the single maximum concentrated visible period. B +t Γ1 or t B +t Γ2 The periodicity of the maximum concentrated visible time is represented by τ1, τ2, or τ3, which represents the long-term drift time. φ n Let θ represent the latitude of the ground target n, θ represent the geocentric angle of the coverage of satellite s, and i represent the orbital inclination of orbital plane p.

[0091] Step S2.3.2: For each satellite on the orbital plane, obtain a list of visible times for each ground target based on the maximum concentrated visible time period of the satellite for each ground target and the initial AoL.

[0092] Step S2.3.2 includes:

[0093] Based on the satellite's surface-to-mass ratio, coverage half-angle, and initial working half-major axis, the orbital attenuation coefficient, RAAN rate of change, and the initial RAAN-AoL visibility conditions of the satellite on the ground target are obtained from the first and second lookup tables.

[0094] By iterating through each moment in the maximum concentrated visible period, the RAAN value and AoL value of the satellite at the current moment are obtained based on the satellite's orbital attenuation coefficient, RAAN change rate, and initial AoL.

[0095] The initial RAAN-AoL visibility conditions of the satellite to the ground target are mapped to the current time to obtain the RAAN-AoL visibility conditions of the satellite to the ground target at the current time.

[0096] Visibility is determined based on the satellite's RAAN value, AoL value, and RAAN-AoL visibility conditions for ground targets at the current moment: if the RAAN value is within the RAAN visibility condition and the AoL value is within the AoL visibility condition, then the satellite is visible to the ground targets at the current moment; otherwise, the satellite is not visible to the ground targets at the current moment.

[0097] By integrating the visibility assessment results of all times within the maximum concentrated visibility period, a list of times when the satellite is visible to ground targets is obtained.

[0098] Step S2.3.3: Merge the lists of visible times of all satellites in the orbital plane to the same ground target to obtain a list of visible time periods of all satellites in the orbital plane to the ground target.

[0099] Step S2.3.4: Based on the list of visible time periods for each ground target by all satellites in the orbital plane, calculate the coverage index and coverage index consistency for each ground target by all satellites in the orbital plane.

[0100] Optionally, the coverage index for each ground target by all satellites in the orbital plane is the average visibility duration or average revisit time of the ground target by all satellites in the orbital plane.

[0101] Step S2.3.5: Obtain the average coverage index and average coverage index consistency of all satellites on the orbital plane for all ground targets. The two together constitute the coverage index of the orbital plane.

[0102] Step S2.4: If the coverage indicators of the orbital plane meet the preset conditions, then the initial working semi-major axis and initial AoL obtained in step S2.1 are used as the correction values ​​for the satellite's working semi-major axis and latitude argument AoL, respectively; otherwise, proceed to step S2.5. In specific implementation, the preset conditions can be determined according to the coverage performance requirements in the mission design phase.

[0103] Step S2.5: Using an optimization algorithm, the initial working semi-major axis and initial AoL of all satellites on the orbital plane are used as variables to be optimized, and the coverage index of the orbital plane within the mission cycle is minimized as the optimization objective. During the optimization process, the coverage index is obtained according to step S2.3, and finally the working semi-major axis correction value and latitude argument AoL correction value of each satellite on the orbital plane are obtained.

[0104] Step S3: For each orbital plane, calculate the right ascension (RAAN) correction value of the ascending node of each satellite based on the working semi-major axis correction value of each satellite in the orbital plane and in conjunction with the first lookup table; the working ascending node right ascension correction can ensure the relative configuration stability between orbital planes under the different working semi-major axes.

[0105] Step S3 includes the following steps:

[0106] For each satellite on the orbital plane, the orbital attenuation parameters are obtained from the first lookup table based on the satellite's working semi-major axis correction value and surface mass ratio, and then the RAAN change of the satellite during the mission cycle is calculated.

[0107] Calculate the average RAAN change of all satellites on the orbital plane during the mission period;

[0108] The difference between the nominal RAAN and the average RAAN change of the orbital plane is obtained to obtain the RAAN correction value of the orbital plane.

[0109] For each satellite on the orbital plane, the difference between the RAAN correction value on the orbital plane and the change in RAAN of the satellite during the mission period is obtained to obtain the RAAN correction value of the right ascension of the ascending node of the satellite.

[0110] In step S3.1, the change in RAAN of the satellite during the mission period is calculated according to the following formula:

[0111] △Ω=-(7Ω'·k1·T 3 ) / (6a)-(7Ω'·k2·T 2 ) / (4a)-[(7Ω'·T) / (2a)]·(a s -a)

[0112] Ω'=-[(3n e ·J2·R e 2 ) / (2a 2 )]·cosi

[0113] In the formula, ΔΩ represents the change in RAAN of the satellite during the mission period, Ω' represents the long-term rate of change of RAAN, k1 represents the semi-major axis attenuation parameter of the quadratic term, k2 represents the semi-major axis attenuation parameter of the linear term, T represents the mission period, and a represents the nominal semi-major axis of the satellite. s n represents the initial working semi-major axis of the satellite. e J2 represents the orbital angular velocity, and R represents the Earth's oblateness perturbation coefficient. e denoted by , where i represents the Earth's radius and i represents the orbital inclination.

[0114] Step S4: Use the working semi-major axis correction value, latitude argument AoL correction value, and ascending node right ascension RAAN correction value of all satellites on all orbital planes as the configuration parameters of the heterogeneous satellite constellation.

[0115] Specific embodiments of the present invention are as follows:

[0116] Example 1

[0117] In this embodiment, the analytical method mainly uses the geometric method proposed by Ulybyshev Y et al. to realize the analytical calculation of the heterogeneous constellation coverage index, and further derives parameters such as the satellite-ground centralized visible period on this basis. The geometric method adopted by the analytical method maps the relationship between the right ascension of the ascending node and the latitude argument of the satellite that satisfies the satellite-ground visibility condition changing with time into a graph, which is called the satellite-ground visibility relationship mapping graph, and then calculates the satellite-ground visible period and other coverage parameters through the graph parameters. For the detailed content, please refer to the following papers of this team:

[0118] Ulybyshev Y. Satellite Constellation Design for Complex Coverage[J]. Journal of Spacecraft and Rockets, 2008, 45(4): 843-849.

[0119] Ulybyshev, Yuri. Geometric Analysis and Design Method for Discontinuous Coverage Satellite Constellations[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(2): 549-557.

[0120] The geometric method adopted by the analytical method in this embodiment is described below:

[0121] Figure 1 is a schematic diagram of the sub-satellite point trajectory and the ground station coverage. For the target ground station M(λ, φ) with the geocentric longitude and latitude being λ and φ respectively, coverage can be achieved when the sub-satellite point trajectory passes through the circular area with M(λ, φ) as the center and θ as the radius. Figure 1 In, i represents the orbital inclination, △λ represents the longitude difference between the boundary orbits that can achieve ground station coverage, and θ represents the coverage geocentric angle.

[0122] As Figure 1 shown in (a) of, when φ < i - θ, both the up-track and the down-track may achieve coverage. In the figure, L A1 , L A2 , L A3 are the up-tracks, L D1 , L D2 , L D3 are the down-tracks, L A2 and L D2 are the longest tracks passing through the center of the circle, and the rest are the tracks tangent to the coverage circle, and the tangent points N A , S A , ND and S D are boundary points. Among them, N A and N D are uniformly denoted as the northern boundary point N, and S A and S D are uniformly denoted as the southern boundary point S. The subscript A represents the upward orbit, and D represents the downward orbit.

[0123] As Figure 1 shown in (b) of, when i - θ < φ < i + θ, the trajectory where the inflection point is within the coverage circle can achieve coverage. In the figure, L2 is the longest trajectory, and the rest are the trajectories tangent to the coverage circle. The tangent points S A and S D are boundary points.

[0124] The above two papers give the conditions that RAAN and AoL need to satisfy when the satellite can cover the ground station in two cases. On this basis, the present embodiment specifically adopts the following geometric method:

[0125] (1) If φ < i - θ, then when the satellite can cover the ground station, the RAANs corresponding to the target ground station M, the northern boundary point N, and the southern boundary point S in the orbit are respectively:

[0126] Ω M = λ + G - arctan(tanu M ·cosi)

[0127] Ω NA = λ + G - △λ MN - arctan(tanu NA ·cosi), Ω ND = 2(λ + G) - Ω NA - π

[0128] Ω SA = λ + G + △λ MS - arctan(tanu SA ·cosi), Ω SD = 2(λ + G) - Ω SA - π

[0129] In the formula, Ω M , Ω NA , Ω ND , Ω SA , Ω SD respectively represent the target ground station M, the northern tangent point N A of the upward orbit and the coverage circle D , the northern tangent point N A of the downward orbit and the coverage circleD The corresponding RAAN in the orbit; △λ MN , △λ MS These represent the longitude differences between the target ground station M and its northern boundary point N and southern boundary point S, respectively; u M u NA u SA These represent the target ground station M, the point of tangency N between the uplink orbit and the northern side of the coverage circle, respectively. A The point of tangency S between the upward trajectory and the south side of the covering circle. A In the orbit, AoL corresponds to λ, which represents the geocentric longitude of the target ground station M, φ represents the geocentric latitude of the target ground station M, θ represents the geocentric angle of coverage, i represents the orbital inclination, and G represents the Greenwich sidereal hour angle.

[0130] The corresponding AoL:u of the target ground station M, the northern boundary point N, and the southern boundary point S in orbit M u NA u ND u SA u SD They are obtained from the following formulas respectively:

[0131] u M =arcsin(sinφ / sini)

[0132] u NA =arcsin(sinφ N / sini), u ND =π-u NA

[0133] u SA =arcsin(sinφ S / sini), u SD =π-u SA

[0134] In the formula, u M u NA u ND u SA u SD These represent the target ground station M, the point of tangency N between the uplink orbit and the northern side of the coverage circle, respectively. A The point of tangency N between the down-track trajectory and the north side of the covering circle. D The point of tangency S between the upward trajectory and the south side of the covering circle. A The point of tangency S between the downward trajectory and the south side of the covering circle D The corresponding AoL in the orbit; φ, φ N φ S Let M, N (north boundary point), and S (south boundary point) represent the latitudes of the target ground station, respectively; φ represents the latitude of the target ground station M; φ NN represents the point of tangency between the upward trajectory and the north side of the covering circle. A The latitude or descending orbit and the point of tangency N on the north side of the covering circle. D The latitudes are the same for both; φ S The point S represents the southern tangent point between the upward trajectory and the covering circle. A The latitude or descending orbit and the point of tangency S on the south side of the covering circle D The latitudes are the same for both.

[0135] The longitude difference Δλ between the target ground station M and the northern boundary point N and the southern boundary point S MN , △λ MS They are obtained from the following formulas respectively:

[0136] △λ MN =arccos[(cosθ-sinφ N ·sinφ) / (cosφ N ·cosφ)]

[0137] △λ MS =arccos[(cosθ-sinφ S ·sinφ) / (cosφ S ·cosφ)]

[0138] sinφ N =sinφ·cosθ+sinθ·cosi

[0139] sinφ S =sinφ·cosθ-sinθ·cosi

[0140] cosφ N =(1-sin 2 φ N ) 1 / 2

[0141] cosφ S =(1-sin 2 φ S ) 1 / 2

[0142] In the formula, △λ MN , △λ MS Represent the longitude difference between the target ground station M and the northern boundary point N and the southern boundary point S, respectively; θ represents the geocentric angle of coverage, i represents the orbital inclination angle, and φ and φ N φ S These represent the latitudes of the target ground station M, the northern boundary point N, and the southern boundary point S, respectively.

[0143] Considering the Earth's rotation, the actual RAAN coverage area △L that can be achieved is:

[0144] △L = △λ N+ △λ MN+ △λ MS -△λ S -(w e T d / (2π))|u N -u S |

[0145] △λ N = arcsin[(cot i)·(sin φ N / cos φ N )]

[0146] △λ S = arcsin[(cot i)·(sin φ S / cos φ S )]

[0147] In the formula, △λ N , △λ S respectively represent the longitude spans of the satellite flying from the equator to the north boundary point N and the south boundary point S along the orbit; △λ MN , △λ MS respectively represent the longitude differences between the target ground station M and the north boundary point N and the south boundary point S; u N , u S respectively represent the AoL corresponding to the north boundary point N and the south boundary point S in the orbit; w e represents the angular velocity of the Earth's rotation, T d represents the orbital period; i represents the orbital inclination, φ N , φ S respectively represent the latitudes of the north boundary point N and the south boundary point S; | | represents the absolute value.

[0148] It should be noted that in the specific implementation, when obtaining the RAAN range of the upward orbit by the geometric method, u N , u S are respectively taken as u NA , u SA , on the contrary, when obtaining the RAAN range of the downward orbit, u N , u S are respectively taken as u ND , u SD .

[0149] (2) If i - θ < φ < i + θ, then when the satellite can cover the ground station, the RAAN corresponding to the boundary point in the orbit is:

[0150] Ω SA = λ + G + △λ MS - arctan(tan uSA ·cosi)

[0151] Ω SD =λ+G-△λ MS -arctan(tanu SD ·cosi)-π

[0152] In the formula, Ω SA Ω SD They represent the tangent point S respectively. A Tangent point S D The corresponding RAAN in the orbit; u SA u SD They represent the tangent point S respectively. A Tangent point S D The corresponding AoL in the orbit; △λ MS λ represents the longitude difference between the target ground station M and the boundary point S; λ represents the geocentric longitude of the target ground station M; i represents the orbital inclination; and G represents the Greenwich Mean Time (GMT).

[0153] Considering Earth's rotation, the actual RAAN coverage area △L that can be achieved is:

[0154] △L=2(△λ MS -△λ S )+π-(w e T d / (2π))|u SD -u SA |

[0155] In the formula, △λ MS Δλ represents the longitude difference between the target ground station M and the boundary point S. S This represents the longitude span of a satellite traveling along its orbit from the equator to the boundary point S; u SA u SD They represent the tangent point S respectively. A Tangent point S D The corresponding AoL in the orbit; w e T represents the angular velocity of Earth's rotation. d represents the orbital period, and | represents the absolute value.

[0156] (3) When the satellite's RAAN is inside the boundary point of (1) or (2) above, and the following AoL condition is met, the satellite and the ground target are just visible:

[0157] Assuming the actual RAAN of the satellite is Ω, and the nadir trajectory intersects the coverage circle at points S1 and S2, solve the following system of equations to obtain the AoL of intersection points S1 and S2: u1, u2:

[0158] α = Ω + arctan(tanu·cosi)

[0159] δ = arcsin(sinu·sini)

[0160] cosθ=(cosφ)·(cosδ)·cos(|α M -α|)+(sinφ C )·(sinδ1)

[0161] In the formula, α represents the right ascension of S1 or S2, δ represents the declination of S1 or S2, φ represents the latitude of the target ground station, and α M Let represent the right ascension of the target ground station, and u represent the AoL to be solved, with two values ​​of u, u1 and u2. In practice, the intersection points S1 and S2 are substituted into this system of equations to obtain the AoL of intersection points S1 and S2: u1 and u2.

[0162] Except for boundary points, each RAAN value corresponds to two AoL values, and the visible range is closed. For example, Figure 2 The RAAN-AoL visibility range of satellites at different inclinations on a 400km orbit is plotted for ground stations with the same longitude and latitudes of 10°, 45°, and 85°, with a coverage half-angle of 30°. The shape enclosed by the solid line in the figure represents the current RAAN-AoL visibility range. If the combination of satellite RAAN-AoL values ​​falls exactly within the solid line, satellite-to-ground coverage is achieved.

[0163] In this embodiment, the orbit and orbital parameters of the approximately uniform constellation to be corrected are referred to as the nominal orbit and nominal orbital parameters; the actual orbits and orbital parameters of each satellite after correction are referred to as the working orbit and working orbital parameters. It is assumed that the constellation to be corrected consists of P orbital planes, with the subscript p representing the orbital plane number. The number of satellites in each orbital plane is S. p Each satellite is arranged according to 1~S p The sequential numbering. The nominal RAAN of orbital plane p is Ω. p The adjustment range of RAAN is [Ω]. min,p Ω max,p The nominal semi-major axis of a constellation is 'a', and the allowable adjustment range of the semi-major axis is [a]. min a max The inclination angle is i, the eccentricity is 0, and the mission duration is T. After time T, constellation reconfiguration requires configuration maintenance. Each satellite in the constellation carries a sensor, and the sensor's coverage half-angle is α. y y∈[1,Y], where Y is the number of different covering half-angles in the constellation. The surface-to-mass ratio of the windward side of each satellite in the constellation is S / m. z , z∈[1,Z], where Z is the number of different face-to-mass ratios in the constellation.

[0164] like Figure 7 and Figure 8 As shown, the specific steps of the method in this embodiment are as follows:

[0165] Step 1: Data Preprocessing

[0166] 1.1 Feature semi-major axis sampling: within the allowable adjustment range of the constellation semi-major axis [a min a max Several semi-major axis values ​​are selected at equal intervals within the [inner region] as the characteristic semi-major axis a r We take X characteristic semi-major axes, denoted as a. r,x , x∈[1,X].

[0167] 1.2 Calculation of the coverage angle:

[0168] 1.2.1 Let the semi-major axis of the feature be x=1.

[0169] 1.2.2 Let the index of the half-angle covered be y=1.

[0170] 1.2.3 By solving the following equation, the value located on the characteristic semi-major axis a is obtained. r,x Above, and carrying a coverage half-angle of α y The geocentric angle θ corresponding to the Earth coverage range of the satellite of the sensor x,y :

[0171] a r,x =(R e sinθ x,y ) / tanα y +R e cosθ x,y

[0172] Among them, a r,x α represents the characteristic semi-major axis of the satellite. y θ represents the half-angle of coverage corresponding to the satellite. x,y R represents the geocentric angle of the satellite's Earth coverage area. e R is the Earth's radius. e =6378.139km. The relationship between the half-angle of cover and the geocentric angle of cover is as follows: Figure 3 As shown.

[0173] 1.2.4 Let the index of the half-angle being covered be y = y + 1. If y ≤ Y, go to step 1.2.3; otherwise, go to step 1.2.5.

[0174] 1.2.5 Let the semi-major axis of the feature be x = x + 1. If x ≤ X, go to step 1.2.2; otherwise, go to step 1.3.

[0175] 1.3 Calculation of satellite orbital drift rate:

[0176] 1.3.1 Let the characteristic semi-major axis number x=1.

[0177] 1.3.2 Let the surface-to-mass ratio number z = 1.

[0178] 1.3.3 Using simulation methods, within the task duration T, the region located on the characteristic semi-major axis a is simulated. r,x Above, and the surface-to-mass ratio is S / m z The satellite's trajectory is recorded, along with data on the semi-major axis variation and the RAAN variation ΔΩ. x,z .

[0179] 1.3.4 Using a quadratic fitting method, the semi-major axis variation data is fitted to the following formula:

[0180] a x,z (t)=a r,x +k1 x,z t 2 +k2 x,z t, t∈[0,T]

[0181] In the formula, k1 is the semi-major axis decay parameter of the quadratic term, k2 is the semi-major axis decay parameter of the linear term, and a(t) represents the semi-major axis at time t; the subscript x represents the characteristic semi-major axis index, and the subscript z represents the surface-to-mass ratio index. Wherein, k2 <k1<0。

[0182] The reason for using a quadratic fitting method to fit the satellite's semi-major axis decay data is that, under the influence of atmospheric drag, the decay of the semi-major axis of a low-Earth orbit satellite is approximately linear in the short term, but the decay rate increases over a longer period. Quadratic fitting can better simulate the medium- to long-term decay of the semi-major axis. After fitting, the semi-major axis during the mission period can be described by the semi-major axis decay parameters and the initial semi-major axis, avoiding the introduction of empirical or exponential models and facilitating subsequent numerical calculations.

[0183] 1.3.5 Calculate the linearized RAAN rate of change γ according to the following formula. Ω,x,z :

[0184] γ Ω,x,z =△Ω x,z / T,x∈[1,X],z∈[1,Z]

[0185] In the formula, γ Ω,x,z This represents the value located on the characteristic semi-major axis a within the simulation task duration T. r,x Above, and the surface-to-mass ratio is S / m z The linearized RAAN rate of change of the satellite, ΔΩ x,z This represents the value located on the characteristic semi-major axis a within the simulation task duration T. r,x Above, and the surface-to-mass ratio is S / m z The RAAN variation of the satellite.

[0186] 1.3.6 Let the surface-to-mass ratio number z = z + 1. If z ≤ Z, go to step 1.3.3; otherwise, go to step 1.3.7.

[0187] 1.3.7 Let the semi-major axis of the feature be x = x + 1. If x ≤ X, go to step 1.3.2; otherwise, go to step 2.

[0188] Step 2: Calculation of Star-to-Ground Visibility Conditions

[0189] 2.1 Ground Target Grid Generation: Assuming the ground targets include n1 regional targets and n2 point targets, the n1 regional targets are divided into several point targets using a grid method, and then merged with the n2 point targets to form N ground point targets, denoted by the subscript n. The longitude and latitude of the ground point target n are represented by (λ... n φ n ).

[0190] 2.2 Calculation of satellite-to-ground visibility conditions:

[0191] 2.2.1 Let the ground target number n=1.

[0192] 2.2.2 Calculation of RAAN visibility conditions for ground target n at initial time:

[0193] Using geometric methods (1) or (2), the initial time t=0 is calculated by iterating through the semi-major axis a. r,x On the orbital plane x∈[1,X], using a covering half-angle of α y For a satellite with sensor y∈[1,Y], when the satellite covers a ground point target n, the RAAN boundary condition Ω of the satellite must satisfy. A Ω D , center point Ω M And the visible range △L of RAAN, where the geocentric angle is calculated by step 1.2; and the boundary condition Ω under either (1) or (2) orbital conditions. A Ω D Use the left boundary Ω uniformly L Ω R This indicates that the initial RAAN visibility condition set for ground target n is obtained:

[0194] A n Ω,t=0 =(Ω n L,x,y Ω n M,x,y Ω n R,x,y , △L n x,y ), x∈[1,X], y∈[1,Y], n∈[1,N]

[0195] Among them, each Ω L Ω R Ω M Combination of △L (Ω) L Ω R Ω M , △L) constitute a set of initial RAAN visibility opportunities; the subscript x indicates the feature semi-major axis number, and the subscript z indicates the surface-to-mass ratio number.

[0196] Meanwhile, the initial RAAN visibility condition set A for ground target n n Ω,t=0 The CCP contains the initial RAAN visibility chance of group H, where H∈[XY, 2XY].

[0197] 2.2.3 Discretization of the initial RAAN visibility condition for ground target n:

[0198] Traverse A n Ω,t=0 For each set of initial RAAN visibility opportunities, perform the following calculations for each set of initial RAAN visibility opportunities: at the boundary point Ω L Ω R Selecting points including boundary points Ω L Ω R and the center point Ω M Several feature points are used to obtain a set of initial RAAN visible feature points.

[0199] Finally, we obtain the discretized initial RAAN visible feature point set R as shown in the following formula. n Ω,t=0 And the initial RAAN visible feature range set △L:

[0200]

[0201] Among them, R n Ω,t=0 In △L, each row represents a set of initial RAAN visibility opportunities, and the two are in one-to-one correspondence, for a total of H sets.

[0202] 2.2.4 Calculation of the initial RAAN-AoL visibility condition feature set for ground target n:

[0203] For R n Ω,t=0 In the initial RAAN visibility opportunities of groups 1 to H, for each visible RAAN value, the AoL boundary points u1 and u2 corresponding to each RAAN value are calculated using the geometric method (3), forming the RAAN-AoL visibility condition feature set R at the initial time. n Ω-u,t=0 :

[0204]

[0205] In the formula, Ω L,1 (u1 L,1 u2 l,1 ) indicates the visible RAAN value Ω L,1 The corresponding AoL boundary points u1, u2, and so on.

[0206] 2.2.5 Initial RAAN-AoL Visibility Condition Fitting and Reconstruction of Ground Target n:

[0207] Using R n Ω,t=0 Polynomial fitting is performed on the feature points in the matrix to reconstruct the feature point set R within the visible range of RAAN at the initial time. n Ω,t=0 The AoL visibility range corresponding to points other than those in the initial time step forms a complete set of RAAN-AoL visibility conditions:

[0208] f Ω-u n (t=0)=f n (Ω t=0 (u1, u2) h h∈[1,H]

[0209] f Ω-u n (t=0) There are H rows in total, each row representing a set of initial RAAN-AoL visibility opportunities. Each set of RAAN-AoL visibility opportunities is represented by the function f. n (Ω t=0 (u1, u2) describes the visible polygon, polygon shape reference. Figure 2 .

[0210] The reason for using polynomial fitting instead of traversal calculation is that the computational cost of calculating the visibility conditions of RAAN is relatively small, but accurately obtaining the visibility conditions of AoL corresponding to each RAAN and drawing the shape of the visible range requires a large amount of computation. Therefore, by sacrificing some feasible solutions as a cost to reduce the computational burden, the complete shape can be preserved under sufficient computing power.

[0211] Furthermore, it is preferable to use polynomials of degree 6 or higher for fitting, thereby further balancing computational accuracy and computational cost.

[0212] 2.2.7 Let the ground target number n = n + 1. If n ≤ N, then go to step 2.2.2; otherwise, go to step 3.

[0213] Step 3: Create the query table

[0214] 3.1 Create Query Table 1:

[0215] Based on the calculation results of step 1, a lookup table 1 is created, which contains X characteristic semi-major axes a. r Z surface mass ratio S / m z The corresponding satellite's semi-major axis attenuation parameters k1, k2, and linearized RAAN rate of change γ Ω Refer to the table for the content and structure of Table 1.

[0216] Table 1. Schematic diagram of the content and structure of query table 1

[0217]

[0218] 3.2 Create Query Table 2:

[0219] Based on the calculation results of steps 1 and 2, query table 2 is created, which includes:

[0220] a. X characteristic semi-major axes r The geocentric angle θ corresponding to each of the Y half-angles of coverage;

[0221] b. X characteristic semi-major axes a r The visible feature range ΔL of RAAN at the initial time corresponding to Y covering half-angle α and N ground targets;

[0222] c. Initial RAAN-AoL visibility condition f Ω-u (t=0).

[0223] Refer to Table 2 for the content and structure of the query table.

[0224] Table 2. Schematic diagram of the content and structure of Table 2

[0225]

[0226] The purpose of creating a lookup table is that the data in the lookup table is data that the algorithm needs to use frequently in subsequent steps, and pre-calculation can save resources and time in subsequent calculations.

[0227] Step 4: Modeling the working semi-major axis and latitude argument correction values

[0228] Traverse orbital planes 1 through P. For each orbital plane p, perform semi-major axis and AoL corrections one by one. Set the initial working semi-major axis a of each satellite within orbital plane p. s , s∈[1,S p ] and initial AoL value u s , s∈[2,S p As the variable to be corrected, an optimization model for the working trajectory parameters is established.

[0229] The nominal initial RAAN value of orbital plane p is Ω0, and the nominal initial semi-major axis is a. This is implemented through the following steps:

[0230] 4.1 Let the orbital plane number p=1.

[0231] 4.2 Let the ground target number n=1.

[0232] 4.3 Set the satellite serial number s=1.

[0233] 4.4 Calculation of the maximum concentrated visible period of ground target n for orbital plane p satellite s:

[0234] The significance of the maximum concentrated visibility period is: assuming there are enough satellites in the orbital plane to ensure that the AoL visibility condition is always met, then if the RAAN visibility condition is also met, there will inevitably be satellites in that orbit capable of providing coverage. For example, Figure 4 The graph plots the relationship between RAAN visibility conditions at a ground station at a latitude of 30° and different orbits over time, retaining only the boundaries of the visible range. As time increases, the boundary points of the RAAN visibility conditions shift with the Earth's rotation rate. The band between the two blue diagonal lines represents the visible RAAN range for a 30° inclination orbit. The band between the red diagonal lines represents the visible RAAN range for a 50° inclination orbit. The dashed line represents the downhill orbit, and the solid line represents the uphill orbit. Satellite 1 is a 30° inclination satellite; when its RAAN trajectory (blue vertical line) intersects with the blue band, Satellite 1 has a chance of being visible to the target ground station. The period from t1 to t2 is the period of maximum concentrated visibility for this orbit. Satellite 2 is a 50° inclination satellite; when its RAAN trajectory (red vertical line) intersects with the red band, Satellite 2 has a chance of being visible to the target ground station. It has a downhill visibility opportunity during the period from t3 to t4 and an uphill visibility opportunity during the period from t5 to t6. t3 to t4 and t5 to t6 are the periods of maximum concentrated visibility for this orbit.

[0235] Visibility between satellites and the ground is only possible during the period of maximum concentrated visibility. The drift differences between heterogeneous satellites will cause continuous changes in their relative positions during this period, thus affecting coverage performance. Therefore, designing the positions and relative relationships of each satellite within each period of maximum concentrated visibility is crucial for improving constellation coverage performance.

[0236] 4.4.1 In the data of ground target n in Table 2, find the sensor half-angle that is closest to the half-angle of the current satellite s, and whose characteristic semi-major axis is closest to the initial working semi-major axis a of the current satellite s. s The closest case is recorded, including its geocentric angle θ and the initial RAAN visible boundary point Ω. R And the initial visible range of RAAN △L.

[0237] Among them, Ω R For O0 or O u0 O d0, depending on the relationship between i-θ and φ. Specifically: O0 is the RAAN value corresponding to the right boundary point in the orbit when i-θ < φ < i+θ, that is, S A the RAAN value corresponding in the orbit; O u0 、O d0 are respectively the RAAN values corresponding to the right boundary points of the ascending orbit and the descending orbit in the orbit when φ < i-θ, that is, S A 、N D the RAAN values corresponding in the orbit.

[0238] Among them, △L is △L1 or △L2, depending on the relationship between i-θ and φ. Specifically: △L1 is the RAAN visible range when i-θ < φ < i+θ, and △L2 is the RAAN visible range when φ < i-θ.

[0239] Note: The initial working semi-major axis a of the current satellite s s will be assigned an initial value by the solution algorithm, and the same applies hereinafter.

[0240] 4.4.2 In the query table 1, find the case with the surface mass ratio closest to that of the current satellite s and the characteristic semi-major axis closest to the initial working semi-major axis a of the current satellite s s and record its RAAN change rate as γ Ω value.

[0241] 4.4.3 The maximum concentrated visible period Γ of the ground target n for the satellite s on the orbital plane p Ω0,p,s :

[0242]

[0243] In the formula, Γ Ω0,p,s represents the maximum concentrated visible period of the satellite s on the orbital plane p for the ground target n; A1 represents the initial offset coefficient, A2 represents the long-term drift proportional coefficient, t Γ1 or t Γ2 represents the duration of a single maximum concentrated visible period, t B +t Γ1 or t B +t Γ2 represents the repetition period of the maximum concentrated visible period, τ1, τ2 or τ3 represents the long-term drift time, φ n represents the latitude of the ground target n, θ represents the earth-centered angle covered by the satellite s, i represents the orbital inclination of the orbital plane p; r is a natural number.

[0244] Among them, the specific values of each parameter are:

[0245] A1 = Ω0tanγ, A2 = tanγ Ω / (tanγ Ω +tanγ)

[0246] t B1 =(2π-△L1) / w e , t B2 =(2π-△L2) / w e γ=arctan(1 / w e )

[0247] t Γ1 =△L1tanγ,t Γ2 =△L2tanγ

[0248] τ1=(2π-O0)tanγ,τ2=(2π-O u0 )tanγ,τ3=(2π-O d0 )tanγ.

[0249] In the formula, Ω0 represents the nominal initial RAAN value of the orbital surface p, and γ Ω w represents the rate of change of RAAN of satellite s on orbital plane p. e The value represents the Earth's rotational angular velocity, and γ is an auxiliary parameter representing the arctangent function value of the Earth's rotational angular velocity.

[0250] 4.5 Calculation of the visible time of orbital plane p satellite s relative to ground target n:

[0251] Suppose that within time T, there are J maximum concentrated visible periods for ground target n relative to orbital plane satellite s, denoted by subscript j, where j∈[1, J]. Perform the following steps:

[0252] 4.5.1 Let the maximum concentrated visible time period number j=1.

[0253] 4.5.2 Let the current time t be the initial time of the current maximum concentrated visible time period j.

[0254] 4.5.3 Calculation of satellite orbital parameters:

[0255] In lookup table 1, find the satellite whose surface-to-mass ratio is closest to that of the current satellite s, and whose characteristic semi-major axis is the same as the initial working semi-major axis a of the current satellite s. s The closest case is recorded, showing its orbital decay parameters k1 and k2, as well as the RAAN rate of change γ. Ω Calculate the RAAN and AoL values ​​of orbital plane p satellite s at the current time t, denoted as Ω. t,p,s u t,p,s :

[0256] Ω t,p,s =Ω0+γ Ω ·t

[0257] u t,p,s =us -[(7λ'+3n e )·k1·t 3 ] / (6a)-[(7λ'+3n e )·k2·t 2 ] / (4a)-[(7λ'+3n e )·(a s -a)·t] / (2a)

[0258] λ'=[(3n e J2R e 2 ) / (2a 2 )]·(4cos 2 i-1)

[0259] In the formula, Ω t,p,s Ω represents the RAAN value of satellite s at orbital plane p at the current time t, Ω0 represents the nominal initial RAAN value of the satellite, and γ Ω Indicates the RAAN rate of change of the satellite; u t,p,s U represents the AoL value of satellite s in orbital plane p at the current time t. s Let λ' represent the initial AoL of the satellite, and λ' represent the long-term rate of change of the mean anomaly angle, n e This represents the orbital angular velocity corresponding to the nominal semi-major axis, k1 represents the satellite's quadratic semi-major axis attenuation parameter, k2 represents the satellite's linear semi-major axis attenuation parameter, and a... s This represents the initial operational semi-major axis of the satellite, i.e., the semi-major axis value currently being optimized. 'a' represents the satellite's orbital semi-major axis before correction, i.e., the nominal semi-major axis. J2 = 1082.63 × 10 -6 .

[0260] 4.5.4 Calculation of RAAN-AoL visibility conditions of satellite s to ground target n at time t:

[0261] In the data of ground target n in Table 2, find the sensor half-angle that is closest to the half-angle of the current satellite s, and whose characteristic semi-major axis is closest to the initial working semi-major axis a of the current satellite s. s The closest case is to record the initial RAAN-AoL visibility condition f. Ω-u n (t=0)=f n (Ω t=0 (u1, u2). And f Ω-u n (t=0) is mapped to the current time t, forming the RAAN-AoL visibility condition of the current satellite s for the ground target n at time t:

[0262] f Ω-u n (t)=fn (Ω t=0 +w e t, u1, u2).

[0263] The reason why the RAAN-AoL visibility condition of satellite s relative to ground target n at time t can be calculated using the above formula is that, compared with the initial RAAN-AoL visibility condition, the RAAN condition at time t shifts with the Earth's rotation, while the AoL visibility condition remains unchanged. For example, Figure 5 This figure shows the visibility range of a 500km, 30° inclination satellite to a ground station at 30° latitude. The visibility range in the figure is plotted discretely with one orbital period, omitting numerous similar visibility range curves between adjacent visible ranges. Due to the Earth's rotation, the visibility range continuously shifts to the right.

[0264] 4.5.5 Visibility assessment of satellite s relative to ground target n at time t:

[0265] If the RAAN value of satellite s at the current time t is Ω t,p,s (Calculated according to step 4.5.3) Within the RAAN visibility condition (calculated according to step 4.5.4), and the AoL value u t,p,s (Calculated according to step 4.5.3) If the AoL visibility condition (u1, u2) (calculated according to step 4.5.4) is within the range, then it is denoted as the visibility of the current satellite s to the ground target n at time t, and denoted as AccFlag. n p,s If (t) = 1, otherwise mark it as invisible, and denote it as AccFlag. n p,s (t)=0.

[0266] 4.5.6 Let the current time t = t + Δt, where Δt is the step size. If t is less than or equal to the end time of the current maximum concentrated visible time period j, then proceed to step 4.5.3; otherwise, proceed to step 4.5.6.

[0267] 4.5.6 Organize the visibility indicators of orbital plane p satellite s for ground target n at each time point within the j-th maximum concentrated visibility period, forming a list of visible times, denoted as l. n acc,p,s (j)=[t1,t2,…t n ],t∈Γ n Ω0,p,s (j).

[0268] 4.5.7 Let the index of the most concentrated visible time period be j = j + 1. If j ≤ J, go to step 4.5.2; otherwise, go to step 4.5.8.

[0269] 4.5.8 Arrange the list of J visible moments of orbital plane p satellite s relative to ground target n in chronological order and integrate them to form a list of visible moments of orbital plane p satellite s relative to ground target n within the mission period T, denoted as l. n acc,p,s =[t1,t2,…t n ],t∈Γ n Ω0,p,s .

[0270] 4.6 Let the satellite number s = s + 1. If s ≤ S p If yes, proceed to step 4.4; otherwise, proceed to step 4.7.

[0271] 4.7 Calculation of the list of visible time periods for all satellites on orbital plane p relative to ground target n:

[0272] List the S visible times of all satellites on orbital plane p relative to ground target n during mission period T. n acc,p,s The data is then merged and arranged chronologically to form a list l of the visible times of all satellites on orbital plane p relative to ground target n within mission period T. n acc,p :

[0273]

[0274] Then l n acc,p The consecutive moments are merged and arranged according to the start time of the visible time period to form a visible time period list L. n acc,p :

[0275]

[0276] Each row represents a visible time period, assuming there are a total of Q. n p The visible time period is divided into three columns: the first column is the start time of the visible time period, the second column is the end time of the visible time period, and the third column is the name of the visible satellite, which may include one or more satellites.

[0277] 4.8 Calculation of coverage index of all satellites on orbital plane p for ground target n:

[0278] Within mission period T, the coverage index of all satellites on orbital plane p for ground target n is expressed as ε. n p Based on the focus of the task, ε n p Indicators such as average visibility duration or average revisit time can be used.

[0279] Among them, the average visibility index τn p Specifically:

[0280]

[0281] In the formula, q represents the sequence number of the visible time period, Q n p EndT represents the total number of visible time periods. q begT represents the end time of the q-th visible time period. q This represents the start time of the qth visible time period.

[0282] Among them, the average revisit time index △t n p Specifically:

[0283]

[0284] In the formula, q represents the sequence number of the visible time period, Q n p EndT represents the total number of visible time periods. q-1 begT represents the end time of the (q-1)th visible time period. q This represents the start time of the qth visible time period.

[0285] Within mission period T, the consistency of coverage indicators for ground target n by all satellites on orbital plane p is expressed as follows:

[0286]

[0287] In the formula, ε n’ p,q This represents the coverage metric value for each visible time period. Specifically, when the coverage metric is selected as the visible duration, When the coverage metric is selected as the revisit time, Var n p This indicates consistency in coverage metrics.

[0288] 4.9 Let the ground target number n = n + 1. If n ≤ N, go to step 4.3; otherwise, go to step 4.10.

[0289] 4.10 Calculation of the coverage index of orbital plane p for all ground targets:

[0290] The average coverage index is: ε p =(1 / N)·∑ N n=1 (w·ε n p );

[0291] The average coverage consistency index is: Varp =(1 / N)·∑ N n=1 (w·Var n p );

[0292] Where w represents the weighting coefficient.

[0293] 4.11 Optimization modeling of working trajectory parameters on orbital plane p:

[0294] The optimization model for the working orbital parameters of all satellites within orbital plane p is as follows: the optimization variable is the initial semi-major axis a of each satellite. s , s∈[1,S p ] and initial AoL, s∈[2,S p ], 2S in total p -1 satellites, with satellite number 1 having an AoL of 0. The optimization objective is the average coverage index ε. p Consistency index Var with average coverage p Minimum, constrained to the semi-major axis adjustment range, is expressed as:

[0295]

[0296] 4.12 Let the orbital plane number p = p + 1. If p ≤ P, go to step 4.2; otherwise, go to step 5.

[0297] Step 5: Solving for the correction values ​​of the working semi-major axis and latitude argument

[0298] Traverse orbital planes 1-P, and solve for the semi-major axis and AoL for each orbital plane p. Calculate the initial working semi-major axis a for each satellite within orbital plane p. s , s∈[1,S p ] and initial AoL, s∈[2,S p As the variable to be corrected, it is solved using analytical methods or optimization algorithms.

[0299] 5.1 Let the orbital plane number p=1.

[0300] 5.2 Solving for the working trajectory parameters of the orbital surface p using analytical methods:

[0301] 5.2.1 Calculate the initial working semi-major axis of all satellites on orbital plane p:

[0302] For each satellite in orbital plane p, perform the following operations: In lookup table 1, find the case whose surface-to-mass ratio is closest to that of the current satellite s, and whose characteristic semi-major axis is closest to the nominal semi-major axis a of the current satellite s, and record its orbital attenuation parameters k1 and k2. Calculate the initial working semi-major axis of each satellite in orbital plane p:

[0303] as =a-(k1 s T 2 ) / 3-(k2 s T) / 2, s∈[1,S] p ]

[0304] In the formula, a s Let represent the initial operating semi-major axis of satellite s, 'a' represent the nominal semi-major axis of satellite s, 'k1' represent the quadratic semi-major axis attenuation parameter of satellite s, 'k2' represent the linear semi-major axis attenuation parameter of satellite s, 'T' represent the mission period of satellite s, and 'S' represent the satellite's mission period. p This represents the total number of satellites on orbital plane p where satellite s is located.

[0305] 5.2.2 Calculate the nominal AoL for all satellites on orbital plane p:

[0306] u s =[2(s-1)π] / S p , s∈[1,S p ]

[0307] In the formula, u s 'Indicates the nominal AoL of satellite s.

[0308] 5.2.3 Calculate the AoL offset for all satellites on orbital plane p:

[0309] When the initial semi-major axis of the satellite is a s At that time, within period T, the maximum AoL offset Δu of satellite s max,s It is obtained through the following formula:

[0310] △u max =[(7λ'+3n e ) / (12a)]·[-2k1·T m 3 -3k2·T m 2 +(2k1·T 2 +3k2·T)·T m ]

[0311] λ'=[(3n e ·J2·R e 2 ) / (2a 2 )]·(4cos 2 i-1)

[0312] T m =-k2 / (2k1)-[3k2 2 +2k1·(2k1·T 2 +3k2·T)] 1 / 2 / (2×31 / 2 ×k1)

[0313] In the formula, △u max λ' represents the maximum AoL offset of the satellite, λ' represents the long-term rate of change of the mean anomaly angle, and n e Let represent the orbital angular velocity, 'a' represent the satellite's nominal semi-major axis, 'k1' represent the satellite's quadratic semi-major axis attenuation parameter, 'k2' represent the satellite's linear semi-major axis attenuation parameter, and 'T' represent the orbital angular velocity. m J2 represents the time required to reach the maximum AoL offset, and R represents the Earth's oblateness perturbation coefficient. e denoted by , where i represents the Earth's radius and i represents the orbital inclination.

[0314] 5.2.4 Perform the following operations on each satellite in orbital plane p: In the data from Table 2, find the case where the sensor's half-angle is closest to the half-angle of the current satellite s, and the characteristic semi-major axis is closest to the nominal semi-major axis a of the current satellite s. Record the coverage geocentric angle θ. s .

[0315] 5.2.5 Calculate the initial AoL for all satellites on orbital plane p:

[0316] u s =u s '-△u max,s / 2, s∈[1,S] p ]

[0317] In the formula, u s Denotes the initial AoL, u of satellite s s 'Indicates the nominal AoL of satellite s, Δu max,s This represents the maximum AoL offset of satellite s.

[0318] 5.2.6 If all satellites on orbital plane p satisfy:

[0319] △u max,s ≤2θ s , s∈[1,S p ]

[0320] In the formula, △u max,s θ represents the maximum AoL offset of satellite s. s This represents the geocentric angle of satellite s' coverage.

[0321] If the result is positive, proceed to step 5.3; otherwise, proceed to step 5.4.

[0322] The reason for using the analytical method to correct the working orbit parameters is as follows: Using the virtual satellite at the nominal AoL position of each satellite as the reference point for each satellite, and setting the relative AoL change of the satellite within a time interval to 0, the semi-major axis of the initial working orbit of each satellite is obtained. Then, within period t, the satellite approximately performs symmetrical motion relative to its respective reference point. Figure 6 The figure shows the relative motion trajectories of three satellites, S1, S2, and S3, within one period T. When the nadir trajectory just passes over the ground target, the satellite's AoL deviation is within ±θ, and the satellite can still complete coverage within the specified time without affecting the drift index. When the nadir trajectory crosses the ground coverage circle but does not pass over the ground target, the redundant deviation is less than ±θ. Therefore, when the maximum AoL offset Δu of satellite s... max,s ≤2θ s When this is the case, the above analytical method can be used to calculate the approximate optimal solution of the working orbit parameters of each satellite, thereby reducing the computational complexity.

[0323] Using this method, the semi-major axis of each satellite needs to be adjusted to a at the end of each mission cycle T. s .

[0324] 5.3 Based on steps 4.1 to 4.10, calculate the coverage index of orbital plane p within the mission period T. The initial working semi-major axis and initial AoL for each satellite are taken from the values ​​calculated in step 5.2. If the coverage index of the orbital plane meets expectations, the initial working semi-major axis obtained in step 5.2.1 and the initial AoL obtained in step 5.2.5 are used as the working semi-major axis correction value and the latitude argument AoL correction value, respectively, and the process proceeds to step 5.5; otherwise, proceed to step 5.4.

[0325] Among them, "coverage indicators meeting expectations" means that the coverage indicators meet the coverage performance requirements of the task design phase.

[0326] 5.4 Solving for the working trajectory parameters of orbital surface p using an optimization algorithm:

[0327] The working orbit parameter optimization model for orbital plane p in step 4.11 is solved using intelligent optimization algorithms such as genetic algorithms. The initial working semi-major axis a of each satellite is then determined. s , s∈[1,S p ] and initial AoL value u s , s∈[2,S p ], 2S in total p -1 parameters are used as the working orbital parameters for orbital plane p, ultimately yielding the working semi-major axis correction value and latitude argument (AoL) correction value for each satellite. In each iteration, for each individual satellite, according to steps 4.1 to 4.10, the coverage index and coverage index consistency index of orbital plane p within the mission period T are calculated.

[0328] The reason for using optimization algorithms to correct working trajectory parameters is that when the conditions for using analytical methods are not met, or when the coverage index obtained by using analytical methods is not ideal, all 2S... pUsing one working trajectory parameter as an optimization variable and employing an intelligent optimization algorithm can greatly increase the solution space, thereby obtaining a higher-quality optimal solution and improving the coverage index.

[0329] Using this method, it is necessary to recalculate and reconstruct the working orbit parameters of each satellite in step 5.4 at the end of each mission cycle T.

[0330] 5.5 Let the orbital plane number p = p + 1. If p ≤ P, go to step 5.2; otherwise, go to step 6.

[0331] Step 6: Correction of right ascension at the ascending node of the working point

[0332] The reason for constellation RAAN correction is to ensure the stability of relative AoL and transit time within the plane, different semi-major axes are designed for satellites in the same orbit. Under perturbation, the differentiated semi-major axes will cause the RAAN of each satellite in the same orbit to disperse. Therefore, the actual working RAAN is corrected based on the nominal RAAN of each orbit to counteract the RAAN drift caused by the semi-major axis.

[0333] Calculate the RAAN correction values ​​for orbital planes 1 to P using the following steps:

[0334] 6.1 Let the orbital plane number p=1.

[0335] 6.2 Set satellite serial number s=1.

[0336] 6.3 In lookup table 1, find the satellite whose surface-to-mass ratio is closest to that of the current satellite s, and whose characteristic semi-major axis is the same as the optimized initial working semi-major axis (working semi-major axis correction value) of the current satellite s. s For the closest case, record its orbital attenuation parameters k1 and k2. Calculate the RAAN change ΔΩ of satellite s at orbital plane p during mission period T. p,s :

[0337] △Ω p,s =-(7Ω'·k1·T 3 ) / (6a)-(7Ω'·k2·T 2 ) / (4a)-[(7Ω'·T) / (2a)]·(a s -a)

[0338] Ω'=-[(3n e ·J2·R e 2 ) / (2a 2 )]·cosi

[0339] In the formula, △Ω p,sLet Ω' represent the change in RAAN of satellite s at orbital plane p during the mission period, Ω' represent the long-term rate of change of RAAN, k1 represent the quadratic term semi-major axis attenuation parameter, k2 represent the linear term semi-major axis attenuation parameter, T represent the mission period, and a represent the nominal semi-major axis of the satellite. s This represents the optimized initial working semi-major axis, i.e., the working semi-major axis correction value, n. e J2 represents the orbital angular velocity, and R represents the Earth's oblateness perturbation coefficient. e denoted by , where i represents the Earth's radius and i represents the orbital inclination.

[0340] 6.4 Let the satellite number s = s + 1. If s ≤ S p If yes, proceed to step 6.3; otherwise, proceed to step 6.5.

[0341] 6.5 Calculate the average RAAN change of all satellites at orbital plane p within mission period T:

[0342]

[0343] In the formula, △Ω p ΔΩ represents the average change in RAAN for all satellites on orbital plane p over mission period T. p,s This represents the change in RAAN of satellite s on orbital plane p during mission period T.

[0344] 6.6 Update the nominal RAAN of orbital plane p: Ω p =Ω0-△Ω p Let the orbital plane number p = p + 1. If p ≤ P, go to step 6.2; otherwise, go to step 6.7.

[0345] 6.7 Let the orbital plane number p=1.

[0346] 6.8 Set satellite serial number s=1.

[0347] 6.9 Update the working RAAN parameters of the orbital plane p-satellite s: Ω p,s =Ω p -△Ω p,s .

[0348] 6.10 Let the satellite number s = s + 1. If s ≤ S p Proceed to step 6.9, otherwise proceed to step 6.11.

[0349] 6.11 Let the orbital plane number p = p + 1. If p ≤ P, go to step 6.8; otherwise, end the process.

[0350] Example 2

[0351] This embodiment verifies the orbital constellation through numerical simulation under conditions of an orbital altitude of 500km and an orbital inclination of 30°. Compared to the traditional Walker configuration, the optimized configuration obtained in this embodiment reduces non-uniform coverage by approximately 17%, and significantly increases the consistency of coverage indicators, such as revisit time, thus resolving the problem of coverage indicator drift in heterogeneous constellations. To achieve this goal, the control frequency of the traditional configuration needs to be increased by 3 to 6 times, or the number of satellites needs to be increased by approximately 30 to 40%.

[0352] The above specific embodiments are used to explain and illustrate the present invention, but not to limit the present invention. Any modifications and changes made to the present invention within the spirit and scope of the claims shall fall within the protection scope of the present invention.

Claims

1. A heterogeneous satellite constellation configuration modification method based on the stability of the coverage performance to the ground, characterized in that, The method comprises the following steps: Step S1: constructing a first query table and a second query table; Step S2: for each orbital plane, using a working orbit parameter optimization method, combining the first query table and the second query table, calculating the dynamic visibility condition and the visibility list, and then correcting the working semi-major axis and the AoL of each satellite on the orbital plane; In the step S2, the working orbit parameter optimization method comprises: Step S2.1: for each satellite on the orbital plane, according to the surface quality ratio, the coverage half-angle and the nominal semi-major axis of the satellite, the initial working semi-major axis, the initial AoL, the coverage central angle and the maximum AoL offset of the satellite are obtained by an analytical method in combination with the first query table and the second query table; Step S2.2: if the maximum AoL offset of each satellite on the orbital plane is less than or equal to twice the coverage central angle, step S2.3 is entered, otherwise, step S2.5 is entered; Step S2.3: according to the initial working semi-major axis and the initial AoL of all satellites on the orbital plane, the dynamic visibility condition and the visibility list are calculated by an analytical method to obtain the coverage index of the orbital plane within a mission period; The step S2.3 comprises: Step S2.3.1: for each satellite on the orbital plane, according to the surface quality ratio, the coverage half-angle and the initial working semi-major axis of the satellite, the RAAN change rate, the coverage central angle, the initial time RAAN-AoL visibility condition and the initial RAAN visibility feature range of each ground target of the satellite are obtained from the first query table and the second query table, and then the maximum concentrated visible period of the satellite to each ground target is obtained; Step S2.3.2: for each satellite on the orbital plane, according to the maximum concentrated visible period of the satellite to each ground target and the initial AoL, the visible time list of the satellite to each ground target is obtained; Step S2.3.3: merging the visible time lists of all satellites on the orbital plane to the same ground target, the visible period list of all satellites on the orbital plane to the ground target is obtained; Step S2.3.4: according to the visible period list of all satellites on the orbital plane to each ground target, the coverage index and the coverage index consistency of all satellites on the orbital plane to each ground target are calculated; Step S2.3.5: the average coverage index and the average coverage index consistency of all satellites on the orbital plane to all ground targets are obtained, and the coverage index of the orbital plane is obtained; Step S2.4: if the coverage index of the orbital plane meets the preset condition, the initial working semi-major axis and the initial AoL obtained in step S2.1 are taken as the working semi-major axis correction value and the latitude angle AoL correction value of the satellite respectively; otherwise, step S2.5 is entered; Step S2.5: using an optimization algorithm, taking the initial working semi-major axis and the initial AoL of all satellites on the orbital plane as optimization variables, and taking the minimum coverage index of the orbital plane within a mission period as an optimization target, in the optimization process, the coverage index of the orbital plane is obtained according to step S2.3, and finally the working semi-major axis correction value and the latitude angle AoL correction value of each satellite on the orbital plane are obtained; Step S3: For each orbital plane, according to the semi-major axis correction value of each satellite on the orbital plane, combined with the first query table, the RAAN correction value is calculated respectively; Step S4: The semi-major axis correction value, the AoL correction value and the RAAN correction value of each satellite on all orbital planes are taken as the configuration parameters of the heterogeneous satellite constellation.

2. The method according to claim 1, wherein the method is characterized by: The first query table contains semi-major axis decay parameters and RAAN change rates of satellites with different face-to-area ratios at each characteristic semi-major axis; the second query table contains coverage central angles of satellites with different coverage half-angles at each characteristic semi-major axis, and initial time RAAN-AoL visibility conditions and initial RAAN visibility feature ranges of each ground target.

3. The method according to claim 1, wherein the method is characterized by: The step S2.3.2 comprises: According to the face-to-area ratio, the coverage half-angle and the initial semi-major axis of the satellite, the orbital decay coefficient, the RAAN change rate and the initial time RAAN-AoL visibility condition of the satellite to the ground target are obtained from the first query table and the second query table; The RAAN value and the AoL value of the satellite at the current time are obtained according to the orbital decay coefficient, the RAAN change rate and the initial AoL of the satellite; the initial time RAAN-AoL visibility condition of the satellite to the ground target is mapped to the current time to obtain the RAAN-AoL visibility condition of the satellite to the ground target at the current time; the visibility is judged according to the RAAN value, the AoL value and the RAAN-AoL visibility condition of the satellite to the ground target at the current time: if the RAAN value is within the RAAN visibility condition, and the AoL value is within the AoL visibility condition, the satellite is visible to the ground target at the current time, otherwise, the satellite is not visible to the ground target at the current time; The visibility judgment results of all times in the maximum concentrated visible period are integrated to obtain the visible time list of the satellite to the ground target.

4. The method according to claim 1, wherein the method is characterized by: The coverage index of all satellites on the orbital plane to a single ground target is the average visibility duration or the average revisit time of all satellites on the orbital plane to the ground target.

5. The method of claim 1, wherein the method is characterized by: The step S3 comprises the following steps: For each satellite on the orbital plane, according to the semi-major axis correction value of the satellite and the face-to-area ratio, the orbital decay parameter is obtained from the first query table, and then the RAAN change amount of the satellite in the mission period is calculated; The average RAAN change amount of all satellites on the orbital plane in the mission period is calculated; The difference between the nominal RAAN of the orbital plane and the average RAAN change amount is obtained to obtain the RAAN correction value of the orbital plane; For each satellite on the orbital plane, according to the RAAN correction value of the orbital plane and the RAAN change amount of the satellite in the mission period, the RAAN correction value of the satellite is obtained.

6. The method of claim 1, wherein the method is characterized by: The step S1 comprises: Firstly, the coverage range data and the orbit drift data under different characteristic semi-major axes are acquired; the coverage range data includes the coverage central angle of the satellite with different coverage half-angle under each characteristic semi-major axis; the orbit drift data includes the semi-major axis attenuation parameter and the RAAN change rate of the satellite with different surface quality ratio under each characteristic semi-major axis, and the semi-major axis attenuation parameter includes the quadratic semi-major axis attenuation parameter and the linear semi-major axis attenuation parameter; Finally, the initial time RAAN-AoL visible condition set and the initial RAAN visible characteristic range set of each ground target of the satellite are obtained in combination with the satellite-ground visible condition.

7. The method according to claim 6, wherein the method is characterized by: In the step S1, the initial time RAAN-AoL visible condition and the initial RAAN visible characteristic range of the satellite to each ground target are calculated by the following process: According to the characteristic semi-major axis and the coverage half-angle of the satellite, one or two groups of initial time RAAN visible conditions of the satellite to the ground target are calculated, and each group of initial time RAAN visible condition is mainly composed of a boundary condition, a center point and a visible range; For each group of initial time RAAN visible condition, a plurality of RAAN characteristic points are generated by discretization according to the boundary condition and the center point, and two AoL boundary points corresponding to each RAAN characteristic point are calculated, and the initial time RAAN-AoL visible condition is obtained by polynomial fitting of all RAAN characteristic points and corresponding AoL boundary points.

8. The method according to claim 7, wherein the method is characterized by: When the satellite satisfies φ n < i - θ, the satellite has two sets of initial RAAN visibility conditions for the ground target n, and the RAAN difference between the two sets of boundary points is the same but the numerical value is different; when the satellite satisfies i - θ < φ n < θ, the satellite has one set of initial RAAN visibility conditions for the ground target n; wherein φ n represents the latitude of the ground target n, θ represents the coverage central angle of the satellite, and i represents the orbital inclination of the orbital plane to which the satellite belongs.

Images