Heterogeneous satellite constellation configuration correction method based on ground coverage performance stability

Abstract

The invention discloses a heterogeneous satellite constellation configuration correction method based on ground coverage performance stability. Comprising the steps of constructing a first look-up table and a second look-up table; for each orbital plane, analyzing and calculating a dynamic visibility condition and a visibility list by adopting a working orbital parameter optimization method in combination with the first lookup table and the second lookup table, and further correcting a working semi-major axis and a latitude argument AoL of each satellite on the orbital plane; for each orbital plane, calculating an ascending node right ascension RAAN correction value of the satellite according to the working semi-major axis correction value of the satellite on the orbital plane in combination with the first query table; the working semi-major axis correction values, the latitude argument AoL correction values and the ascending node right ascension RAAN correction values of all satellites on all orbit planes form the configuration parameters of the heterogeneous satellite constellation. According to the method, the earth coverage index drift of the earth coverage heterogeneous constellation is restrained, the constellation control frequency is reduced, and extra increase of the number of satellites due to the coverage index drift and multi-time artificial iteration design are avoided.

Description

TECHNICAL FIELD

[0001] The present application relates to a heterogeneous constellation configuration correction method, in particular to a heterogeneous constellation configuration correction method. BACKGROUND

[0002] Heterogeneous constellation refers to a constellation composed of satellites or formations with different functions, appearances, physical structures, and carried payloads. Such a constellation can integrate different types of satellites such as communication, navigation, and remote sensing, improve the working efficiency, functional density, cost-effectiveness ratio, and inter-satellite cooperation possibility of the constellation through information interconnection and resource integration, enhance the system risk resistance, and save the orbital resources. It is an inevitable trend of current and future constellation development.

[0003] Compared with traditional homogeneous constellations, the differences in the windward surface mass ratio of each satellite and the coverage range of the sensor in the heterogeneous constellation are the main factors affecting the constellation configuration design scheme and the stability of the configuration. The differences in the coverage range of each satellite make the analysis of the coverage capability of the constellation and the design of the configuration unable to rely on the traditional ideas of homogenization, uniformization, and linear superposition. At the same time, the differences in the drift characteristics of heterogeneous satellites cause the drift of the inter-satellite configuration and the coverage characteristics and coverage indicators. The design and control method of traditional homogeneous constellations will increase the number of satellites and the control cost. Among them, the coverage drift of the coplanar follow-up configuration caused by the Argument of Latitude (AoL) drift is the most serious, and the semi-major axis drift and the Right Ascension of the Ascending Node (RAAN) drift also have adverse effects on the stability of the constellation configuration and the coverage indicators.

[0004] The constellation design method can be roughly divided into design methods based on classical configuration, analytical methods based on geometric principles, and design methods based on modern optimization methods and intelligent algorithms according to the design principles. In the design method based on classical configuration, the common classical configurations include Walker constellation and its derivative constellation, Iridium constellation, etc. The analytical method based on geometric principles includes the coverage band method, ground track method, etc., which has strict mathematical proof and can be customized flexibly according to the needs, and is usually used for special configurations such as retrograde orbit. The design method based on modern optimization methods and intelligent algorithms focuses on the establishment and solution of the optimization model, rather than directly arranging the orbits using the projection of the subsatellite point. However, the current extensive research only solves the specific coverage problem of homogeneous constellations, does not solve the coverage performance analysis and evaluation of heterogeneous constellations under the difference in satellite coverage capability, and the configuration and coverage indicator drift caused by the difference in satellite drift trajectory, and does not solve the design problem of heterogeneous constellations from the aspects of coverage efficiency, configuration stability, coverage indicator stability, etc. At the same time, the performance analysis and design of heterogeneous constellations rely on simulation iteration, which greatly increases the complexity of manual design.

[0005] In the field of heterogeneous constellation design, some general methods have been proposed to design heterogeneous constellation configuration, but the heterogeneity involved is single, which is applicable to various heterogeneous satellite payloads, but the stability drift of coverage index caused by the difference in satellite face-mass ratio is not considered, and the constellation design based on the stability of the coverage performance to the earth is not considered when the heterogeneous payload and the off-plane mass ratio act simultaneously.

[0006] In summary, the traditional constellation configuration design method usually takes minimizing cost or maximizing coverage index as the goal, and since the configuration stability and coverage capacity stability of the heterogeneous constellation are not fully considered, directly applying it to the heterogeneous constellation will lead to the problem of coverage index drift over time, which requires additional satellites or constellation configuration maintenance, thereby increasing the manufacturing and control cost of the constellation. SUMMARY

[0007] For the effective optimization problem of the configuration of the heterogeneous constellation for earth coverage, the present application provides a heterogeneous satellite constellation configuration correction method based on the stability of the coverage performance to the earth. The method of the present application not only realizes the analytical calculation of the coverage index of the heterogeneous constellation through geometric methods, but also solves the problem of coverage index drift of the heterogeneous constellation for earth coverage in medium and low orbits, and the problem of high constellation control frequency or additional increase in the number of satellites caused by coverage index drift. The present application provides an effective analytical method for the configuration optimization design of the heterogeneous constellation for earth coverage in medium and low orbits, avoiding multiple manual iterative designs.

[0008] The technical scheme adopted by the present application is: The method of the present application specifically comprises the following steps: Step S1: constructing a first query table and a second query table. The first query table contains the semi-major axis decay parameter and the RAAN change rate of satellites with different face-mass ratios at each characteristic semi-major axis; the second query table contains the coverage central angle of satellites with different coverage half-opening angles at each characteristic semi-major axis, and the initial time RAAN-AoL visibility condition and the initial RAAN visibility feature range for each ground target.

[0009] The step S1 comprises: Firstly, obtain the coverage range data and orbit drift data at different characteristic semi-major axes; the coverage range data includes the coverage central angle of satellites with different coverage half-opening angles at each characteristic semi-major axis; the orbit drift data includes the semi-major axis decay parameter and the RAAN change rate of satellites with different face-mass ratios at each characteristic semi-major axis, the semi-major axis decay parameter includes the quadratic term semi-major axis decay parameter and the linear term semi-major axis decay parameter; Finally, the initial RAAN-AoL visibility condition set and the initial RAAN visibility feature range set of each ground target of the satellite are obtained in combination with the satellite-ground visibility condition.

[0010] The initial RAAN-AoL visibility condition and the initial RAAN visibility feature range of each ground target of the satellite are calculated by the following process: first, according to the characteristic semi-major axis and the coverage half-angle of the satellite, a group or two groups of initial RAAN visibility conditions of the satellite to the ground target are calculated, and each group of initial RAAN visibility conditions mainly consists of boundary conditions (two boundary points), a center point and a visible range. Then, for each group of initial RAAN visibility conditions, a plurality of RAAN feature points including the boundary conditions Ω L , Ω R and the center point Ω M are generated by discretization according to the boundary conditions and the center point, and two AoL boundary points corresponding to each RAAN feature point are calculated. The polynomial fitting is performed on all RAAN feature points and corresponding AoL boundary points to obtain the initial RAAN-AoL visibility condition.

[0011] The process of a group or two groups of initial RAAN visibility conditions of the satellite to the ground target is specifically: when the satellite to the ground target n satisfies φ n < i-θ, the satellite to the ground target n has two groups of initial RAAN visibility conditions, and the RAAN difference between the two groups of boundary points is the same, but the numerical values of the two groups of boundary points are different; when the satellite to the ground target n satisfies i-θ < φ n < θ, the satellite to the ground target n has one group of initial RAAN visibility conditions. 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.

[0012] Step S2: for each orbital plane, the working orbit parameter optimization method is adopted to analyze and calculate the dynamic visibility condition and the visibility list in combination with the first query table and the second query table, and then the working semi-major axis and the latitude argument AoL of each satellite on the orbital plane are corrected, and finally the working semi-major axis correction value and the latitude argument AoL correction value of each satellite on the orbital plane are obtained.

[0013] In the step S2, the working orbit parameter optimization method includes: Step S2.1: for each satellite on the orbital plane, the initial working semi-major axis, the initial AoL, the coverage central angle and the maximum AoL offset of the satellite are obtained by the analytical method in combination with the first query table and the second query table.

[0014] The analytical method includes: According to the face quality ratio and the nominal semi-major axis of the satellite, a semi-major axis attenuation parameter of the satellite is obtained from a first query table; according to the coverage half-angle and the nominal semi-major axis of the satellite, a coverage central angle of the satellite is obtained from a second query table; According to the nominal semi-major axis and the semi-major axis attenuation parameter of the satellite, an initial working semi-major axis is calculated; According to the total number of satellites on the orbital plane, a nominal AoL of each satellite is obtained; According to the initial working semi-major axis, the nominal semi-major axis and the semi-major axis attenuation parameter of the satellite, a maximum AoL offset is obtained; According to the nominal AoL and the maximum AoL offset of the satellite, an initial AoL of the satellite is obtained.

[0015] 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, then step S2.3 is entered, otherwise, step S2.5 is entered.

[0016] 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 analytically to obtain the coverage index of the orbital plane within the mission period.

[0017] The step S2.3 includes: Step S2.3.1: For each satellite on the orbital plane, according to the face quality ratio, the coverage half-angle and the initial working semi-major axis of the satellite, the RAAN change rate, the coverage central angle of the satellite are obtained from the first query table and the second query table, and the initial time RAAN-AoL visibility condition and the initial RAAN visibility feature range of each ground target are obtained, and then the maximum concentrated visible period of the satellite to each ground target is obtained.

[0018] The maximum concentrated visible period of the satellite to each ground target is represented as:

[0019] In the formula, Γ Ω0,p,s represents the maximum concentrated visible period of the satellite s to the ground target n on the orbital plane p; A1 represents the initial offset coefficient, A2 represents the long-term drift proportion coefficient, t Γ1 or t Γ2 represents the single maximum concentrated visible period duration, t B +t Γ1 or t B +t Γ2 represents the maximum concentrated visible period repetition period, τ1, τ2 or τ3 represents the long-term drift time, φ n represents the latitude of the ground target n, θ represents the coverage central angle of the satellite s, and i represents the orbital inclination of the orbital plane p.

[0020] Step S2.3.2: For each satellite on the orbital plane, according to the maximum concentrated visible time period of the satellite to each ground target and the initial AoL, obtain the list of visible time instants of the satellite to each ground target.

[0021] The step S2.3.2 comprises: According to the surface-to-mass ratio of the satellite, the coverage half-angle and the initial working semi-major axis, obtain the orbital decay coefficient, the RAAN change rate and the initial time RAAN-AoL visibility condition of the satellite to the ground target from the first query table and the second query table; Traverse each time in the maximum concentrated visible time period, and according to the orbital decay coefficient, the RAAN change rate and the initial AoL of the satellite, obtain the RAAN value and the AoL value of the satellite at the current time; map the initial time RAAN-AoL visibility condition of the satellite to the ground target to the current time to obtain the RAAN-AoL visibility condition of the satellite to the ground target at the current time; and 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, judge the visibility: 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; Integrate the visibility judgment results of all times in the maximum concentrated visible time period to obtain the list of visible time instants of the satellite to the ground target.

[0022] Step S2.3.3: Merge the visible time instant lists of all satellites on the orbital plane to the same ground target to obtain the visible time period list of all satellites on the orbital plane to the ground target.

[0023] Step S2.3.4: According to the visible time period list of all satellites on the orbital plane to each ground target, calculate to obtain the coverage index and coverage index consistency of all satellites on the orbital plane to each ground target.

[0024] The coverage index of all satellites on the orbital plane to each ground target is the average visible time length or the average revisit time of all satellites on the orbital plane to the ground target.

[0025] Step S2.3.5: Obtain the average coverage index and the average coverage index consistency of all satellites on the orbital plane to all ground targets to form the coverage index of the orbital plane.

[0026] Step S2.4: If the coverage index of the orbital plane meets the preset condition, the initial working semi-major axis obtained in step S2.1 is taken as the working semi-major axis correction value of the satellite, and the initial AoL is taken as the latitude argument AoL correction value; otherwise, step S2.5 is entered. In specific implementation, the preset condition can be determined according to the coverage performance requirement in the task design stage.

[0027] Step S2.5: using an optimization algorithm, taking the initial working semi-major axis and initial AoL of all satellites on the orbital plane as optimization variables, and taking the minimum coverage index of the orbital plane in the mission cycle as the optimization target, in the optimization process, the coverage index is obtained according to step S2.3, and finally the working semi-major axis correction value and the latitude argument AoL correction value of each satellite on the orbital plane are obtained.

[0028] Step S3: for each orbital plane, according to the working semi-major axis correction value of each satellite on the orbital plane, the first query table is combined to calculate the RAAN correction value of each satellite.

[0029] The step S3 comprises the following steps: For each satellite on the orbital plane, the orbital decay parameter is obtained from the first query table according to the working semi-major axis correction value and the face quality ratio, and then the RAAN change amount of the satellite in the mission cycle is calculated; The average RAAN change amount of all satellites on the orbital plane in the mission cycle is calculated; The difference between the nominal RAAN of the orbital plane and the average RAAN change amount is obtained, and the RAAN correction value of the orbital plane is obtained; For each satellite on the orbital plane, the difference between the RAAN correction value of the orbital plane and the RAAN change amount of the satellite in the mission cycle is obtained, and the RAAN correction value of the ascending node of the satellite is obtained.

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

[0031] The beneficial effects of the present application are: 1. On the basis of the traditional uniform constellation such as walker constellation, the present application obtains the differentiated working orbit parameters of each satellite through the compensation and correction of the three types of orbit parameters of semi-major axis, AoL and RAAN, and increases the stability of the coverage index of the low earth orbit circular heterogeneous constellation composed of satellites with different coverage half-angle and different coverage value. The stability is represented by the consistency of the coverage index, that is, the coverage index in the mission cycle is always close to the design value. Common coverage indexes include revisit time, transit time, transit times, coverage number, regional coverage imaging time, etc.

[0032] 2. The present application realizes the beneficial effects of suppressing the coverage index drift of the heterogeneous constellation, reducing the constellation control frequency, avoiding the additional increase of the number of satellites due to the coverage index drift, and avoiding multiple manual iterative design. BRIEF DESCRIPTION OF DRAWINGS

[0033] Figure 1 Figure 1 is a schematic diagram of a subsatellite point trajectory and ground station coverage.

[0034] Figure 2 Figure 2 is a schematic diagram of the relationship between the ascending node right ascension-ecliptic latitude-ecliptic longitude range and the orbit inclination and the ground station latitude.

[0035] Figure 3 Figure 3 is a schematic diagram of the relationship between the satellite coverage half-angle and the coverage central angle.

[0036] Figure 4 Figure 4 is a schematic diagram of the ascending node right ascension range and the maximum concentrated visible time period.

[0037] Figure 5 Figure 5 is a schematic diagram of the ascending node right ascension-ecliptic latitude-ecliptic longitude range over time.

[0038] Figure 6 Figure 6 is a satellite motion trajectory obtained by solving satellite orbit parameters by an analytical method.

[0039] Figure 7 Figure 7 is a schematic diagram of the overall process of the method of the present application.

[0040] Figure 8 Figure 8 is a schematic diagram of the process of the method of the present application. DETAILED DESCRIPTION

[0041] In order to make the purpose, technical solutions and advantages of the present application clearer and more apparent, the present application is further described in detail below in combination with embodiments. The specific embodiments described herein are only used to explain the present application and do not constitute any limitation on the present application.

[0042] The method of the present application specifically comprises the following steps: Step S1: constructing a first query table and a second query table; the purpose of this step is to save subsequent computing resources. The first query table contains the semi-major axis decay parameters and the RAAN change rate of satellites with different surface-to-mass ratios at each characteristic semi-major axis. The second query table contains the coverage central angle of satellites with different coverage half-angles at each characteristic semi-major axis, and the initial time RAAN-AoL visibility condition and the initial RAAN visibility characteristic range of each satellite for each ground target.

[0043] Specifically, step S1 comprises: First, obtain coverage range data and orbit drift data at different characteristic semi-major axes; the coverage range data includes the coverage central angle of satellites with different coverage half-angles at each characteristic semi-major axis; the orbit drift data includes the semi-major axis decay parameters and the RAAN change rate of satellites with different surface-to-mass ratios at each characteristic semi-major axis, the semi-major axis decay parameters including the quadratic semi-major axis decay parameters and the linear semi-major axis decay parameters; Finally, the initial RAAN-AoL visibility condition set and the initial RAAN visibility feature range set of each ground target of the satellite are obtained in combination with the satellite-ground visibility condition. The RAAN-AoL visibility condition refers to the RAAN and AoL conditions that need to be met for the satellite and the ground target to be visible.

[0044] Further, the initial RAAN-AoL visibility condition and the initial RAAN visibility feature range of each ground target of the satellite can be calculated by the following process: First, according to the characteristic semi-major axis a r,x and the covering half-angle a y of the satellite at the initial time, a set or two sets of initial RAAN visibility conditions of the satellite to the ground target are calculated, each set of initial RAAN visibility conditions (Ω L , Ω R , Ω M , △L) is mainly composed of boundary conditions Ω A , Ω D , a center point Ω M and a visible range △L. Specifically: when the satellite satisfies φ n < i-θ to the ground target n, the satellite has two sets of initial RAAN visibility conditions to the ground target n, the RAAN difference between the two sets of boundary points is the same, but the numerical values of the two sets of boundary points are different; when the satellite satisfies i-θ < φ n < θ to the ground target n, the satellite has one set of initial RAAN visibility conditions to the ground target n. Wherein, φ n represents the latitude of the ground target n, θ represents the covering central angle of the satellite, and i represents the orbital inclination of the orbital plane to which the satellite belongs.

[0045] Finally, for each set of initial RAAN visibility conditions (Ω L , Ω R , Ω M , △L), according to the boundary conditions and the center point, a plurality of RAAN feature points including the boundary conditions Ω L , Ω R and the center point Ω M are generated by discretization, and two AoL boundary points u1, u2 corresponding to each RAAN feature point are calculated, all RAAN feature points and corresponding AoL boundary points are polynomial fitted to obtain an initial RAAN-AoL visibility condition f n (Ω t=0 , u1, u2), and the initial RAAN-AoL visibility condition is a function f n (Ω t=0 , u1, u2) of the visibility polygon.

[0046] Further, the satellite semi-major axis decay data can be obtained by the following process: simulating the motion trajectory of a satellite with a specific characteristic semi-major axis and a specific area-mass ratio within a specific mission duration T by using a simulation method, extracting semi-major axis variation data and RAAN variation amount ΔΩ x,z . The semi-major axis variation data is fitted by using a quadratic fitting method to obtain a quadratic term semi-major axis decay parameter and a linear term semi-major axis decay parameter. The RAAN variation amount ΔΩ x,z is divided by the simulation mission duration T to obtain a RAAN variation rate.

[0047] Step S2: for each orbital plane, a working orbit parameter optimization method is used to analyze and calculate the dynamic visibility condition and the visibility list in combination with the first query table and the second query table, and then the working semi-major axis and the latitude angle AoL of each satellite on the orbital plane are corrected.

[0048] In step S2, the working orbit parameter optimization method includes: Step S2.1: for each satellite on the orbital plane, the initial working semi-major axis, the initial AoL, the coverage central angle, and the maximum AoL offset of the satellite are obtained by using an analytical method in combination with the first query table and the second query table.

[0049] The analytical method specifically includes: For each satellite on the orbital plane, the semi-major axis decay parameter of the satellite is obtained from the first query table according to the area-mass ratio and the nominal semi-major axis of the satellite, and the coverage central angle of the satellite is obtained from the second query table according to the coverage half-angle and the nominal semi-major axis of the satellite; The initial working semi-major axis is calculated according to the nominal semi-major axis and the semi-major axis decay parameter of the satellite; The maximum AoL offset is obtained according to the initial working semi-major axis, the nominal semi-major axis, and the semi-major axis decay parameter of the satellite; The nominal AoL of each satellite is calculated according to the total number of satellites on the orbital plane according to the following formula: u s '=[2(s-1)π] / S p , s∈[1,S p ] In the formula, u s ' represents the nominal AoL of the satellite s, and S p represents the total number of satellites on the orbital plane where the satellite s is located; The initial AoL of each satellite is calculated according to the nominal AoL and the maximum AoL offset of each satellite according to the following formula: u s =u s '-△u max,s / 2 In the formula, u sdenotes the initial AoL of satellite s, u s denotes the nominal AoL of satellite s, Δu max,s denotes the maximum AoL offset of satellite s.

[0050] 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 of itself, that is, the use condition of the analytical method is met, then go to step S2.3, otherwise, go to step S2.5.

[0051] 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 analytically to obtain the coverage index of the orbital plane within the mission period. The coverage index of the orbital plane includes the average coverage index of all satellites on the orbital plane to all ground targets and the consistency of the average coverage index.

[0052] Step S2.3 includes: 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 the satellite to each ground target 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. The maximum concentrated visible period of the satellite to each ground target is expressed as:

[0053] In the formula, Γ Ω0,p,s denotes the maximum concentrated visible period of satellite s to ground target n on the orbital plane p; A1 denotes the initial offset coefficient, A2 denotes the long-term drift proportion coefficient, t Γ1 or t Γ2 denotes the single maximum concentrated visible period duration, t B +t Γ1 or t B +t Γ2 denotes the maximum concentrated visible period repetition period, τ1, τ2 or τ3 denotes the long-term drift time, φ n denotes the latitude of ground target n, θ denotes the coverage central angle of satellite s, i denotes the orbital inclination of the orbital plane p.

[0054] 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 visibility time list of the satellite to each ground target is obtained.

[0055] Step S2.3.2 includes: According to the area-to-mass ratio, the half opening angle of coverage and the initial working semi-major axis of the satellite, the orbit decay coefficient, the RAAN change rate and the initial RAAN-AoL visibility condition of the satellite to the ground target are obtained from the first query table and the second query table; According to the orbit decay coefficient, the RAAN change rate and the initial AoL of the satellite, the RAAN value and the AoL value of the satellite at the current time are obtained by traversing each time in the maximum concentrated visible period. The initial 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. 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, the visibility is judged: 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 a visible time list of the satellite to the ground target.

[0056] Step S2.3.3: merging the visible time lists of all satellites on the orbit plane to the same ground target to obtain a visible period list of all satellites on the orbit plane to the ground target.

[0057] Step S2.3.4: according to the visible period list of all satellites on the orbit plane to each ground target, the coverage index and the coverage index consistency of all satellites on the orbit plane to each ground target are calculated.

[0058] Optionally, the coverage index of all satellites on the orbit plane to each ground target is the average visible duration or the average revisit time of all satellites on the orbit plane to the ground target.

[0059] Step S2.3.5: obtaining the average coverage index and the average coverage index consistency of all satellites on the orbit plane to all ground targets, which constitute the coverage index of the orbit plane.

[0060] Step S2.4: if the coverage index of the orbit plane meets the preset condition, the initial working semi-major axis obtained in step S2.1 and the initial AoL are taken as the working semi-major axis correction value and the latitude argument AoL correction value of the satellite respectively; otherwise, step S2.5 is entered. In specific implementation, the preset condition can be determined according to the coverage performance requirement in the task design stage.

[0061] Step S2.5: Using an optimization algorithm, taking the initial working semi-major axes and initial AoLs of all satellites on the orbital plane as optimization variables, and taking the minimum coverage index of the orbital plane within the mission period as the optimization target, in the optimization process, the coverage index is obtained according to step S2.3, and finally the working semi-major axis correction value and the latitude argument AoL correction value of each satellite on the orbital plane are obtained.

[0062] Step S3: For each orbital plane, the working semi-major axis correction value of each satellite on the orbital plane is combined with the first query table to calculate the RAAN correction value of each satellite respectively; the working RAAN correction can ensure the stability of the relative configuration between the orbital planes under the differentiated working semi-major axis.

[0063] Step S3 includes the following steps: For each satellite on the orbital plane, the orbital decay parameter is obtained from the first query table according to the working semi-major axis correction value and the specific ratio of the satellite, and then the RAAN change amount of the satellite within the mission period is calculated; The average RAAN change amount of all satellites on the orbital plane within 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, the difference between the RAAN correction value of the orbital plane and the RAAN change amount of the satellite within the mission period is obtained to obtain the RAAN correction value of the satellite.

[0064] In step S3.1, the RAAN change amount of the satellite within the mission period is calculated according to the following formula: △Ω=-(7Ω’·k1·T 3 ) / (6a)-(7Ω’·k2·T 2 ) / (4a)-[(7Ω’·T) / (2a)]·(a s -a) Ω’=-[(3n e ·J2·R e 2 ) / (2a 2 )]·cosi In the formula, △Ω represents the RAAN change amount of the satellite within the mission period, Ω' represents the long-term change rate of RAAN, k1 represents the quadratic semi-major axis decay parameter, k2 represents the linear semi-major axis decay parameter, T represents the mission period, a represents the nominal semi-major axis of the satellite, a s represents the initial working semi-major axis of the satellite, n e represents the orbital angular velocity, J2 represents the Earth's oblateness perturbation coefficient, R e represents the Earth's radius, and i represents the orbital inclination.

[0065] Step S4: The working 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.

[0066] The specific implementation of the present application is as follows:

[0067] Embodiment 1 In this embodiment, the analytic method mainly uses the geometric method proposed by Ulybyshev Y et al. to realize the analytic calculation of the coverage index of the heterogeneous constellation, and further deduces the parameters such as the satellite-ground concentrated visible period. The geometric method used by the analytic method maps the relationship between the RAAN and the AoL of the satellite satisfying the satellite-ground visible condition with time into a graph, which is called satellite-ground visible relationship mapping graph, and then calculates the satellite-ground visible period and other coverage parameters through the graph parameters. For details, please refer to the following papers of the team: Ulybyshev Y. Satellite Constellation Design for Complex Coverage[J]. Journal of Spacecraft and Rockets, 2008, 45(4): 843-849. Ulybyshev, Yuri. Geometric Analysis and Design Method for Discontinuous Coverage Satellite Constellations[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(2): 549-557. The geometric method used by the analytic method in this embodiment is described as follows: Figure 1 is a schematic diagram of the subsatellite point trajectory and the ground station coverage. For a target ground station M(λ, φ) with geocentric longitude and latitude of λ and φ, when the subsatellite point trajectory passes through the circular region with M(λ, φ) as the center and θ as the radius, the coverage can be realized. Figure 1 In the formula, i represents the orbital inclination, △λ represents the longitude difference between the boundary orbits that can realize the ground station coverage, and θ represents the coverage geocentric angle.

[0068] As shown in (a) of FIG. 1, Figure 1 When φ < i-θ, both the ascending trajectory and the descending trajectory can realize the coverage. In the figure, L A1 , L A2 , L A3 is the ascending trajectory, L D1 , L D2 , LD3 is the downward trajectory, L A2 and L D2 is the longest trajectory passing through the center of the circle, and the rest are trajectories tangent to the covering circle, with the tangent point N A 、S A 、N D 、S D is the boundary point. A With N D Unified as the north boundary point N, S A With S D They are uniformly denoted as the south boundary point S, with the subscript A indicating the upward track and D indicating the downward track.

[0069] like Figure 1 As shown in (b), when i-θ<φ <i+θ时,拐点位于覆盖圆之内的轨迹可以实现覆盖。图中,L2为最长轨迹,其余为与覆盖圆相切的轨迹,切点S A 、S D is the boundary point.

[0070] The two aforementioned papers provide two conditions for RAAN and AoL to be met when satellites can cover ground stations. Based on these conditions, this embodiment specifically adopts the following geometric method: (1) If φ <i-θ,则当卫星能够覆盖到地面站时,目标地面站M、北侧边界点N和南侧边界点S在轨道中对应的RAAN分别为: Ω M =λ+G-arctan(tanu M ·cosi) Ω NA =λ+G-△λ MN -arctan(tanu NA ·cosi), Ω ND =2(λ+G)-Ω NA -π Ω SA =λ+G+△λ MS -arctan(tanu SA ·cosi), Ω SD =2(λ+G)-Ω SA -π Where, Ω M ,Ω NA ,Ω ND ,Ω SA ,Ω SD They represent the target ground station M, the north tangent point N between the uplink track and the coverage circle, and A , the north tangent point N of the descending track and the covering circle D, the south side tangent point S of the uplink orbit and the coverage circle A , the south side tangent point S of the downlink orbit and the coverage circle D the corresponding RAAN in the orbit; Δλ MN , Δλ MS respectively represent the longitude difference of the target ground station M and the north side boundary point N, the south side boundary point S; u M , u NA , u SA respectively represent the target ground station M, the north side tangent point N of the uplink orbit and the coverage circle A , the south side tangent point S of the uplink orbit and the coverage circle A the corresponding AoL in the orbit; λ represents the geocentric longitude of the target ground station M, φ represents the geocentric latitude of the target ground station M, θ represents the coverage geocentric angle, i represents the orbit inclination, and G represents the Greenwich hour angle.

[0071] the corresponding AoL in the orbit of the target ground station M, the north side boundary point N and the south side boundary point S: u M , u NA , u ND , u SA , u SD are obtained by the following formulas respectively: u M = arcsin (sinφ / sini) u NA = arcsin (sinφ N / sini), u ND = π - u NA u SA = arcsin (sinφ S / sini), u SD = π - u SA In the formula, u M , u NA , u ND , u SA , u SD respectively represent the target ground station M, the north side tangent point N of the uplink orbit and the coverage circle A , the north side tangent point N of the downlink orbit and the coverage circle D , the south side tangent point S of the uplink orbit and the coverage circle A , the south side tangent point S of the downlink orbit and the coverage circle D the corresponding AoL in the orbit; φ, φ N , φ S respectively represent the latitude of the target ground station M, the north side boundary point N, the south side boundary point S; φ represents the latitude of the target ground station M; φ N represents the north side tangent point N of the uplink orbit and the coverage circleA The latitude or the point where the descending track touches the north side of the covered circle N D The latitude of the two is the same; φ S The south tangent point S of the upward track and the coverage circle A The latitude or the point of tangency of the descending track with the south side of the covered circle S D The latitude of , the two values ​​are the same.

[0072] 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 by the following formulas: △λ MN =arccos[(cosθ-sinφ N ·sinφ) / (cosφ N ·cosφ)] △λ MS =arccos[(cosθ-sinφ S ·sinφ) / (cosφ S ·cosφ)] sinφ N =sinφ·cosθ+sinθ·cosi sinφ S =sinφ·cosθ-sinθ·cosi cosφ N =(1-sin 2 φ N ) 1 / 2 cosφ S =(1-sin 2 φ S ) 1 / 2 Where, △λ MN , △λ MS They represent the longitude differences between the target ground station M and the northern boundary point N and the southern boundary point S respectively; θ represents the coverage geocentric angle, i represents the orbital inclination, φ, φ N 、φ S They represent the latitudes of the target ground station M, the northern boundary point N, and the southern boundary point S respectively.

[0073] Taking into account the rotation of the earth, the actual RAAN range that can be covered is: △L=△λ N+ △λ MN+ △λ MS -△λ S -(w e T d / (2π))|u N-u S | △λ N =arcsin[(coti)·(sinφ N / cosφ N )] △λ S =arcsin[(coti)·(sinφ S / cosφ S )] wherein, △λ N , △λ S respectively represent longitude span of the satellite along the orbit from the equator to the northern boundary point N and the southern boundary point S; △λ MN , △λ MS respectively represent longitude difference between the target ground station M and the northern boundary point N and the southern boundary point S; u N , u S respectively represent AoL corresponding to the northern boundary point N and the southern boundary point S in the orbit; w e represents the earth rotation angular velocity, T d represents the orbit period; i represents the orbit inclination, φ N , φ S respectively represent latitude of the northern boundary point N and the southern boundary point S; | | represents absolute value.

[0074] It should be noted that in the specific implementation, when the geometric method is used to obtain the RAAN range of the uplink orbit, u N , u S respectively take u NA , u SA , and vice versa, when the RAAN range of the downlink orbit is obtained, u N , u S respectively take u ND , u SD .

[0075] (2) If i-θ<φ<i+θ, when the satellite can cover the ground station, the RAAN corresponding to the boundary point in the orbit is: Ω SA =λ+G+△λ MS -arctan(tanu SA ·cosi) Ω SD =λ+G-△λ MS -arctan(tanu SD ·cosi)-π wherein, Ω SA , Ω SD respectively represent the tangent point S A , the tangent point S Dthe corresponding RAAN in the orbit; u SA , SD respectively represent the tangent point S A , the tangent point S D the corresponding AoL in the orbit; Δλ MS represents the difference in longitude 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 hour angle.

[0076] Considering the Earth's rotation, the actual RAAN range ΔL that can be covered is: ΔL = 2(Δλ MS - Δλ S ) + π - (w e T d / (2π)) |u SD - u SA | In the formula, Δλ MS represents the difference in longitude between the target ground station M and the boundary point S; Δλ S represents the longitude span of the satellite along the orbit from the equator to the boundary point S; u SA , u SD respectively represent the tangent point S A , the tangent point S D the corresponding AoL in the orbit; w e represents the Earth's rotation angular velocity, T d represents the orbital period, and | | represents the absolute value.

[0077] (3) When the RAAN of the satellite is within the boundary points in the above (1) or (2), and the following AoL conditions are met, the satellite and the ground target are just visible: Assuming that the actual RAAN of the satellite is Ω, the sub-satellite point trajectory intersects the coverage circle at two points S1 and S2, and the following equation set is solved to obtain the AoL of the intersection point S1 and the intersection point S2: u1, u2: α = Ω + arctan(tan u · cos i) δ = arcsin(sin u · sin i) cos θ = (cos φ) · (cos δ) · cos(|α M - α |) + (sin φ C ) · (sin δ1) 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, α M represents the right ascension of the target ground station, and u represents the AoL to be solved, with two values of u being u1 and u2. In specific implementation, the intersection point S1 and the intersection point S2 are substituted into this equation set to further obtain the AoL of the intersection point S1 and the intersection point S2: u1 and u2.

[0078] Each RAAN value corresponds to two AoL values except for the boundary points, and the visible range is in a closed shape. For example, Figure 2 The RAAN-AoL visible ranges of satellites with different inclinations on a 400 km orbit to ground stations with the same longitude and latitude of 10°, 45°, 85° are plotted, covering a half opening angle of 30°. The shape enclosed by the solid line in the figure is the current RAAN-AoL visible range, and if the RAAN-AoL value combination of the satellite falls exactly within the solid line, the satellite-ground coverage can be achieved.

[0079] In the method of the embodiment, the orbit and orbit parameters of the approximate uniform constellation to be corrected are referred to as nominal orbit and nominal orbit parameters, and the actual orbit and orbit parameters of each satellite after correction are referred to as working orbit and working orbit parameters. It is assumed that the constellation to be corrected consists of P orbital planes, and the orbital plane number is represented by subscript p. The number of satellites in each orbital plane is S p , and each satellite is numbered in the order of 1~S p . The nominal RAAN of the orbital plane p is Ω p , and the adjustment range of the RAAN is [Ω min,p , Ω max,p ]. The nominal semi-major axis of the constellation is a, the adjustment range of the semi-major axis is [a min , a max ], the inclination is i, the eccentricity is 0, the mission time is T, and the constellation needs to be reconstructed by configuration maintenance beyond T time. Each satellite in the constellation carries a sensor, and the coverage half opening angle of the sensor is α y , y∈[1,Y],where Y is the number of different coverage half opening angles in the constellation. The face quality ratio of the windward face of each satellite in the constellation is S / m z , z∈[1,Z],where Z is the number of different face quality ratios in the constellation.

[0080] As shown in Figure 7 and Figure 8 , the specific steps of the method of the embodiment are as follows: Step 1: Data preprocessing 1.1 Feature semi-major axis sampling: select several semi-major axis values as feature semi-major axes a r in the adjustment range [a min , a max ] of the semi-major axis of the constellation at equal intervals. A total of X feature semi-major axes are taken, denoted as a r,x , x∈[1,X]。

[0081] 1.2 Coverage central angle calculation: 1.2.1 Let the feature semi-major axis number x=1.

[0082] 1.2.2 Let the coverage half opening angle number y=1.

[0083] 1.2.3 Obtain the semi-major axis a r,x of the satellite located on the semi-major axis a y and carrying the sensor with the coverage half-angle α x,y : a r,x =(R e sinθ x,y ) / tanα y +R e cosθ x,y where a r,x denotes the semi-major axis of the satellite, α y denotes the coverage half-angle of the satellite, θ x,y denotes the coverage central angle corresponding to the coverage range of the satellite, R e is the radius of the earth, and R e = 6378.139 km. The relationship between the coverage half-angle and the coverage central angle is shown in FIG. 1. Figure 3

[0084] 1.2.4 Let the coverage half-angle sequence number y = y + 1. If y ≤ Y, go to step 1.2.3, otherwise go to step 1.2.5.

[0085] 1.2.5 Let the semi-major axis sequence number x = x + 1. If x ≤ X, go to step 1.2.2, otherwise go to step 1.3.

[0086] 1.3 Satellite orbit drift rate calculation: 1.3.1 Let the semi-major axis sequence number x = 1.

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

[0088] 1.3.3 By simulation method, simulate the motion trajectory of the satellite located on the semi-major axis a r,x and with the surface-to-mass ratio S / m z in the task duration T, record the semi-major axis change data and the RAAN change amount ΔΩ x,z .

[0089] 1.3.4 Adopt the quadratic fitting method to fit the semi-major axis change data into the following formula: a x,z (t)=a r,x +k1 x,z t 2 +k2 x,z t, t ∈ [0, T] ​In the formula, k1 is a quadratic term semi-major axis decay parameter, k2 is a linear term semi-major axis decay parameter, a(t) represents the semi-major axis at time t; subscript x represents the characteristic semi-major axis sequence number, subscript z represents the surface quality ratio sequence number. Wherein, k2 < k1 < 0.

[0090] The reason for using the quadratic fitting method to fit the satellite semi-major axis decay data is that under the action of atmospheric resistance, the decay of the semi-major axis of the low earth satellite is approximately linear in a short time, and the decay rate increases when the time is longer. The quadratic fitting can better simulate the medium and long term decay of the semi-major axis. The semi-major axis in the task period can be described by the semi-major axis decay parameter and the initial semi-major axis after fitting, avoiding the introduction of empirical model or exponential model, and facilitating subsequent numerical calculation.

[0091] 1.3.5 Calculate the linearized RAAN rate γ according to the following formula Ω,x,z : γ Ω,x,z =△Ω x,z / T, x∈[1,X],z∈[1,Z] In the formula, γ Ω,x,z represents the linearized RAAN rate of the satellite located on the characteristic semi-major axis a r,x , and the surface quality ratio is S / m z , in the simulation task length T, △Ω x,z represents the RAAN change of the satellite located on the characteristic semi-major axis a r,x , and the surface quality ratio is S / m z , in the simulation task length T.

[0092] 1.3.6 Let the surface quality ratio sequence number z=z+1. If z≤Z, go to step 1.3.3, otherwise go to step 1.3.7.

[0093] 1.3.7 Let the characteristic semi-major axis sequence number x=x+1. If x≤X, go to step 1.3.2, otherwise go to step 2.

[0094] Step 2: Star-ground visible condition calculation 2.1 Ground target grid division: assuming that the ground target contains n1 area targets and n2 point targets, the n1 area targets are divided into several point targets by grid method, combined with the n2 point targets, forming N ground point targets, represented by subscript n, the longitude and latitude of ground point target n are represented as (λ n , φ n ).

[0095] 2.2 Star-ground visible condition calculation: 2.2.1 Let the ground target sequence number n=1.

[0096] 2.2.2 Initial time RAAN visibility condition calculation of ground target n: Using geometric method (1) or (2), the boundary conditions Ω r,x , Ω y , central point Ω A and RAAN visibility range ΔL that the satellite's RAAN should satisfy when the satellite covers the ground point target n are calculated, where the geocentric angle is calculated in step 1.2; and the boundary conditions Ω D , Ω M in the two orbit conditions (1) or (2) are unified as left boundary Ω A , Ω D , to obtain the initial time RAAN visibility condition set of ground target n: A L R = (Ω 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] where each combination of Ω L , Ω R , Ω M , ΔL (Ω L , Ω R , Ω M , ΔL) constitutes a set of initial RAAN visibility opportunities; subscript x represents the characteristic semi-major axis sequence number, and subscript z represents the surface quality ratio sequence number.

[0097] At the same time, the initial time RAAN visibility condition set A n Ω,t=0 of ground target n contains H sets of initial RAAN visibility opportunities, H∈[XY,2XY].

[0098] 2.2.3 Discretization of initial time RAAN visibility condition of ground target n: Each set of initial RAAN visibility opportunities in A n Ω,t=0 is traversed, and the following calculation is performed for each set of initial RAAN visibility opportunities: the boundary points Ω L , Ω R , the central point Ω L , Ω R and the center point Ω​M a set of initial RAAN visible feature points.

[0099] Finally, the discrete initial RAAN visible feature point set R n Ω,t=0 and the initial RAAN visible range set ΔL:

[0100] wherein R n Ω,t=0 and ΔL, each row is a set of initial RAAN visible opportunities, one-to-one correspondence, a total of H sets.

[0101] 2.2.4 Initial time RAAN-AoL visible condition feature set calculation of ground target n: For R n Ω,t=0 1~H sets of initial RAAN visible opportunities in R n Ω-u,t=0 , for each visible RAAN value, the geometric method (3) is used to calculate the AoL boundary points u1, u2 corresponding to each RAAN value, forming the initial time RAAN-AoL visible condition feature set R

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

[0103] 2.2.5 Initial time RAAN-AoL visible condition fitting restoration of ground target n: Polynomial fitting is performed on the feature points in R n Ω,t=0 to restore the AoL visible range corresponding to the points in the initial time RAAN visible range other than the feature point set R n Ω,t=0 , forming the complete initial time RAAN-AoL visible condition set: f Ω-u n (t=0)=f n (Ω t=0 , u1, u2) h , h∈[1, H] f Ω-u n (t=0) has H rows, each row is a set of initial RAAN-AoL visible opportunities. Each set of RAAN-AoL visible opportunities is represented by the function fn (Ω t=0 , u1, u2) describes the visible polygon, the polygon shape is referenced Figure 2 .

[0104] The reason for using polynomial fitting instead of traversal calculation is that the calculation amount of calculating the RAAN visible condition is small, but a large amount of calculation is required to accurately obtain the AoL visible condition corresponding to each RAAN and draw the visible range shape. Therefore, reduce the calculation burden at the cost of losing some feasible solutions, and under the condition of sufficient computing power, the complete shape can be preserved.

[0105] Further, it is preferred to use a polynomial of 6 times or more for fitting, so as to further balance the calculation accuracy and the calculation amount.

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

[0107] Step 3: Query table establishment 3.1 Query table 1 establishment: According to the calculation results of step 1, a query table 1 is established, which contains X characteristic semi-major axes a r , Z surface quality ratios S / m z Corresponding semi-major axis attenuation parameters k1, k2 and linearized RAAN change rate γ Ω of the satellite. The content and structure of the query table 1 are referenced to the table.

[0108] Table 1. Content and structure of query table 1

[0109] 3.2 Query table 2 establishment: According to the calculation results of step 1 and step 2, a query table 2 is established, which contains: a, X characteristic semi-major axes a r , Y corresponding coverage half-angle coverage central angle θ; b, X characteristic semi-major axes a r , Y coverage half-angle α, N ground targets corresponding initial time RAAN visible characteristic range △L; c, initial time RAAN-AoL visible condition f Ω-u (t = 0).

[0110] The content and structure of the query table 2 are referenced to table 2.

[0111] Table 2. Content and structure of query table 2

[0112] The purpose of establishing the query table is that the data in the query table is the data frequently used in the subsequent steps of the algorithm, and the pre-computation can save the resources and time of subsequent computation.

[0113] Step 4: Modeling of the semi-major axis and the latitude argument correction value Traverse the orbit planes 1~P, and for each orbit plane p, the semi-major axis and AoL correction is performed one by one. The initial working semi-major axis a s , s∈[1,S p ] and the initial AoL value u s , s∈[2,S p ] of each satellite in the orbit plane p are taken as variables to be corrected, and a working orbit parameter optimization model is established.

[0114] The nominal initial RAAN value of the orbit plane p is Ω0, and the nominal initial semi-major axis is a. The following steps are specifically performed: 4.1 Let the orbit plane number p=1.

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

[0116] 4.3 Let the satellite number s=1.

[0117] 4.4 The maximum concentrated visible period of the ground target n for the satellite s in the orbit plane p is calculated: The meaning of the maximum concentrated visible period is that if there are enough satellites in the orbit plane to satisfy the AoL visible condition at all times, then only the RAAN visible condition needs to be satisfied, and there must be a satellite in the orbit that can achieve coverage. For example, Figure 4 The RAAN visible condition of the ground station with a latitude of 30° for different orbits is plotted with respect to time, and only the boundaries of the visible range are retained. With the increase of time, the boundary points of the RAAN visible condition are translated at the rate of the earth rotation. The band-shaped range between the two blue oblique lines in the figure is the visible RAAN range of the 30° inclination orbit. The band-shaped range between the red oblique lines is the visible RAAN range of the 50° inclination orbit, and the dashed line is the descending orbit and the solid line is the ascending orbit. When the RAAN trajectory (blue vertical line) of satellite 1 intersects with the blue band-shaped range, satellite 1 has a visible opportunity to the target ground station, and the t1~t2 period is the maximum concentrated visible period of the orbit. When the RAAN trajectory (red vertical line) of satellite 2 intersects with the red band-shaped range, satellite 2 has a visible opportunity to the target ground station, and the t3~t4 period has a descending visible opportunity and the t5~t6 period has an ascending visible opportunity, and t3~t4, t5~t6 are the maximum concentrated visible periods Γ of the orbit.

[0118] Only in the maximum concentrated visible period, the satellite has the visible opportunity. The drift difference of heterogeneous satellites will cause the relative position of the satellite in the maximum concentrated visible period to change constantly, and then affect the coverage performance. Therefore, the design of the position of each satellite in each maximum concentrated visible period and the relative position relationship plays a key role in improving the constellation coverage performance.

[0119] 4.4.1 In the query table 2-ground target n data, find the sensor half-angle closest to the half-angle of the current satellite s, and the characteristic semi-major axis is closest to the initial working semi-major axis a s of the current satellite s. The closest case records its coverage central angle θ, initial RAAN visible boundary point Ω R and initial RAAN visible range ΔL.

[0120] Where Ω R is O0or O u0 , O d0 , depending on the relationship between i-θ and φ. Specifically: O0is the right boundary point in the orbit corresponding to the RAAN value of S A in the orbit; O u0 , O d0 are the right boundary points in the orbit corresponding to the RAAN value of S A , N D in the orbit, respectively.

[0121] Where Δ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-θ.

[0122] Note: The initial working semi-major axis a s of the current satellite s will be given by the solving algorithm, the same below.

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

[0124] 4.4.3 The maximum concentrated visible period Γ Ω0,p,s of the ground target n to the orbit surface p satellite s:

[0125] Where Γ Ω0,p,srepresents the maximum concentrated visible time period of the satellite s on the orbital plane p to the ground target n; A1 represents an initial offset coefficient, A2 represents a long-term drift proportional coefficient, t Γ1 or t Γ2 represents the single maximum concentrated visible time period duration, t B +t Γ1 or t B +t Γ2 represents the maximum concentrated visible time period repetition period, τ1, τ2 or τ3 represents a long-term drift time, φ n represents the latitude of the ground target n, θ represents the coverage central angle of the satellite s, i represents the orbital inclination of the orbital plane p; r is a natural number.

[0126] Wherein, each parameter is specifically: A1=Ω0tanγ, A2=tanγ Ω / (tanγ Ω +tanγ) t B1 =(2π-△L1) / w e , t B2 =(2π-△L2) / w e , γ=arctan(1 / w e ) t Γ1 =△L1tanγ, t Γ2 =△L2tanγ τ1=(2π-O0)tanγ, τ2=(2π-O u0 )tanγ, τ3=(2π-O d0 )tanγ.

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

[0128] 4.5 Calculation of the visible time of the satellite s on the orbital plane p to the ground target n: Suppose that the maximum concentrated visible time period of the ground target n to the satellite s on the orbital plane p is J, which is described by subscript j, j∈[1, J] within T time. The following steps are performed: 4.5.1 Let the maximum concentrated visible time period number j=1.

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

[0130] 4.5.3 Satellite s orbit parameter calculation: In the query table 1, find the closest case to the face quality ratio of the current satellite s, and the characteristic semi-major axis and the initial working semi-major axis a of the current satellite s s The closest case, record its orbit decay parameters k1 and k2 and RAAN rate γ Ω ; Calculate the RAAN value and AoL value of the orbit plane p of the satellite s at the current time t, denoted as Ω t,p,s , u t,p,s : Ω t,p,s =Ω0+γ Ω ·t u t,p,s =u s -[(7λ’+3n e )·k1·t 3 ] / (6a)-[(7λ’+3n e )·k2·t 2 ] / (4a)-[(7λ’+3n e )·(a s -a)·t] / (2a) λ’=[(3n e J2R e 2 ) / (2a 2 )]·(4cos 2 i-1) In the formula, Ω t,p,s represents the RAAN value of the orbit plane p of the satellite s at the current time t, Ω0 represents the nominal initial RAAN value of the satellite, γ Ω represents the RAAN rate of the satellite; u t,p,s represents the AoL value of the orbit plane p of the satellite s at the current time t, u s represents the initial AoL of the satellite, λ’ represents the long-term change rate of the mean anomaly, n e represents the nominal semi-major axis corresponding to the orbit angular velocity, k1 represents the quadratic semi-major axis decay parameter of the satellite, k2 represents the linear semi-major axis decay parameter of the satellite, a s represents the initial working semi-major axis of the satellite, i.e. the semi-major axis value currently being optimized, a represents the orbit semi-major axis of the satellite before correction, i.e. the nominal semi-major axis, J2=1082.63×10 -6 .

[0131] 4.5.4 RAAN-AoL visibility condition calculation of the current satellite s to the ground target n at time t: In the query table 2-ground target n data, find the closest case to the sensor half-angle and the half-angle of the current satellite s, and the characteristic semi-major axis and the initial working semi-major axis a of the current satellite s sThe closest case, 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: f Ω-u n (t)=f n (Ω t=0 +w e t, u1, u2).

[0132] The reason why the RAAN-AoL visibility condition of the current satellite s for the ground target n at time t can be calculated according to the above formula is that, compared with the initial RAAN-AoL visibility condition, the RAAN condition in the RAAN-AoL visibility condition at time t is shifted with the Earth rotation, and the AoL visibility condition is unchanged. For example, Figure 5 The figure shows the five-orbit visibility range of a 500km, 30°inclination orbit satellite for a ground station with latitude 30°. The visibility range in the figure is discretely plotted with a period of one orbit, and the innumerable similar visibility range curves between the adjacent two visibility ranges are omitted. Due to the influence of the Earth rotation, the visibility range is constantly shifted to the right.

[0133] 4.5.5 Visibility judgment of the current satellite s for the ground target n at time t: If the RAAN value Ω t,p,s (computed according to step 4.5.3) of the satellite s at the current time t is within the RAAN visibility condition (computed according to step 4.5.4), and the AoL value u t,p,s (computed according to step 4.5.3) is within the AoL visibility condition (u1, u2) (computed according to step 4.5.4), then it is recorded that the current satellite s is visible for the ground target n at time t, and AccFlag n p,s (t) = 1, otherwise it is recorded as invisible, and AccFlag n p,s (t) = 0.

[0134] 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 visibility period j, go to step 4.5.3, otherwise go to step 4.5.6.

[0135] 4.5.6 Arrange the visibility flags of the satellite s in the orbit plane p to the ground target n at each time instant within the jth maximum concentrated visible time period, form a list of visible time instants, denoted as l n acc,p,s (j) = [t1, t2,... t n ], t e Γ n Ω0,p,s (j).

[0136] 4.5.7 Let the maximum concentrated visible time period index j = j + 1. If j < J, go to step 4.5.2, otherwise go to step 4.5.8.

[0137] 4.5.8 Arrange the J lists of visible time instants of the satellite s in the orbit plane p to the ground target n in time order, and integrate to form a list of visible time instants of the satellite s in the orbit plane p to the ground target n within the mission period T, denoted as l n acc,p,s = [t1, t2,... t n ], t e Γ n Ω0,p,s .

[0138] 4.6 Let the satellite index s = s + 1. If s < S p , go to step 4.4, otherwise go to step 4.7.

[0139] 4.7 The list of visible time periods of all satellites in the orbit plane p to the ground target n is calculated: The S lists of visible time instants of all satellites in the orbit plane p to the ground target n within the mission period T, denoted as l n acc,p,s , are combined, where all the visible time instants are arranged in time order, to form a list of visible time instants of all satellites in the orbit plane p to the ground target n within the mission period T, denoted as l n acc,p :

[0140] The consecutive time instants in l n acc,p are combined and arranged in the order of the start time of the visible time period, to form a list of visible time periods L n acc,p :

[0141] where each row represents a visible time period, and there are Q n p visible time periods in total, 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 can contain one or more satellites.

[0142] 4.8 The coverage index of all satellites of orbit plane p to ground target n is calculated: The coverage index of all satellites of orbit plane p to ground target n in mission cycle T is denoted as ε n p According to the focus of the mission, ε n p The average visibility duration or average revisit time can be used as the index.

[0143] The average visibility duration index τ n p Specifically:

[0144] In the formula, q represents the sequence number of the visibility period, Q n p represents the total number of visibility periods, endT q represents the end time of the qth visibility period, begT q represents the start time of the qth visibility period.

[0145] The average revisit time index Δt n p Specifically:

[0146] In the formula, q represents the sequence number of the visibility period, Q n p represents the total number of visibility periods, endT q-1 represents the end time of the qth visibility period, begT q represents the start time of the qth visibility period.

[0147] The coverage index consistency of all satellites of orbit plane p to ground target n in mission cycle T is denoted as:

[0148] In the formula, ε n’ p,q represents the coverage index value of each visibility period, specifically: when the coverage index is selected as the visibility duration, ; when the coverage index is selected as the revisit time, Var n p represents the coverage index consistency.

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

[0150] 4.10 The coverage index of orbit plane p to all ground targets is calculated: The average coverage index is: ε p = (1 / N) ·∑ N n=1 (w · ε n p ) ; The average coverage index consistency index is: Var p = (1 / N) ·∑ N n=1 (w · Var n p ) ; Wherein, w represents a weight coefficient.

[0151] 4.11 Working orbit parameter optimization modeling of orbital plane p: The optimization model of working orbit parameters of all satellites in the orbital plane p is: the optimization variables are the initial semi-major axis a s , s ∈ [1, S p ] and initial AoL, s ∈ [2, S p ], a total of 2S p -1, wherein the AoL of the satellite numbered 1 is defined as 0. The optimization objective is to minimize the average coverage index ε p and the average coverage index consistency index Var p , and the constraint is the semi-major axis adjustment range, which is represented as:

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

[0153] Step 5: Working semi-major axis and latitude argument correction value solving Iterate through the orbital planes 1-P, and solve the semi-major axis and AoL one by one for each orbital plane p. The initial working semi-major axis a s , s ∈ [1, S p ] and initial AoL, s ∈ [2, S p ] of each satellite in the orbital plane p are taken as variables to be corrected, and an analytical method or an optimization algorithm is used to solve.

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

[0155] 5.2 Solve the working orbit parameters of the orbital plane p using an analytical method: 5.2.1 Calculate the initial working semi-major axis of all satellites in the orbital plane p: For each satellite on the orbit plane p, the following operations are performed: in the query table 1, find the case where the face quality ratio is closest to the face quality ratio 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 its orbit decay parameters k1, k2. Calculate the initial working semi-major axis of each satellite on the orbit plane p: a s s T 2 ) / 3-(k2 s T) / 2, s∈[1,S p ] In the formula, a s represents the initial working semi-major axis of the satellite s, a represents the nominal semi-major axis of the satellite s, k1 represents the quadratic semi-major axis decay parameter of the satellite s, k2 represents the linear semi-major axis decay parameter of the satellite s, T represents the mission period of the satellite s, and S p represents the total number of satellites on the orbit plane p where the satellite s is located.

[0156] 5.2.2 Calculate the nominal AoL of all satellites on the orbit plane p: u s '=[2(s-1)π] / S p , s∈[1,S p ] In the formula, u s ' represents the nominal AoL of the satellite s.

[0157] 5.2.3 Calculate the AoL offset of all satellites on the orbit plane p: When the initial semi-major axis of the satellite is a s , the maximum AoL offset △u max,s of the satellite s in T period is obtained by the following formula: △u max =[(7λ’+3n e ) / (12a)]·[-2k1·T m 3 -3k2·T m 2 +(2k1·T 2 +3k2·T)·T m ] λ’=[(3n e ·J2·R e 2 ) / (2a 2 )]·(4cos 2 i-1) T m =-k2 / (2k1)-[3k2 2 +2k1·(2k1·T 2 ​+ 3k2-T) 1 / 2 (2 x 3 1 / 2 x k1) wherein, △u max represents the maximum AoL offset of the satellite, λ' represents the long-term change rate of the mean anomaly, n e represents the orbital angular velocity, a represents the nominal semi-major axis of the satellite, k1 represents the quadratic semi-major axis decay parameter of the satellite, k2 represents the linear semi-major axis decay parameter of the satellite, T m represents the time required to reach the maximum AoL offset, J2 represents the Earth oblateness perturbation coefficient, R e represents the Earth radius, and i represents the orbital inclination.

[0158] 5.2.4 For each satellite s in the orbital plane p, the following operation is performed: in the data of the query table 2, the case in which the sensor 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 is found, and the covered central angle θ s is recorded.

[0159] 5.2.5 The initial AoL of each satellite in the orbital plane p is calculated as follows: u s = u s ' - △u max,s / 2, s ∈ [1, S p ] wherein, u s represents the initial AoL of the satellite s, u s ' represents the nominal AoL of the satellite s, and △u max,s represents the maximum AoL offset of the satellite s.

[0160] 5.2.6 If all satellites in the orbital plane p satisfy: △u max,s ≤ 2θ s , s ∈ [1, S p ] wherein, △u max,s represents the maximum AoL offset of the satellite s, and θ s represents the covered central angle of the satellite s.

[0161] Step 5.3 is entered, otherwise Step 5.4 is entered.

[0162] The reason for using the analytical method to correct the working orbit parameters is that the virtual satellite at the nominal AoL position of each satellite is taken as the reference point of each satellite, the relative AoL change of the satellite in time is set to 0, and the initial working orbit semi-major axis of each satellite is obtained. In the t period, the satellite moves approximately symmetrically relative to the reference point, and the reference Figure 6The figure shows the relative motion trajectories of the three satellites S1, S2, and S3 within one period T. When the sub-satellite point trajectory just passes through the ground target, the satellite's AoL deviation is ±θ, and the satellite can still complete coverage at the specified time, with no impact on the drift index. When the sub-satellite point trajectory crosses the ground coverage circle but does not pass through the ground target, the redundant deviation is less than ±θ. Therefore, when the maximum AoL offset of satellite s is △u max,s ≤2θ s When , the above analytical method can be used to calculate the approximate optimal solution of each satellite's orbital parameters to reduce the computational complexity.

[0163] Using this method, it is necessary to adjust the semi-major axis of each satellite to a at the end of each mission cycle T. s .

[0164] 5.3 Calculate the coverage index for orbital plane p during mission period T according to steps 4.1 through 4.10. The initial operating semi-major axis and initial AoL for each satellite are the values ​​calculated in step 5.2. If the coverage index for the orbital plane meets the expectations, use the initial operating semi-major axis obtained in step 5.2.1 and the initial AoL obtained in step 5.2.5 as the correction values ​​for the operating semi-major axis and the argument of latitude AoL, respectively, and proceed to step 5.5. Otherwise, proceed to step 5.4.

[0165] Among them, coverage indicators meet expectations means that the coverage indicators meet the coverage performance requirements in the task design phase.

[0166] 5.4 Use the optimization algorithm to solve the working orbit parameters of the orbital surface p: Use genetic algorithm or other intelligent optimization algorithms to solve the orbital parameter optimization model of the orbital plane p in step 4.11. s , s∈[1,S p ] and the initial AoL value u s , s∈[2,S p ], a total of 2S p -1 parameter as the working orbital parameter of orbital plane p, and ultimately obtain the working semi-major axis correction value and latitude argument AoL correction value of each satellite. In each iteration, for each individual, according to steps 4.1 to 4.10, calculate the coverage index and coverage index consistency index of orbital plane p within the mission period T.

[0167] The reason for using the optimization algorithm to modify the working orbit parameters is that when the conditions for using the analytical method are not met, or the coverage index obtained by using the analytical method is not ideal, all 2S p -1 working orbit parameter as the optimization variable and the intelligent optimization algorithm can greatly increase the solution space, thereby obtaining a better optimal solution and improving the coverage index.

[0168] With this method, it is required to recalculate the working orbit parameters of each satellite and reconstruct at the end of each task period T with step 5.4.

[0169] 5.5 Let the orbit plane index p = p + 1. If p ≤ P, go to step 5.2, otherwise go to step 6.

[0170] Step 6: Working RAAN correction The reason for the constellation RAAN correction is that different semi-major axes are designed for satellites in the same orbit to ensure the stability of the relative AoL and transit time in the plane. Under the action of perturbation, the differentiated semi-major axes will cause the dispersion of the RAAN of each satellite in the same orbit. Therefore, the actual working RAAN is corrected on the basis of the nominal RAAN of each orbit to resist the RAAN drift caused by the semi-major axis. The RAAN correction value of the orbit plane 1~P is calculated according to the following steps:

[0171] 6.1 Let the orbit plane index p = 1. 6.2 Let the satellite index s = 1.

[0172] 6.3 In the lookup table 1, find the plane quality ratio closest to the plane quality ratio of the current satellite s, and the characteristic semi-major axis and the initial working semi-major axis (working semi-major axis correction value) a s of the current satellite s after optimization.

[0173] The closest case, record its orbit decay parameters k1, k2. Calculate the RAAN change amount of the orbit plane p satellite s in the task period T p,s : △Ω p,s =-(7Ω’·k1·T 3 ) / (6a)-(7Ω’·k2·T 2 ) / (4a)-[(7Ω’·T) / (2a)]·(a s -a) Ω’=-[(3n e ·J2·R e 2 ) / (2a 2 )]·cosi In the formula, △Ω p,s represents the RAAN change amount of the orbit plane p satellite s in the task period, Ω’ represents the long-term change rate of RAAN, k1 represents the quadratic term semi-major axis decay parameter, k2 represents the linear term semi-major axis decay parameter, T represents the task period, a represents the nominal semi-major axis of the satellite, a s represents the initial working semi-major axis after optimization, i.e. the working semi-major axis correction value, n e represents the orbital angular velocity, J2 represents the Earth oblateness perturbation coefficient, R eR represents the radius of the earth, i represents the orbit inclination.

[0174] 6.4 Let the satellite number s = s + 1. If s ≤ S p , go to step 6.3, otherwise go to step 6.5.

[0175] 6.5 Calculate the average RAAN change of all satellites in orbit plane p over the mission period T:

[0176] where ΔΩ p represents the average RAAN change of all satellites in orbit plane p over the mission period T, and ΔΩ p,s represents the RAAN change of satellite s in orbit plane p over the mission period T.

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

[0178] 6.7 Let the orbit plane number p = 1.

[0179] 6.8 Let the satellite number s = 1.

[0180] 6.9 Update the working RAAN parameter of satellite s in orbit plane p: Ω p,s = Ω p - ΔΩ p,s .

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

[0182] 6.11 Let the orbit plane number p = p + 1. If p ≤ P, go to step 6.8, otherwise end the flow.

[0183] Embodiment 2 This embodiment carries out numerical simulation verification of the orbit constellation under the condition of an orbit height of 500 km and an orbit inclination of 30°. Compared with the traditional walker configuration, the optimized configuration obtained by the embodiment reduces the non-uniform coverage by about 17%, greatly increases the consistency of the coverage index represented by revisit time, and solves the problem of coverage index drift in the heterogeneous constellation. In order to achieve this goal, the control frequency of the traditional configuration needs to be increased by 3-6 times, or the number of satellites needs to be increased by about 30-40%.

[0184] The above detailed description is intended to explain and describe the application, but not to limit the application. Any modification and change within the spirit and scope of the application will be included in the scope of the application.

Claims

1. A method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability, characterized in that: The following steps are involved: Step S1: constructing a first query table and a second query table; Step S2: For each orbital plane, the operating orbit parameter optimization method is used, combined with the first lookup table and the second lookup table, to analyze and calculate the dynamic visibility conditions and visibility list, and then correct the operating semi-major axis and latitude argument AoL of each satellite on the orbital plane; Step S3: For each orbital plane, the RAAN correction value of the ascending node is calculated based on the semi-major axis correction value of each satellite on the orbital plane and the first lookup table; Step S4: The working semi-major axis correction value, the argument of latitude AoL correction value and the right ascension of ascending node RAAN correction value of each satellite in all orbital planes are used as the configuration parameters of the heterogeneous satellite constellation.

2. The method for correcting a heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 1, characterized in that: The first lookup table contains the semi-major axis attenuation parameters and RAAN change rates for satellites with different area-to-mass ratios at each characteristic semi-major axis. The second lookup table contains the geocentric coverage angles for satellites with different coverage half angles at each characteristic semi-major axis, as well as the initial RAAN-AoL visibility conditions and initial RAAN visibility characteristic range for each ground target.

3. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 1, characterized in that: In step S2, the working track parameter optimization method includes: Step S2.1: For each satellite on the orbital plane, the initial operating semi-major axis, initial AoL, geocentric angle of coverage, and maximum AoL offset of the satellite are obtained by analytically combining the first and second lookup tables based on the satellite's area-to-mass ratio, half-angle of coverage, and nominal semi-major axis. Step S2.2: If the maximum AoL offset of each satellite on the orbital plane is less than or equal to twice the coverage angle of the Earth's center, proceed to step S2.3; otherwise, proceed to step S2.5; Step S2.3: Based on the initial semi-major axes and initial AoLs of all satellites on the orbital plane, the coverage index of the orbital plane during the mission period is obtained by analytically calculating the dynamic visibility conditions and the visibility list; Step S2.4: If the orbital plane coverage index meets the preset conditions, the initial working semi-major axis and initial AoL obtained in step S2.1 are used as the satellite's working semi-major axis correction value and latitude argument AoL correction value, respectively; otherwise, proceed to step S2.5; 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 the variables to be optimized, and minimizing the coverage index of the orbital plane within the mission cycle is used as the optimization goal. During 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 latitude angle AoL correction value of each satellite on the orbital plane are obtained.

4. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 3, characterized in that: The step S2.3 includes: Step S2.3.1: For each satellite on the orbital plane, based on the satellite's area-to-mass ratio, half-angle of coverage, and initial operating semi-major axis, obtain the satellite's RAAN change rate, geocentric angle of coverage, and initial RAAN visibility conditions and initial RAAN visibility feature range for each ground target from the first and second lookup tables, thereby obtaining the maximum concentrated visibility period of the satellite for each ground target. The maximum concentrated visible period of the satellite to each ground target is expressed as: ； Where, Γ Ω0,p,s It represents the maximum concentrated visible period of satellite s to ground target n on orbital plane p; A1 represents the initial offset coefficient, A2 represents the long-term drift proportional coefficient, t Γ1 or t Γ2 Indicates 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 coverage geocentric angle of the satellite s, and i represents the orbital inclination of the orbital plane p; Step S2.3.2: For each satellite on the orbital plane, obtain a list of satellite visibility times for each ground target based on the satellite's maximum concentrated visibility period for each ground target and the initial AoL; Step S2.3.3: Merge the visibility time lists of all satellites in the orbital plane for the same ground target to obtain a list of visibility time periods of all satellites in the orbital plane for the ground target; Step S2.3.4: Based on the list of visible periods of all satellites in the orbital plane for each ground target, calculate the coverage index and coverage index consistency of all satellites in the orbital plane for each ground target; 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 to obtain the coverage index of the orbital plane.

5. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 4, characterized in that: The step S2.3.2 includes: Obtaining the orbital attenuation coefficient, the RAAN change rate, and the initial RAAN-AoL visibility condition of the satellite to the ground target from the first lookup table and the second lookup table based on the satellite's area-to-mass ratio, coverage half angle, and initial working semi-major axis; Traverse each moment in the maximum concentrated visibility period and obtain the satellite's RAAN value and AoL value at the current moment based on the satellite's orbital attenuation coefficient, RAAN change rate, and initial AoL. Map the satellite's initial RAAN-AoL visibility condition for ground targets to the current moment to obtain the satellite's RAAN-AoL visibility condition for ground targets at the current moment. Determine visibility based on the satellite's RAAN value, AoL value, and RAAN-AoL visibility condition 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 target at the current moment; otherwise, the satellite is not visible to the ground target at the current moment. The visibility judgment results of all moments in the maximum concentrated visible period are integrated to obtain a list of visible moments of the satellite to ground targets.

6. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 4, characterized in that: The coverage index of a single ground target by all satellites in the orbital plane is the average visibility time or average revisit time of all satellites in the orbital plane to the ground target.

7. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 1, characterized in that: The step S3 comprises the following steps: For each satellite on the orbital plane, the orbital attenuation parameter is obtained from the first lookup table based on the satellite's semi-major axis correction value and the surface mass ratio, and the RAAN change of the satellite during the mission cycle is calculated. Calculate the average RAAN change of all satellites in the orbital plane during the mission cycle; Obtain the difference between the nominal RAAN and the average RAAN change of the orbital surface to obtain the RAAN correction value of the orbital surface; For each satellite on the orbital plane, the RAAN correction value of the satellite's ascending node right ascension is obtained based on the RAAN correction value of the orbital plane and the RAAN change of the satellite during the mission cycle.

8. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 1, characterized in that: The step S1 comprises: First, coverage data and orbital drift data under different characteristic semi-major axes are obtained; the coverage data includes the coverage geocentric angle of satellites with different coverage half-angles under each characteristic semi-major axis; the orbital drift data includes the semi-major axis attenuation parameter and RAAN change rate under each characteristic semi-major axis for satellites with different surface-to-mass ratios, wherein the semi-major axis attenuation parameter includes a quadratic semi-major axis attenuation parameter and a linear semi-major axis attenuation parameter; Finally, combined with the satellite-ground visibility conditions, the initial RAAN-AoL visibility condition set and the initial RAAN visibility feature range set of each ground target on the satellite are obtained.

9. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 8, characterized in that: In step S1, the satellite's initial RAAN-AoL visibility condition and initial RAAN visibility feature range for each ground target are calculated by the following process: Based on the characteristic semi-major axis and coverage half angle of the satellite, one or two sets of initial RAAN visibility conditions for the satellite to the ground target are calculated. Each set of initial RAAN visibility conditions mainly consists of boundary conditions, center point and visibility range. For each set of initial RAAN visibility conditions, several RAAN feature points are generated through discretization based on the boundary conditions and the center point. The two AoL boundary points corresponding to each RAAN feature point are calculated. A polynomial fit is performed on all RAAN feature points and the corresponding AoL boundary points to obtain the initial RAAN-AoL visibility conditions.

10. The method for correcting heterogeneous satellite constellation configuration based on ground coverage performance stability according to claim 9, characterized in that: When the satellite satisfies φ for the ground target n n <i - θ, the satellite has two sets of RAAN visibility conditions at the initial time for the ground target n. The RAAN difference between the two boundary points in the two sets of RAAN visibility conditions at the initial time is the same, but the values are different; when the satellite satisfies i - θ < φ for the ground target n n <θ, the satellite has one set of RAAN visibility conditions at the initial time.

Images