Low-performance chip-based satellite orbit extrapolation rapid calculation method and system

Abstract

The application relates to the technical field of satellite orbit extrapolation rapid calculation, and discloses a satellite orbit extrapolation rapid calculation method and system based on a low-performance chip, wherein the satellite orbit extrapolation rapid calculation method based on the low-performance chip comprises the following steps: constructing a simplified satellite orbit calculation model based on main perturbation forces; adopting an optimized fourth-order Runge-Kutta integrator to perform state vector integration; calculating main perturbation accelerations; converting the main perturbation accelerations into the J2000 coordinate system through coordinate system conversion, synthesizing total perturbation accelerations in the J2000 coordinate system, constructing a six-dimensional state derivative function; cyclically integrating the six-dimensional state derivative function to update a state vector until a target extrapolation time is reached; and outputting a state vector of the satellite at the target time as a prediction result; the optimized fourth-order Runge-Kutta integrator is adopted, the calculation complexity is reduced, and the orbit extrapolation calculation can be efficiently performed on the low-performance chip.

Description

Technical Field

[0001] This invention relates to the field of rapid calculation technology for satellite orbit extrapolation, and more specifically, to a method and system for rapid calculation of satellite orbit extrapolation based on low-performance chips. Background Technology

[0002] With the development of aerospace technology, miniaturized and low-cost satellite terminal products are increasingly widely used in communication, navigation and remote sensing. These terminal devices usually need to perform satellite orbit prediction calculations to determine key information such as satellite visibility time windows and communication link establishment conditions. However, traditional satellite orbit extrapolation calculation methods usually rely on high-performance computing equipment, using complex numerical integration algorithms and comprehensive perturbation models, such as variable step size Runge-Kutta-Fehlberg7(8) integrators, which have a large computational load and high resource consumption. At present, most miniaturized satellite terminal products on the market use low-performance chips (such as ARM Cortex-M4 series, main frequency 72MHz, no floating-point arithmetic unit or weak floating-point arithmetic capability) as core processors to reduce costs, size and power consumption. Such chips face serious performance bottlenecks when processing complex satellite orbit extrapolation calculations, specifically manifested as follows:

[0003] The traditional RKF7(8) integrator runs slowly on low-performance chips without floating-point units, making it difficult to meet real-time computing requirements;

[0004] A complete model of Earth's non-spherical gravitational field (typically requiring 70th order or higher) has high computational complexity, exceeding the processing capabilities of low-performance chips;

[0005] Precise calculations of perturbations such as the gravitational pull of the sun and moon and atmospheric drag require a large number of floating-point operations, which are difficult to achieve efficiently in resource-constrained environments.

[0006] These issues have led existing miniaturized satellite terminal products to either rely on external computing resources for orbit prediction, increasing system complexity and communication overhead, or employ extremely simplified orbit models, sacrificing prediction accuracy and affecting the practical application of the terminal products. Therefore, how to achieve efficient and accurate satellite orbit extrapolation calculations in a low-performance chip environment has become a key technical problem restricting the development of miniaturized satellite terminal products. Summary of the Invention

[0007] This invention provides a method and system for rapid calculation of satellite orbit extrapolation based on low-performance chips, solving the technical problem of how to achieve rapid calculation of satellite orbit extrapolation on low-performance hardware chips in related technologies.

[0008] This invention provides a fast method for satellite orbit extrapolation based on low-performance chips, including:

[0009] A simplified satellite orbit calculation model based on the main perturbation force is constructed and the initial state vector of the satellite is set;

[0010] An optimized fourth-order Runge-Kutta integrator is used for state vector integration, and the initial state of the optimized fourth-order Runge-Kutta integrator is determined based on the satellite's initial state vector.

[0011] When the optimized fourth-order Runge-Kutta integrator performs integration calculations, it calculates the main perturbation accelerations based on both the ecliptic coordinate system and the Earth-fixed coordinate system.

[0012] The main perturbation accelerations are uniformly transformed to the J2000 coordinate system through coordinate system transformation. The total perturbation acceleration is synthesized in the J2000 coordinate system. A six-dimensional state derivative function is constructed based on the total perturbation acceleration and the central gravitational acceleration.

[0013] The optimized fourth-order Runge-Kutta integrator is used to iteratively integrate the six-dimensional state derivative function to update the state vector until the target extrapolation time is reached.

[0014] The satellite's state vector at the target time is output as the prediction result.

[0015] Furthermore, the main perturbation forces include the perturbation force of the Earth's non-spherical gravitational field, the perturbation force of the Sun and Moon's gravity, and the perturbation force of atmospheric drag.

[0016] Furthermore, the optimized fourth-order Runge-Kutta integrator is expressed as:

[0017]

[0018] in This represents the state vector at step n3; h3 represents the state vector at step n3+1; h3 represents the integration step size; k1 represents the first slope estimate; k2 represents the second slope estimate; k3 represents the third slope estimate; k4 represents the fourth slope estimate; f represents the right-hand side function of the differential equation. This represents the time for the n3rd step.

[0019] Furthermore, the initial state vector of the satellite is represented as:

[0020]

[0021] in Represents the initial state vector. and These are the initial position vector and velocity vector of the satellite in the J2000 coordinate system;

[0022] The initial state of the optimized fourth-order Runge-Kutta integrator, determined based on the satellite's initial state vector, refers to the optimized fourth-order Runge-Kutta integrator calculating the satellite state sequentially from the initial time t0 at times t0+h3, t0+2h3, ..., t0+n3h3 according to the stated integration step size, until the target extrapolated time t is reached. f The state vector of the optimized fourth-order Runge-Kutta integrator at initial time t0 is:

[0023] Furthermore, the calculation of the main perturbation accelerations based on the ecliptic coordinate system and the Earth-fixed coordinate system respectively specifically includes determining the gravitational position vector of the Sun and Moon based on the ecliptic coordinate system and then calculating the gravitational perturbation acceleration of the Sun and Moon based on the gravitational position vector of the Sun and Moon, and calculating the perturbation acceleration of the Earth's non-spherical gravitational field and the atmospheric drag perturbation acceleration based on the Earth-fixed coordinate system.

[0024] Furthermore, determining the gravitational position vector of the Sun and Moon based on the ecliptic coordinate system, and then calculating the gravitational perturbation acceleration of the Sun and Moon based on the gravitational position vector, specifically includes:

[0025] Calculate the orbital elements of the Sun in the geocentric ecliptic coordinate system;

[0026] Calculate the orbital parameters of the Moon in the ecliptic coordinate system;

[0027] Calculate the rectangular coordinates of the Sun and Moon relative to the vernal equinox at an epoch;

[0028] Transform the positions of the Sun and Moon from the ecliptic coordinate system to the J2000 coordinate system;

[0029] Calculation of lunar and solar gravitational perturbation accelerations based on the J2000 coordinate system.

[0030] Furthermore, the calculation of the perturbation acceleration of the Earth's non-spherical gravitational field based on the Earth-fixed coordinate system specifically includes:

[0031] Perform coordinate transformation between the J2000 coordinate system and the Earth-fixed coordinate system;

[0032] Calculate the spherical coordinates of the satellite in the Earth-fixed coordinate system;

[0033] Calculation of the acceleration perturbed by the Earth's non-spherical gravitational field based on the expansion of spherical harmonic functions.

[0034] Furthermore, the calculation of atmospheric drag perturbation acceleration based on the Earth-fixed coordinate system specifically includes:

[0035] Calculate the satellite's velocity relative to the Earth's atmosphere;

[0036] Calculate satellite altitude based on Earth's oblateness;

[0037] Atmospheric density is calculated based on the exponential decay model.

[0038] Calculate the atmospheric drag perturbation acceleration.

[0039] Furthermore, the extrapolated time to reach the target refers to t0 + nh3 ≥ t f And when t0+nh3≥t f The state vector at the target time is directly output, represented as:

[0040]

[0041] in Represents the satellite's final state vector; and The target extrapolation time t is respectively f The satellite's position and velocity vectors in the J2000 coordinate system; t f The value represents the target extrapolation time; n represents the current integration step number; h4 represents the integration step size; and t0 represents the initial time.

[0042] This invention provides a fast satellite orbit extrapolation calculation system based on a low-performance chip, used to execute the aforementioned fast satellite orbit extrapolation calculation method based on a low-performance chip, comprising:

[0043] The orbit calculation model simplification module is used to construct a simplified satellite orbit calculation model based on the main perturbation force and set the satellite's initial state vector;

[0044] The integrator integration module is used to perform state vector integration using an optimized fourth-order Runge-Kutta integrator, and to determine the initial state of the optimized fourth-order Runge-Kutta integrator based on the satellite's initial state vector.

[0045] The perturbation force fast calculation module, when the optimized fourth-order Runge-Kutta integrator performs integration calculation, calculates the main perturbation acceleration based on the ecliptic coordinate system and the Earth-fixed coordinate system respectively;

[0046] The coordinate transformation optimization module is used to uniformly transform the main perturbation acceleration to the J2000 coordinate system through coordinate system transformation, synthesize the total perturbation acceleration in the J2000 coordinate system, and construct a six-dimensional state derivative function based on the total perturbation acceleration and the central gravitational acceleration.

[0047] The six-dimensional state derivative function construction module is used to perform cyclic integration on the six-dimensional state derivative function using the optimized fourth-order Runge-Kutta integrator to update the state vector until the target extrapolation time is reached.

[0048] The result output module is used to output the satellite's state vector at the target time as the prediction result.

[0049] The beneficial effects of this invention are as follows:

[0050] The invention improves computational speed by using an optimized fourth-order Runge-Kutta integrator: The invention replaces the traditional RKF7(8) integrator with a specially optimized fourth-order Runge-Kutta integrator, which reduces computational complexity and enables orbit extrapolation calculations to be executed efficiently on low-performance chips, meeting the needs of real-time applications.

[0051] By simplifying the Earth's non-spherical gravitational field model to maintain a balance between accuracy and efficiency, the computational load is reduced while maintaining the accuracy required for engineering applications. This allows the orbit extrapolation accuracy of near-Earth, medium-Earth orbit, and high-Earth orbit satellites to meet the practical application needs of communication, navigation, and remote sensing satellite terminals.

[0052] By improving the calculation method of solar and lunar gravitational perturbation force to reduce resource consumption, adopting a simplified method for calculating the positions of the sun and moon, optimizing the calculation process of solar and lunar gravitational perturbation acceleration, reducing complex mathematical operations, and reducing dependence on the floating-point computing power of the chip, the calculation of solar and lunar gravitational perturbation force can be efficiently implemented on low-performance chips.

[0053] By optimizing the calculation method for atmospheric drag perturbation, computational efficiency is improved. A simplified exponential atmospheric density model and an optimized method for calculating satellite relative atmospheric velocity are adopted, which reduces the complexity of atmospheric drag perturbation calculation, improves computational efficiency, and makes atmospheric drag perturbation calculation adaptable to low-performance chip environments. Attached Figure Description

[0054] Figure 1 This is a flowchart of the fast satellite orbit extrapolation calculation method based on low-performance chips in this invention;

[0055] Figure 2 This is a simplified satellite orbit extrapolation calculation model diagram in this invention. Detailed Implementation

[0056] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, some features described in the examples may be combined in other examples.

[0057] At least one embodiment of the present invention discloses a method for fast calculation of satellite orbit extrapolation based on low-performance chips, such as... Figures 1 to 2 As shown, it includes:

[0058] Step 1: Construct a simplified satellite orbit calculation model based on the main perturbation force and set the initial state vector of the satellite;

[0059] Step 1.1: Determine the calculation model based on engineering accuracy requirements;

[0060] In this embodiment, due to the resource constraints of low-performance chips, a simplified satellite orbit calculation model is adopted, considering only the perturbations that have the most significant impact on satellite orbits: perturbations of the Earth's non-spherical gravitational field (mainly the J2 term, i.e., the Earth's oblateness term, which represents the main component of the Earth's shape deviating from a sphere), perturbations of the Sun and Moon's gravity, and perturbations of atmospheric drag.

[0061] The basic dynamic equations of satellite motion in the simplified satellite orbit calculation model are expressed as follows:

[0062]

[0063] in This represents the satellite's position vector in the J2000 coordinate system; This represents the satellite's acceleration vector in the J2000 coordinate system; μ represents the Earth's gravitational constant, with a value of 3.986004418 × 10⁻⁶. 14 m 3 / s 2 ; r represents the geocentric distance of the satellite, i.e. This represents the sum of accelerations generated by the main perturbation forces in the J2000 coordinate system. All perturbation forces refer to the perturbation of the Earth's non-spherical gravitational field, the gravitational perturbation of the Sun and Moon, and the perturbation of atmospheric drag.

[0064] In this embodiment, the simplified model is implemented in the following ways: the Earth's non-spherical gravitational field model is only retained up to order 20; the positions of the Sun and Moon are calculated using a simplified method; and the atmospheric density model adopts an exponential decay simplified model.

[0065] Step 1.2, Preparation of the satellite's initial state vector;

[0066] Define the initial state vector of the satellite:

[0067]

[0068] in Represents the initial state vector. and These are the initial position vector and velocity vector of the satellite in the J2000 coordinate system, respectively.

[0069] Step 2: Use an optimized fourth-order Runge-Kutta integrator to integrate the state vector, and determine the initial state of the optimized fourth-order Runge-Kutta integrator based on the satellite's initial state vector.

[0070] Conventional high-precision orbit extrapolation calculations typically employ the variable-step Runge-Kutta-Fehlberg 7(8) integrator, hereinafter referred to as RKF7(8). However, the RKF7(8) integrator operates relatively slowly on low-performance hardware chips, especially for chips lacking hardware floating-point units. The calculation speed of the RKF7(8) integrator with low-performance hardware chips is difficult to meet the requirements of engineering applications. The calculation form of the RKF7(8) integrator is as follows:

[0071]

[0072] in Represents the state variable at step n1; This represents the state variable at step n1+1; represents the higher-order state variable at step n1+1; h1 represents the integration step size; ∑ represents the summation operator; i2 represents the summation index; This represents the weight coefficient of the i2th function; These represent the higher-order weight coefficients of the i-th 3rd degree function; This indicates the evaluation of the function at the i2th degree; This indicates the evaluation of the function at the i-th degree; This represents the 8th-order error term for step size h1; This represents the 9th-order error term for step size h1.

[0073] The calculation requires multiple function evaluations, resulting in a large computational load.

[0074]

[0075] Where f0 represents the initial function evaluation; f represents the right-hand side function of the differential equation; This represents the time for step n2; Represents the state variable at step n2; This indicates the evaluation of the function at the i2th degree; Indicates the coefficient at the time node; m represents the coefficients of the weight matrix; m4 represents the order of the integrator.

[0076] However, this integrator requires a large amount of computation and is not suitable for low-performance chips.

[0077] Through experimental comparison and analysis, this embodiment selects to use an optimized fourth-order Runge-Kutta (RK4) integrator, which reduces computational complexity while ensuring computational accuracy. Specifically, this is expressed as follows:

[0078]

[0079] in This represents the state vector at step n3; h3 represents the state vector at step n3+1; h3 represents the integration step size; k1 represents the first slope estimate; k2 represents the second slope estimate; k3 represents the third slope estimate; k4 represents the fourth slope estimate; f represents the right-hand side function of the differential equation. This represents the time for the n3rd step.

[0080] and The six-dimensional state vectors for steps n3 and n3+1 are respectively. The first three components represent the satellite's position coordinates in the J2000 coordinate system, and the last three components represent the satellite's velocity components in the J2000 coordinate system. T represents the matrix transpose symbol. Through this state vector representation, the second-order differential equation can be transformed into a system of first-order differential equations for numerical solution, which facilitates the application of the RK4 integrator.

[0081] The initial state of the integrator is determined by the satellite initial state vector prepared in step 1.2. supply.

[0082] The integration process begins, and the satellite state at times t0+h3, t0+2h3, ..., t0+n3h3 is calculated sequentially according to the set integration step size h3 (60 seconds in this case), until the target extrapolation time t is reached. f .

[0083] in The initial state vector of the satellite is represented by t; t0 represents the initial time; h3 represents the integration step size; t0+h3 represents the time of the first integration step; t0+2h3 represents the time of the second integration step; t0+n3h3 represents the time of the n3rd integration step; n3 represents the number of integration steps; t f Indicates the target extrapolation time.

[0084] Step 3: Calculate the main perturbation accelerations based on the ecliptic coordinate system and the Earth-fixed coordinate system, respectively;

[0085] In each integration step, the optimized fourth-order Runge-Kutta integrator needs to calculate the state derivatives k1, k2, k3, and k4 four times. Each calculation requires evaluating all forces (including central gravity and main perturbation forces) acting on the satellite in a specific state, and then updating the satellite state based on a weighted average of these derivatives.

[0086] Step 3.1: Determine the gravitational position vector of the Sun and Moon based on the ecliptic coordinate system, and then calculate the gravitational perturbation acceleration of the Sun and Moon based on the gravitational position vector of the Sun and Moon;

[0087] The gravitational perturbation accelerations of the Sun and Moon include the gravitational perturbation accelerations of the Sun and the Moon.

[0088] Calculate the Julian century number T1 starting from the standard epoch (January 1, 2000):

[0089]

[0090] Where T1 represents the Julian century number; JD represents the Julian day at the current moment.

[0091] The orbital elements of the Sun in the geocentric ecliptic coordinate system are:

[0092] a s =23454.803844a e ;

[0093]

[0094] Ω s =0;

[0095]

[0096] Solving Kepler's equations yields the aperitoneal angle E. s :

[0097] E s =M s +e s sinE s ;

[0098] Where a s Indicates the semi-major axis of the sun; a e Indicates the Earth's radius; e s i represents the eccentricity of the solar orbit; i6 represents the inclination of the solar orbit; Ω s ω represents the right ascension of the ascending node of the Sun's orbit. s Indicates the perigee argument of the solar orbit; M s E represents the mean aperitoneal angle of the sun. s T1 indicates the angle of approach to the sun; T1 indicates the Julian century number.

[0099] Calculate the rectangular coordinates of the Sun relative to the mean vernal equinox at that epoch:

[0100]

[0101] Where x s1 The x-axis coordinates of the sun relative to the mean vernal equinox of the epoch; y s1 The z-axis coordinates represent the sun's position relative to the mean vernal equinox; s1 i represents the z-axis coordinate of the Sun relative to the mean vernal equinox of the epoch; i6 represents the inclination of the Sun's orbit; ω s E represents the perigee argument of the solar orbit. sθ represents the angle of approach to the sun; cos represents the cosine function; sin represents the sine function.

[0102] Consider the angle of translation of the vernal equinox:

[0103]

[0104] Calculate the rectangular coordinates of the Sun relative to the orbital coordinate system:

[0105]

[0106] Where x s The x-axis coordinate of the sun in the orbital coordinate system; y s Represents the y-axis coordinate of the Sun in the orbital coordinate system; z s θ represents the z-axis coordinate of the Sun in the orbital coordinate system. s T1 represents the angle of translation of the vernal equinox; T1 represents the Julian century number; cos represents the cosine function; sin represents the sine function.

[0107] First, calculate the orbital parameters of the Sun and Moon in the ecliptic coordinate system. These parameters are the basis for calculating the Moon's position:

[0108]

[0109] Where ω m The perigee argument of the lunar orbit; λ m Indicates the ecliptic longitude of the moon; Ω m T1 represents the ecliptic longitude of the ascending node of the moon's orbit; T1 represents the Julian century number.

[0110] Auxiliary parameters for calculating the relative position of the moon:

[0111] α m =Ω m -0.0023sin2(λ m -λ s )+0.00213sin2(λ m -Ω m )+0.02724sin2(λ s -Ω m -0.00283sinM s +0.00047sinM m ;

[0112] β m1 =λ m +0.10976sinM m +0.00373sin2M m +0.02223sin[2(λ m -λ s )-Mm ]+0.01161sin2(λ m -λ s -0.00324sinM s ;

[0113] β m =β m1 -α m ;

[0114] β m2 =0.00241cos[2(λ m -Ω m )-2(λ m -λ s )]+0.089532511;

[0115] Δr1=cosβ m2 ;

[0116] Δr2=sinβ m2 ;

[0117] Where α m Auxiliary parameters representing the lunar orbit; λ m λ represents the ecliptic longitude of the moon. s Indicates the sun's longitude; Ω m M represents the ecliptic longitude of the ascending node of the moon's orbit; m β represents the mean anomaly of the Moon. m1 The first auxiliary parameter indicating the position of the moon; 2M m β represents twice the mean apogee of the moon; m The correction parameter representing the moon's position; β m2 The second auxiliary parameter represents the position of the moon; Δr1 represents the first component of the radial position correction of the moon; Δr2 represents the second component of the radial position correction of the moon; cos represents the cosine function; sin represents the sine function.

[0118] Calculate the rectangular coordinates of the Moon relative to the mean vernal equinox at that epoch:

[0119] x m1 =cosβ m cosα m -Δr1sinβ m sinα m ;

[0120] y m1 =cosi6(cosβ) m sinα m +Δr1sinβ m cosα m )-Δr2sini6sinβm ;

[0121] z m1 =sini6(cosβ) m sinα m +Δr1sinβ m cosα m )+Δr2cosi6sinβ m ;

[0122] Where x m1 The x-axis coordinates of the Moon relative to the mean vernal equinox; y m1 The z-axis coordinates represent the Moon's position relative to the mean vernal equinox; m1 The z-axis coordinates of the Moon relative to the mean vernal equinox are represented by: Δr1; Δr2; cos; sin; β. m Correction parameters for the moon's position; α m i represents the auxiliary parameters of the lunar orbit; i6 represents the inclination of the solar orbit.

[0123] Considering the translation angle of the vernal equinox, calculate the rectangular coordinates of the Moon relative to the orbital coordinate system:

[0124]

[0125] Where x m Represents the x-axis coordinate of the Moon in the orbital coordinate system; y m Represents the y-axis coordinate of the Moon in the orbital coordinate system; z m This represents the z-axis coordinate of the Moon in the orbital coordinate system; x m1 The x-axis coordinates of the Moon relative to the mean vernal equinox; y m1 The z-axis coordinates represent the Moon's position relative to the mean vernal equinox; m1 The z-axis coordinates of the Moon relative to the mean vernal equinox; cos represents the cosine function; sin represents the sine function; θ s This represents the angle of translation of the vernal equinox.

[0126] Transform the positions of the Sun and Moon from the orbital coordinate system to the J2000 coordinate system:

[0127]

[0128] in R represents the position vector of the sun in the J2000 coordinate system. OJ This represents the rotation matrix from the orbital coordinate system to the J2000 coordinate system; x s The x-axis coordinate of the sun in the orbital coordinate system; y sRepresents the y-axis coordinate of the Sun in the orbital coordinate system; z s This represents the z-axis coordinate of the Sun in the orbital coordinate system.

[0129]

[0130] in The x-axis represents the position vector of the Moon in the J2000 coordinate system. m Represents the x-axis coordinate of the Moon in the orbital coordinate system; y m Represents the y-axis coordinate of the Moon in the orbital coordinate system; z m R represents the z-axis coordinate of the Moon in the orbital coordinate system. OJ This represents the rotation matrix from the orbital coordinate system to the J2000 coordinate system.

[0131] Calculate the perturbation accelerations of the satellite by the Sun and Moon:

[0132]

[0133]

[0134] in μ represents the solar gravitational perturbation acceleration in the J2000 coordinate system. s Represents the solar gravitational constant; Represents the satellite's position vector; The vector represents the position of the sun; || represents the vector magnitude operator. μ represents the lunar gravitational perturbation acceleration in the J2000 coordinate system. m Represents the lunar gravitational constant; This represents the lunar position vector.

[0135] Step 3.2: Calculate the perturbation acceleration of the Earth's non-spherical gravitational field based on the Earth-fixed coordinate system;

[0136] In calculating the perturbation acceleration of the Earth's non-spherical gravitational field, the satellite's position in the Earth-fixed coordinate system is first calculated based on its position in the J2000 coordinate system (the J2000 coordinate system is the reference system for satellite orbit calculation; GPS, BeiDou, and other navigation systems provide J2000 coordinate system data for satellite positions). The satellite's position in the J2000 coordinate system is represented as follows:

[0137]

[0138] in y5 represents the satellite position vector; z5 represents the satellite's x-axis coordinate component in the J2000 coordinate system; y5 represents the satellite's y-axis coordinate component in the J2000 coordinate system; z5 represents the satellite's z-axis coordinate component in the J2000 coordinate system; T represents the matrix transpose symbol.

[0139] The formula for calculating the position of a satellite in the Earth-fixed coordinate system is as follows:

[0140]

[0141] in R represents the satellite's position vector in the Earth-fixed coordinate system. J2000→ECEF This represents the rotation matrix from the J2000 coordinate system to the Earth-fixed coordinate system; This represents the satellite's position vector in the J2000 coordinate system.

[0142] Then, calculate the spherical coordinates of the satellite in the Earth-fixed coordinate system:

[0143]

[0144] Where r6 represents the satellite's distance from the Earth's center; x6 represents the x-axis coordinate component in the Earth-fixed coordinate system; y6 represents the y-axis coordinate component in the Earth-fixed coordinate system; and z6 represents the z-axis coordinate component in the Earth-fixed coordinate system.

[0145] λ6 = atan2(y6, x6);

[0146] Where λ6 represents longitude; atan2 represents the arctangent function in the four quadrants; y6 represents the y-axis coordinate component in the Earth-fixed coordinate system; and x6 represents the x-axis coordinate component in the Earth-fixed coordinate system.

[0147]

[0148] Where φ6 represents latitude; atan represents the arctangent function; and z6 represents the z-axis coordinate component in the Earth-fixed coordinate system.

[0149] When calculating the perturbation force of Earth's non-spherical gravitational field, Legendre polynomials are needed to represent Earth's gravitational potential function. Earth's non-spherical gravitational field is usually represented by a spherical harmonic function expansion, where the normalized Legendre polynomial is the foundation for constructing the spherical harmonic function. Calculating the normalized Legendre polynomial... Using a recursive formula:

[0150]

[0151] in φ represents a normalized Legendre polynomial of degree 0; sin represents the sine function; φ6 represents latitude.

[0152]

[0153] in Let represent a normalized Legendre polynomial of degree 1.

[0154] For l5≥2:

[0155]

[0156] in Represents a normalized Legendre polynomial of degree l5; Represents a normalized Legendre polynomial of degree l5-1; This represents a normalized Legendre polynomial of degree l5-2.

[0157] For the associated Legendre function:

[0158]

[0159] in Let denote the first-order normalized correlation Legendre function.

[0160] For l6≥2:

[0161]

[0162] in This represents the l6th order l6 normalized correlation Legendre function; This represents the l6-1 order l6-1 normalized correlation Legendre function.

[0163] For l7≥2:

[0164]

[0165] in This represents the l7th order 1st normalized correlation Legendre function; This represents the l7-1th order first-order normalized correlation Legendre function; Let l7-2 be the first-order normalized correlation Legendre function.

[0166] For l8≥3, m7=2, 3, ..., l8-1:

[0167]

[0168] in This represents the l8th order m7 normalized correlation Legendre function; m7 represents the order. This represents the l8-1 order m7-1 normalized correlation Legendre function; This represents the l8-2 m7th order normalized correlation Legendre function.

[0169] In calculations of Earth's non-spherical gravitational field, the perturbation acceleration of the spacecraft under Earth's non-spherical gravitational field is:

[0170]

[0171] Where aN This represents the acceleration caused by the perturbation of Earth's non-spherical gravitational field; Represents the gradient operator; U represents the Earth's gravitational potential;

[0172] The spherical coordinate components projected in the Earth-fixed coordinate system are:

[0173]

[0174] Where μ represents the Earth's gravitational constant; r7 represents the satellite's geocentric distance; Represents the radial unit vector; cos represents the cosine function; φ7 represents latitude; ∑ represents the summation operator; n6 represents the order of the spherical harmonic function; m8 represents the order of the spherical harmonic function; R e π represents the Earth's radius; π represents pi. Represents the normalized Legendre function of order n6 m8; sin represents the sine function; The normalized gravitational coefficient C is represented by λ7; longitude is represented by λ7. The normalized gravitational coefficient S term; Represents a unit vector in the latitudinal direction; This represents the normalized Legendre function of order m8 and degree n6-1; It represents the unit vector of longitude direction.

[0175] If φ7 is When the above expression is odd, use the following expression to replace it:

[0176]

[0177] Where a N ′ indicates that φ7 is The acceleration due to the perturbation of the Earth's non-spherical gravitational field; μ represents the Earth's gravitational constant; r8 represents the distance from the satellite to the Earth's center; Represents the radial unit vector; ∑ represents the summation operator; n7 represents the degree of the spherical harmonic function; R e Indicates the Earth's radius; The normalized gravitational coefficient C is represented by m9; m9 represents the order of the spherical harmonic function.

[0178] Given the spacecraft's rectangular coordinates in the Earth-fixed rectangular coordinate system as (x7, y7, z7), and the corresponding unit vectors in the three directions as (i5, j4, k2), then:

[0179]

[0180] in x7 represents the radial unit vector; i5 represents the x-axis coordinate in the Earth-fixed coordinate system; y7 represents the y-axis coordinate in the Earth-fixed coordinate system; j4 represents the y-axis unit vector; z7 represents the z-axis coordinate in the Earth-fixed coordinate system; and k2 represents the z-axis unit vector.

[0181]

[0182] in y7 represents the unit vector in the longitude direction; y7 represents the y-axis coordinate in the Earth-fixed coordinate system; i5 represents the unit vector in the x-direction; x7 represents the x-axis coordinate in the Earth-fixed coordinate system; j4 represents the unit vector in the y-direction.

[0183]

[0184] in Represents a unit vector in the latitudinal direction; Represents a radial unit vector; It represents the unit vector of longitude direction.

[0185] In the calculation of the non-spherical gravitational field of the Earth, the perturbation acceleration of the non-spherical gravitational field of the Earth is first calculated in the spherical coordinate system because the spherical harmonic expansion of the Earth's gravitational potential function is most naturally expressed in the spherical coordinate system. Then, these accelerations need to be transformed to the Earth-fixed coordinate system, and finally to the J2000 coordinate system, so as to be synthesized with other perturbation accelerations.

[0186] Convert the acceleration in the spherical coordinate system to the acceleration in the Earth-fixed coordinate system:

[0187]

[0188] in This represents the x-axis acceleration component in the Earth-fixed coordinate system. This represents the y-direction acceleration component in the Earth-fixed coordinate system. This represents the z-axis acceleration component in the Earth-fixed coordinate system. Indicates the radial acceleration component; Indicates the latitudinal acceleration component; φ8 represents the acceleration component along the longitude direction; φ8 represents latitude; λ8 represents longitude. denoted as the Earth's non-spherical gravitational acceleration vector in the Earth-fixed coordinate system; cos represents the cosine function; sin represents the sine function.

[0189] Step 3.3: Calculate the atmospheric drag perturbation acceleration based on the Earth-fixed coordinate system;

[0190] Calculate the satellite's velocity relative to Earth's atmosphere:

[0191]

[0192] in This represents the satellite's velocity vector relative to the atmosphere; This represents the satellite's velocity vector in the J2000 coordinate system; This represents the Earth's rotational angular velocity vector; This represents the satellite's position vector in the J2000 coordinate system.

[0193] The magnitude of the Earth's rotational angular velocity vector is 7.29211515 × 10⁻⁶. -5 rad / s, direction along the Earth's axis of rotation.

[0194] A simplified exponential atmospheric density model is used:

[0195]

[0196] Where ρ1 represents the atmospheric density at the satellite's altitude; ρ0 represents the atmospheric density at the reference altitude; exp represents the exponential function; h5 represents the satellite's altitude; h0 represents the reference altitude; and H represents the atmospheric density elevation.

[0197] Satellite altitude h5 calculation:

[0198] h5 = r9 - R e -Δh(φ9);

[0199] Where Δh(φ9) represents the Earth's surface height correction term considering the Earth's oblateness; R e φ9 represents the Earth's radius; r9 represents the distance from the Earth's center to the satellite; φ9 represents latitude.

[0200] Calculate atmospheric drag perturbation acceleration:

[0201]

[0202] in C represents the atmospheric drag perturbation acceleration in the Earth-fixed coordinate system. D The vector represents the drag coefficient; S represents the satellite's cross-sectional area; m represents the satellite's mass; ρ1 represents the atmospheric density; and V represents the satellite's relative atmospheric velocity. This represents the satellite's velocity vector relative to the atmosphere.

[0203] Step 4: Transform the main perturbation accelerations to the J2000 coordinate system, and synthesize the total perturbation accelerations in the J2000 coordinate system;

[0204] The gravitational perturbation acceleration caused by the Sun and Moon is calculated using a model in the ecliptic coordinate system, and then the acceleration is calculated directly in the J2000 coordinate system, thus eliminating the need for additional coordinate transformation.

[0205]

[0206] in and These represent the solar gravitational perturbation acceleration and the lunar gravitational perturbation acceleration in the J2000 coordinate system, respectively.

[0207] Convert the acceleration in the Earth-fixed coordinate system to the J2000 coordinate system:

[0208]

[0209] in R represents the acceleration caused by the perturbation of the Earth's non-spherical gravitational field in the J2000 coordinate system. ECEF→J2000 This represents the rotation matrix from the Earth-fixed coordinate system to the J2000 coordinate system; This represents the acceleration vector of the Earth's non-spherical gravitational field in the Earth-fixed coordinate system.

[0210] Transform the atmospheric drag perturbation acceleration from the Earth-fixed coordinate system to the J2000 coordinate system:

[0211]

[0212] in R represents the atmospheric drag perturbation acceleration in the J2000 coordinate system. ECEF→J2000 This represents the rotation matrix from the Earth-fixed coordinate system to the J2000 coordinate system; This represents the atmospheric drag perturbation acceleration in the Earth-fixed coordinate system; and These represent the solar gravitational perturbation acceleration and the lunar gravitational perturbation acceleration in the J2000 coordinate system, respectively.

[0213]

[0214] in This represents the total perturbation acceleration in the J2000 coordinate system; This represents the atmospheric drag perturbation acceleration in the J2000 coordinate system; This represents the acceleration caused by the perturbation of the Earth's non-spherical gravitational field in the J2000 coordinate system.

[0215] Step 5: Construct a six-dimensional state derivative function based on the total perturbation acceleration and the central gravitational acceleration;

[0216] By combining the central gravitational acceleration and the total perturbation acceleration, a complete six-dimensional state derivative function is constructed:

[0217]

[0218] in This indicates that the satellite is at time t and in state t. The rate of change of the six-dimensional state at time, Represents a six-dimensional state vector, and Represents the velocity vector. This represents the satellite's position vector in the J2000 coordinate system; This represents the total perturbation acceleration in the J2000 coordinate system; This represents the satellite's acceleration vector in the J2000 coordinate system; μ represents the Earth's gravitational constant, with a value of 3.986004418 × 10⁻⁶. 14 m 3 / s 2 .

[0219] Constructed six-dimensional state derivative function The optimized fourth-order Runge-Kutta integrator used directly in step 2 is used to integrate the state vector and update the state vector:

[0220]

[0221] in Let represent the six-dimensional state vector of the nth integration step, where n represents the current integration step number; k1 represents the first slope estimate; k2 represents the second slope estimate; k3 represents the third slope estimate; k4 represents the fourth slope estimate; f represents the right-hand side function of the differential equation; t n This represents the time of the nth integration step.

[0222] Each evaluation of the six-dimensional state derivative function requires calling the perturbation force calculation method in step 3; performing coordinate transformation and acceleration synthesis in step 4, and then constructing the complete six-dimensional state derivative; determining whether the target has been reached until the extrapolated time t is reached. f If the target extrapolation time t is not reached f Then continue the iterative integration.

[0223] When t0+nh3≥t f When the integration process ends, the integration result of step n is directly output as the state vector at the target time.

[0224] This design enables a complete technology chain from perturbation calculation to integrator application, ensuring the continuity and efficiency of the calculation process.

[0225] Step 6: Output the satellite's state vector at the target time as the prediction result;

[0226] After completing the integration calculation, the system outputs the satellite's state vector at the target time:

[0227]

[0228] in Represents the satellite's final state vector; and The target extrapolation time t is respectively f The satellite's position and velocity vectors in the J2000 coordinate system.

[0229] The system can convert the extrapolation results into the following forms according to user needs:

[0230] The satellite's position and velocity at a specific point in time;

[0231] Satellite visibility time window;

[0232] Observation parameters such as satellite transit time, azimuth angle, and elevation angle;

[0233] Conditions for establishing a communication link between the satellite and the ground station.

[0234] Through the above six main steps, this embodiment realizes rapid calculation of satellite orbit extrapolation in a low-performance chip environment, meeting the real-time calculation needs of application scenarios such as satellite communication, navigation and remote sensing.

[0235] A fast satellite orbit extrapolation calculation system based on low-performance chips is used to execute the aforementioned fast satellite orbit extrapolation calculation method based on low-performance chips, including:

[0236] The orbit calculation model simplification module is used to construct a simplified satellite orbit calculation model based on the main perturbation force and set the satellite's initial state vector;

[0237] The integrator integration module is used to perform state vector integration using an optimized fourth-order Runge-Kutta integrator, and to determine the initial state of the optimized fourth-order Runge-Kutta integrator based on the satellite's initial state vector.

[0238] The perturbation force fast calculation module, when the optimized fourth-order Runge-Kutta integrator performs integration calculation, calculates the main perturbation acceleration based on the ecliptic coordinate system and the Earth-fixed coordinate system respectively;

[0239] The coordinate transformation optimization module is used to uniformly transform the main perturbation acceleration to the J2000 coordinate system through coordinate system transformation, synthesize the total perturbation acceleration in the J2000 coordinate system, and construct a six-dimensional state derivative function based on the total perturbation acceleration and the central gravitational acceleration.

[0240] The six-dimensional state derivative function construction module is used to perform cyclic integration on the six-dimensional state derivative function using the optimized fourth-order Runge-Kutta integrator to update the state vector until the target extrapolation time is reached.

[0241] The result output module is used to output the satellite's state vector at the target time as the prediction result.

[0242] Here, the present invention provides an implementation example:

[0243] This embodiment simulates a microsatellite communication terminal device developed by an aerospace technology company. The device uses an ARM Cortex-M4 processor (72MHz clock speed, no hardware floating-point unit) and has limited memory resources (64KB Flash storage, 20KB RAM). The device needs to calculate the satellite's orbital position in real time to establish a communication link, but traditional orbit calculation methods cannot meet the real-time requirements under these hardware conditions.

[0244] For a certain communication satellite constellation (orbital altitude 1200km), the model is simplified as follows based on the performance limitations of the terminal equipment:

[0245] Earth's non-spherical gravitational field model: simplified from the original 70th order to a 20th order;

[0246] Calculation of the positions of the Sun and Moon: Analytical formulas are used instead of the DE405 ephemeris.

[0247] Atmospheric density model: The exponential decay model is used instead of the NRLMSISE-00 model.

[0248] Example of initial state vector (J2000 coordinate system): Position vector The velocity vector is [-3862.725, 3278.941, 4834.279] km. The speed is [-5.7631, -4.0257, 2.3088] km / s; the initial time is 2023-06-15 12:00:00 UTC.

[0249] Through experimental comparison and analysis, we screened and tested various RK4 integrator coefficients. The performance comparison of different combinations of RK4 integrator coefficients is shown in Table 1:

[0250] Table 1: Performance Comparison of Different RK4 Integrator Coefficient Combinations

[0251] Coefficient combination 24-hour extrapolation position error (m) Calculation time (ms) Memory usage (KB) Classic RK4 285.6 620 14.2 Gill variant 263.8 645 14.5 Kutta variant 312.4 590 14.0 Optimize RK4 198.3 510 13.8 Fehlberg variant 176.2 890 16.3

[0252] Experiments show that the optimized fourth-order Runge-Kutta integrator achieves the best balance between computational accuracy and efficiency, improving accuracy by 30.6% compared to the classic RK4, while increasing computational speed by 17.7%.

[0253] Taking the calculation of perturbation forces in the non-spherical gravitational field of Earth as an example, the optimization process is as follows:

[0254] Pre-compute and store Legendre polynomial values;

[0255] Simplify the coefficient matrix of the Earth's non-spherical gravitational field (retaining only the significant influence terms);

[0256] Optimize the coordinate transformation algorithm.

[0257] Test results show that the optimized calculation speed for Earth's non-spherical gravitational field is increased by 4.8 times, while the accuracy loss is only 2.3%.

[0258] Implement a complete orbit extrapolation system on the target hardware platform:

[0259] Static memory allocation: 12KB of RAM space is pre-allocated;

[0260] Calculation process optimization: The calculation order of perturbation forces has been adjusted to prioritize effects from largest to smallest;

[0261] Adaptive parameter configuration: Dynamically adjust the complexity of the calculation model according to accuracy requirements.

[0262] The results of the technical effectiveness verification are as follows:

[0263] Table 2 shows the extrapolation accuracy and computational performance at different orbital altitudes:

[0264] Table 2: Extrapolation accuracy and computational performance at different orbital altitudes

[0265]

[0266] The system resource usage comparison is shown in Table 3:

[0267] Table 3: Comparison of System Resource Usage

[0268] index Optimization Algorithm Traditional Algorithm Improved proportions Code storage space (KB) 38.5 125.3 69.3％ Runtime RAM (KB) 11.7 42.8 72.7％ CPU utilization (%) 12.5 56.2 77.8％ Battery consumption (mAh / day) 18.3 83.6 78.1％

[0269] The actual application effects of the terminal devices are shown in Table 4:

[0270] Table 4: Actual Application Effects of Terminal Equipment

[0271] Application metrics Before optimization After optimization Improvement effect <![CDATA[Terminal device volume (cm 3 )]]> 320 224 Reduce by 30% Equipment cost (RMB) 1250 980 Decreased by 21.6% Battery life (h) 36 72 100% increase Satellite acquisition time (s) 12.5 3.2 Shortened by 74.4%

[0272] Experimental results show that the proposed fast satellite orbit extrapolation calculation method based on low-performance chips improves computational efficiency while maintaining accuracy, enabling efficient execution of satellite orbit calculations on resource-constrained low-performance chips. Application in a satellite communication terminal product demonstrates that this method can meet the needs of practical engineering applications, providing effective technical support for the miniaturization and cost reduction of satellite communication, navigation, and remote sensing terminals.

[0273] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A fast calculation method for satellite orbit extrapolation based on low-performance chips, characterized in that, include: A simplified satellite orbit calculation model based on the main perturbation force is constructed and the initial state vector of the satellite is set; An optimized fourth-order Runge-Kutta integrator is used for state vector integration, and the initial state of the optimized fourth-order Runge-Kutta integrator is determined based on the satellite's initial state vector. When the optimized fourth-order Runge-Kutta integrator performs integration calculations, it calculates the main perturbation accelerations based on both the ecliptic coordinate system and the Earth-fixed coordinate system. The main perturbation accelerations are uniformly transformed to the J2000 coordinate system through coordinate system transformation. The total perturbation acceleration is synthesized in the J2000 coordinate system. A six-dimensional state derivative function is constructed based on the total perturbation acceleration and the central gravitational acceleration. The optimized fourth-order Runge-Kutta integrator is used to iteratively integrate the six-dimensional state derivative function to update the state vector until the target extrapolation time is reached. The satellite's state vector at the target time is output as the prediction result; The optimized fourth-order Runge-Kutta integrator is expressed as follows: ； in Indicates the first The state vector of the step; Indicates the first The state vector of the step; Indicates the integration step size; This represents the first slope estimate; This indicates the second slope estimate; This represents the third slope estimate; This indicates the fourth slope estimate; Represent the function on the right-hand side of the differential equation; Indicates the first The time it takes to walk; The initial state vector of the satellite is represented as follows: ； in Represents the initial state vector. and These are the initial position vector and velocity vector of the satellite in the J2000 coordinate system; The initial state of the optimized fourth-order Runge-Kutta integrator, determined based on the satellite's initial state vector, refers to the state of the optimized fourth-order Runge-Kutta integrator from its initial time. Calculate sequentially according to the aforementioned integration step size. , , ..., The satellite status at any given moment, until the target extrapolation time is reached. The initial time of the optimized fourth-order Runge-Kutta integrator The state vector is .

2. The method for rapid calculation of satellite orbit extrapolation based on low-performance chips according to claim 1, characterized in that, The main perturbation forces include the perturbation force of the Earth's non-spherical gravitational field, the perturbation force of the Sun and Moon's gravity, and the perturbation force of atmospheric drag.

3. The method for rapid calculation of satellite orbit extrapolation based on low-performance chips according to claim 1, characterized in that, The calculation of the main perturbation accelerations based on the ecliptic coordinate system and the Earth-fixed coordinate system specifically includes determining the gravitational position vector of the Sun and Moon based on the ecliptic coordinate system and then calculating the gravitational perturbation acceleration of the Sun and Moon based on the gravitational position vector of the Sun and Moon, and calculating the perturbation acceleration of the Earth's non-spherical gravitational field and the atmospheric drag perturbation acceleration based on the Earth-fixed coordinate system.

4. The method for rapid calculation of satellite orbit extrapolation based on low-performance chips according to claim 3, characterized in that, Determining the gravitational position vector of the Sun and Moon based on the ecliptic coordinate system, and then calculating the gravitational perturbation acceleration of the Sun and Moon based on the gravitational position vector, specifically includes: Calculate the orbital elements of the Sun in the geocentric ecliptic coordinate system; Calculate the orbital parameters of the Moon in the ecliptic coordinate system; Calculate the rectangular coordinates of the Sun and Moon relative to the vernal equinox at an epoch; Transform the positions of the Sun and Moon from the ecliptic coordinate system to the J2000 coordinate system; Calculation of lunar and solar gravitational perturbation accelerations based on the J2000 coordinate system.

5. The method for rapid calculation of satellite orbit extrapolation based on low-performance chips according to claim 3, characterized in that, The calculation of the perturbation acceleration of the Earth's non-spherical gravitational field based on the Earth-fixed coordinate system specifically includes: Perform coordinate transformation between the J2000 coordinate system and the Earth-fixed coordinate system; Calculate the spherical coordinates of the satellite in the Earth-fixed coordinate system; Calculation of the acceleration perturbed by the Earth's non-spherical gravitational field based on the expansion of spherical harmonic functions.

6. The method for rapid calculation of satellite orbit extrapolation based on low-performance chips according to claim 3, characterized in that, The calculation of atmospheric drag perturbation acceleration based on the Earth-fixed coordinate system specifically includes: Calculate the satellite's velocity relative to the Earth's atmosphere; Calculate satellite altitude based on Earth's oblateness; Atmospheric density is calculated based on the exponential decay model. Calculate the atmospheric drag perturbation acceleration.

7. The method for rapid calculation of satellite orbit extrapolation based on low-performance chips according to claim 1, characterized in that, The time to reach the target extrapolation refers to , and when The state vector at the target time is directly output, represented as: ； in Represents the satellite's final state vector; and The target extrapolation time is respectively The satellite's position and velocity vectors in the J2000 coordinate system; Indicates the target extrapolation time; Indicates the current number of integration steps; Indicates the integration step size; Indicates the initial time.

8. A fast calculation system for satellite orbit extrapolation based on low-performance chips, characterized in that, A method for performing rapid satellite orbit extrapolation calculation based on low-performance chips as described in any one of claims 1-7, comprising: The orbit calculation model simplification module is used to construct a simplified satellite orbit calculation model based on the main perturbation force and set the satellite initial state vector; The integrator integration module is used to perform state vector integration using an optimized fourth-order Runge-Kutta integrator, and to determine the initial state of the optimized fourth-order Runge-Kutta integrator based on the satellite's initial state vector. The perturbation force fast calculation module, when the optimized fourth-order Runge-Kutta integrator performs integration calculation, calculates the main perturbation acceleration based on the ecliptic coordinate system and the Earth-fixed coordinate system respectively; The coordinate transformation optimization module is used to uniformly transform the main perturbation acceleration to the J2000 coordinate system through coordinate system transformation, synthesize the total perturbation acceleration in the J2000 coordinate system, and construct a six-dimensional state derivative function based on the total perturbation acceleration and the central gravitational acceleration. The six-dimensional state derivative function construction module is used to perform cyclic integration on the six-dimensional state derivative function using the optimized fourth-order Runge-Kutta integrator to update the state vector until the target extrapolation time is reached. The result output module is used to output the satellite's state vector at the target time as the prediction result.

Images