A calculation method for near-circular low- and medium-orbit satellite orbits

Abstract

The present invention discloses a method for calculating the orbit of a near-circular medium-low orbit satellite, and relates to the technical field of satellite orbit calculation. The present invention realizes real-time collection and classification processing of disturbance data through parallel computing and a multi-disturbance factor preprocessing mechanism, so that the orbit calculation can quickly respond to environmental changes. The parallel computing framework independently processes the influence of different disturbance factors at each node and quickly forms an orbit deviation trend, so that the satellite can make adaptive adjustments more quickly when the orbit deviates. At the same time, an adaptive compensation model is adopted to dynamically adjust the compensation coefficient according to the actual orbit deviation. By adaptively adjusting different disturbance coefficients and automatically optimizing compensation according to real-time orbit data, the orbit prediction can adapt to rapidly changing disturbance conditions. When the short-term disturbance is large, the satellite orbit calculation accuracy can still be guaranteed, and the deviation accumulation is greatly reduced. The extended Kalman filter (EKF) is executed at each node, thereby improving the orbit state estimation accuracy.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite orbit calculation, and in particular to a method for calculating the orbit of a near-circular medium-low orbit satellite. Background Art

[0002] In the orbit calculation of low-orbit satellites, disturbances from factors such as Earth's gravity, the Earth's non-spherical effect, atmospheric drag, and solar radiation pressure can cause significant and subtle changes in the satellite orbit in a short period of time. Therefore, in order to achieve high real-time performance and calculation accuracy in a high-frequency changing environment, orbit calculations must have strong responsiveness. However, traditional orbit calculation methods mostly use serial calculations and fixed-step iterative solutions, which have a lag in capturing orbital perturbations. The commonly used finite difference method, multi-step method, or RungeKutta method have low orbit prediction efficiency due to their strong serial dependence, making it difficult to meet the rapidly changing orbit calculation needs of low-orbit satellites.

[0003] In addition, there is uncertainty about the impact of orbital perturbations on satellite positions, so orbital compensation algorithms need to have strong adaptability. However, traditional compensation algorithms are mostly based on static models, which make it difficult to dynamically adapt to rapid changes in environmental perturbations. For example, although the static Earth non-spherical effect and atmospheric drag models can provide relatively accurate corrections in a certain orbital segment, when the satellite enters different altitudes or areas, the model may lag. In addition, compensation algorithms are mostly based on periodically collected observation data for orbit correction, which makes it difficult to cope with the rapidly changing orbital environment of low-orbit satellites.

[0004] Some solutions use parallel computing to improve orbit calculation efficiency, but traditional orbit calculation methods have difficulty fully utilizing the advantages of parallelization in dynamic compensation and correction. Serial calculation methods such as the finite difference method update point by point, making it difficult to effectively utilize parallel computing resources, resulting in low computational efficiency. At the same time, orbit compensation is often performed in fixed cycles rather than adaptively adjusting according to real-time conditions, further limiting the compensation effect and computational efficiency. For example, during periods of large changes in atmospheric drag, satellite orbit deviations increase rapidly, while static models are unable to adaptively compensate in real time, resulting in error accumulation and affecting the accuracy of orbit prediction. Therefore, a calculation method for near-circular low- and medium-orbit satellite orbits is urgently needed to solve this problem. Summary of the Invention

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

[0006] The present invention provides a method for calculating the orbit of a near-circular low-orbit satellite to solve the problem that the current low-orbit satellite orbit calculation lacks a sufficiently flexible and fast dynamic compensation mechanism, and it is also difficult to efficiently utilize parallel computing resources.

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

[0008] The present invention provides a method for calculating the orbit of a near-circular low- and medium-orbit satellite, which includes:

[0009] Step S1, parallel preprocessing of multiple disturbance factors,

[0010] Through sensors on the ground station and satellite, the initial position and velocity of the satellite are collected in real time, and disturbance data is obtained at the same time; all types of collected data are pre-processed in parallel;

[0011] Step S2, parallel prediction of track status,

[0012] Taking the initial orbit parameters and the multi-source disturbance data generated in step S1 as input, each computing node independently runs the orbit state iterative prediction algorithm based on the parallel computing framework; outputs the time series orbit state data and forms the orbit deviation trend;

[0013] Step S3, dynamic adjustment of the adaptive compensation model,

[0014] Based on the track status data output in step S2, each node analyzes the current track deviation trend, builds an adaptive compensation model, and adaptively adjusts the disturbance compensation coefficient based on the real-time position and speed;

[0015] Step S4, parallel error correction and data fusion,

[0016] In the orbit prediction after the adaptive compensation in step S3, each node uses the extended Kalman filter to perform error correction and integrate the orbit state data;

[0017] Step S5, orbital parameter output and control input,

[0018] The orbital state parameters corrected and fused in step S4 are transmitted to the satellite control system for use as orbit correction and control input.

[0019] Furthermore, in step S1, the disturbance data includes atmospheric density, earth non-spherical term and solar radiation pressure;

[0020] In step S1, the collected data are classified according to the disturbance type and assigned to different computing nodes for parallel processing. The preliminary impact of each disturbance factor on the orbit is estimated, and the classified multi-source disturbance data is output.

[0021] Furthermore, in step S1, the disturbance data is classified and allocated:

[0022] Based on the collected atmospheric density data ρ, the influence of atmospheric drag on satellite orbit is estimated and the atmospheric drag acceleration is calculated. Among them, a d represents the atmospheric drag acceleration, that is, the effect of atmospheric drag on the satellite, Cd represents the drag coefficient, A represents the surface area of ​​the satellite facing the airflow, v represents the speed of the satellite in orbit, and m represents the mass of the satellite;

[0023] Use J2 and J3 terms to represent the non-spherical gravitational term of the Earth and estimate the gravitational perturbation of the non-spherical effect of the Earth:

[0024] Among them, a g represents the acceleration of the Earth's non-spherical gravitational term, representing the impact of the Earth's non-spherical effect on the satellite orbit, G represents the gravitational constant, M represents the mass of the Earth, r represents the distance between the satellite and the center of the Earth, J2 and J3 represent the second-order and third-order coefficients of the Earth's non-spherical effect, R represents the radius of the Earth, and φ represents the geocentric latitude of the satellite;

[0025] Preliminary estimate of solar radiation pressure disturbance, Among them, a s Represents the solar radiation pressure acceleration, representing the disturbance of solar radiation on the satellite orbit, P s represents the solar radiation constant, A * represents the effective area of ​​the satellite facing the sun's radiation, q represents the satellite's reflection coefficient, m represents the satellite's mass, and d represents the distance between the satellite and the sun;

[0026] The classified perturbation data (a d 、a g 、a s ) are distributed to different computing nodes, and the disturbance acceleration data are processed in parallel on each node, and the multi-source disturbance data matrix D is output:

[0027] Furthermore, in step S2, each node periodically synchronizes the time step and the prediction state.

[0028] Furthermore, in step S2, the track state parallel prediction method is:

[0029] Determine the initial state based on position and velocity, and combine the multi-source disturbance data D generated in step S1 as the initial input for parallel prediction. in, represents the initial state vector, represents the initial position vector, represents the initial velocity vector;

[0030] At each time step Δt, the node makes a track state prediction:

[0031] in, represents the state vector at time step k+1, represents the state vector at time step k, Δt represents the time step, represents the acceleration disturbance caused by the non-spherical gravity of the Earth, represents the acceleration disturbance caused by atmospheric drag, represents the acceleration disturbance caused by solar radiation pressure, where:

[0032] Among them, GM represents the Earth's gravitational parameter, To express the distance between the satellite and the center of the Earth, use The modulus value calculation, J2, J3 represent the second-order and third-order coefficients of the earth's non-spherical effect, R represents the earth's radius, φ k represents the geocentric latitude of the satellite at time step k;

[0033] Among them, C d represents the drag coefficient, A represents the satellite's frontal area, ρ k represents the atmospheric density at time step k, represents the velocity modulus of the satellite, and m represents the mass of the satellite;

[0034] Among them, P s represents the solar radiation constant, q represents the satellite reflection coefficient, d k represents the distance between the satellite and the sun at time step k, Represents the position vector of the sun at time step k.

[0035] Furthermore, in step S3, the model automatically optimizes the step size and calculation accuracy to generate compensation parameters that adapt to the current disturbance characteristics.

[0036] Furthermore, in step S3, each computing node constructs an adaptive compensation model based on the track state data output in step S2:

[0037] The time series track status data output from step S2 , calculate the orbit deviation at each time step in, represents the orbit deviation at time step k, represents the actual state vector at time step k, represents the reference orbit state vector at time step k;

[0038] The deviation change rate is used to evaluate the trend of orbit disturbance, and the deviation change rate is defined as:

[0039] in, represents the orbit deviation change rate at time step k, and Δt represents the time step length;

[0040] Based on the rate of change of deviation Adaptively adjust the compensation coefficient in the compensation term of each disturbance factor:

[0041] Among them, C d,adj Indicates the adjusted atmospheric drag compensation coefficient, C d represents the initial atmospheric drag coefficient, α d The adjustment factor of the atmospheric drag compensation coefficient is determined according to the amplitude of the orbit deviation change rate. The modulus value represents the rate of change of the deviation. Other disturbance factors are processed in a similar way.

[0042] The adaptive model dynamically optimizes the time step according to the deviation change trend, and the new step size is updated as follows:

[0043] Where Δt new represents the new time step, Δt represents the original time step, β represents the step adjustment coefficient, The modulus value representing the rate of change of the deviation;

[0044] Dynamically determine the accuracy adjustment coefficient based on the deviation change trend: Among them, γ represents the calculation accuracy adjustment coefficient, δ represents the accuracy adjustment factor, The modulus of the deviation change rate is expressed as follows: d,adj , optimization step length Δt new The compensation parameters at the current moment are generated based on the and precision coefficient γ.

[0045] Furthermore, in step S3, error correction data is exchanged between nodes to further improve the calculation consistency; the errors of position, velocity and various disturbance factors are comprehensively corrected to form the final orbit state.

[0046] Furthermore, in step S4, parallel error correction and data fusion are performed:

[0047] At time step k, each node predicts the state at the next moment:

[0048] in, represents the predicted state vector at time step k+1, based on the state at time step k, Represents the system state transfer function, describing the state At the control input The changes below, represents the adaptive compensation parameter output from step S3, as the control input, represents the process noise, which is assumed to obey a Gaussian distribution with a mean of zero;

[0049] Update the forecast covariance matrix: in, represents the prediction covariance matrix at time step k+1, F k Represents the state transfer matrix, which is equal to the state transfer function about The Jacobian matrix, P k represents the covariance matrix at time step k, Q k represents the process noise covariance matrix;

[0050] Each node is based on sensor observations Update the state estimate and the observation forecast as:

[0051] in, represents the predicted observation value at time step k+1, Represents the observation model and describes the state The mapping relationship from t to observation space is: represents the measurement noise, which is assumed to obey a Gaussian distribution with mean zero;

[0052] Measures the residuals, which are the differences between the predictions and the actual observations: in, represents the measurement residual at time step k+1, represents the actual observed value.

[0053] Furthermore, in step S4, the Kalman gain is calculated:

[0054] Among them, K k+1 represents the Kalman gain matrix at time step k+1, H k+1 Represents the Jacobian matrix of the observation model, R k+1 represents the measurement noise covariance matrix;

[0055] Update the state vector using the Kalman gain, in, represents the updated state vector at time step k+1;

[0056] The covariance matrix is ​​updated to Among them, P k+1 represents the updated covariance matrix at time step k+1, and I represents the identity matrix;

[0057] State estimation of each node and covariance P k+1 After exchange, the comprehensive state after fusion is expressed as:

[0058] in, represents the fusion state vector of time step k+1, N represents the number of parallel computing nodes, Represents the updated state vector of the i-th node. This fused state is used for comprehensive error correction of the final position, velocity, and disturbance factor.

[0059] The beneficial effects of the present invention are:

[0060] The present invention uses parallel computing and a multi-disturbance factor preprocessing mechanism to achieve real-time collection and classification processing of disturbance data, enabling orbit calculations to quickly respond to environmental changes. The parallel computing framework independently processes the impact of different disturbance factors at each node and quickly forms an orbit deviation trend, allowing the satellite to make adaptive adjustments more quickly when the orbit shifts.

[0061] The present invention adopts an adaptive compensation model to dynamically adjust the compensation coefficient according to the actual orbit deviation. By adaptively adjusting different disturbance coefficients and automatically optimizing compensation according to real-time orbit data, the orbit prediction can adapt to rapidly changing disturbance conditions. Even when the short-term disturbance is large, the satellite orbit calculation accuracy can still be guaranteed, greatly reducing the accumulation of deviations.

[0062] The present invention executes an extended Kalman filter (EKF) at each node, thereby improving the accuracy of orbital state estimation. By performing multi-node fusion and comprehensive correction on the errors of position, velocity and various disturbance factors, the accuracy of orbital state calculation is greatly improved. Error correction data is exchanged and fused between nodes to ensure consistency of orbital state data at each node, thus solving the problem of error accumulation caused by independent calculation of each node in traditional serial methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0064] Figure 1 The figure is a flow chart of the method for calculating the orbit of a near-circular low- and medium-orbit satellite according to the present invention. DETAILED DESCRIPTION

[0065] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0066] In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those described herein. Those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0067] Secondly, the term "one embodiment" or "embodiment" herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in various places throughout this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive of other embodiments.

[0068] Example 1, reference Figure 1 This embodiment provides a method for calculating the orbit of a near-circular low- and medium-orbit satellite, comprising the following steps:

[0069] Step S1, parallel preprocessing of multiple disturbance factors,

[0070] Through sensors on the ground station and satellite, the initial position and velocity of the satellite are collected in real time, and disturbance data is obtained at the same time; all types of collected data are pre-processed in parallel;

[0071] In step S1, the disturbance data includes atmospheric density, earth non-spherical term and solar radiation pressure;

[0072] In step S1, the collected data are classified according to the disturbance type and distributed to different computing nodes for parallel processing. The preliminary impact of each disturbance factor on the orbit is estimated, and the classified multi-source disturbance data is output.

[0073] In step S1, disturbance data classification and distribution are performed:

[0074] Based on the collected atmospheric density data ρ, the influence of atmospheric drag on satellite orbit is estimated and the atmospheric drag acceleration is calculated. Among them, a d represents the atmospheric drag acceleration, that is, the effect of atmospheric drag on the satellite, C d represents the drag coefficient, A represents the surface area of ​​the satellite facing the airflow, v represents the speed of the satellite in orbit, and m represents the mass of the satellite;

[0075] Use J2 and J3 terms to represent the non-spherical gravitational term of the Earth and estimate the gravitational perturbation of the non-spherical effect of the Earth:

[0076] Among them, a grepresents the acceleration of the Earth's non-spherical gravitational term, representing the impact of the Earth's non-spherical effect on the satellite orbit, G represents the gravitational constant, M represents the mass of the Earth, r represents the distance between the satellite and the center of the Earth, J2 and J3 represent the second-order and third-order coefficients of the Earth's non-spherical effect, R represents the radius of the Earth, and φ represents the geocentric latitude of the satellite;

[0077] Preliminary estimate of solar radiation pressure disturbance, Among them, a s Represents the solar radiation pressure acceleration, representing the disturbance of solar radiation on the satellite orbit, P s represents the solar radiation constant, A represents the effective area of ​​the satellite facing the solar radiation, q represents the satellite reflection coefficient, m represents the satellite mass, and d represents the distance between the satellite and the sun;

[0078] The classified perturbation data (a d 、a g 、a s ) are distributed to different computing nodes, and the disturbance acceleration data are processed in parallel on each node, and the multi-source disturbance data matrix D is output:

[0079] Step S2, parallel prediction of track status,

[0080] Taking the initial orbit parameters and the multi-source disturbance data generated in step S1 as input, each computing node independently runs the orbit state iterative prediction algorithm based on the parallel computing framework; outputs the time series orbit state data and forms the orbit deviation trend;

[0081] In step S2, each node periodically synchronizes the time step and prediction status;

[0082] In step S2, the track state parallel prediction method is:

[0083] Determine the initial state based on position and velocity, and combine the multi-source disturbance data D generated in step S1 as the initial input for parallel prediction. in, represents the initial state vector, represents the initial position vector, represents the initial velocity vector;

[0084] At each time step Δt, the node makes a track state prediction:

[0085] in, represents the state vector at time step k+1, represents the state vector at time step k, Δt represents the time step, represents the acceleration disturbance caused by the non-spherical gravity of the Earth, represents the acceleration disturbance caused by atmospheric drag, represents the acceleration disturbance caused by solar radiation pressure, where:

[0086] Among them, GM represents the Earth's gravitational parameter, To express the distance between the satellite and the center of the Earth, use The modulus value calculation, J2, J3 represent the second-order and third-order coefficients of the earth's non-spherical effect, R represents the earth's radius, φ k represents the geocentric latitude of the satellite at time step k;

[0087] Among them, C d represents the drag coefficient, A represents the satellite's frontal area, ρ k represents the atmospheric density at time step k, represents the velocity modulus of the satellite, and m represents the mass of the satellite;

[0088] Among them, P s represents the solar radiation constant, q represents the satellite reflection coefficient, d k represents the distance between the satellite and the sun at time step k, Represents the position vector of the sun at time step k;

[0089] Specifically, in the parallel computing framework, each node autonomously iterates to achieve independent prediction of the orbital status, and regularly synchronizes time steps and status, so that the entire system maintains a high degree of consistency; the parallel framework greatly increases the update frequency of orbital status data, and can capture orbital deviation trends more quickly, and can respond quickly when encountering sudden disturbances.

[0090] Step S3, dynamic adjustment of the adaptive compensation model,

[0091] Based on the track status data output in step S2, each node analyzes the current track deviation trend, builds an adaptive compensation model, and adaptively adjusts the disturbance compensation coefficient based on the real-time position and speed;

[0092] In step S3, the model automatically optimizes the step size and calculation accuracy to generate compensation parameters that adapt to the current disturbance characteristics;

[0093] In step S3, each computing node constructs an adaptive compensation model based on the track state data output in step S2:

[0094] The time series track status data output from step S2 , calculate the orbit deviation at each time step in, represents the orbit deviation at time step k, represents the actual state vector at time step k, represents the reference orbit state vector at time step k;

[0095] The deviation change rate is used to evaluate the trend of orbit disturbance, and the deviation change rate is defined as:

[0096] in, represents the orbit deviation change rate at time step k, and Δt represents the time step length;

[0097] Based on the rate of change of deviation Adaptively adjust the compensation coefficient in the compensation term of each disturbance factor:

[0098] Among them, C d,adj Indicates the adjusted atmospheric drag compensation coefficient, C d represents the initial atmospheric drag coefficient, α d The adjustment factor of the atmospheric drag compensation coefficient is determined according to the amplitude of the orbit deviation change rate. The modulus value represents the rate of change of the deviation. Other disturbance factors are processed in a similar way.

[0099] The adaptive model dynamically optimizes the time step according to the deviation change trend, and the new step update method is:

[0100] Where Δt new represents the new time step, Δt represents the original time step, β represents the step adjustment coefficient, The modulus value representing the rate of change of the deviation;

[0101] Dynamically determine the accuracy adjustment coefficient based on the deviation change trend: Among them, γ represents the calculation accuracy adjustment coefficient, δ represents the accuracy adjustment factor, The modulus of the deviation change rate is expressed as follows: d,adj , optimization step length Δt new Generate the compensation parameters at the current moment based on the and precision coefficient γ;

[0102] Specifically, in the adaptive compensation model, each node dynamically adjusts the disturbance compensation coefficient according to the real-time orbit deviation, while optimizing the calculation step size and accuracy to enhance the adaptability of the model in changing environments. For different disturbance factors, the model automatically generates compensation parameters that adapt to the current disturbance characteristics through deviation trend analysis, which not only improves the robustness of the orbit calculation, but also achieves a balance between the model convergence speed and accuracy, and can still effectively reduce the orbit prediction error in complex disturbance environments.

[0103] Step S4, parallel error correction and data fusion,

[0104] In the orbit prediction after adaptive compensation in step S3, each node uses the extended Kalman filter to perform error correction and integrate the orbit state data. In step S3, error correction data is exchanged between nodes to further improve calculation consistency. The errors of position, velocity, and various disturbance factors are comprehensively corrected to form the final orbit state.

[0105] In step S4, parallel error correction and data fusion are performed:

[0106] At time step k, each node predicts the state at the next moment:

[0107] in, represents the predicted state vector at time step k+1, based on the state at time step k, Represents the system state transfer function, describing the state At the control input The changes below, represents the adaptive compensation parameter output from step S3, as the control input, represents the process noise, which is assumed to obey a Gaussian distribution with a mean of zero;

[0108] Update the forecast covariance matrix: in, represents the prediction covariance matrix at time step k+1, F k Represents the state transfer matrix, which is equal to the state transfer function about The Jacobian matrix, P k represents the covariance matrix at time step k, Q k represents the process noise covariance matrix;

[0109] Each node is based on sensor observations Update the state estimate and the observation forecast as:

[0110] in, represents the predicted observation value at time step k+1, Represents the observation model and describes the state The mapping relationship from t to observation space is: represents the measurement noise, which is assumed to obey a Gaussian distribution with mean zero;

[0111] Measures the residuals, which are the differences between the predictions and the actual observations: in, represents the measurement residual at time step k+1, represents the actual observed value;

[0112] In step S4, the Kalman gain is calculated:

[0113] Among them, K k+1 represents the Kalman gain matrix at time step k+1, H k+1 Represents the Jacobian matrix of the observation model, R k+1 represents the measurement noise covariance matrix;

[0114] Update the state vector using the Kalman gain, in, represents the updated state vector at time step k+1;

[0115] The covariance matrix is ​​updated to Among them, P k+1 represents the updated covariance matrix at time step k+1, and I represents the identity matrix;

[0116] State estimation of each node and covariance P k+1 After exchange, the comprehensive state after fusion is expressed as: in, represents the fusion state vector of time step k+1, N represents the number of parallel computing nodes, Represents the updated state vector of the i-th node. This fused state is used for comprehensive error correction of the final position, velocity, and disturbance factor;

[0117] Specifically, the orbital state data output by each node is corrected and integrated through the extended Kalman filtering technology, which significantly improves the accuracy of orbital state prediction; the EKF algorithm can effectively filter noise in state estimation, making the orbit calculation under multiple disturbance conditions more accurate; at the same time, the error correction data between nodes is exchanged to make the system more consistent in overall calculation, which can effectively avoid the error accumulation caused by independent calculation of each node, and finally form a comprehensive orbital state, which greatly reduces the cumulative error of various disturbance factors on orbital position and velocity.

[0118] Step S5, orbital parameter output and control input,

[0119] The orbital state parameters corrected and fused in step S4 are transmitted to the satellite control system for use as orbit correction and control input.

[0120] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A method for calculating the orbit of a near-circular low- and medium-orbit satellite, characterized by: include, Step S1, parallel preprocessing of multiple disturbance factors, Through sensors on the ground station and satellite, the initial position and velocity of the satellite are collected in real time, and disturbance data is obtained at the same time; all types of collected data are pre-processed in parallel; Step S2, parallel prediction of track status, Taking the initial orbit parameters and the multi-source disturbance data generated in step S1 as input, each computing node independently runs the orbit state iterative prediction algorithm based on the parallel computing framework; outputs the time series orbit state data and forms the orbit deviation trend; Step S3, dynamic adjustment of the adaptive compensation model, Based on the track status data output in step S2, each node analyzes the current track deviation trend, builds an adaptive compensation model, and adaptively adjusts the disturbance compensation coefficient based on the real-time position and speed; Step S4, parallel error correction and data fusion, In the orbit prediction after the adaptive compensation in step S3, each node uses the extended Kalman filter to perform error correction and integrate the orbit state data; Step S5, orbital parameter output and control input, The orbital state parameters corrected and fused in step S4 are transmitted to the satellite control system for use as orbit correction and control input; In step S1, disturbance data classification and distribution are performed: Based on the collected atmospheric density data ρ, the influence of atmospheric drag on satellite orbit is estimated and the atmospheric drag acceleration is calculated. Among them, a d represents the atmospheric drag acceleration, that is, the effect of atmospheric drag on the satellite, C d represents the drag coefficient, A represents the surface area of ​​the satellite facing the airflow, v represents the speed of the satellite in orbit, and m represents the mass of the satellite; Use J2 and J3 terms to represent the non-spherical gravitational term of the Earth and estimate the gravitational perturbation of the non-spherical effect of the Earth: Among them, a g represents the acceleration of the Earth's non-spherical gravitational term, representing the impact of the Earth's non-spherical effect on the satellite orbit, G represents the gravitational constant, M represents the mass of the Earth, r represents the distance between the satellite and the center of the Earth, J2 and J3 represent the second-order and third-order coefficients of the Earth's non-spherical effect, R represents the radius of the Earth, and Φ represents the geocentric latitude of the satellite; Preliminary estimate of solar radiation pressure disturbance, Among them, a s Represents the solar radiation pressure acceleration, representing the disturbance of solar radiation on the satellite orbit, P s represents the solar radiation constant, A * represents the effective area of ​​the satellite facing the sun's radiation, q represents the satellite's reflection coefficient, m represents the satellite's mass, and d represents the distance between the satellite and the sun; The classified perturbation data (a d 、a g 、a s ) are distributed to different computing nodes, and the disturbance acceleration data are processed in parallel on each node, and the multi-source disturbance data matrix D is output:

2. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 1, wherein: In step S1, the disturbance data includes atmospheric density, earth non-spherical term and solar radiation pressure; In step S1, the collected data are classified according to the disturbance type and assigned to different computing nodes for parallel processing. The preliminary impact of each disturbance factor on the orbit is estimated, and the classified multi-source disturbance data is output.

3. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 2, wherein: In step S2, each node periodically synchronizes the time step and prediction status.

4. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 3, wherein: In step S2, the track state parallel prediction method is: Determine the initial state based on position and velocity, and combine the multi-source disturbance data D generated in step S1 as the initial input for parallel prediction. in, represents the initial state vector, represents the initial position vector, represents the initial velocity vector; At each time step Δt, the node makes a track state prediction: in, represents the state vector at time step k+1, represents the state vector at time step k, Δt represents the time step, represents the acceleration disturbance caused by the non-spherical gravity of the Earth, represents the acceleration disturbance caused by atmospheric drag, represents the acceleration disturbance caused by solar radiation pressure, where: Among them, GM represents the Earth's gravitational parameter, To express the distance between the satellite and the center of the Earth, use The modulus value calculation, J2, J3 represents the second-order and third-order coefficients of the earth's non-spherical effect, R represents the earth's radius, Φ k represents the geocentric latitude of the satellite at time step k; Among them, C d represents the drag coefficient, A represents the surface area of ​​the satellite facing the airflow, and ρ k represents the atmospheric density at time step k, represents the velocity modulus of the satellite, and m represents the mass of the satellite; Among them, P s represents the solar radiation constant, q represents the satellite reflection coefficient, d k represents the distance between the satellite and the sun at time step k, Represents the position vector of the sun at time step k.

5. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 4, characterized in that: In step S3, the model automatically optimizes the step size and calculation accuracy to generate compensation parameters that adapt to the current disturbance characteristics.

6. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 5, characterized in that: In step S3, each computing node constructs an adaptive compensation model based on the track state data output in step S2: The time series track status data output from step S2 , calculate the orbit deviation at each time step in, represents the orbit deviation at time step k, represents the actual state vector at time step k, represents the reference orbit state vector at time step k; The deviation change rate is used to evaluate the trend of orbit disturbance, and the deviation change rate is defined as: in, represents the orbit deviation change rate at time step k, and Δt represents the time step length; Based on the rate of change of deviation Adaptively adjust the compensation coefficient in the compensation term of each disturbance factor: Among them, C d,adj Indicates the adjusted atmospheric drag compensation coefficient, C d represents the initial atmospheric drag coefficient, α d The adjustment factor of the atmospheric drag compensation coefficient is determined according to the amplitude of the orbit deviation change rate. The modulus value representing the rate of change of the deviation; The adaptive model dynamically optimizes the time step according to the deviation change trend, and the new step update method is: Where Δt new represents the new time step, Δt represents the original time step, β represents the step adjustment coefficient, The modulus value representing the rate of change of the deviation; Dynamically determine the accuracy adjustment coefficient based on the deviation change trend: Among them, γ represents the calculation accuracy adjustment coefficient, δ represents the accuracy adjustment factor, The modulus of the deviation change rate is expressed as follows: d,adj , optimization step length Δt new The compensation parameters at the current moment are generated based on the and precision coefficient γ.

7. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 6, characterized in that: In step S3, error correction data is exchanged between nodes; errors in position, velocity, and various disturbance factors are comprehensively corrected to form the final orbital state.

8. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 7, characterized in that: In step S4, parallel error correction and data fusion are performed: At time step k, each node predicts the state at the next moment: in, represents the predicted state vector at time step k+1, based on the state at time step k, Represents the system state transfer function, describing the state At the control input The changes below, represents the adaptive compensation parameter output from step S3, as the control input, represents the process noise, which is assumed to obey a Gaussian distribution with a mean of zero; Update the forecast covariance matrix: in, represents the prediction covariance matrix at time step k+1, F k Represents the state transfer matrix, which is equal to the state transfer function about The Jacobian matrix, P k represents the covariance matrix at time step k, Q k represents the process noise covariance matrix; Each node is based on sensor observations Update the state estimate and the observation forecast as: in, represents the predicted observation value at time step k+1, Represents the observation model and describes the state The mapping relationship from t to observation space is: represents the measurement noise, which is assumed to obey a Gaussian distribution with mean zero; Measures the residuals, which are the differences between the predictions and the actual observations: in, represents the measurement residual at time step k+1, represents the actual observed value.

9. The method for calculating the orbit of a near-circular low- and medium-orbit satellite according to claim 8, characterized in that: In step S4, the Kalman gain is calculated: Among them, K k+1 represents the Kalman gain matrix at time step k+1, H k+1 Represents the Jacobian matrix of the observation model, R k+1 represents the measurement noise covariance matrix; Update the state vector using the Kalman gain, in, represents the updated state vector at time step k+1; The covariance matrix is ​​updated to Among them, P k+1 represents the updated covariance matrix at time step k+1, and I represents the identity matrix; State estimation of each node and covariance P k+1 After exchange, the comprehensive state after fusion is expressed as: in, represents the fusion state vector of time step k+1, N represents the number of parallel computing nodes, Represents the updated state vector of the i-th node. This fused state is used for comprehensive error correction of the final position, velocity, and disturbance factor.

Images