Calculation method for orbit of near-circle medium-and-low-orbit satellite

Abstract

The invention discloses a calculation method for an orbit of a near-circle medium-and-low-orbit satellite, and relates to the technical field of satellite orbit calculation. Real-time acquisition and classification processing of disturbance data are realized through parallel calculation and a multi-disturbance-factor preprocessing mechanism, so that orbit calculation can quickly respond to environmental changes, and the calculation accuracy is improved. The parallel computing framework independently processes the influence of different disturbance factors on each node, and quickly forms an orbit deviation trend, so that the satellite can perform adaptive adjustment more quickly when the orbit deviates, and meanwhile, an adaptive compensation model is adopted to dynamically adjust a compensation coefficient according to the actual orbit deviation, so that the adaptive compensation is realized. Adaptive adjustment is carried out on different disturbance coefficients, and automatic optimization compensation is carried out according to real-time orbit data, so that orbit prediction can adapt to rapidly changing disturbance conditions; when short-term disturbance is large, satellite orbit calculation precision can still be guaranteed, deviation accumulation is greatly reduced, extended Kalman filtering (EKF) is executed on each node, and orbit state estimation precision is improved.

Description

Technical Field

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

[0002] In the orbit calculation of low-earth orbit satellites, perturbations caused by factors such as the earth's gravity, the earth's non-spherical effect, atmospheric drag, and solar radiation pressure will cause significant 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 calculation needs to have strong response capabilities. However, traditional orbit calculation methods mostly use serial calculation and fixed-step iterative solution, resulting in a lag in capturing orbit perturbations. Commonly used finite difference methods, multi-step methods, or Runge-Kutta methods have strong serial dependencies, leading to low orbit prediction efficiency and difficulty in meeting the orbit calculation requirements of rapidly changing low-earth orbit satellites.

[0003] In addition, the impact of orbit perturbations on the satellite position is uncertain, so the orbit compensation algorithm needs to have strong adaptability. However, most traditional compensation algorithms are based on static models and are difficult to dynamically adapt to the rapid changes in environmental perturbations. For example, static models of the earth's non-spherical effect and atmospheric drag can provide relatively accurate corrections in a certain orbit segment, but when the satellite enters different altitudes or regions, the model may lag. In addition, compensation algorithms mostly perform orbit corrections based on periodically collected observation data and are difficult to cope with the rapidly changing orbit environment of low-earth orbit satellites.

[0004] Some solutions improve orbit calculation efficiency through parallel calculation, but traditional orbit calculation methods are difficult to fully apply the advantages of parallelization in dynamic compensation and correction. Serial calculation methods such as finite difference methods update point by point and are difficult to effectively utilize parallel computing resources, resulting in low calculation efficiency. At the same time, orbit compensation is mostly executed at fixed intervals rather than adaptively adjusted according to real-time status, further limiting the compensation effect and calculation efficiency. For example, during periods of large changes in atmospheric drag, the satellite orbit deviation increases rapidly, while the static model cannot adaptively compensate in real time, resulting in error accumulation and affecting the accuracy of orbit prediction. Therefore, there is an urgent need for a method for calculating the orbit of a near-circular medium-low earth orbit satellite to solve such problems. 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 medium-low earth orbit satellite to solve the problem that the current low-earth orbit satellite orbit calculation lacks a sufficiently flexible and fast dynamic compensation mechanism and is also difficult to efficiently utilize parallel computing resources.

[0007] 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 medium and low Earth orbit satellite, which includes,

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

[0010] Through the sensors on the ground station and the satellite, the initial position and velocity of the satellite are collected in real time, and at the same time, perturbation data is obtained; various types of collected data are preprocessed in parallel;

[0011] Step S2, parallel prediction of orbit state,

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

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

[0014] According to the orbit state data output in Step S2, each node analyzes the current orbit deviation trend, constructs an adaptive compensation model, and adaptively adjusts the perturbation compensation coefficient based on the real-time position and velocity;

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

[0016] In the orbit prediction after adaptive compensation in Step S3, each node uses the extended Kalman filter for error correction and fuses the orbit state data;

[0017] Step S5, output of orbit parameters and control input,

[0018] The orbit state parameters after correction and fusion in Step S4 are transmitted to the satellite control system and used as orbit correction and control input.

[0019] Further, in Step S1, the perturbation data includes atmospheric density, the Earth's non-spherical terms, and solar radiation pressure;

[0020] In Step S1, various types of collected data are classified according to the perturbation type and allocated to different computing nodes for parallel processing, and the preliminary influence of each perturbation factor on the orbit is estimated, and the classified multi-source perturbation data is output.

[0021] Further, in Step S1, the classification and allocation of perturbation data are as follows:

[0022] Based on the collected atmospheric density data ρ, the influence of atmospheric drag on the satellite orbit is estimated, and the atmospheric drag acceleration is calculated, where, a d represents the atmospheric drag acceleration, that is, the influence of atmospheric drag on the satellite, Cd where \(C_D\) represents the drag coefficient, \(A\) represents the surface area of the satellite facing the airflow, \(v\) represents the velocity of the satellite in orbit, and \(m\) represents the mass of the satellite;

[0023] Use \(J\) 2 and \(J\) 3 terms to represent the non - spherical gravitational terms of the Earth, and estimate the gravitational perturbation of the non - spherical effect of the Earth:

[0024] where, \(a\) g represents the acceleration of the non - spherical gravitational term of the Earth, representing the influence of the non - spherical effect of the Earth 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, \(J\) 2 、\(J\) 3 represent the second - order and third - order term coefficients of the non - spherical effect of the Earth, \(R\) represents the radius of the Earth, and \(\varphi\) represents the geocentric latitude of the satellite;

[0025] Preliminarily estimate the solar radiation pressure perturbation, where, \(a\) s represents the solar radiation pressure acceleration, representing the perturbation of the 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 mass of the satellite, and \(d\) represents the distance between the satellite and the sun;

[0026] Allocate the classified perturbation data (\(a\) d 、\(a\) g 、\(a\) s ) to different computing nodes, process the perturbation acceleration data in parallel on each node, and output the multi - source perturbation data matrix \(D\):

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

[0028] Furthermore, in step S2, the parallel prediction method of the orbit state is as follows:

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

[0030] At each time step \(\Delta t\), the node performs orbit state prediction:

[0031] where, represents the state vector at time step \(k + 1\), represents the state vector at time step k, and Δt represents the time step size. represents the acceleration perturbation caused by the non-spherical gravitational terms of the Earth. represents the acceleration perturbation caused by atmospheric drag. represents the acceleration perturbation caused by solar radiation pressure, where:

[0032]

[0033] Among them, GM represents the Earth's gravitational parameter. represents the distance between the satellite and the Earth's center, calculated using the modulus of, and J 2 and J 3 represent the second- and third-order term coefficients of the Earth's non-spherical effects, R represents the Earth's radius, and φ k represents the geocentric latitude of the satellite at time step k.

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

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

[0036] Furthermore, in step S3, the model automatically optimizes the step size and calculation accuracy to generate compensation parameters adapted to the current perturbation characteristics.

[0037] Furthermore, in step S3, according to the orbital state data output in step S2, each calculation node constructs an adaptive compensation model:

[0038] From the time series orbital state data output in step S2 calculate the orbital deviation at each time step Among them, represents the orbital deviation at time step k. represents the actual state vector at time step k. represents the reference orbital state vector at time step k.

[0039] Adopt the deviation change rate to evaluate the trend of orbital perturbation, and the deviation change rate is defined as:

[0040] Among them, represents the rate of change of the orbital deviation at time step k, and Δt represents the time step size;

[0041] Based on the rate of change of the deviation adaptively adjust the compensation coefficient in the compensation term of each perturbation factor:

[0042] Among them, C d，adj represents the adjusted atmospheric drag compensation coefficient, C d represents the initial atmospheric drag coefficient, α d represents the adjustment factor of the atmospheric drag compensation coefficient, which is determined according to the amplitude of the rate of change of the orbital deviation, represents the modulus value of the rate of change of the deviation, and other perturbation factors are processed using a similar method;

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

[0044] Among them, Δt new represents the new time step size, Δt represents the original time step size, and β represents the time step size adjustment coefficient, represents the modulus value of the rate of change of the deviation;

[0045] Dynamically determine the precision adjustment coefficient according to the deviation change trend: Among them, γ represents the calculation precision adjustment coefficient, and δ represents the precision adjustment factor, represents the modulus value of the rate of change of the deviation, and finally, the model generates the compensation parameters at the current moment based on the adjusted perturbation compensation coefficient C d，adj 、optimized time step size Δt new and precision coefficient γ.

[0046] Furthermore, in step S3, error correction data exchange is performed between nodes to further improve the calculation consistency; the errors of position, velocity, and each perturbation factor are comprehensively corrected to form the final orbital state.

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

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

[0049] Among them, represents the predicted state vector at time step k + 1, based on the state at time step k, represents the system state transition function, which describes the state under the control input and its change, Denote the adaptive compensation parameter output by step S3, which is used as the control input. Denote the process noise, which is assumed to follow a Gaussian distribution with zero mean.

[0050] Update the predicted covariance matrix: where denotes the predicted covariance matrix at time step k + 1, F k denotes the state transition matrix, which is equal to the Jacobian matrix of the state transition function with respect to , P k denotes the covariance matrix at time step k, Q k denotes the process noise covariance matrix;

[0051] Each node updates the state estimate based on the sensor observation value The observation prediction is:

[0052] where denotes the predicted observation value at time step k + 1, denotes the observation model, which describes the mapping relationship from the state to the observation space, denotes the measurement noise, which is assumed to follow a Gaussian distribution with zero mean;

[0053] The measurement residual, i.e., the difference between the prediction and the actual observation value: where denotes the measurement residual at time step k + 1, denotes the actual observation value.

[0054] Furthermore, in step S4, calculate the Kalman gain:

[0055] where K k+1 denotes the Kalman gain matrix at time step k + 1, H k+1 denotes the Jacobian matrix of the observation model, R k+1 denotes the measurement noise covariance matrix;

[0056] Update the state vector using the Kalman gain, where denotes the updated state vector at time step k + 1;

[0057] The covariance matrix is updated to where P k+1 denotes the updated covariance matrix at time step k + 1, and I denotes the identity matrix;

[0058] The state estimate of each node and the covariance P k+1After the exchange, the combined state after fusion is expressed as:

[0059] Among them, represents the fusion state vector at 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 the comprehensive error correction of the final position, velocity, and perturbation factor.

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

[0061] In the present invention, through parallel computing and multi-perturbation factor preprocessing mechanism, the real-time acquisition and classification processing of perturbation data are realized, enabling the orbit calculation to quickly respond to environmental changes. The parallel computing framework independently processes the influence of different perturbation factors on each node and quickly forms the trend of orbit deviation, enabling the satellite to make a more rapid adaptive adjustment when the orbit deviates.

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

[0063] In the present invention, the extended Kalman filter EKF is executed on each node, improving the accuracy of orbit state estimation. Through multi-node fusion and comprehensive correction of the errors of position, velocity, and each perturbation factor, the accuracy of orbit state calculation is greatly improved; the error correction data exchange and fusion are carried out among nodes to make the orbit state data of each node consistent, solving the problem of error accumulation caused by independent calculation of each node in the traditional serial method. Description of the Drawings

[0064] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0065] Figure 1 It is a schematic flow chart of the calculation method for the near-circular medium and low-earth orbit satellite orbit of the present invention. Detailed Embodiments

[0066] In order to make the above objects, features, and advantages of the present invention more obvious and understandable, the following will describe the detailed embodiments of the present invention with reference to the drawings in the specification.

[0067] In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention. However, the present invention may be practiced in other ways than those specifically described herein, and those skilled in the art can make similar extensions without departing from the spirit of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed below.

[0068] Secondly, as used herein, an "embodiment" or "embodiments" refer to specific features, structures, or characteristics that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it an individual or alternative embodiment that is mutually exclusive with other embodiments.

[0069] Embodiment 1, referring to Figure 1 , this embodiment provides a method for calculating the orbit of a near-circular medium and low Earth orbit satellite, including the following steps:

[0070] Step S1, parallel preprocessing of multiple perturbation factors,

[0071] Through sensors on the ground station and the satellite, the initial position and velocity of the satellite are collected in real time, and at the same time, perturbation data is obtained; various types of collected data are preprocessed in parallel;

[0072] In step S1, the perturbation data includes atmospheric density, Earth's non-spherical terms, and solar radiation pressure;

[0073] In step S1, various types of collected data are classified according to the perturbation type and distributed to different computing nodes for parallel processing to estimate the preliminary impact of each perturbation factor on the orbit, and multi-source perturbation data after classification is output;

[0074] In step S1, perturbation data classification and distribution are performed:

[0075] Based on the collected atmospheric density data ρ, the impact of atmospheric drag on the satellite orbit is estimated, and the atmospheric drag acceleration is calculated, where a d represents the atmospheric drag acceleration, that is, the impact 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 velocity of the satellite in the orbit, and m represents the mass of the satellite;

[0076] Use J 2 and J 3 terms to represent the non-spherical gravitational terms of the Earth, and estimate the gravitational perturbation of the non-spherical effect of the Earth:

[0077] where a grepresents the acceleration of the Earth's non-spherical gravitational terms, which represents the impact of the Earth's non-spherical effect on the satellite orbit. G represents the universal gravitational constant, M represents the mass of the Earth, r represents the distance between the satellite and the Earth's center, and J 2 and J 3 represent the second- and third-order term coefficients of the Earth's non-spherical effect. R represents the radius of the Earth, and φ represents the geocentric latitude of the satellite;

[0078] Preliminarily estimate the solar radiation pressure perturbation, where a s represents the solar radiation pressure acceleration, which represents the perturbation of the 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 mass of the satellite, and d represents the distance between the satellite and the sun;

[0079] Assign the classified perturbation data (a d , a g , a s ) to different computing nodes, process the perturbation acceleration data in parallel on each node, and output the multi-source perturbation data matrix D:

[0080] Step S2, parallel prediction of the orbit state,

[0081] Using the initial orbit parameters and the multi-source perturbation data generated in step S1 as input, based on the parallel computing framework, each computing node independently runs the orbit state iterative prediction algorithm; output the time series orbit state data and form the orbit deviation trend;

[0082] In step S2, each node regularly synchronizes the time step and the predicted state;

[0083] In step S2, the parallel prediction method of the orbit state is as follows:

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

[0085] At each time step Δt, the node performs orbit state prediction:

[0086] where, represents the state vector at time step k + 1, represents the state vector at time step k, and Δt represents the time step, represents the acceleration perturbation caused by the non-spherical gravitational terms of the Earth, represents the acceleration perturbation caused by atmospheric drag, represents the acceleration perturbation caused by solar radiation pressure, where:

[0087] where, GM represents the Earth's gravitational parameter, represents the distance between the satellite and the Earth's center, calculated using the modulus of J 2 、J 3 represent the second and third order term coefficients of the Earth's non-spherical effects, R represents the Earth's radius, φ k represents the geocentric latitude of the satellite at time step k;

[0088] where, C d represents the drag coefficient, A represents the satellite's windward area, ρ k represents the atmospheric density at time step k, represents the modulus of the satellite's velocity, m represents the satellite's mass;

[0089] where, P s represents the solar radiation constant, q represents the satellite's reflection coefficient, d k represents the distance between the satellite and the sun, represents the position vector of the sun at time step k;

[0090] Specifically, in the parallel computing framework, each node independently iterates to achieve independent prediction of the orbital state, and periodically synchronizes the time step and state. The entire system maintains a high degree of consistency; the parallel framework greatly improves the update frequency of the orbital state data and can capture the orbital deviation trend more quickly. When encountering sudden perturbations, it can respond quickly.

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

[0092] According to the orbital state data output in step S2, each node analyzes the current orbital deviation trend, constructs an adaptive compensation model, and adaptively adjusts the perturbation compensation coefficient based on the real-time position and velocity;

[0093] In step S3, the model automatically optimizes the step size and calculation accuracy to generate compensation parameters suitable for the current perturbation characteristics;

[0094] In step S3, according to the orbital state data output in step S2, each computing node constructs an adaptive compensation model:

[0095] From the time series orbital state data output in step S2 Among them, calculate the orbital deviation at each time step Among them, represents the orbital deviation at time step k, represents the actual state vector at time step k, represents the reference orbital state vector at time step k;

[0096] Adopt the deviation change rate to evaluate the trend of orbital perturbation, and the deviation change rate is defined as:

[0097] Among them, represents the orbital deviation change rate at time step k, and Δt represents the time step size;

[0098] Based on the deviation change rate Adaptive adjustment of the compensation coefficient in the compensation term of each perturbation factor:

[0099] Among them, C d,adj represents the adjusted atmospheric drag compensation coefficient, C d represents the initial atmospheric drag coefficient, α d represents the adjustment factor of the atmospheric drag compensation coefficient, determined according to the amplitude of the orbital deviation change rate, represents the modulus value of the deviation change rate, and other perturbation factors are processed using a similar method;

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

[0101] Among them, Δt new represents the new time step size, Δt represents the original time step size, and β represents the step size adjustment coefficient, represents the modulus value of the deviation change rate;

[0102] Dynamically determine the precision adjustment coefficient according to the deviation change trend: Among them, γ represents the calculation precision adjustment coefficient, and δ represents the precision adjustment factor, represents the modulus value of the deviation change rate. Finally, the model generates the compensation parameters at the current moment based on the adjusted perturbation compensation coefficient C d,adj , the optimized step size Δt new and the precision coefficient γ;

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

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

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

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

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

[0108] where, represents the predicted state vector at time step k + 1, based on the state at time step k, represents the system state transition function, describing the state under the control input and its change, represents the adaptive compensation parameter output by step S3 and serves as the control input, represents the process noise, assumed to follow a Gaussian distribution with a mean of zero;

[0109] Update the predicted covariance matrix: where, represents the predicted covariance matrix at time step k + 1, F k represents the state transition matrix, equal to the Jacobian matrix of the state transition function with respect to P k represents the covariance matrix at time step k, Q k represents the process noise covariance matrix;

[0110] Each node updates the state estimate based on the sensor observation value and the observation prediction is:

[0111] where, represents the predicted observation value at time step k + 1, Represents the observation model, which describes the mapping relationship from the state to the observation space, represents the measurement noise, assumed to follow a Gaussian distribution with zero mean;

[0112] The measurement residual, that is, the difference between the prediction and the actual observation value: where, represents the measurement residual at time step k + 1, represents the actual observation value;

[0113] In step S4, calculate the Kalman gain:

[0114] where k k+1 represents the Kalman gain matrix at time step k + 1, and H k+1 represents the Jacobian matrix of the observation model, and R k+1 represents the measurement noise covariance matrix;

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

[0116] The covariance matrix is updated to where P k+1 represents the updated covariance matrix at time step k + 1, and I represents the identity matrix;

[0117] The state estimates of each node and the covariance P k+1 are exchanged, and the integrated state after fusion is represented as:

[0118] where, represents the fused state vector at time step k + 1, N represents the number of parallel computing nodes, represents the updated state vector of the i-th node, and this fused state is used for the comprehensive error correction of the final position, velocity, and perturbation factor;

[0119] Specifically, the orbit state data output by each node is corrected and fused through the extended Kalman filtering technology, significantly improving the accuracy of orbit state prediction; the EKF algorithm can effectively filter noise in state estimation, making the orbit calculation under multi-perturbation conditions more accurate; at the same time, the error correction data between nodes is exchanged, making the system more consistent in overall calculation, effectively avoiding the error accumulation caused by independent calculation of each node, and finally forming a comprehensive orbit state, greatly reducing the cumulative error of various perturbation factors on the orbit position and velocity.

[0120] Step S5, orbit parameter output and control input,

[0121] Transmit the orbit state parameters after calibration and fusion in step S4 to the satellite control system for use as inputs for orbit correction and control.

[0122] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that the technical solutions of the present invention can be modified or equivalently replaced without departing from the spirit and scope of the technical solutions of the present invention, and they should all be covered within the scope of the claims of the present invention.

Claims

1. A method for calculating the orbit of a near-circular medium-low orbit satellite, characterized in that: include, Step S1, parallel preprocessing of multiple disturbance factors, Through sensors on the ground station and satellite, the initial position and speed of the satellite are collected in real time, and disturbance data is obtained at the same time; various 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, According to the track state 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 orbit state parameters corrected and fused in step S4 are transmitted to the satellite control system for use as orbit correction and control input.

2. The method for calculating the orbit of a near-circular low-mid-orbit satellite according to claim 1, characterized in that: In step S1, the disturbance data includes atmospheric density, earth non-spherical terms and solar radiation pressure; In step S1, the various types of 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-mid-orbit satellite according to claim 2, characterized in that: In step S1, the disturbance data is classified and allocated: 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 disturbance of the non-spherical effect of the Earth: Among them, a g represents the acceleration of the non-spherical gravitational term of the earth, representing the influence of the non-spherical effect of the earth 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 non-spherical effect of the earth, 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, which represents the disturbance of solar radiation to 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; The classified perturbation data (a d 、a g 、a s ) are distributed to different computing nodes, the disturbance acceleration data are processed in parallel on each node, and the multi-source disturbance data matrix D is output:

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

5. The method for calculating the orbit of a near-circular medium-low orbit satellite according to claim 4, characterized in that: 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 is calculated, J2 and 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; Among them, C d represents the drag coefficient, A represents the satellite's windward 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; 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, Represents the position vector of the sun at time step k.

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

7. The method for calculating the orbit of a near-circular medium-low orbit satellite according to claim 6, 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; 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 for 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: Among them, Δ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 accuracy coefficient γ.

8. The method for calculating the orbit of a near-circular low-mid-orbit satellite according to claim 7, 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 orbit state.

9. The method for calculating the orbit of a near-circular medium-low orbit satellite according to claim 8, 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 control input The changes below represents the adaptive compensation parameter outputted from step S3, as the control input, represents 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 observe the prediction 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 the observation space is: represents the measurement noise, which is assumed to follow a Gaussian distribution with a mean of zero; Measures residuals, which are the differences between predictions and actual observations: in, represents the measurement residual at time step k+1, represents the actual observed value.

10. The method for calculating the orbit of a near-circular medium-low orbit satellite according to claim 9, 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; The state vector is updated 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 the exchange, the comprehensive state after fusion is expressed as: in, represents the fusion state vector at 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