A near-circular orbit analytical prediction method suitable for ultra-low earth orbit satellites

Abstract

The application provides a near-circular orbit analytical prediction method suitable for ultra-low earth orbit satellites, and comprises the following steps: considering the orbit environment of the ultra-low earth orbit satellite under the action of atmospheric resistance and J2 perturbation, constructing a nonlinear dynamics equation of the ultra-low earth orbit satellite orbit by taking the geocentric distance, the speed module value, the speed inclination angle, the orbit inclination angle, the ascending node right ascension and the latitude amplitude angle as state variables; after the current state variable and time are dimensionless, linearizing and expanding the nonlinear dynamics equation at the reference circular orbit, obtaining a linearized dynamics equation and solving the linearized dynamics equation to obtain a near-circular orbit analytical solution of the ultra-low earth orbit satellite; if the analytical solution satisfies a preset iteration termination condition, generating a prediction result of the near-circular orbit based on the analytical solution. The application fully considers the orbit characteristics of the ultra-low earth orbit satellite under the coupling action of atmospheric resistance and non-spherical gravity perturbation, and the scheme is simple and effective, and has a wide application prospect in the deployment of future ultra-low earth orbit satellites.

Description

Technical Field

[0001] This invention belongs to the field of orbital dynamics of ultra-low orbit satellites in aerospace technology, and specifically relates to an analytical prediction method for near-circular orbits of ultra-low orbit satellites. Background Technology

[0002] With the continuous evolution and deep integration of aerospace technology, human exploration of the aerospace field is constantly reaching new heights. Currently, a trend is emerging of "aviation reaching increasingly higher altitudes and spaceflight gradually moving lower altitudes": In aviation, near space (20-100km airspace), as a key region connecting aviation and aerospace, has become an international research hotspot, with its unique space environment and strategic value attracting widespread attention; in aerospace, ultra-low Earth orbit (usually referring to the space region with an orbital altitude below 300 km), as a special type of near-Earth space with high application potential, is rapidly emerging as a new research focus. This orbital region not only possesses significant advantages in remote sensing observation and communication enhancement, but also has the potential to develop into the "fifth airspace" after traditional aviation, aerospace, near space, and low altitude, providing crucial support for the integrated aerospace technology system and the expansion of future space applications.

[0003] Ultra-low Earth orbit (ULE) satellites orbit at significantly lower altitudes than conventional low Earth orbit (LEO) satellites, giving them superior spatial resolution and signal transmission performance for Earth observation and communication, thus possessing significant civilian and strategic value. Furthermore, developing ULE satellite systems also holds significant national security implications for enhancing space situational awareness and addressing potential competition for frequency and orbital resources and security threats from foreign mega-constellations (such as the US Starlink).

[0004] Although research on ultra-low Earth orbit (ULE) satellites is increasing, current work mainly focuses on engine design and aerodynamic-thermodynamic environment analysis, while research on their orbital dynamics and control remains relatively weak. ULE satellites experience significant atmospheric drag during their operation and are also subject to strong gravitational perturbations from the Earth's non-spherical shape (especially the J2 term). Some studies have attempted to model the orbital evolution under the coupled effects of atmospheric drag and J2 perturbations and derive corresponding analytical solutions, but these solutions are often complex in structure, have numerous parameters, and limited engineering applicability. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing technologies and propose an analytical prediction method for near-circular orbits of ultra-low Earth orbit (ULE) satellites, addressing the problem of accurate and rapid orbit prediction under the tight coupling of aerodynamic forces and non-spherical J2 gravitational perturbations. This invention is simple in form, efficient in its solution process, and combines good prediction accuracy with engineering practicality, enabling direct application to future ULE satellite orbit deployment, mission planning, and on-orbit management.

[0006] This invention proposes a near-circular orbit analytical prediction method suitable for ultra-low Earth orbit satellites, comprising:

[0007] Considering the atmospheric drag and J2 perturbation effect of the ultra-low orbit satellite orbit environment, a nonlinear dynamic equation for the ultra-low orbit satellite orbit is constructed using the geocentric distance, velocity modulus, velocity inclination angle, orbital inclination angle, right ascension of the ascending node, and latitudinal argument as state variables.

[0008] After dimensionless processing of the current state variables and time, the nonlinear dynamic equations are linearized and expanded at the reference circular orbit to obtain linearized dynamic equations.

[0009] Solving the linearized dynamic equations yields an analytical solution for the near-circular orbit of the ultra-low Earth orbit satellite. If the analytical solution satisfies a preset iteration termination condition, a prediction result for the near-circular orbit is generated based on the analytical solution.

[0010] In one specific embodiment of the present invention, the expression of the nonlinear dynamic equation is as follows:

[0011] (1)

[0012] In the formula, The distance from the Earth's center. For velocity modulus, The velocity tilt angle, For the track inclination angle, Right ascension of the ascending node, Argument of latitude; The gravitational constant of Earth, Atmospheric density, It is the product of the satellite's surface mass ratio and drag coefficient. , , and These represent the long-term term, first-harmonic sine coefficient, second-harmonic cosine coefficient, and second-harmonic sine coefficient corresponding to the non-spherical gravitational perturbation of Earth J2, respectively. This is a correction term for the Earth's rotation, expressed as the square of the velocity. This is a fundamental correction term for Earth's rotation.

[0013] parameter , , , , , and The calculation expression is:

[0014] (2)

[0015] (3)

[0016] (4)

[0017] In the formula, For the satellite's windward area, For satellite quality, The drag coefficient, The value is the J2 perturbation value. For the Earth's radius, The angular velocity of Earth's rotation. The relative velocity of the satellite with respect to the rotating Earth;

[0018] (5).

[0019] In one specific embodiment of the present invention, it further includes:

[0020] Let the number of iterations Set the initial value of the reference radius for the reference circular orbit during iteration. Initial values ​​for reference orbit inclination angle iteration Set reference speed tilt angle It is always zero;

[0021] Among them, the initial time is recorded. The initial value of the distance from the Earth's center is The initial value of the velocity modulus is The initial value of the velocity tilt angle is The initial value of the orbital inclination angle is The initial value of the right ascension of the ascending node is The initial value of the latitude argument is ;

[0022] The initial value of the reference radius for the reference circular orbit is calculated using the following formula. Initial values ​​for the inclination angle of the reference circular orbit. :

[0023] (6).

[0024] In a specific embodiment of the present invention, the dimensionless processing of the current state quantity and time includes:

[0025] Reference radius based on the reference circle orbit in the current iteration round and reference orbital inclination Calculate the reference velocity modulus :

[0026] (7)

[0027] In the formula, This is a reference value for the long-term term corresponding to the non-spherical gravitational perturbation of Earth J2;

[0028] The distance from the Earth's center, the velocity magnitude, the velocity tilt angle, and time are dimensionless and expressed as follows:

[0029] (8)

[0030] in, , and Let A, B, and C be the dimensionless distance from the geocenter, velocity modulus, and velocity tilt angle, respectively, with initial values ​​of [missing values]. , and ; The time is dimensionless, and the initial value is... .

[0031] In one specific embodiment of the present invention, obtaining the linearized dynamic equation includes:

[0032] Linearize equation (1) around the reference circular orbit:

[0033] (9)

[0034] (10)

[0035] In the formula, This is the dimensionless atmospheric density term. This is the dimensionless atmospheric density gradient term. This is a dimensionless atmospheric drag correction term. The J2 perturbation correction term is dimensionless. This is a dimensionless correction term for Earth's rotation. The coefficients of the second harmonic cosine term, which are dimensionless derivatives of the velocity modulus, are... The coefficients of the second harmonic sine term, which is a dimensionless derivative of the velocity modulus, are... For the dimensionless velocity tilt angle derivative, the second harmonic cosine term coefficient;

[0036] (11)

[0037] In the formula, For reference density; Reference density For reference radius The partial derivative;

[0038] parameter , , The definition is as follows:

[0039] (12)

[0040] In the formula, and These are the Earth rotation correction terms for the square of the velocity. The corresponding long-term and second-harmonic cosine coefficients, and These are the fundamental correction terms for Earth's rotation. The corresponding long-term and second-harmonic cosine coefficients;

[0041] (13)

[0042] In the formula, , and These are the coefficients of the long-term term, the second harmonic cosine term, and the second harmonic sine term of the dimensionless orbital inclination derivative, respectively. , and These are the long-term term, second-harmonic cosine term coefficient, and second-harmonic sine term coefficient of the dimensionless right ascension derivative of the ascending node, respectively. , and These are the long-term term, second-harmonic cosine term coefficient, and second-harmonic sine term coefficient of the dimensionless latitudinal argument derivative, respectively.

[0043] (14)

[0044] In the formula, For the fundamental correction of Earth's rotation The reference value satisfies:

[0045] (15).

[0046] In one specific embodiment of the present invention, solving the linearized dynamic equation includes:

[0047] The solutions to differential equations (9) and (10) are in the following forms:

[0048] (16)

[0049] (17)

[0050] In the formula, coordinates The constant term in the analytical solution, coordinates The coefficient of the long-term variation term in the analytical solution. coordinates The coefficients of the first harmonic cosine term in the analytical solution, coordinates The coefficients of the first harmonic sine term in the analytical solution, coordinates The coefficients of the second harmonic cosine term in the analytical solution, coordinates The coefficients of the second harmonic sine term in the analytical solution; The orbital altitude attenuation coefficient is... This is the orbital ellipticity attenuation coefficient; It is the average latitude argument;

[0051] in:

[0052] parameter and satisfy:

[0053] (18)

[0054] For coordinates Pick , and In this case, the corresponding coefficient and The solution expression is as follows:

[0055] The intermediate variable is obtained by solving equation (19). and Then, using equation (20) to obtain and Values;

[0056] (19)

[0057] (20)

[0058] Then, the dimensionless average geocentric distance is obtained using equation (21). initial value Dimensionless average velocity modulus initial value and dimensionless average velocity tilt angle initial value ;

[0059] (twenty one)

[0060] Further, intermediate variables are obtained by solving equation (22). , , and Then, using equation (23) to obtain , , and ;

[0061] (twenty two)

[0062] (twenty three)

[0063] For coordinates Pick , and In this case, the coefficient , , and The calculation expression is as follows:

[0064] (twenty four)

[0065] In the formula, , and These are the average orbital inclination angles. Mean ascending node right ascension and average latitudinal arc Initial values:

[0066] (25).

[0067] In one specific embodiment of the present invention, it further includes:

[0068] Determine whether the analytical solution for the near-circular orbit of a very low Earth orbit satellite satisfies the following equation:

[0069] (27)

[0070] In the formula, The set threshold;

[0071] If equation (27) is not satisfied, then the reference radius and inclination angle of the reference circular orbit are updated using the following equation:

[0072] (28)

[0073] In the formula, Let be the reference radius of the reference circular orbit in the (j+1)th iteration. Let be the inclination angle of the reference circular orbit for the (j+1)th iteration; then, let the iteration number be... Then proceed to the next round of iteration to solve;

[0074] If criterion (27) is satisfied, then the analytical solution of the trajectory is output according to the results of equations (16) and (17) of the current iteration round. When outputting, the state variables are included. , , The dimensionless process is performed as follows:

[0075] (29)

[0076] Then, orbital analysis prediction is performed based on the results of equation (29).

[0077] Features and beneficial effects of the present invention:

[0078] 1. This invention abandons the classic six-element description method and innovatively uses geocentric distance, velocity modulus, velocity inclination angle, orbital inclination angle, right ascension of the ascending node, and latitudinal argument as state parameters. This description method effectively avoids the singularity problem of near-circular orbits, and because it directly includes geocentric distance and velocity modulus, it is more direct and convenient to calculate the effect of atmospheric drag, especially suitable for constructing dynamic models of ultra-low orbit satellites.

[0079] 2. This invention linearizes the nonlinear dynamic equations near a reference circular orbit, resulting in a clearly structured linearized equation. This equation has well-defined parameters and intuitive physical meaning, facilitating the analysis of the impact of various perturbations on orbit evolution. Simultaneously, the linear model significantly simplifies the controller design process, providing convenience for the design of ultra-low Earth orbit satellite control systems.

[0080] 3. The analytical solution derived in this invention is concise in form and has a clear physical meaning. It has high accuracy in short-term orbit prediction and can clearly reveal the orbital evolution law of ultra-low orbit satellites. It has good engineering application prospects and can provide an efficient and reliable theoretical tool for orbit prediction and control tasks. Attached Figure Description

[0081] Figure 1 This is a flowchart of a near-circular orbit analytical prediction method for ultra-low Earth orbit satellites according to an embodiment of the present invention.

[0082] Figure 2 This is a comparison chart of the results of the analytical prediction method and the numerical integration method in a specific embodiment of the present invention.

[0083] Figure 3 This is a graph showing the error variation between the analytical prediction method and the numerical integration method in a specific embodiment of the present invention. Detailed Implementation

[0084] This invention proposes an analytical prediction method for near-circular orbits of ultra-low Earth orbit satellites, which is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0085] This invention proposes a near-circular orbit analytical prediction method suitable for ultra-low Earth orbit satellites, comprising:

[0086] Considering the atmospheric drag and J2 perturbation effect of the ultra-low orbit satellite orbit environment, a nonlinear dynamic equation for the ultra-low orbit satellite orbit is constructed using the geocentric distance, velocity modulus, velocity inclination angle, orbital inclination angle, right ascension of the ascending node, and latitudinal argument as state variables.

[0087] After dimensionless processing of the current state variables and time, the nonlinear dynamic equations are linearized and expanded at the reference circular orbit to obtain linearized dynamic equations.

[0088] Solving the linearized dynamic equations yields an analytical solution for the near-circular orbit of the ultra-low Earth orbit satellite. If the analytical solution satisfies a preset iteration termination condition, a prediction result for the near-circular orbit is generated based on the analytical solution.

[0089] In a specific embodiment of the present invention, the overall process of the near-circular orbit analytical prediction method suitable for ultra-low orbit satellites is as follows: Figure 1 As shown, it includes the following steps:

[0090] 1) Considering the orbital environment of ultra-low orbit satellites under atmospheric drag and J2 perturbation, a nonlinear dynamic equation is constructed with the geocentric distance, velocity modulus, velocity inclination angle, orbital inclination angle, right ascension of the ascending node, and latitudinal argument as orbital state variables.

[0091] In this embodiment, atmospheric drag is taken into account the effect of the rotation of the central celestial body.

[0092] This embodiment uses the distance from the center of the earth. Velocity modulus Velocity tilt angle Track inclination Right ascension of ascending node and latitude angle To describe the orbital state of a satellite, the nonlinear dynamic equation for the orbit of a very low Earth orbit satellite is:

[0093] (1)

[0094] In the formula, The gravitational constant of Earth, Atmospheric density, It is the product of the satellite's surface mass ratio and drag coefficient. , , and These represent the long-term term, first-harmonic sine coefficient, second-harmonic cosine coefficient, and second-harmonic sine coefficient corresponding to the non-spherical gravitational perturbation of Earth J2, respectively. This is a correction term for the Earth's rotation, expressed as the square of the velocity. This is a basic correction term for Earth's rotation.

[0095] parameter , , , , , and The specific calculation expression is as follows:

[0096] (2)

[0097] (3)

[0098] (4)

[0099] In the formula, For the satellite's windward area, For satellite quality, The drag coefficient, The value is the J2 perturbation value. For the Earth's radius, The angular velocity of Earth's rotation. This represents the satellite's relative velocity to the rotating Earth.

[0100] (5)

[0101] 2) Let the number of iterations be... Set the initial value of the reference radius for the reference circular orbit during iteration. Initial values ​​for reference orbit inclination angle iteration Set reference speed tilt angle It is always zero.

[0102] In this embodiment, the initial time is recorded. The initial value of the distance from the Earth's center is The initial value of the velocity modulus is The initial value of the velocity tilt angle is The initial value of the orbital inclination angle is The initial value of the right ascension of the ascending node is The initial value of the latitude argument is The above values ​​can be obtained by converting orbit determination data or numerical simulation data. The reference radius of the reference circular orbit is calculated using the following formula. initial values ​​of iteration and reference orbital inclination initial value of iteration :

[0103] (6)

[0104] 3) Based on the results of step 2), calculate the reference velocity modulus. :

[0105] (7)

[0106] In the formula, This is a reference value for the long-term term corresponding to the non-spherical gravitational perturbation of Earth J2. In this embodiment, the value in equation (3) is used. The calculation expression, when calculated, needs to include the expression. and Replace with reference values ​​respectively and It can be calculated .

[0107] 4) The initial values ​​of distance from the center of the earth, velocity modulus, velocity tilt angle and time are dimensionless.

[0108] In this embodiment, the state variables , , and time After dimensionless processing, the expression is as follows:

[0109] (8)

[0110] in, , and Let A, B, and C be the dimensionless distance from the geocenter, velocity modulus, and velocity tilt angle, respectively, with initial values ​​of [missing values]. , and . For dimensionless time, the initial value is... .

[0111] 5) Based on the results of step 4), the nonlinear dynamic equations obtained in step 1) are linearized and expanded at the reference circular orbit to obtain the linearized dynamic equations.

[0112] In this embodiment, equation (1) is linearized and expanded near the reference circular orbit, resulting in the following form:

[0113] (9)

[0114] (10)

[0115] In the formula, This is the dimensionless atmospheric density term. This is the dimensionless atmospheric density gradient term. This is a dimensionless atmospheric drag correction term. The J2 perturbation correction term is dimensionless. This is a dimensionless correction term for Earth's rotation. The coefficients of the second harmonic cosine term, which are dimensionless derivatives of the velocity modulus, are... The coefficients of the second harmonic sine term, which is a dimensionless derivative of the velocity modulus, are... For the dimensionless velocity tilt angle derivative, the second harmonic cosine term coefficient.

[0116] (11)

[0117] In the formula, Reference density, i.e., reference radius The corresponding density value; Reference density For reference radius The partial derivatives, and all other physical quantities marked with "*" represent their reference values. The parameters appearing within... , , The definition is as follows:

[0118] (12)

[0119] In the formula, and These are the Earth rotation correction terms for the square of the velocity. The corresponding long-term and second-harmonic cosine coefficients, and These are the fundamental correction terms for Earth's rotation. The corresponding long-term and second-harmonic cosine coefficients.

[0120] (13)

[0121] In the formula, , and These are the coefficients of the long-term term, the second harmonic cosine term, and the second harmonic sine term of the dimensionless orbital inclination derivative, respectively. , and These are the long-term term, second-harmonic cosine term coefficient, and second-harmonic sine term coefficient of the dimensionless right ascension derivative of the ascending node, respectively. , and These are the long-term term, second harmonic cosine term coefficient, and second harmonic sine term coefficient of the dimensionless latitudinal argument derivative, respectively.

[0122] (14)

[0123] In the formula, For the fundamental correction of Earth's rotation The reference value satisfies:

[0124] (15)

[0125] 6) Solve the linearized dynamic equations obtained in step 5) to obtain the analytical solution for the near-circular orbit of the ultra-low orbit satellite.

[0126] In this embodiment, the solutions to differential equations (9) and (10) are in the following forms:

[0127] (16)

[0128] (17)

[0129] In the formula, coordinates The constant term in the analytical solution, coordinates The coefficient of the long-term variation term in the analytical solution. coordinates The coefficients of the first harmonic cosine term in the analytical solution, coordinates The coefficients of the first harmonic sine term in the analytical solution, coordinates The coefficients of the second harmonic cosine term in the analytical solution, coordinates The coefficients of the second harmonic sine term in the analytical solution, where the coordinates... Desirable ; The orbital altitude attenuation coefficient is... This is the orbital ellipticity attenuation coefficient; This represents the average latitudinal argument. Wherein:

[0130] parameter and satisfy:

[0131] (18)

[0132] For coordinates Pick , and In this case, the corresponding coefficient and The solution expression is as follows:

[0133] The intermediate variable is obtained by solving equation (19). and Then, using equation (20) to obtain and Values.

[0134] (19)

[0135] (20)

[0136] Then, the dimensionless average geocentric distance is obtained using equation (21). initial value Dimensionless average velocity modulus initial value and dimensionless average velocity tilt angle initial value .

[0137] (twenty one)

[0138] Further, intermediate variables are obtained by solving equation (22). , , and Then, using equation (23) to obtain , , and .

[0139] (twenty two)

[0140] (twenty three)

[0141] For coordinates Pick , and In this case, the coefficient , , and The calculation expression is as follows:

[0142] (twenty four)

[0143] In the formula, , and These are the average orbital inclination angles. Mean ascending node right ascension and average latitudinal arc Initial values:

[0144] (25)

[0145] 7) Judgment:

[0146] If the initial reference circular orbit is chosen appropriately, it should satisfy the following:

[0147] (26)

[0148] Therefore, the following iterative criterion can be adopted:

[0149] (27)

[0150] In the formula, The threshold value is a small amount, for example, it can be taken as... .

[0151] If criterion (27) is not satisfied, then the reference radius and reference orbit inclination are updated using the following formula:

[0152] (28)

[0153] In the formula, Let be the reference radius of the reference circular orbit in the (j+1)th iteration. Let be the inclination angle of the reference circular orbit in the (j+1)th iteration. Then, let the iteration number... Return to step 3).

[0154] If criterion (27) is satisfied, then the analytical solution of the trajectory is output according to the results of equations (16) and (17) of the current iteration. When outputting, the state variables need to be included. , , The dimensionless process is performed as follows:

[0155] (29)

[0156] This allows for subsequent trajectory analysis and prediction.

[0157] Furthermore, in a specific embodiment of the present invention, the following simulation example is considered, with relevant parameter values ​​as shown in Table 1:

[0158] Table 1. Parameter value table in a specific embodiment of the present invention

[0159] Parameter name Value Parameter name Value <![CDATA[3.986012E14 m 3 / s 2 ]]> 2.2 6378145 m <![CDATA[0.5 m 2 ]]> 7.29211585E-5 rad / s <![CDATA[5.74E-7 kg / m 3 ]]> 0.0010826359 <![CDATA[4.51E-11 m -2 ]]> 20 kg <![CDATA[-4.78E-5 m -1 ]]>

[0160] Where parameters , and Used to calculate atmospheric density The expression is as follows:

[0161]

[0162] In the formula, The atmospheric density at sea level. The coefficient of the quadratic term for height in the density function. The coefficient of the first-order term of the density function is the height. Current altitude (using geocentric distance) Subtract Earth's radius get).

[0163] In this embodiment, the initial orbital values ​​are shown in Table 2:

[0164] Table 2. Initial orbital values ​​in a specific embodiment of the present invention.

[0165] Parameter name Value Parameter name Value 240 km 60° 7754.85 m / s 30° 0° 0°

[0166] The numerical integration results (using the Runge-Kutta method, with an integration step size of 1 s and a total integration time of 24 h) are compared with the analytical prediction method proposed in this patent, as follows: Figures 2 to 3 As shown.

[0167] Figure 2 This is a comparison chart of the results of the analytical prediction method and the numerical integration method in a specific embodiment of the present invention. Figure 2 In the diagram, the horizontal axis of all six subplots represents time. (Unit: h), the vertical axis corresponds to altitude. (Unit: km), velocity modulus (Unit: m / s), velocity tilt angle (Unit: deg) Track inclination (Unit: deg), Right Ascension of Ascending Node (Unit: deg) and latitudinal argument (Unit: deg), dotted lines represent numerical solutions, and dashed lines represent analytical solutions proposed in this patent. Sub-graph a shows the altitude. Over time The changes show an overall downward trend. The analytical solution agrees well in the first 18 hours, but a clear separation occurs after 18 hours. Subplot b shows the velocity magnitude. Over time The changes show an overall upward trend. The analytical solution matches well for the first 18 hours, but a clear separation occurs after 18 hours. Subplot c shows the velocity tilt angle. Over time The change in θ oscillated around 0°, with the amplitude showing a slow contraction trend. The analytical solution matched well throughout the entire process. The d-subplot shows the orbital inclination. Over time The changes show a slow downward trend overall, and the analytical solution matches well throughout the process; the e-subgraph shows the right ascension of the ascending node. Over time The changes show an overall downward trend, and the analytical solution matches well throughout; the f subgraph shows the latitudinal argument. Over time The change in θ, from 0° to 360°, is cyclical, and the analytical solution matches it well throughout the entire process.

[0168] Figure 3This is a graph showing the error variation between the analytical prediction method and the numerical integration method in a specific embodiment of the present invention. Figure 3 In the middle, the horizontal and vertical axes of the 6 subplots are all... Figure 2 The meaning is the same. Sub-graph a shows the altitude. Error over time The absolute value of the error showed an overall oscillating increasing trend, with the maximum error reaching 6 km within 24 hours; subplot b shows the velocity magnitude. Error over time The absolute value of the error showed an overall oscillating increasing trend, with the maximum error reaching -4 m / s within 24 hours; subplot c shows the velocity tilt angle. Error over time The error fluctuated around 0°, but the amplitude of the fluctuation gradually increased, reaching a maximum of 0.027° within 24 hours; the d-subplot shows the orbital inclination. Error over time The error fluctuates around 0° overall, but the amplitude of the fluctuation gradually increases, reaching a maximum of 0.014° within 24 hours; the e-subplot shows the right ascension of the ascending node. Error over time The absolute value of the error showed an overall oscillating increasing trend, with the maximum error reaching 0.05° within 24 hours; the f subplot shows the latitude argument. Error over time The absolute value of the error showed an overall increasing trend, with the maximum error reaching -18° within 24 hours.

[0169] Figure 2 and Figure 3 This indicates that although the method described in this embodiment can predict the changing trend of the actual orbital state relatively well in a short period of time and analytically explain the orbital change law of ultra-low orbit satellites, error divergence is inevitable. In practical applications, error divergence can be effectively controlled by limiting the single prediction duration of the analytical solution to a specific time and periodically re-updating it.

Claims

1. A method for analytical prediction of near-circular orbits suitable for ultra-low Earth orbit satellites, characterized in that, include: Considering the atmospheric drag and J2 perturbation effect of the ultra-low orbit satellite orbit environment, a nonlinear dynamic equation for the ultra-low orbit satellite orbit is constructed using the geocentric distance, velocity modulus, velocity inclination angle, orbital inclination angle, right ascension of the ascending node, and latitudinal argument as state variables. After dimensionless processing of the current state variables and time, the nonlinear dynamic equations are linearized and expanded at the reference circular orbit to obtain linearized dynamic equations. Solving the linearized dynamic equations yields an analytical solution for the near-circular orbit of the ultra-low Earth orbit satellite. If the analytical solution satisfies a preset iteration termination condition, a prediction result for the near-circular orbit is generated based on the analytical solution.

2. The method according to claim 1, characterized in that, The expression for the nonlinear dynamic equation is as follows: （1） In the formula, The distance from the Earth's center. For velocity modulus, The velocity tilt angle, For the track inclination angle, Right ascension of the ascending node, Argument of latitude; The gravitational constant of Earth, Atmospheric density, It is the product of the satellite's surface mass ratio and drag coefficient. , , and These represent the long-term term, first-harmonic sine coefficient, second-harmonic cosine coefficient, and second-harmonic sine coefficient corresponding to the non-spherical gravitational perturbation of Earth J2, respectively. This is a correction term for the Earth's rotation, expressed as the square of the velocity. This is a fundamental correction term for Earth's rotation. parameter , , , , , and The calculation expression is: （2） （3） （4） In the formula, For the satellite's windward area, For satellite quality, The drag coefficient, The value is the J2 perturbation value. For the Earth's radius, The angular velocity of Earth's rotation. The relative velocity of the satellite with respect to the rotating Earth; （5）。 3. The method according to claim 2, characterized in that, Also includes: Let the number of iterations Set the initial value of the reference radius for the reference circular orbit during iteration. Initial values ​​for reference orbit inclination angle iteration Set reference speed tilt angle It is always zero; Among them, the initial time is recorded. The initial value of the distance from the Earth's center is The initial value of the velocity modulus is The initial value of the velocity tilt angle is The initial value of the orbital inclination angle is The initial value of the right ascension of the ascending node is The initial value of the latitude argument is ; The initial value of the reference radius for the reference circular orbit is calculated using the following formula. Initial values ​​for the inclination angle of the reference circular orbit. : （6）。 4. The method according to claim 3, characterized in that, The dimensionless processing of the current state quantity and time includes: Reference radius based on the reference circle orbit in the current iteration round and reference orbital inclination Calculate the reference velocity modulus : （7） In the formula, This is a reference value for the long-term term corresponding to the non-spherical gravitational perturbation of Earth J2; The distance from the Earth's center, the velocity magnitude, the velocity tilt angle, and time are dimensionless and expressed as follows: （8） in, , and Let A, B, and C be the dimensionless distance from the geocenter, velocity modulus, and velocity tilt angle, respectively, with initial values ​​of [missing values]. , and ; The time is dimensionless, and the initial value is... .

5. The method according to claim 4, characterized in that, The obtained linearized dynamic equations include: Linearize equation (1) around the reference circular orbit: （9） （10） In the formula, This is the dimensionless atmospheric density term. This is the dimensionless atmospheric density gradient term. This is a dimensionless atmospheric drag correction term. The J2 perturbation correction term is dimensionless. This is a dimensionless correction term for Earth's rotation. The coefficients of the second harmonic cosine term, which are dimensionless derivatives of the velocity modulus, are... The coefficients of the second harmonic sine term, which is a dimensionless derivative of the velocity modulus, are... For the dimensionless velocity tilt angle derivative, the second harmonic cosine term coefficient; （11） In the formula, For reference density; Reference density For reference radius The partial derivative; parameter , , The definition is as follows: （12） In the formula, and These are the Earth rotation correction terms for the square of the velocity. The corresponding long-term and second-harmonic cosine coefficients, and These are the fundamental correction terms for Earth's rotation. The corresponding long-term and second-harmonic cosine coefficients; （13） In the formula, , and These are the coefficients of the long-term term, the second harmonic cosine term, and the second harmonic sine term of the dimensionless orbital inclination derivative, respectively. , and These are the long-term term, second-harmonic cosine term coefficient, and second-harmonic sine term coefficient of the dimensionless right ascension derivative of the ascending node, respectively. , and These are the long-term term, second-harmonic cosine term coefficient, and second-harmonic sine term coefficient of the dimensionless latitudinal argument derivative, respectively. （14） In the formula, For the fundamental correction of Earth's rotation The reference value satisfies: （15）。 6. The method according to claim 5, characterized in that, The solution of the linearized dynamic equations includes: The solutions to differential equations (9) and (10) are in the following forms: （16） （17） In the formula, coordinates The constant term in the analytical solution, coordinates The coefficient of the long-term variation term in the analytical solution. coordinates The coefficients of the first harmonic cosine term in the analytical solution, coordinates The coefficients of the first harmonic sine term in the analytical solution, coordinates The coefficients of the second harmonic cosine term in the analytical solution, coordinates The coefficients of the second harmonic sine term in the analytical solution; The orbital altitude attenuation coefficient is... This is the orbital ellipticity attenuation coefficient; It is the average latitude argument; in: parameter and satisfy: （18） For coordinates Pick , and In this case, the corresponding coefficient and The solution expression is as follows: The intermediate variable is obtained by solving equation (19). and Then, using equation (20) to obtain and Values; （19） （20） Then, the dimensionless average geocentric distance is obtained using equation (21). initial value Dimensionless average velocity modulus initial value and dimensionless average velocity tilt angle initial value ; （21） Further, intermediate variables are obtained by solving equation (22). , , and Then, using equation (23) to obtain , , and ; （22） （23） For coordinates Pick , and In this case, the coefficient , , and The calculation expression is as follows: （24） In the formula, , and These are the average orbital inclination angles. Mean ascending node right ascension and average latitudinal arc Initial values: （25）。 7. The method according to claim 6, characterized in that, Also includes: Determine whether the analytical solution for the near-circular orbit of a very low Earth orbit satellite satisfies the following equation: （27） In the formula, The set threshold; If equation (27) is not satisfied, then the reference radius and inclination angle of the reference circular orbit are updated using the following equation: （28） In the formula, Let be the reference radius of the reference circular orbit in the (j+1)th iteration. The reference orbit inclination angle is the reference circular orbit for the (j+1)th iteration. Then, let the number of iterations... Then proceed to the next round of iterations to solve the problem; If criterion (27) is satisfied, then the analytical solution of the trajectory is output according to the results of equations (16) and (17) of the current iteration round. When outputting, the state variables are included. , , The dimensionless process is performed as follows: （29） Then, orbital analysis prediction is performed based on the results of equation (29).

Images