Space target orbit calculation method based on SGP4 model

Abstract

The invention discloses a space target orbit calculation method based on an SGP4 model, and relates to the technical field of space target orbit calculation. Dynamic weighting of a high-order perturbation item is achieved through coupling of a continuous window function and height / eccentricity ratio double factors, and the moment calculation precision of a perigee window is remarkably improved; based on a dynamic threshold mechanism of track statistical characteristics, the high-efficiency characteristic of an analytical solution is maintained to the maximum extent on the premise of ensuring the precision; in addition, double timestamp output is adopted, so that the time axis dislocation risk caused by leap second adjustment is fundamentally avoided; the exponential decay design in the weight model guarantees the stable performance of the algorithm at different orbit heights.

Description

Technical Field

[0001] The present invention relates to the technical field of space target orbit calculation, in particular to a space target orbit calculation method based on an SGP4 model. Background Art

[0002] Due to its unique coverage characteristics, highly elliptical orbits have become the preferred configuration for space missions such as regional communication enhancement and polar monitoring. The perigee altitude of such orbits is usually less than 2,000 kilometers, and the apogee can reach tens of thousands of kilometers. The targets in orbit experience drastic changes from the dense atmosphere of low orbit to the weak perturbation environment of deep space.

[0003] With the deployment of new quantum communication satellite constellations, the use of special orbital configurations such as Molniya orbits has increased dramatically, and the requirements for orbit prediction technology have become increasingly stringent. The traditional SGP4 model is based on the analytical solution of the average orbital elements. Although it can quickly generate orbit predictions, its perturbation force modeling method reveals inherent limitations under extreme orbital conditions.

[0004] Existing improvement schemes generally adopt dynamic parameter compensation strategies, such as integrating real-time atmospheric density data correction methods, and dynamically updating atmospheric model parameters by introducing GNSS occultation observation data; in view of the characteristics of highly elliptical orbits, the latest research attempts to divide the perigee sensitive area during the orbit integration process, and adopts a hybrid numerical-analytical method to deal with the critical point perturbation mutation; however, if these methods are used in engineering practice, although the accuracy can be improved by using full numerical solutions, the calculation time will increase sharply, which is difficult to meet the real-time requirements of satellite-borne equipment; if the analytical method framework is maintained, the coupling effect of high-order perturbation terms (such as J3 and J4 terms) near the perigee is insufficiently modeled, resulting in systematic deviations in the perigee altitude prediction. Summary of the Invention

[0005] In view of the above existing problems, the present invention is proposed.

[0006] The present invention provides a space target orbit calculation method based on the SGP4 model to solve the problems that the existing SGP4 improvement scheme faces the dilemma of choosing between computational efficiency and modeling accuracy when dealing with sudden perturbations of highly elliptical orbits, and lacks an effective time system deviation correction mechanism.

[0007] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0008] The embodiment of the present invention provides a method for calculating a space target orbit based on the SGP4 model, which includes:

[0009] Step S1, receiving two lines of orbital element data TLE and UTC time input of the space target;

[0010] Step S2, in the orbit extrapolation cycle of the SGP4 model, the perturbation term calculation strategy is dynamically selected based on the orbit phase characteristic parameter;

[0011] Step S3, when the orbit angular velocity change rate is detected to exceed a dynamic threshold, a variable step numerical integration mode is switched, and the dynamic threshold is automatically adapted based on the target orbit type to generate;

[0012] Step S4, in the orbit extrapolation process, the time reference deviation is corrected through a double timestamp mechanism;

[0013] Step S5, outputting the orbit prediction result containing the position vector, the velocity vector and the corrected timestamp.

[0014] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, the orbit phase characteristic parameter comprises:

[0015] The normalized ratio of the instantaneous orbit height to the standard earth radius;

[0016] The difference parameter of the orbit eccentricity and the critical threshold;

[0017] The derivative parameter of the orbit angular velocity change rate.

[0018] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, the perturbation term calculation strategy in step S2 comprises:

[0019] The weight proportion of the high-order earth non-spherical perturbation term is increased at the perigee window;

[0020] The calculation priority of the solar and lunar gravity perturbation term is enhanced and the weight of the atmospheric resistance term is reduced at the apogee window.

[0021] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, the weight adjustment of the high-order earth non-spherical perturbation term comprises:

[0022] The harmonic compensation term is introduced in the longitude keeping calculation of the geosynchronous orbit target;

[0023] The harmonic compensation term realizes the J22 term effect compensation through the mapping relationship between the perturbation term and the orbit element.

[0024] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, in step S2, the weight adjustment mode of the high-order earth non-spherical perturbation term is:

[0025] Firstly, the perigee window function is defined, and the perigee window function is expressed as:

[0026]

[0027] where σ p represents the perigee determination output, k is the window steepness coefficient, Δr is the difference between the instantaneous radius and the perigee radius, δ r is the window radius threshold value;

[0028] The height normalized and eccentricity difference value is defined as:

[0029] η = e - e c ,

[0030] where r is the instantaneous orbit radius, R e is the average radius of the Earth, is the height normalized ratio, ξ is the height deviation, e is the orbit eccentricity, e c is the critical eccentricity threshold value, and η is the eccentricity difference value;

[0031] Under the coupling of the window function and the two types of difference values, the harmonic order i takes 3, 4, 5, or 22, and the dynamic weight of the harmonic order is represented as:

[0032] ω i = α i (1 + β i ξ + γ i η)σ p ,

[0033] where ω i is the order i weight, α i is the order inherent reference weight, β i is the height sensitivity coefficient, and γ i is the eccentricity sensitivity coefficient;

[0034] In the formula, the reference weight and the sensitivity coefficient are:

[0035] α i = α0i p ,

[0036] where α0 is a reference constant, p is an order attenuation index, β0 and γ0 are initial values of the sensitivity coefficients, λ β , λ γ are exponential attenuation rates, and ∈ is a natural constant;

[0037] When the target is in the geosynchronous resonance region, a compensation term is constructed as:

[0038]

[0039] where ρ is the resonance compensation amplitude, a and b are mapping coefficients, is the difference between the geocentric meridian angle and the satellite longitude;

[0040] The final adjusted high-order perturbation coefficient is written as:

[0041]

[0042] wherein, is the adjusted harmonic coefficient, J i is the original high-order aspherical term, δ i,22 is the Kronecker delta, 1 if i = 22, otherwise 0, 22 belongs to the compensation order, and 3-5 are processed respectively.

[0043] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, the dynamic threshold determination method of step S3 comprises:

[0044] The threshold interval is set according to the angular velocity change rate statistical characteristics in the orbit history data;

[0045] When the target is in the geosynchronous orbit, a threshold correction factor based on the orbit inclination is adopted.

[0046] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, the variable step numerical integration mode comprises:

[0047] The step length control parameter is generated according to the weighted combination of the orbit semi-major axis and the inclination change rate;

[0048] The high-order numerical integration method is used for intensive calculation at the perigee window, and the analytical solution method is restored at the apogee window.

[0049] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, in step S3, in the process of adaptively generating a dynamic threshold, whether the angular velocity change rate triggers the variable step integration in the orbit extrapolation cycle,

[0050] First, the historical samples ω j ′ are statistically estimated, ω j ′ is the first-order derivative of the angular velocity, and the statistical formula is:

[0051]

[0052] wherein, μ ω′ represents the sample mean of the angular velocity derivative, σ ω′ represents the sample standard deviation, L is the sample length, and ω′ j is the angular velocity derivative at the jth observation time;

[0053] The reference threshold is given based on the statistical results:

[0054] θ0= μ ω′ + λ σ σ ω′ ,

[0055] Wherein, θ0is the reference threshold, λ σ is the standard deviation amplification factor;

[0056] According to the target orbit type, adaptive correction is added, and an orbit category index is first constructed:

[0057]

[0058] Wherein, χis the category index, a s is the orbit semi-major axis, R e is the average radius of the earth;

[0059] The inclination adjustment factor is defined:

[0060] υ=1+ κ ψ |sinψ|,

[0061] Wherein, υis the inclination adjustment factor, κ ψ is the inclination sensitivity coefficient, and ψ is the orbit inclination;

[0062] The final dynamic threshold is obtained by introducing the category index and the inclination adjustment:

[0063] Θ dyn = θ0(1+ κ χ χ)υ,

[0064] Wherein, Θ dyn is the dynamic threshold, and κ χ is the category index sensitivity coefficient;

[0065] When ω ′ > Θ dyn , switch to variable step integration, otherwise maintain the analytical step size.

[0066] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, the double time stamp mechanism of step S4 comprises:

[0067] Embedding the cumulative second offset mark of international atomic time TAI in the UTC time output;

[0068] Eliminate the time system deviation of UT1 and UTC by pre-storing a leap second compensation table.

[0069] As a preferred scheme of the space target orbit calculation method based on the SGP4 model, it further comprises:

[0070] Constructing an atmospheric parameter mutation detection model based on real-time space weather data;

[0071] When the abnormal fluctuation of the solar radiation flux index is detected, the atmospheric resistance parameter in the TLE is dynamically scaled.

[0072] The present application has the beneficial effects that: the present application realizes dynamic weighting of high-order perturbation terms by coupling a continuous window function with height / eccentricity double factors, significantly improves the moment calculation accuracy of the perigee window; and based on a dynamic threshold mechanism of the orbital statistical characteristics, the high efficiency characteristics of the analytical solution method are maximally maintained under the premise of ensuring accuracy; in addition, the present application uses double timestamp output to fundamentally avoid the time axis misplacement risk caused by leap second adjustment; the exponential decay design in the weight model guarantees the stable performance of the algorithm at different orbital altitudes. BRIEF DESCRIPTION OF DRAWINGS

[0073] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed to be used in the embodiment description will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0074] Figure 1 A flowchart of a space target orbit calculation method based on the SGP4 model in embodiment 1. DETAILED DESCRIPTION

[0075] In order to make the above-mentioned purposes, features and advantages of the present application more apparent and easy to understand, the specific embodiments of the present application will be described in detail below with reference to the drawings of the specification.

[0076] In the following description, many specific details are set forth in order to provide a thorough understanding of the present application, but the present application can also be implemented in other ways different from those described herein, and those skilled in the art can make similar generalizations without departing from the connotation of the present application, therefore the present application is not limited to the specific embodiments disclosed below.

[0077] Secondly, the "one embodiment" or "embodiment" referred to herein means that the specific features, structures or characteristics can be included in at least one implementation of the present application. "In one embodiment" appearing in different places in the specification does not mean the same embodiment, nor is it an independent or alternative embodiment that excludes other embodiments.

[0078] The key challenge in the field of current space target orbit calculation is to balance the contradiction between calculation efficiency and modeling accuracy, especially for the perturbation mutation problem of special orbit configurations such as high elliptical orbit and geosynchronous orbit. The traditional SGP4 model is based on the analytical solution of the average orbit element, which can quickly generate orbit prediction, but has significant limitations in dealing with high-order perturbation coupling of perigee window, sudden space weather disturbance and time reference drift. The existing improvement scheme mostly uses static parameter correction or full numerical integration strategy, the former is difficult to dynamically adapt to the dramatic changes of the orbit environment, and the latter is too large in the consumption of computing resources and cannot meet the real-time demand of the on-board equipment. In this background, an innovative method is needed to inherit the high efficiency of the SGP4 model and break through the bottleneck of perturbation modeling. In this background, the present application is proposed.

[0079] Embodiment 1, refer to Figure 1 The embodiment provides a space target orbit calculation method based on SGP4 model, comprising the following steps:

[0080] Step S1, receiving two-line orbit element TLE data of a space target and UTC time input;

[0081] Step S2, in the orbit extrapolation cycle of SGP4 model, dynamically selecting perturbation term calculation strategy based on orbit phase characteristic parameters;

[0082] The orbit phase characteristic parameters include:

[0083] The normalized ratio of instantaneous orbit height to standard earth radius;

[0084] The difference parameter of orbit eccentricity and critical threshold;

[0085] The derivative parameter of orbit angular velocity change rate;

[0086] The perturbation term calculation strategy in step S2 includes:

[0087] Increasing the weight proportion of high-order earth non-spherical perturbation terms (J3-J5 terms) in the perigee window;

[0088] Enhancing the calculation priority of the sun and moon gravity perturbation terms in the apogee window and reducing the weight of the atmospheric drag term;

[0089] The weight adjustment of high-order earth non-spherical perturbation terms includes:

[0090] Introducing a resonance compensation term in the longitude keeping calculation of geosynchronous orbit target;

[0091] The resonance compensation term realizes J22 term effect compensation through the mapping relationship between perturbation terms and orbit elements;

[0092] In step S2, the weight adjustment mode of high-order earth non-spherical perturbation term is:

[0093] Firstly, the perigee window function is defined, which is expressed as:

[0094]

[0095] Wherein, σ p represents the perigee determination output, k is the window steepness coefficient, Δr is the difference between the instantaneous radius and the perigee radius, δ r is the window radius threshold;

[0096] The height normalization and eccentricity difference are defined, which are expressed as:

[0097] η=e-e c ,

[0098] Wherein, r is the instantaneous orbit radius, R e is the average radius of the earth, is the height normalization ratio, ξ is the height deviation, e is the orbit eccentricity, e c is the critical eccentricity threshold, and η is the eccentricity difference;

[0099] Under the coupling of the window function and the two types of differences, the harmonic order i is 3, 4, 5 or 22, and the dynamic weight of the harmonic order is expressed as:

[0100] ω i =α i (1+β i ξ+γ i η)σ p ,

[0101] Wherein, ω i is the order i weight, α i is the order inherent reference weight, β i is the height sensitivity coefficient, and γ i is the eccentricity sensitivity coefficient;

[0102] In the formula, the reference weight and the sensitivity coefficient are:

[0103] α i =α0i p ,

[0104] Wherein, α0 is the reference constant, p is the order attenuation index, β0 and γ0 are the sensitivity coefficient initial values, λ β , λ γ is the exponential attenuation rate, and ∈ is the natural constant;

[0105] When the target is in the geosynchronous resonance region, the compensation term is constructed:

[0106]

[0107] wherein, ρ is the resonance compensation amplitude, a, b are mapping coefficients, is the difference between the prime vertical and the satellite longitude;

[0108] The final adjusted high-order perturbation coefficient is written as:

[0109]

[0110] wherein, is the adjusted harmonic coefficient, J i is the original high-order spheroidal term, δ i,22 is the Kronecker delta, 1 if i = 22, otherwise 0, 22 belongs to the compensation order, and 3-5 are processed respectively;

[0111] Specifically, the weight model is σ p A continuous window is generated to make the weight rise exponentially with the perigee distance and then decay smoothly with the apogee distance, avoiding the shock caused by the hard threshold. The high and eccentricity factors are adjusted by β i and γ i respectively, so that high and low orbits can obtain adaptive scaling in the same architecture. The reference weight α i decays according to the power law, and the exponential sensitivity coefficient is used to make the weight monotonically decrease with the order and remain differentiable, improving the analytical stability. The resonance compensation ρ directs the dynamic energy of J4 and J5 to J 22 , providing additional torque for longitude maintenance in the geosynchronous zone and reducing long-term drift. The overall method introduces a smooth and adjustable high-order perturbation weight curve on the premise of maintaining analytical efficiency, suppresses high-frequency errors that occur during integral analytical switching, and improves long-term prediction consistency and short-term window reliability.

[0112] Step S3, when the orbit angular velocity change rate is detected to exceed the dynamic threshold, switch to the variable step size numerical integration mode, and the dynamic threshold is automatically adapted based on the target orbit type to generate;

[0113] The dynamic threshold determination method of step S3 includes:

[0114] Set the threshold interval according to the angular velocity change rate statistical characteristics in the orbit historical data;

[0115] When the target is in geosynchronous orbit, use the threshold correction factor based on the orbit inclination;

[0116] The variable step size numerical integration mode includes:

[0117] Generate a step size control parameter according to the weighted combination of the orbit semi-major axis and the inclination rate of change;

[0118] High order numerical integration method is used for intensive calculation in perigee window, and analytical solution is used in apogee window;

[0119] In step S3, in the process of adaptively generating dynamic threshold, whether the angular velocity change rate triggers variable step integration in the orbit extrapolation cycle is determined in real time,

[0120] First, the historical sample ω j ′ is statistically estimated, ω j ′ is the first derivative of angular velocity, and the statistical formula is:

[0121]

[0122] wherein μ ω′ represents the sample mean of angular velocity derivative, σ ω′ represents the sample standard deviation, L is the sample length, and ω′ j is the angular velocity derivative at the jth observation time;

[0123] The reference threshold is given based on the statistical results:

[0124] θ0=μ ω′ +λ σ σ ω′ ,

[0125] wherein θ0 is the reference threshold, and λ σ is the standard deviation amplification coefficient;

[0126] According to the target orbit type, adaptive correction is added, and first, the orbit category index is constructed:

[0127]

[0128] wherein χ is the category index, a s is the semi-major axis of the orbit, and R e is the average radius of the earth;

[0129] The inclination adjustment factor is defined:

[0130] υ=1+κ ψ |sinψ|,

[0131] wherein υ is the inclination adjustment factor, κ ψ is the inclination sensitivity coefficient, and ψ is the orbit inclination;

[0132] The final dynamic threshold is obtained by introducing the category index and the inclination adjustment:

[0133] Θ dyn =θ0(1+κ χ χ)υ,

[0134] wherein Θ dyn is a dynamic threshold, κ χ is a category index sensitivity coefficient;

[0135] When ω' > Θ dyn , switch to variable step-size integration, otherwise maintain the analytical step-size;

[0136] Specifically, the threshold model captures the statistical fluctuations of the historical angular velocity derivative with θ0, and continuously maps the orbital distance information with the logarithmic χ, so that near-Earth and far-Earth targets are amplified or contracted differently in the same framework. The inclination factor υ introduces the |sinψ| form, so that polar and inclined orbits obtain more stringent thresholds under the same semi-major axis. The two types of sensitivity coefficients κ χ , κ ψ provide adjustable margins, making it easy to adjust offline according to new samples. The overall threshold increases logarithmically with the semi-major axis and linearly with the absolute value of the inclination. The curve is continuous and differentiable, and can be directly embedded into the error control module of the numerical integrator, avoiding the frequent mode jitter caused by hard thresholds and improving the time step stability and computational efficiency of high eccentricity and high inclination orbits.

[0137] Step S4, correcting the time reference deviation through a double timestamp mechanism during orbit extrapolation;

[0138] The double timestamp mechanism of step S4 includes:

[0139] Embedding the cumulative second offset mark of international atomic time TAI in UTC time output;

[0140] Eliminating the time system deviation between UT1 and UTC by pre-storing a leap second compensation table;

[0141] Step S5, outputting the orbit prediction results containing position vectors, velocity vectors and corrected timestamps;

[0142] The space target orbit calculation method further includes:

[0143] Based on real-time space weather data, an atmospheric parameter mutation detection model is constructed;

[0144] When abnormal fluctuations of the solar radiation flux index are detected, the atmospheric drag parameters in the TLE are dynamically scaled.

[0145] In summary, the present scheme proposes a hybrid calculation method that combines dynamic perturbation weight adjustment and adaptive integration control to reconstruct the core architecture of the SGP4 model from three dimensions:

[0146] 1. Dynamic perturbation modeling system

[0147] Continuous Window Function Mechanism: By introducing a perigee window determination model based on a Sigmoid function, the orbital phase is divided into a perigee-sensitive region and an apogee-stable region, achieving a smooth transition in perturbation force calculation. The exponential response of the window function accurately captures orbital parameter mutation points, avoiding the numerical oscillations caused by traditional hard threshold segmentation.

[0148] A multi-factor weighted coupling model: This model constructs a joint control mechanism combining height normalization ratios, eccentricity differences, and perturbation orders, dynamically adjusting the weights of the J3-J5 terms with the orbital state. This model adaptively suppresses high-order perturbations through exponential decay parameters, ensuring computational accuracy in the perigee window while avoiding unnecessary computational overhead in the apogee region.

[0149] Geosynchronous resonance compensation: To address the longitude drift problem of GEO targets, a compensation term based on the directional superposition of J4 / J5 energy is designed. Through the dynamic mapping of the geocentric meridian angle and orbital inclination, high-order perturbation effects are converted into longitude-maintaining control torque, significantly reducing fuel consumption.

[0150] 2. Adaptive integral control strategy

[0151] Statistics-driven dynamic threshold: Based on the statistical characteristics of the angular velocity change rate of historical orbit data, a dynamic threshold model is constructed that includes logarithmic scaling of the orbit semi-major axis and an inclination adjustment factor. This model automatically adapts the trigger conditions based on the target orbit type, ensuring that variable-step integration is initiated only during the necessary window.

[0152] Hybrid integration mode switching: The fourth-order Runge-Kutta method is used for short, intensive computations in the perigee window, while the SGP4 analytical solution is restored in the apogee window. A step-size control parameter is generated through a weighted combination of the semi-major axis and the rate of change of the inclination angle, achieving an optimal balance between computational accuracy and efficiency.

[0153] 3. Enhanced robustness of the time system

[0154] Dual timestamp output mechanism: Synchronously outputs UTC standard time and International Atomic Time (TAI) mark, combined with a pre-stored leap second compensation table, to fundamentally eliminate time axis drift errors caused by changes in the Earth's rotation rate.

[0155] Real-time space weather response: A sliding window detection model for the solar radiation flux index F10.7 is constructed. When sudden ionospheric disturbances are detected, the atmospheric drag parameters in the TLE are dynamically scaled to improve the model's adaptability to sudden changes in the space environment.

[0156] This invention not only solves the inherent defects of the SGP4 model in highly dynamic track scenarios, but also creates a new paradigm for integrating classic analytical models with intelligent adaptive algorithms. Its technical framework can be extended to the following fields:

[0157] Deep space navigation: improve the orbit determination accuracy of Mars probe and asteroid probe by extending the weight model of the sun and moon perturbation.

[0158] Space situation awareness: build a high-precision space target catalog database combined with inter-satellite networking observation data.

[0159] Commercial space application: provide a low-cost and high-reliability autonomous orbit management solution for low-orbit Internet constellation.

[0160] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present application and are not limiting. Although the present application has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solutions of the present application can be modified or replaced by equivalents without departing from the spirit and scope of the present application, which should be covered in the scope of the claims of the present application.

Claims

1. A space target orbit calculation method based on the SGP4 model, characterized in that: include, Step S1, receiving two lines of orbital element data TLE and UTC time input of the space target; Step S2, in the orbital extrapolation cycle of the SGP4 model, dynamically selecting the perturbation term calculation strategy based on the orbital phase characteristic parameters; Step S3: When it is detected that the rate of change of the orbital angular velocity exceeds a dynamic threshold, the system switches to a variable step-size numerical integration mode, where the dynamic threshold is automatically generated based on the target orbit type. Step S4, correcting the time base deviation through the dual timestamp mechanism during the orbit extrapolation process; Step S5: Output the orbit prediction result including the position vector, velocity vector and correction time stamp.

2. A method for calculating a space target orbit based on the SGP4 model according to claim 1, characterized in that: The orbital phase characteristic parameters include: Normalized ratio of instantaneous orbital altitude to standard Earth radius; The difference parameter between orbital eccentricity and critical threshold; The derivative parameter of the rate of change of orbital angular velocity.

3. A method for calculating a space target orbit based on the SGP4 model as claimed in claim 2, characterized in that: The perturbation term calculation strategy in step S2 includes: Increase the weight ratio of high-order Earth non-spherical perturbation terms in the perigee window; Increase the calculation priority of the solar and lunar gravitational perturbation terms and reduce the weight of the atmospheric drag terms in the apogee window.

4. A method for calculating a space target orbit based on the SGP4 model as claimed in claim 3, characterized in that: The weight adjustment of the high-order non-spherical perturbation term of the earth includes: Introducing a resonance compensation term into the longitude hold calculation for geosynchronous orbit targets; The resonance compensation term realizes the compensation of the J22 term effect through the mapping relationship between the perturbation term and the number of orbital elements.

5. The method for calculating the space target orbit based on the SGP4 model according to claim 4, characterized in that: In step S2, the weight adjustment method of the high-order earth non-spherical perturbation term is: First define the perigee window function, which is expressed as: Among them, σ p represents the perigee determination output, k is the window steepness coefficient, Δr is the difference between the instantaneous radius and the perigee radius, δ r is the window radius threshold; Define the height normalization and eccentricity difference as: Where r is the instantaneous orbit radius, R e is the average radius of the Earth, is the normalized ratio of height, ξ is the height deviation, e is the orbital eccentricity, e c is the critical eccentricity threshold, η is the eccentricity difference; Under the coupling of the window function and the two types of differences, the harmonic order i takes 3, 4, 5 or 22, and the dynamic weight of the harmonic order is expressed as: oh i =a i (1+b i ξ+γ i h)s p , Among them, ω i is the weight of order i, α i is the order inherent benchmark weight, β i is the height sensitivity coefficient, γ i is the eccentricity sensitivity coefficient; In the formula, the benchmark weight and sensitivity coefficient are: Among them, α0 is the reference constant, p is the order attenuation index, β0,γ0 are the initial values ​​of the sensitivity coefficient, λ β ,λ γ is the exponential decay rate, ∈ is a natural constant; When the target is in the geosynchronous resonance region, the compensation term is constructed: Among them, ρ is the resonance compensation amplitude, a, b are mapping coefficients, is the difference between the geocentric meridian and the satellite longitude; The final adjusted high-order perturbation coefficients are written as: in, is the adjusted harmonic coefficient, J i is the original high-order non-spherical term, δ i,22 is Kroneckerδ, which takes the value 1 if i=22 and 0 otherwise. 22 belongs to the compensation order and is processed separately from 3-5.

6. A method for calculating a space target orbit based on the SGP4 model according to claim 1, characterized in that: The dynamic threshold determination method of step S3 includes: Setting a threshold interval based on the statistical characteristics of the angular velocity change rate in the historical orbit data; When the target is in geosynchronous orbit, a threshold correction factor based on orbital inclination is used.

7. A method for calculating a space target orbit based on the SGP4 model according to claim 6, characterized in that: The variable step size numerical integration mode includes: The step size control parameter is generated based on a weighted combination of the orbit semi-major axis and the rate of change of the inclination angle; High-order numerical integration methods are used for intensive calculations in the perigee window, and the analytical solution is restored in the apogee window.

8. A method for calculating a space target orbit based on the SGP4 model according to claim 7, characterized in that: In step S3, during the process of adaptively generating the dynamic threshold, in order to determine in real time whether the angular velocity change rate triggers the variable step-size integration in the orbit extrapolation cycle, First, the historical sample ω j ′ Perform statistical estimation, ω j ′ is the first-order derivative of angular velocity, statistical formula: Among them, μ ω′ represents the sample mean of the angular velocity derivative, σ ω′ represents the sample standard deviation, L is the sample length, ω′ j is the angular velocity derivative at the jth observation moment; Give a benchmark threshold based on statistical results: θ0=μ ω′ +λ σ s ω′ , Among them, θ0 is the reference threshold, λ σ is the standard deviation magnification factor; Adaptive correction is added according to the target orbit type. The orbit category index is first constructed: Among them, χ is the category index, a s is the semi-major axis of the orbit, R e is the average radius of the Earth; Define the tilt adjustment factor: u=1+k ψ |sinψ|, Among them, υ is the tilt adjustment factor, κ ψ is the inclination sensitivity coefficient, ψ is the orbit inclination; The final dynamic threshold is obtained by introducing the category index and tilt adjustment: I dyn =θ0(1+κ χ x)y, Among them, Θ dyn is the dynamic threshold, κ χ is the category index sensitivity coefficient; When ′ >Θ dyn When , switch to variable step size integration, otherwise maintain the analytical step size.

9. The method for calculating a space target orbit based on the SGP4 model according to claim 1, wherein: The dual timestamp mechanism of step S4 includes: Embed the cumulative seconds offset mark of the International Atomic Time TAI in the UTC time output; The time system deviation between UT1 and UTC is eliminated by pre-storing leap second compensation tables.

10. The method for calculating the space target orbit based on the SGP4 model according to claim 1, wherein: Also includes: Build an atmospheric parameter mutation detection model based on real-time space weather data; When abnormal fluctuations in the solar radiation flux index are detected, the atmospheric drag parameters in the TLE are dynamically scaled.

Images