A method for adaptive calibration of orbit control thrust based on inter-satellite distance and error compensation

Abstract

The application provides a kind of orbit control thrust self-adapting calibration method based on inter-satellite distance and error compensation, belongs to spacecraft orbit control thrust calibration technical field.According to satellite relative distance and velocity and relative dynamics equation, establish thrust calibration system;With satellite inter-satellite distance and relative velocity information as observation, establish state equation, evolution equation and observation equation;Use symmetric sampling method to select sample point set and its weight value, state update equation and observation update equation;According to residual covariance, calculate noise covariance adaptive matrix, the observation noise covariance and system noise covariance are adaptively updated, and the satellite relative orbit state estimation value is obtained;According to satellite relative orbit state estimation value, the least square solution of satellite orbit control thrust vector is calculated, and real-time calibration of orbit control thrust is realized.The method of the application overcomes the limitation of the source of measurement and control information, makes up for the error of engine and other hardware defects, and improves the thrust calibration accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of spacecraft orbit control thrust calibration, and relates to a method for adaptive calibration of orbit control thrust based on inter-satellite distance and error compensation. BACKGROUND

[0002] The increasing diversification of military and civilian application requirements makes the world's space powers pay more and more attention to the development and construction of large-scale constellations. In recent years, domestic and foreign aerospace research institutions and enterprises have successively launched large-scale constellation plans composed of tens of thousands of small satellites under the support of the state, and have accelerated the implementation and construction, which are applied in various fields such as national defense, military, global communication, Internet and the like. With the increasing scale of constellation missions, the number of satellites in orbit is increasing, and the satellite orbit control is facing many challenges in resource scheduling, implementation accuracy and the like.

[0003] Accurate orbit control thrust calibration is one of the key technologies for improving the control accuracy and autonomy of satellites. The commonly used orbit control thrust calibration methods are divided into two types: ground calibration test and on-orbit iterative correction method. The ground calibration test method usually carries out a large number of full-coverage experiments in a ground vacuum chamber to obtain the performance parameters of the thruster, and the test time is usually long and the development cost is high. Due to the great difference between the ground environment and the space environment, the reliability and accuracy of the calibration results cannot be guaranteed due to the influence of complex factors such as temperature, radiation and pressure. The use of high specific impulse electric thrusters for high-precision orbit control and maintenance makes it difficult to strictly complete the high-precision thrust calibration through ground tests due to the small level of thrust provided. The on-orbit iterative correction method is based on ground-based tracking and control information or satellite GPS information, and the orbit control is iterated through a dynamic model until the control orbit calculated by iteration is consistent with the control orbit obtained by precise tracking and control, and then the control effect and related parameters are obtained. The method is usually restricted by the location and number of ground stations, and requires continuous provision of high-precision orbit determination data from the ground, which cannot meet the demand for high-precision and full-real-time orbit parameter measurement.

[0004] At the same time, the environmental perturbation force as a system error is one of the important factors that interfere with the orbit evolution of near-earth orbit satellites, and its magnitude is comparable to that of the engine thrust. If the model of the perturbation force received by the satellite is not accurate, it will also cause a large thrust estimation error, leading to a decrease in filtering accuracy and even filter divergence, affecting the thrust calibration effect and having limitations.

[0005] Therefore, in order to enable the satellite to better perform complex orbit control, it is necessary to establish an accurate dynamic model, to compensate for the multi-source errors caused by the defects of the measurement equipment and the engine and the like, to estimate in real time on orbit through a filtering algorithm, and to realize accurate thrust calibration. How to get rid of the dependence on ground tracking and control stations, satellite GPS and the like for orbit determination data, to utilize the inter-satellite state information, to compensate for errors adaptively under limited measurement conditions, and to obtain accurate thrust estimation values is a bottleneck that needs to be broken through. SUMMARY

[0006] The present application solves the problem that the existing rail control thrust self-adaptive calibration method has many environmental uncertainties, large sample requirements, poor real-time performance, uncertainties caused by engine hardware defects, and low calibration accuracy.

[0007] The present application provides a satellite on-orbit autonomous, high-precision, multi-directional rail control thrust calibration method based on inter-satellite real-time ranging information and using adaptive estimation method for thrust error compensation.

[0008] The technical solution of the present application is as follows.

[0009] A rail control thrust self-adaptive calibration method based on inter-satellite distance and error compensation, comprising:

[0010] Step one, establishing a thrust calibration system according to the relative distance and velocity of the satellite in the reference star orbit coordinate system;

[0011] Step two, taking the inter-satellite distance and relative velocity information of the satellite as the observation, establishing the state equation of the unscented Kalman filter process, the evolution equation of the unscented Kalman filter process, and the observation equation of the unscented Kalman filter process;

[0012] Step three, based on the state equation of the unscented Kalman filter process, the evolution equation of the unscented Kalman filter process, and the observation equation of the unscented Kalman filter process of step two, according to the unscented change, using the symmetric sampling method to select the sample point set and its weight value, determining the state update equation and the observation update equation of the thrust calibration system;

[0013] Step four, based on the state update equation and the observation update equation of the thrust calibration system of step three, according to the residual covariance estimation value, calculating the noise covariance adaptive matrix, and through the noise covariance adaptive matrix, adaptively updating the observation noise covariance and the system noise covariance, and then obtaining the satellite relative orbit state estimation value;

[0014] Step five, according to the satellite relative orbit state estimation value of step four, calculating the least squares solution of the satellite rail control thrust vector, which is the on-orbit real-time calibration result of the satellite rail control thrust, so as to realize the adaptive calibration of the satellite rail control thrust.

[0015] Note: The thrust calibration system refers to the dynamics system used for satellite thrust calibration; it does not specifically refer to a certain equation, but rather the dynamics nature followed in the calibration process, mainly represented by the satellite motion equation and the relative dynamics equation based on the Gauss perturbation equation.

[0016] As an aspect of the present application, in step one, the relative distance and velocity of the satellite in the reference star orbit coordinate system include:

[0017] The relative distance and velocity of a satellite within the reference star's orbital coordinate system are defined as (x, y, z, v). x ,v y ,v z ), where x is the relative distance component along the x-direction of the coordinate system, y is the relative distance component along the y-direction of the coordinate system, z is the relative distance component along the z-direction of the coordinate system, and V x V is the component of relative velocity along the x-direction of the coordinate system. y V represents the relative velocity component along the y-direction of the coordinate system. z The relative velocity component along the z-direction of the coordinate system is represented as the difference between the orbital elements of the satellite and the reference star, i.e., the relative orbital elements, which are expressed as (Δa, Δe, Δi, ΔΩ, Δω, Δu), where Δa, Δe, Δi, ΔΩ, Δω, Δu are the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative argument of latitude, respectively.

[0018] In step one, a thrust calibration system is established, which includes satellite motion equations based on relative orbital elements and relative dynamic equations based on Gaussian perturbation equations.

[0019] The equations of motion for a satellite based on relative orbital elements are:

[0020]

[0021] Where x is the relative distance component along the x-direction of the coordinate system, y is the relative distance component along the y-direction of the coordinate system, z is the relative distance component along the z-direction of the coordinate system, and V x V is the component of relative velocity along the x-direction of the coordinate system. y V represents the relative velocity component along the y-direction of the coordinate system. z The relative velocity component along the z-direction of the coordinate system is given by: a, i, u, and V, which represent the satellite's semi-major axis, orbital inclination, latitudinal argument, and orbital velocity, respectively; u0 is the latitudinal argument of the satellite's initial position; Δa, Δe, Δi, ΔΩ, Δω, and Δu represent the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative latitudinal argument, respectively; the relative eccentricity vector between the satellite and the reference star is Δe = [Δe...]. x ,Δe y ] T Replace Δe, Δω; Δe x and Δe y These are the components of the relative eccentricity Δe along the x-direction and y-direction of the coordinate system, respectively.

[0022] The relative dynamic equation based on the Gaussian perturbation equation is:

[0023]

[0024] Where t is time, Let f be the derivative of a physical quantity with time, n be the angular velocity of the satellite orbit, and f be the angular velocity of the satellite orbit. u f r and f h Ω, ω, and u represent the satellite's orbital control thrust along the track, radial, and normal directions, respectively; a, e, i, Ω, ω, and u represent the satellite's semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and argument of latitude, respectively; V is the satellite's orbital velocity; e x and e y These are the components of the satellite eccentricity along the x-direction and y-direction of the coordinate system, respectively.

[0025] As one aspect of the present invention, in step two, the inter-satellite distance and relative velocity information of the satellites includes:

[0026] Using the start time of the inter-satellite measurement as the initial time t0, obtain t k The inter-satellite distance and relative velocity at time (k = 0, 1, ..., n) are used as the observations, and the satellite relative orbital elements are used as the parameters to be estimated by the filter.

[0027] Set the sampling time to Δt, and let X... k For t k The system state estimate at time X, i.e., the satellite's relative orbital elements. k =[Δa k ,Δi k ,ΔΩ k ,Δe x,k ,Δe y,k ,Δu k ] T X k+1 For t k+1 The system state estimate at time Z; k For t k System observations at any given time, namely the inter-satellite distances and relative velocities of the satellites. f(·) and h(·) are the system state transition function and the system measurement function, respectively; ω k ,υ k t k The system noise and observation noise at time t are independent zero-mean Gaussian white noises, and both are at time t. k The noise covariance matrices at time points Q and Q' are respectively k and R k .

[0028] As one aspect of the present invention, in step two, establishing the state equation, evolution equation, and observation equation of the unscented Kalman filtering process includes:

[0029] The relative dynamic equation based on the Gaussian perturbation equation is used as the state equation for the unscented Kalman filtering process;

[0030] The evolution equations for the unscented Kalman filtering process are established as follows:

[0031]

[0032] The observation equations for the unscented Kalman filtering process are established as follows:

[0033]

[0034] Among them, X k For t k The system state estimate at time X k+1 For t k+1 The system state estimate at time Z k For t k The system observation at time ω k ,υ k t k The system noise and observation noise at time points are given, f(·) and h(·) are the system state transition functions, Δa, Δe, Δi, ΔΩ, Δω, and Δu are the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative argument of latitude, respectively; a, i, u, and V are the satellite's semi-major axis, orbital inclination, argument of latitude, and orbital velocity, respectively; u0 is the argument of latitude of the satellite's initial position; Δe x and Δe y These are the components of the relative eccentricity Δe along the x-direction and the y-direction of the coordinate system, respectively.

[0035] Based on the state equation, evolution equation, and observation equation of the unscented Kalman filter process in step two, and according to the unscented changes, a symmetric sampling method is used to select a set of sample points and their weights to determine the state update equation and observation update equation of the thrust calibration system. This includes the following steps:

[0036] S031: Define t k The filtered prediction state vector of the satellite's relative orbital elements at time (k = 0, 1, ..., n) is: and its estimated variance is P k Initialize system parameters t0 and state vector X0, and set the system state estimate. and its estimated variance P0;

[0037] S032: Based on the state equation, evolution equation, and observation equation of the unscented Kalman filtering process in step two, and according to the unscented variation, the Sigma point set and its weight values ​​are selected using the symmetric sampling method; t is calculated. k-1 Sigma point at time:

[0038]

[0039] in, Indicates t k-1 The Sigma points are symmetrically sampled at each time step, with the superscript i indicating the i-th sampling point; λ is the primary scaling factor. and P k-1 t k-1 The filtered predicted state vector and its estimated variance relative to the orbital elements at time points, where n is the number of samples. This represents the Sigma point symmetrically sampled at time t0;

[0040] S033: Calculate t k The predicted state variables at time t and the corresponding prediction error covariance; define the subscript k|k-1 to indicate the state variables predicted at time t. k-1 Time prediction t k The time estimate; the transformation result of the Sigma point is calculated based on the system state transition function f(·).

[0041]

[0042] Establish the state update equations for the thrust calibration system and calculate the one-step predicted state variables. and the corresponding estimated covariance matrix P k|k-1 :

[0043]

[0044]

[0045] in, To predict the state variables in one step, P k|k-1 To estimate the covariance matrix, where i is the i-th sampling point, n is the number of samples, and Q... k-1 For t k-1 The system noise covariance matrix at time t. To calculate the transformation result of the Sigma point based on the system state transition function f(·), These are the weighted values ​​of the mean and covariance at the Sigma points, respectively.

[0046]

[0047] Where i is the i-th sampling point, n is the number of samples, α is the primary scaling factor, β is the secondary scaling factor, and λ is the primary scaling factor;

[0048] S034: Calculate t k One-step prediction of the Sigma point at time t; based on t k One-step prediction of state quantity at time t and prediction error covariance P k|k-1 Calculate the one-step prediction Sigma point:

[0049]

[0050] in, t represents the state quantity k The Sigma point is predicted at time step i, where the superscript i indicates the i-th sampling point; λ is the primary scaling factor, and n is the number of samples;

[0051] S035: Calculate t k One-step predicted observations at time t and their autocovariance and crosscovariance matrices; calculate the transformation results at the Sigma point based on the system measurement function h(·):

[0052]

[0053] in, t is the observed quantity k Predict the Sigma point one step at a time, where the superscript i represents the i-th sampling point and n represents the number of samples;

[0054] Establish the observation update equations for the thrust calibration system. To calculate the transformation result of the Sigma point based on the system state transition function f(·), the one-step predictive observation is calculated. One-step predictive observations The calculation formula is:

[0055]

[0056] in, t is the observed quantity k Predict the Sigma point one step at a time, where the superscript i represents the i-th sampling point and n represents the number of samples. The weighted average of the Sigma point means; the corresponding autocovariance matrix. for:

[0057]

[0058] Where i is the i-th sampling point, and n is the number of samples. This is the weighted value of the Sigma point covariance. t is the observed quantityk Predict the Sigma point in real time. For t k Time-predicted observations, R k For t k The observation noise covariance matrix at time step; the cross-covariance matrix between the predicted observations and the one-step predicted state variables. for:

[0059]

[0060] As one aspect of the present invention, based on the state update equation and observation update equation of the thrust calibration system in step three, a noise covariance adaptive matrix is ​​calculated according to the residual covariance estimate. The observation noise covariance and system noise covariance are adaptively updated using the noise covariance adaptive matrix, thereby obtaining the satellite relative orbital state estimate, including the following steps:

[0061] S041: Calculate the residual covariance estimate, and define the residual ε. k The difference between the observed actual value and the predicted value in one step:

[0062]

[0063] Among them, Z k For t k Observe the actual value at all times. For t k Predicted observations at specific times; residual covariance estimates The calculation formula is:

[0064]

[0065] Where, ε i For i i Time residual, The mean of the residuals, This is the difference between the residual and the mean of the residuals;

[0066] S042: Based on the state update equation and observation update equation of the thrust calibration system in step three, calculate the estimated value of the adaptive update of the observation noise covariance according to the residual covariance estimate; and set the adaptive matrix S of the observation noise covariance... k Represented as a diagonal matrix, the estimated value of the observation noise covariance is adaptively updated. for:

[0067]

[0068] Among them, R k For t k The observation noise covariance matrix at time step;

[0069] S043: Based on the residual covariance estimate, calculate the adaptive update estimate of the system noise covariance; the system noise covariance adaptive matrix Λ k It can be represented as a diagonal matrix, and the estimated value of the system noise covariance is updated adaptively. for:

[0070]

[0071] Among them, Q k-1 For t k-1 The system noise covariance matrix at time t;

[0072] S044: Update the observation self-covariance matrix and cross-covariance matrix based on the adaptively estimated system noise covariance and observation co-noise covariance; and The calculation formula is the same as step S035;

[0073] S045: Calculate the Kalman gain, state estimation results and their covariance matrix to obtain the satellite's relative orbital state estimate;

[0074] Kalman gain K k The calculation formula is:

[0075]

[0076] in, The cross-covariance matrix between the predicted observations and the one-step predicted state variables; For one step of predicting observations Autocovariance matrix; State estimation results and its covariance matrix P k The calculation formula is:

[0077]

[0078]

[0079] Among them, K k For Kalman gain, For t k The one-step predicted state quantity at time Z k For t k Observe the actual value at all times. For t k Predicted observations at any time, P k|k-1 for The estimated covariance matrix, for The autocovariance matrix.

[0080] As one aspect of the present invention, in step S042, calculating the adaptively updated estimate of the observation noise covariance based on the residual covariance estimate includes:

[0081] If m is the dimension of the observations, then the adaptive matrix S of the observation noise covariance is... k Represented as a diagonal matrix:

[0082]

[0083] Among them, S k (1) S represents the diagonal element of the first row of the adaptive matrix of observation noise covariance. k (2) S represents the diagonal element of the second row of the adaptive matrix of observation noise covariance. k (m) represents the diagonal element of the m-th row of the observation noise covariance adaptive matrix; the diagonal element S k The formula for calculating (i), i = 1, 2, ..., m, is:

[0084]

[0085] in, is the residual covariance estimate; μ is an adjustable parameter, and 1≤μ≤1000; This is the weighted value of the Sigma point covariance; t represents the observed quantity k Predict Sigma points in real time; For t k Time-predicted observations; N k (i,i) and R k (i,i) are matrices N and N respectively. k and R k The diagonal elements in the i-th row and i-th column; the estimated value of the adaptively updated observation noise covariance is:

[0086]

[0087] As one aspect of the present invention, in step S043, calculating the adaptively updated estimate of the system noise covariance based on the residual covariance estimate includes:

[0088] The adaptive matrix Λ of the system noise covariance k Represented as a diagonal matrix:

[0089]

[0090] Among them, Λ k (1) represents the diagonal element of the first row of the adaptive noise covariance matrix, Λ k (2) represents the diagonal element of the second row of the adaptive matrix of system noise covariance, Λk (n) represents the diagonal element of the nth row of the system noise covariance adaptive matrix; the diagonal element Λ k The formula for calculating (i), i = 1, 2, ..., n, is:

[0091]

[0092] in, and For each of the matrices and The diagonal elements in the i-th row and i-th column, where m is the observation dimension and n is the number of samples; the estimated value of the adaptively updated system noise covariance is:

[0093]

[0094] As one aspect of the present invention, based on the estimated satellite relative orbital state value from step four, the least squares solution of the satellite orbital control thrust vector is calculated, which is the on-orbit real-time calibration result of the satellite orbital control thrust, thereby realizing the on-orbit real-time calibration of the satellite orbital control thrust, including:

[0095] Solve for the estimated satellite orbit control thrust. The least squares solution form:

[0096]

[0097] in, For t k The components of the satellite along the track, radial, and normal directions at time (k = 0, 1, ..., n), V k It is t k Constantly refer to the orbital velocity of the star. For the state estimation results, the coefficient matrix A k The calculation formula is:

[0098]

[0099] Among them, a k i k and u k t k By referencing the semi-major axis, orbital inclination, and latitude argument of the satellite at all times, and with Δt representing the sampling time, the least squares solution of the satellite orbit control thrust vector is calculated, which is the real-time on-orbit calibration result of the satellite orbit control thrust, thus realizing the real-time on-orbit calibration of the satellite orbit control thrust.

[0100] Thus, adaptive on-orbit real-time calibration of satellite orbit control thrust based on inter-satellite distance information and multi-source error compensation has been achieved.

[0101] The advantages of this invention compared to the prior art are:

[0102] (1) The method of the present invention is different from the traditional method. It utilizes the relative distance and velocity information measured autonomously between satellites, without the need for ground control stations, on-board GPS, etc. to provide orbit determination data, thus overcoming the limitations of control information sources and improving the autonomy and reliability of thrust calibration.

[0103] (2) The method of the present invention is based on small sample inter-satellite distance data and uses statistical parameters to construct sampling points and corresponding weights to estimate the satellite orbit state, which effectively saves the workload of a large number of sample experiments in the early stage and meets the timeliness requirements of on-orbit thrust calibration.

[0104] (3) The method of the present invention addresses the problem of amplified iteration error in the traditional filtering estimation process. It adaptively matches the process noise covariance matrix based on the error between the actual value and the predicted value to achieve real-time on-orbit estimation. This improves system stability and significantly enhances thrust calibration accuracy. Attached Figure Description

[0105] Figure 1 This is a flowchart of the adaptive calibration method for track control thrust of the present invention;

[0106] Figure 2 This is a flowchart of the adaptive unscented Kalman filter algorithm of the present invention. Detailed Implementation

[0107] The present invention will now be described in further detail with reference to the accompanying drawings.

[0108] This invention utilizes the relative state information measured autonomously between satellites, eliminating the need for orbit determination data support from ground control stations or onboard GPS. It combines filtering methods to estimate satellite thrust. Based on standard unscented Kalman filtering, it incorporates the residual covariance matching principle and the filter divergence criterion, and introduces an adaptive factor matrix to further compensate for errors caused by hardware defects such as engine thrust, reduce the impact of instability in the filtering process, and improve thrust calibration accuracy.

[0109] The present invention provides an adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation, the flowchart of which is shown below. Figure 1 As shown, it includes the following steps:

[0110] Step 1: Establish a thrust calibration system based on the satellite's relative distance and velocity in the reference star's orbital coordinate system. The thrust calibration system does not refer to a specific equation, but rather to the dynamic system used for satellite thrust calibration, which is the dynamic essence followed during the calibration process.

[0111] In step one, the satellite's relative distance and velocity within the reference star's orbital coordinate system include:

[0112] A reference satellite is selected for thrust calibration, and its orbital state is required to be known. The relative distance and relative velocity between the satellite and the reference satellite are obtained through measuring instruments such as optical rangefinders and are used as observations. The changes in the relative distance and relative velocity between the satellite and the reference satellite are directly determined by the difference in orbital elements between the satellite and the reference satellite.

[0113] The relative distance and velocity of a satellite within the reference star's orbital coordinate system are defined as (x, y, z, v). x ,v y ,v z ), where x is the relative distance component along the x-direction of the coordinate system, y is the relative distance component along the y-direction of the coordinate system, z is the relative distance component along the z-direction of the coordinate system, and V x V is the component of relative velocity along the x-direction of the coordinate system. y V represents the relative velocity component along the y-direction of the coordinate system. z The relative velocity component along the z-direction of the coordinate system is represented as the difference between the orbital elements of the satellite and the reference star, i.e., the relative orbital elements, which are expressed as (Δa, Δe, Δi, ΔΩ, Δω, Δu), where Δa, Δe, Δi, ΔΩ, Δω, Δu are the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative argument of latitude, respectively.

[0114] In step one, the thrust calibration system includes satellite motion equations based on relative orbital elements and relative dynamic equations based on Gaussian perturbation equations;

[0115] The equations of motion for a satellite based on relative orbital elements are:

[0116]

[0117] Where x is the relative distance component along the x-direction of the coordinate system, y is the relative distance component along the y-direction of the coordinate system, z is the relative distance component along the z-direction of the coordinate system, and V x V is the component of relative velocity along the x-direction of the coordinate system. y V represents the relative velocity component along the y-direction of the coordinate system. z The relative velocity component along the z-direction of the coordinate system is given by: a, i, u, and V, which represent the satellite's semi-major axis, orbital inclination, latitudinal argument, and orbital velocity, respectively; u0 is the latitudinal argument of the satellite's initial position; Δa, Δe, Δi, ΔΩ, Δω, and Δu represent the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative latitudinal argument, respectively; the relative eccentricity vector between the satellite and the reference star is Δe = [Δe...]. x ,Δe y ] T Instead of Δe and Δω; for ease of calculation, the relative eccentricity vector between the satellite and the reference star is used: Δe = [Δe x ,Δe y] T Replace Δe, Δω, Δe x and Δe y These are the components of the relative eccentricity Δe along the x-direction and y-direction of the coordinate system, respectively; the calculation formula is:

[0118]

[0119] Where e and ω are the satellite's eccentricity and perigee depression angle, respectively. c and ω c These are the eccentricity and perigee depression of the reference star, respectively.

[0120] The orbital control thrust vector applied by the satellite will directly cause a change in the satellite's orbital elements relative to the reference satellite; for orbital control of a satellite in a near-circular orbit, the relative dynamic equation based on the Gaussian perturbation equation is:

[0121]

[0122] Where t is time, Let f be the derivative of a physical quantity with time, n be the angular velocity of the satellite orbit, and f be the angular velocity of the satellite orbit. u f r and f h Ω, ω, and u represent the satellite's orbital control thrust along the track, radial, and normal directions, respectively; a, e, i, Ω, ω, and u represent the satellite's semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and argument of latitude, respectively; V is the satellite's orbital velocity; e x and e y These are the components of the satellite eccentricity along the x-direction and y-direction of the coordinate system, respectively.

[0123] Step 2: Using the inter-satellite distance and relative velocity information of the satellites as observations, establish the state equation, evolution equation, and observation equation of the unscented Kalman filter process.

[0124] In step two, the inter-satellite distance and relative velocity information of the satellites includes:

[0125] Using the start time of the inter-satellite measurement as the initial time t0, obtain t k The inter-satellite distance and relative velocity at time (k = 0, 1, ..., n) are used as the observations, and the satellite relative orbital elements are used as the parameters to be estimated by the filter.

[0126] If the sampling time is set to Δt, then the thrust calibration system in step one can be discretized as follows:

[0127]

[0128] Where, let X k For tk The system state estimate at time X, i.e., the satellite's relative orbital elements. k =[Δa k ,Δi k ,ΔΩ k ,Δe x,k ,Δe y,k ,Δu k ] T X k+1 For t k+1 The system state estimate at time Z; k For t k System observations at any given time, namely the inter-satellite distances and relative velocities of the satellites. f(·) and h(·) are the system state transition function and the system measurement function, respectively; ω k ,υ k t k The system noise and observation noise at time t are independent zero-mean Gaussian white noises, and both are at time t. k The noise covariance matrices at time points Q and Q' are respectively k and R k ;

[0129] The state equation of the unscented Kalman filtering process is the relative dynamic equation based on Gaussian perturbation in step one;

[0130] Step two involves establishing the state equation, evolution equation, and observation equation for the unscented Kalman filter process, including:

[0131] The relative dynamic equation based on the Gaussian perturbation equation is used as the state equation for the unscented Kalman filtering process;

[0132] The evolution equations for the unscented Kalman filtering process are established as follows:

[0133]

[0134] The observation equations for the unscented Kalman filtering process are established as follows:

[0135]

[0136] Among them, X k For t k The system state estimate at time X k+1 For t k+1 The system state estimate at time Z k For t k The system observation at time ω k ,υ k t kThe system noise and observation noise at time points are given, f(·) and h(·) are the system state transition functions, Δa, Δe, Δi, ΔΩ, Δω, and Δu are the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative argument of latitude, respectively; a, i, u, and V are the satellite's semi-major axis, orbital inclination, argument of latitude, and orbital velocity, respectively; u0 is the argument of latitude of the satellite's initial position; Δe x and Δe y These are the components of the relative eccentricity Δe along the x-direction and the y-direction of the coordinate system, respectively;

[0137] Step 3: Based on the state equation, evolution equation, and observation equation of the unscented Kalman filter process in Step 2, and according to the unscented changes, the symmetric sampling method is used to select the sample point set and its weight values ​​to determine the state update equation and observation update equation of the thrust calibration system.

[0138] Unscented Kalman filtering employs unscented transformation, which calculates a series of Sigma sample points and performs a nonlinear function transformation. The transformation result and the corresponding weights are used to calculate the Gaussian distribution, thereby handling the nonlinear propagation problem of mean and covariance.

[0139] Includes the following steps:

[0140] S031: Define t k The filtered prediction state vector of the satellite's relative orbital elements at time (k = 0, 1, ..., n) is: and its estimated variance P k ,in, It consists of the relative orbital element filtered estimates of the satellite. P k The system is a 6×6 matrix; the system parameters t0 and the state vector X0 are initialized, and the system state estimate is set. and its estimated variance P0; when t k At time t0, the filtered predicted state vector is And its estimated variance P0 are:

[0141]

[0142]

[0143] Where E(·) is the mathematical expectation;

[0144] S032: Based on the state equation, evolution equation, and observation equation of the unscented Kalman filtering process in step two, and according to the unscented variation, the Sigma point set and its weight values ​​are selected using the symmetric sampling method; t is calculated. k-1Sigma point at time:

[0145]

[0146] in, Indicates t k-1 The Sigma points are symmetrically sampled at each time step, with the superscript i indicating the i-th sampling point; λ is the primary scaling factor. and P k-1 t k-1 The filtered predicted state vector and its estimated variance relative to the orbital elements at time points, where n is the number of samples. This represents the Sigma point symmetrically sampled at time t0;

[0147] λ=α 2 (n+κ)-n,

[0148] Where α is the primary scaling factor, which determines the distribution of Sigma points, and is typically taken as 10. -4 ≤α≤1; κ is the third-level scaling factor, with a value of κ=3-n;

[0149] S033: Calculate t k The predicted state variables at time t and the corresponding prediction error covariance; define the subscript k|k-1 to indicate the state variables predicted at time t. k-1 Time prediction t k The time estimate; the transformation result of the Sigma point is calculated based on the system state transition function f(·).

[0150]

[0151] Establish the state update equations for the thrust calibration system and calculate the one-step predicted state variables. and the corresponding estimated covariance matrix P k|k-1 :

[0152]

[0153]

[0154] in, To predict the state variables in one step, P k|k-1 To estimate the covariance matrix, where i is the i-th sampling point, n is the number of samples, and Q... k-1 For t k-1 The system noise covariance matrix at time t. To calculate the transformation result of the Sigma point based on the system state transition function f(·), These are the weighted values ​​of the mean and covariance at the Sigma points, respectively.

[0155]

[0156] Where i is the i-th sampling point, n is the number of samples, α is the primary scaling factor, β is the secondary scaling factor, and λ is the primary scaling factor;

[0157] S034: Calculate t k One-step prediction of the Sigma point at time t; based on t k One-step prediction of state quantity at time t and prediction error covariance P k|k-1 Calculate the one-step prediction Sigma point:

[0158]

[0159] in, t represents the state quantity k The Sigma point is predicted at time step i, where the superscript i indicates the i-th sampling point; λ is the primary scaling factor, and n is the number of samples;

[0160] S035: Calculate t k One-step predicted observations at time t and their autocovariance and crosscovariance matrices; calculate the transformation results at the Sigma point based on the system measurement function h(·):

[0161]

[0162] in, t is the observed quantity k Predict the Sigma point one step at a time, where the superscript i represents the i-th sampling point and n represents the number of samples;

[0163] Establish the observation update equations for the thrust calibration system. To calculate the transformation result of the Sigma point based on the system state transition function f(·), the one-step predictive observation is calculated. One-step predictive observations The calculation formula is:

[0164]

[0165] in, t is the observed quantity k Predict the Sigma point one step at a time, where the superscript i represents the i-th sampling point and n represents the number of samples. The weighted average of the Sigma point means; the corresponding autocovariance matrix. for:

[0166]

[0167] Where i is the i-th sampling point, and n is the number of samples. This is the weighted value of the Sigma point covariance. t is the observed quantity k Predict the Sigma point in real time. For t k Time-predicted observations, R k For t k The observation noise covariance matrix at time step; the cross-covariance matrix between the predicted observations and the one-step predicted state variables. for:

[0168]

[0169] Step 4: Based on the state update equation and observation update equation of the thrust calibration system in Step 3, calculate the noise covariance adaptive matrix according to the residual covariance estimate, and adaptively update the observation noise covariance and system noise covariance through the noise covariance adaptive matrix to obtain the satellite relative orbit state estimate.

[0170] Due to inaccuracies in the system statistical model estimation and various interference factors during the measurement process, excessive errors between prior information and the true values ​​of state variables can increase tracking errors and cause filter divergence. Therefore, an adaptive calculation strategy based on covariance matching is adopted to estimate the system noise covariance and observation conoise variance by maintaining consistency between the estimated residuals of system variables and their theoretical values. This includes the following steps:

[0171] S041: Calculate the residual covariance estimate; define the residual ε. k The difference between the observed actual value and the predicted value at each step:

[0172]

[0173] Among them, Z k For t k Observe the actual value at all times. For t k Predicted observations at specific times; theoretical value of residual covariance C k satisfy:

[0174]

[0175] Residual covariance estimate The calculation formula is:

[0176]

[0177] Where, ε i For i i Time residual, The mean of the residuals, This is the difference between the residual and the mean of the residuals;

[0178] S042: Based on the state update equation and observation update equation of the thrust calibration system in step three, calculate the estimated value of the adaptive update of the observation noise covariance according to the residual covariance estimate, including:

[0179] If m is the dimension of the observations, then the adaptive matrix S of the observation noise covariance is... k Represented as a diagonal matrix:

[0180]

[0181] Among them, S k (1) S represents the diagonal element of the first row of the adaptive matrix of observation noise covariance. k (2) S represents the diagonal element of the second row of the adaptive matrix of observation noise covariance. k (m) represents the diagonal element of the m-th row of the observation noise covariance adaptive matrix; the diagonal element S k The formula for calculating (i), i = 1, 2, ..., m, is:

[0182]

[0183] in, is the residual covariance estimate; μ is an adjustable parameter, and 1≤μ≤1000; This is the weighted value of the Sigma point covariance; t represents the observed quantity k Predict Sigma points in real time; For t k Time-predicted observations; N k (i,i) and R k (i,i) are matrices N and N respectively. k and R k The diagonal elements in the i-th row and i-th column; the estimated value of the observation noise covariance adaptively updated. for:

[0184]

[0185] Among them, R k For t k The observation noise covariance matrix at time step;

[0186] S043: Based on the residual covariance estimate, calculate the adaptively updated estimate of the system noise covariance, including:

[0187] The adaptive matrix Λ of the system noise covariance k Represented as a diagonal matrix:

[0188]

[0189] Among them, Λ k (1) represents the diagonal element of the first row of the adaptive noise covariance matrix, Λ k (2) represents the diagonal element of the second row of the adaptive matrix of system noise covariance, Λ k (n) represents the diagonal element of the nth row of the system noise covariance adaptive matrix; the diagonal element Λ k The formula for calculating (i), i = 1, 2, ..., n, is:

[0190]

[0191] in, and For each of the matrices and The diagonal elements in the i-th row and i-th column, where m is the observation dimension and n is the number of samples; the estimated value of the system noise covariance adaptively updated. for:

[0192]

[0193] Among them, Q k-1 For t k-1 The system noise covariance matrix at time t;

[0194] S044: Update the observation self-covariance matrix and cross-covariance matrix based on the adaptively estimated system noise covariance and observation co-noise covariance; and The calculation formula is the same as step S035;

[0195] S045: Calculate the Kalman gain, state estimation results, and their covariance matrix to obtain the estimated value of the satellite's relative orbital state; Kalman gain K k The calculation formula is:

[0196]

[0197] in, The cross-covariance matrix between the predicted observations and the one-step predicted state variables; For one step of predicting observations Autocovariance matrix; State estimation results and its covariance matrix P k The calculation formula is:

[0198]

[0199]

[0200] Among them, K k For Kalman gain, For t k The one-step predicted state quantity at time Z k For t k Observe the actual value at all times. For t k Predicted observations at any time, P k|k-1 for The estimated covariance matrix, for The autocovariance matrix;

[0201] The algorithm flow for satellite relative orbit state estimation using adaptive unscented Kalman filtering in this embodiment of the invention is as follows: Figure 2 As shown; Therefore, based on the interstellar distance and relative velocity, t can be calculated sequentially following the above steps. k The estimated value of the satellite's relative orbital state at time (k = 0, 1, ..., n);

[0202] Step 5: Based on the estimated satellite relative orbital state from Step 4, calculate the least-squares solution of the satellite orbital control thrust vector. This solution provides the real-time on-orbit calibration result of the satellite orbital control thrust, thereby achieving adaptive calibration of the satellite orbital control thrust. This includes:

[0203] t is obtained based on the adaptive unscented Kalman filtering method. k Satellite relative orbital state estimate at time (k = 0, 1, ..., n) Solve for the estimated satellite orbit control thrust. The least squares solution form:

[0204]

[0205] in, For t k The components of the satellite along the track, radial, and normal directions at time (k = 0, 1, ..., n), where V is the value at time t. k Constantly refer to the orbital velocity of the star. For the state estimation results, the coefficient matrix A k The calculation formula is:

[0206]

[0207] Among them, a k i k and u k t k By referencing the semi-major axis, orbital inclination, and latitude argument of the satellite at all times, and with Δt being the sampling time, the least squares solution of the satellite orbit control thrust vector is calculated, which is the real-time on-orbit calibration result of the satellite orbit control thrust, thereby realizing the real-time on-orbit calibration of the satellite orbit control thrust.

[0208] Thus, adaptive on-orbit real-time calibration of satellite orbit control thrust based on inter-satellite distance information and multi-source error compensation has been achieved.

[0209] The parts of this invention not described in detail are well-known to those skilled in the art.

Claims

1. A method for adaptive calibration of orbit control thrust based on inter-satellite distance and error compensation, characterized in that, include: Step 1: Establish a thrust calibration system based on the satellite's relative distance and velocity within the reference star's orbital coordinate system; Step 2: Using the inter-satellite distance and relative velocity information of the satellites as observations, establish the state equation, evolution equation, and observation equation of the unscented Kalman filter process. Step 3: Based on the state equation, evolution equation, and observation equation of the unscented Kalman filter process in Step 2, and according to the unscented changes, a symmetric sampling method is used to select a set of sample points and their weight values ​​to determine the state update equation and observation update equation of the thrust calibration system. Step 4: Based on the state update equation and observation update equation of the thrust calibration system in Step 3, calculate the noise covariance adaptive matrix according to the residual covariance estimate, and adaptively update the observation noise covariance and system noise covariance through the noise covariance adaptive matrix to obtain the satellite relative orbit state estimate. Step 5: Based on the estimated satellite relative orbital state value from Step 4, calculate the least squares solution of the satellite orbital control thrust vector, which is the real-time on-orbit calibration result of the satellite orbital control thrust, thereby realizing the adaptive calibration of the satellite orbital control thrust.

2. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 1, characterized in that, In step one, the relative distance and velocity of the satellite within the reference star orbit coordinate system include: The relative distance and velocity of a satellite within the reference star's orbital coordinate system are defined as (x, y, z, v). x ,v y ,v z ), where x is the relative distance component along the x-direction of the coordinate system, y is the relative distance component along the y-direction of the coordinate system, z is the relative distance component along the z-direction of the coordinate system, and V x V is the component of relative velocity along the x-direction of the coordinate system. y V represents the relative velocity component along the y-direction of the coordinate system. z The relative velocity component along the z-direction of the coordinate system is represented as the difference between the orbital elements of the satellite and the reference star, i.e., the relative orbital elements, which are expressed as (Δa, Δe, Δi, ΔΩ, Δω, Δu), where Δa, Δe, Δi, ΔΩ, Δω, Δu are the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative argument of latitude, respectively. In step one, the thrust calibration system includes satellite motion equations based on relative orbital elements and relative dynamic equations based on Gaussian perturbation equations; The satellite motion equations based on relative orbital elements are as follows: Where x is the relative distance component along the x-direction of the coordinate system, y is the relative distance component along the y-direction of the coordinate system, z is the relative distance component along the z-direction of the coordinate system, and V x V is the component of relative velocity along the x-direction of the coordinate system. y V represents the relative velocity component along the y-direction of the coordinate system. z The relative velocity component along the z-direction of the coordinate system is given by: a, i, u, and V, which represent the satellite's semi-major axis, orbital inclination, latitudinal argument, and orbital velocity, respectively; u0 is the latitudinal argument of the satellite's initial position; Δa, Δe, Δi, ΔΩ, Δω, and Δu represent the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative latitudinal argument, respectively; the relative eccentricity vector between the satellite and the reference star is Δe = [Δe...]. x ,Δe y ] T Replace Δe, Δω; Δe x and Δe y These are the components of the relative eccentricity Δe along the x-direction and y-direction of the coordinate system, respectively. The relative dynamic equation based on the Gaussian perturbation equation is as follows: Where t is time, Let f be the derivative of a physical quantity with time, n be the angular velocity of the satellite orbit, and f be the angular velocity of the satellite orbit. u f r and f h Ω, ω, and u represent the satellite's orbital control thrust along the track, radial, and normal directions, respectively; a, e, i, Ω, ω, and u represent the satellite's semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and argument of latitude, respectively; V is the satellite's orbital velocity; e x and e y These are the components of the satellite eccentricity along the x-direction and y-direction of the coordinate system, respectively.

3. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 2, characterized in that, In step two, the inter-satellite distance and relative velocity information of the satellites includes: Using the start time of the inter-satellite measurement as the initial time t0, obtain t k The inter-satellite distance and relative velocity at time (k = 0, 1, ..., n) are used as the observations, and the satellite relative orbital elements are used as the parameters to be estimated by the filter. Set the sampling time to Δt, and let X... k For t k The system state estimate at time X, i.e., the satellite's relative orbital elements. k =[Δa k ,Δi k ,ΔΩ k ,Δe x,k ,Δe y,k ,Δu k ] T X k+1 For t k+1 The system state estimate at time Z; k For t k System observations at any given time, namely, the inter-satellite distances and relative velocities of the satellites. f(·) and h(·) are the system state transition function and the system measurement function, respectively; ω k ,υ k t k The system noise and observation noise at time t are independent zero-mean Gaussian white noises, and both are at time t. k The noise covariance matrices at time points are Q and Q, respectively. k and R k .

4. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 3, characterized in that, In step two, the state equation, evolution equation, and observation equation of the unscented Kalman filtering process are established, including: The relative dynamic equation based on the Gaussian perturbation equation is used as the state equation for the unscented Kalman filtering process; The evolution equations for the unscented Kalman filtering process are established as follows: The observation equations for the unscented Kalman filtering process are established as follows: Among them, X k For t k The system state estimate at time X k+1 For t k+1 The system state estimate at time Z k For t k The system observation at time ω k ,υ k t k The system noise and observation noise at time points are given, f(·) and h(·) are the system state transition functions, Δa, Δe, Δi, ΔΩ, Δω, and Δu are the relative semi-major axis, relative eccentricity, relative orbital inclination, relative right ascension of the ascending node, relative argument of perigee, and relative argument of latitude, respectively; a, i, u, and V are the satellite's semi-major axis, orbital inclination, argument of latitude, and orbital velocity, respectively; u0 is the argument of latitude of the satellite's initial position; Δe x and Δe y These are the components of the relative eccentricity Δe along the x-direction and the y-direction of the coordinate system, respectively.

5. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 3, characterized in that, Based on the state equation, evolution equation, and observation equation of the unscented Kalman filter process in step two, and according to the unscented changes, a symmetric sampling method is used to select a set of sample points and their weights to determine the state update equation and observation update equation of the thrust calibration system. This includes the following steps: S031: Define t k The filtered prediction state vector of the satellite's relative orbital elements at time (k = 0, 1, ..., n) is: and its estimated variance is P k Initialize system parameters t0 and state vector X0, and set the system state estimate. and its estimated variance P0; S032: Based on the state equation, evolution equation, and observation equation of the unscented Kalman filtering process in step two, and according to the unscented variation, the Sigma point set and its weight values ​​are selected using the symmetric sampling method; t is calculated. k-1 Sigma point at time: in, Indicates t k-1 The Sigma points are symmetrically sampled at each time step, with the superscript i indicating the i-th sampling point; λ is the primary scaling factor. and P k-1 t k-1 The filtered predicted state vector and its estimated variance relative to the orbital elements at time points, where n is the number of samples. This represents the Sigma point symmetrically sampled at time t0; S033: Calculate t k The predicted state variables at time t and the corresponding prediction error covariance; define the subscript k|k-1 to indicate the state variables predicted at time t. k-1 Time prediction t k The time estimate; the transformation result of the Sigma point is calculated based on the system state transition function f(·). Establish the state update equations for the thrust calibration system and calculate the one-step predicted state variables. and the corresponding estimated covariance matrix P k|k-1 : in, To predict the state variables in one step, P k|k-1 To estimate the covariance matrix, where i is the i-th sampling point, n is the number of samples, and Q... k-1 For t k-1 The system noise covariance matrix at time t. To calculate the transformation result of the Sigma point based on the system state transition function f(·), These are the weighted values ​​of the mean and covariance at the Sigma points, respectively. Where i is the i-th sampling point, n is the number of samples, α is the primary scaling factor, β is the secondary scaling factor, and λ is the primary scaling factor; S034: Calculate t k Predicting the Sigma point in one step at time t; based on t k One-step prediction of state quantity at time t and prediction error covariance P k|k-1 Calculate the one-step prediction Sigma point: in, t represents the state quantity k The Sigma point is predicted at time step i, where the superscript i indicates the i-th sampling point; λ is the primary scaling factor, and n is the number of samples; S035: Calculate t k One-step predicted observations at time t and their autocovariance and crosscovariance matrices; calculate the transformation results at the Sigma point based on the system measurement function h(·): in, t is the observed quantity k Predict the Sigma point one step at a time, where the superscript i represents the i-th sampling point and n represents the number of samples; Establish the observation update equations for the thrust calibration system. To calculate the transformation result of the Sigma point based on the system state transition function f(·), the one-step predictive observation is calculated. One-step predictive observations The calculation formula is: in, t is the observed quantity k Predict the Sigma point one step at a time, where the superscript i represents the i-th sampling point and n represents the number of samples. The weighted average of the Sigma point means; the corresponding autocovariance matrix. for: Where i is the i-th sampling point, and n is the number of samples. This is the weighted value of the Sigma point covariance. t is the observed quantity k Predict the Sigma point in real time. For t k Time-predicted observations, R k For t k The observation noise covariance matrix at time step; the cross-covariance matrix between the predicted observations and the one-step predicted state variables. for:

6. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 1, characterized in that, Based on the state update equation and observation update equation of the thrust calibration system obtained in step three, and according to the residual covariance estimate, the noise covariance adaptive matrix is ​​calculated. The observation noise covariance and system noise covariance are then adaptively updated using the noise covariance adaptive matrix to obtain the satellite's relative orbital state estimate. This process includes the following steps: S041: Calculate the residual covariance estimate, and define the residual ε. k The difference between the observed actual value and the predicted value in one step: Among them, Z k For t k Observe the actual value at all times. For t k Predicted observations at specific times; residual covariance estimates The calculation formula is: Where, ε i For i i Time residual, The mean of the residuals, This is the difference between the residual and the mean of the residuals; S042: Based on the state update equation and observation update equation of the thrust calibration system in step three, calculate the estimated value of the adaptive update of the observation noise covariance according to the residual covariance estimate; and set the adaptive matrix S of the observation noise covariance... k Represented as a diagonal matrix, the estimated value of the observation noise covariance is adaptively updated. for: Among them, R k For t k The observation noise covariance matrix at time step; S043: Based on the residual covariance estimate, calculate the adaptive update estimate of the system noise covariance; the system noise covariance adaptive matrix Λ k It can be represented as a diagonal matrix, and the estimated value of the system noise covariance is updated adaptively. for: Among them, Q k-1 For t k-1 The system noise covariance matrix at time t; S044: Update the observation self-covariance matrix and cross-covariance matrix based on the adaptively estimated system noise covariance and observation co-noise covariance; and The calculation formula is the same as step S035; S045: Calculate the Kalman gain, state estimation results and their covariance matrix to obtain the satellite's relative orbital state estimate; The Kalman gain K k The calculation formula is: in, This is the cross-covariance matrix between the predicted observations and the one-step predicted state variables; For one step of predicting observations The autocovariance matrix; the state estimation result and its covariance matrix P k The calculation formula is: Among them, K k For Kalman gain, For t k The one-step predicted state quantity at time Z k For t k Observe the actual value at all times. For t k Predicted observations at any time, P k|k-1 for The estimated covariance matrix, for The autocovariance matrix.

7. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 6, characterized in that, In step S042, the estimated value of the observation noise covariance is calculated adaptively based on the residual covariance estimate, including: If m is the dimension of the observations, then the adaptive matrix S of the observation noise covariance is... k Represented as a diagonal matrix: Among them, S k (1) S represents the diagonal element of the first row of the adaptive matrix of observation noise covariance. k (2) S represents the diagonal element of the second row of the adaptive matrix of observation noise covariance. k (m) represents the diagonal element of the m-th row of the observation noise covariance adaptive matrix; the diagonal element S k The formula for calculating (i), i = 1, 2, ..., m, is: in, is the residual covariance estimate; μ is an adjustable parameter, and 1≤μ≤1000; This is the weighted value of the Sigma point covariance; t represents the observed quantity k Predict Sigma points in real time; For t k Time-predicted observations; N k (i,i) and R k (i,i) are matrices N and N respectively. k and R k The diagonal elements in the i-th row and i-th column; the estimated value of the adaptively updated observation noise covariance is:

8. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 6, characterized in that, In step S043, the estimated value of the system noise covariance is calculated adaptively based on the residual covariance estimate, including: The adaptive matrix Λ of the system noise covariance k Represented as a diagonal matrix: Among them, Λ k (1) represents the diagonal element of the first row of the adaptive noise covariance matrix, Λ k (2) represents the diagonal element of the second row of the adaptive matrix of system noise covariance, Λ k (n) represents the diagonal element of the nth row of the system noise covariance adaptive matrix; the diagonal element Λ k The formula for calculating (i), i = 1, 2, ..., n, is: in, and For each of the matrices and The diagonal elements in the i-th row and i-th column, where m is the observation dimension and n is the number of samples; the estimated value of the adaptively updated system noise covariance is:

9. The adaptive calibration method for orbit control thrust based on inter-satellite distance and error compensation according to claim 6, characterized in that, Based on the estimated satellite relative orbital state value from step four, the least-squares solution of the satellite orbit control thrust vector is calculated, which is the real-time on-orbit calibration result of the satellite orbit control thrust, thereby realizing the real-time on-orbit calibration of the satellite orbit control thrust, including: Solve for the estimated satellite orbit control thrust. The least squares solution form: in, For t k The components of the satellite along the track, radial, and normal directions at time (k = 0, 1, ..., n), V k It is t k Constantly refer to the orbital velocity of the star. For the state estimation results, the coefficient matrix A k The calculation formula is: Among them, a k i k and u k t k By referencing the semi-major axis, orbital inclination, and latitude argument of the satellite at all times, and with Δt representing the sampling time, the least squares solution of the satellite orbit control thrust vector is calculated, which is the real-time on-orbit calibration result of the satellite orbit control thrust, thus realizing the real-time on-orbit calibration of the satellite orbit control thrust.

Images