Satellite formation flat root state continuous estimation method and system based on Kalman filter

Abstract

The present invention discloses a method and system for continuously estimating the flat root state of a satellite formation based on Kalman filtering. The method steps include: obtaining the state equation of the relative orbital elements of the satellite formation based on the relative orbital elements of the satellite formation; obtaining the observation equation of the relative position based on the state equation; predicting the intersatellite state of the satellite formation based on the observation equation; calculating the prediction error variance and the filter gain based on the state equation and the observation equation; obtaining the intersatellite state of the satellite formation based on the intersatellite prediction state, the filter gain and the error variance of the satellite formation. The system includes: a first calculation module, a second calculation module, a first state module, a third calculation module and a second state module; the present invention establishes the intersatellite relative orbital element equation, and then uses the Kalman filter estimation technology to recursively solve the filter gain and the error variance matrix in real time, thereby obtaining the intersatellite state of the satellite formation and overcoming the jump point after the instantaneous flat conversion.

Description

Technical Field

[0001] The invention belongs to the field of satellite formations, and in particular relates to a method and system for continuously estimating the flat-root state of a satellite formation based on Kalman filtering. Background Art

[0002] Satellite formation is one of the key issues in aerospace technology research at home and abroad in recent years. The use of multiple satellites in formation flight can achieve functions that traditional single large satellites do not have. With the continuous development of micro-spacecraft technology, flight formations composed of several or even dozens of low-cost micro-spacecraft have received widespread attention. It can achieve stereoscopic detection of the same target at the same time, and can provide a large aperture and a long measurement baseline. It has extraordinary significance in the fields of communication, remote sensing, navigation, electronic reconnaissance, stereo imaging, precise positioning, atmospheric astronomy and geophysical observation.

[0003] At present, there are several methods for measuring the state of satellite formations, including laser measurement, infrared measurement, visual relative navigation, and GPS measurement. Due to the narrow laser beam, an additional guidance system is required, which increases the complexity of the system and is not conducive to the miniaturization of the measurement system; the infrared measurement method has high accuracy, but can only provide angle measurement information, not distance measurement information, and cannot be used alone; the visual relative navigation method has high measurement accuracy, but its working distance is short, which is not suitable for medium and long-range intersatellite relative navigation; the satellite-based GPS measurement system has high accuracy, but a complex processing process is required to obtain high real-time determination accuracy. Summary of the invention

[0004] The present invention proposes a method and system for continuous estimation of the flat root state of a satellite formation based on Kalman filtering. By establishing the inter-satellite relative orbit element equation and then using the Kalman filter estimation technology, the filter gain and error variance matrix are recursively solved in real time to obtain the inter-satellite state of the satellite formation.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A method for continuous estimation of the flat-root state of a satellite formation based on Kalman filtering, the steps comprising:

[0007] S1. Based on the relative orbit elements of the satellite formation, a state equation of the relative orbit elements of the satellite formation is obtained;

[0008] S2. Based on the state equation, an observation equation of relative position is obtained;

[0009] S3. Predicting the intersatellite state of the satellite formation based on the observation equation;

[0010] S4. Calculate the prediction error variance and the filter gain based on the state equation and the observation equation;

[0011] S5. Obtain an estimate of the inter-satellite state of the satellite formation based on the inter-satellite prediction state of the satellite formation, the filter gain and the error variance.

[0012] Preferably, the method for obtaining the state method in S1 includes:

[0013] S11. Based on the perturbation theory, the expression of the relative orbital elements of the satellite formation changing with time under the first-order long-term perturbation is obtained;

[0014] S12. Based on the expression of the relative orbital elements of the satellite formation changing with time under the first-order long-term perturbation, a state equation of the relative orbital elements of the satellite formation is obtained.

[0015] Preferably, the relative orbital element number of the satellite formation in S11 is defined as δα=[δa,δu,δe x ,δe y ,δi x ,δi y ] T According to the perturbation theory, we can get the relative orbital element number of the satellite formation δα=[δa,δu,δe x ,δe y ,δi x ,δi y ] T In J 2 The expression of the time variation under the first-order long-term perturbation is:

[0016]

[0017] Where a is the semi-major axis, u is the satellite latitude argument, e is the eccentricity, ω is the perigee argument, and e is the x =ecosω,e y =e sinω, i represents the orbital inclination, i x =i cosω,i y =i sinω,u 0 represents the initial value of latitude argument, γ represents a constant that depends on the J2 perturbation, represents the relative inclination vector phase angle, represents the differential of the phase angle of the E vector with respect to the latitude argument u of the formation master star,

[0018] Preferably, the method for obtaining the state equation of the satellite formation relative orbital elements in S12 includes:

[0019] S121, calculating the state transfer matrix;

[0020]

[0021] In the formula, k represents the time, Φ k,k-1 Indicates t k-1 Time to t k Time state transfer matrix; u k Indicates t k The latitude argument at the time, u k-1 Indicates t k-1 The latitude argument at the time;

[0022] S122, calculating the state equation of the relative orbital elements of the satellite formation;

[0023] δα k =Φ k,k-1 δα k-1 +W k-1

[0024] In the formula, δα k is the state vector of the system, W k-1 is the process noise vector.

[0025] Preferably, the method for calculating the state equation in S2 includes:

[0026] S21. Based on the kinematic equations described by the orbital elements, the position formula of the slave satellite relative to the master satellite is obtained:

[0027]

[0028] In the formula, r represents the orbital vector diameter of the slave star relative to the master star, r r , r t , r n Respectively represent the radial, along-track and normal components;

[0029] S22, based on the position formula of the slave satellite relative to the master satellite, obtain the observation matrix H k ;

[0030]

[0031] In the formula, H k Indicates t k The observation matrix at time;

[0032] S23, based on the observation matrix, obtaining an observation equation of the relative position δr;

[0033] δr k =H k δα k +V k

[0034] Where V k represents the observation noise vector; δα k is the state vector of the system.

[0035] Preferably, the prediction error variance and filter gain calculation method in S4 includes:

[0036] The prediction error variance calculation formula is as follows:

[0037]

[0038] Where P k,k-1 represents the prediction error variance matrix; P k-1 Indicates t k-1 The forecast error variance matrix at time Φ k,k-1 Indicates t k-1 Time to t k The n×n dimensional state transfer matrix at the time; Γ k,k-1 represents the n×p dimensional noise input matrix, Q k-1 Indicates t k-1 The process noise W at time k The p×p dimensional symmetric non-negative definite variance matrix of ;

[0039] The filter gain calculation formula is as follows:

[0040]

[0041] In the formula, K k is the gain matrix, represents the m×n dimensional observation matrix, R k represents the observation noise V k The m-by-m symmetric positive definite variance matrix of .

[0042] The present invention also provides a satellite formation flat root state continuous estimation system based on Kalman filtering, comprising: a first calculation module, a second calculation module, a first state module, a third calculation module and a second state module;

[0043] The first calculation module is used to obtain a state equation of the relative orbit elements of the satellite formation based on the relative orbit elements of the satellite formation;

[0044] The second calculation module is used to obtain an observation equation of the relative position based on the state equation;

[0045] The first state module is used to predict the inter-satellite state of the satellite formation based on the observation equation;

[0046] The third calculation module is used to calculate the prediction error variance and the filter gain based on the state equation and the observation equation;

[0047] The second state module is used to obtain an estimate of the inter-satellite state of the satellite formation based on the inter-satellite prediction state of the satellite formation, the filter gain and the error variance.

[0048] Preferably, the third calculation module includes an error variance calculation unit and a filter gain calculation unit;

[0049] The error variance calculation unit is used to calculate the prediction error variance;

[0050] The filter gain calculation unit is used to calculate the filter gain.

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

[0052] (1) By establishing a satellite formation flat-root state perceptual divergence model, the satellite formation state and observation equations are given, which is more conducive to the continuous estimation of the flat-root state and the solution of the relative state between satellite formations;

[0053] (2) By adopting the method of the present invention, the estimated satellite formation tracking state is continuous and real-time, and the estimated tracking number is accurate, highly real-time, short in delay, and continuous. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0055] Figure 1 1 is a flow chart of a method for continuous estimation of the flat-root state of a satellite formation based on Kalman filtering according to Embodiment 1 of the present invention;

[0056] Figure 2 It is a structural diagram of a satellite formation flat root state continuous estimation system based on Kalman filtering according to Embodiment 2 of the present invention. DETAILED DESCRIPTION

[0057] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

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

[0059] Embodiment 1

[0060] like Figure 1 FIG. 1 is a flow chart of a method for continuous estimation of the flat-root state of a satellite formation based on Kalman filtering in accordance with the first embodiment of the present invention. The method steps include:

[0061] S1. Based on the relative orbit elements of the satellite formation, a state equation of the relative orbit elements of the satellite formation is obtained;

[0062] Among them, the methods for obtaining the status method include:

[0063] S11. Based on the perturbation theory, the expression of the relative orbital elements of the satellite formation changing with time under the first-order long-term perturbation is obtained;

[0064] The relative orbital element of the satellite formation is defined as δα = [δa, δu, δe x ,δe y ,δi x ,δi y ] T According to the perturbation theory, we can get the relative orbital elements of the satellite formation δα=[δa,δu,δe x ,δe y ,δi x ,δi y ] T In J 2 The expression of the time variation under the first-order long-term perturbation is:

[0065]

[0066] Where a is the semi-major axis, u is the satellite latitude argument, e is the eccentricity, ω is the perigee argument, and e is the x =ecosω,e y =e sinω, i represents the orbital inclination, i x =i cosω,i y =i sinω,u 0 represents the initial value of latitude argument, γ is a constant that depends on the J2 perturbation, represents the relative inclination vector phase angle, represents the differential of the phase angle of the E vector with respect to the latitude argument u of the formation master star,

[0067] S12. Based on the expression of the relative orbit elements of the satellite formation changing with time under the first-order long-term perturbation, the state equation of the relative orbit elements of the satellite formation is obtained;

[0068] S121, calculating the state transfer matrix;

[0069]

[0070] In the formula, k represents the time, Φ k,k-1 Indicates t k-1 Time to t k Time state transfer matrix; u k Indicates t k The latitude argument at the time, u k-1 Indicates t k-1 The latitude argument at the time;

[0071] S122, calculating the state equation of the relative orbital elements of the satellite formation;

[0072] δα k =Φ k,k-1 δα k-1 +W k-1

[0073] In the formula, δα k is the state vector of the system, W k-1 is the process noise vector.

[0074] S2. Based on the state equation, the observation equation of relative position is obtained;

[0075] S21. Based on the kinematic equations described by the orbital elements, the position formula of the slave satellite relative to the master satellite is obtained:

[0076]

[0077] Where r is the orbital vector diameter of the slave star relative to the master star, r r r t r n Represent the radial, along-track and normal components respectively;

[0078] S22. Based on the position formula of the slave satellite relative to the master satellite, the observation matrix H is obtained k ;

[0079]

[0080] In the formula, H k Indicates t k The observation matrix at time;

[0081] S23, based on the observation matrix, obtaining the observation equation of the relative position δr;

[0082] δr k =H k δα k +V k

[0083] Where V k represents the observation noise vector; δα k is the state vector of the system.

[0084] Among them, the process noise vector W k-1 and the observation noise vector V k Kalman filtering requires that all noise interference is Gaussian white noise, and process noise and observation noise are independent of each other. Therefore, for the statistical characteristics of noise, the following assumptions are made:

[0085]

[0086] In the formula, Q k is the process noise W k The p×p dimensional symmetric non-negative definite variance matrix, R k is the observation noise V k The m×m symmetric positive definite variance matrix, δ kj is the Kronrcker-δ function, that is

[0087] S3. Based on the observation equation, the intersatellite state of the satellite formation is predicted;

[0088] S4. Calculate the prediction error variance and filter gain based on the state equation and the observation equation;

[0089] The forecast error variance is calculated as follows:

[0090]

[0091] Where P k,k-1 represents the prediction error variance matrix; P k-1 Indicates t k-1 The forecast error variance matrix at time Φ k,k-1 Indicates t k-1 Time to t k The n×n dimensional state transfer matrix at the time, Γ k,k-1 represents the n×p dimensional noise input matrix, Q k-1 Indicates t k-1 The process noise W at time k The p×p dimensional symmetric non-negative definite variance matrix of ;

[0092] The filter gain calculation formula is as follows:

[0093]

[0094] In the formula, K k is the gain matrix, represents the m×n dimensional observation matrix, R k represents the observation noise V k The m-by-m symmetric positive definite variance matrix of .

[0095] S5, based on the satellite formation inter-satellite prediction state, filter gain and error variance, using the filter equation

[0096]

[0097] Get an estimate of the intersatellite state of the satellite formation

[0098] Embodiment 2

[0099] The present invention also proposes a satellite formation flat root state continuous estimation system based on Kalman filtering, comprising: a first calculation module, a second calculation module, a first state module, a third calculation module and a second state module;

[0100] The first calculation module is used to obtain the state equation of the relative orbit elements of the satellite formation based on the relative orbit elements of the satellite formation:

[0101] Its working process includes:

[0102] Define the relative orbital elements of the satellite formation as δα=[δa,δu,δe x ,δe y ,δi x ,δi y ] T According to the perturbation theory, the relative orbital elements of the satellite formation are obtained in J 2 The expression of the time-varying first-order long-term perturbation is:

[0103]

[0104] Where a is the semi-major axis, u is the satellite latitude argument, e is the eccentricity, ω is the perigee argument, and e is the x =ecosω,e y =e sinω, i represents the orbital inclination, i x =i cosω,i y =i sinω,u 0 represents the initial value of latitude argument, γ is a constant that depends on the J2 perturbation, represents the relative inclination vector phase angle, represents the differential of the phase angle of the E vector with respect to the latitude argument u of the formation master star,

[0105] Calculate the state transfer matrix of the relative orbital elements of the satellite formation:

[0106]

[0107] In the formula, k represents the time, Φk ,k-1 Indicates t k-1 Time to t kTime state transfer matrix; u k Indicates t k The latitude argument at the time, u k-1 Indicates t k-1 The latitude argument at the time;

[0108] Then calculate the state equation of the relative orbital elements of the satellite formation:

[0109] δα k =Φ k,k-1 δα k-1 +W k-1

[0110] In the formula, δα k is the state vector of the system, W k-1 is the process noise vector.

[0111] The second calculation module is used to obtain the observation equation of the relative position based on the state equation;

[0112] Its working process includes:

[0113] First, based on the kinematic equations described by the orbital elements, the position formula of the slave star relative to the master star is obtained:

[0114]

[0115] Where r is the orbital vector diameter of the slave star relative to the master star, r r r t r n Represent the radial, along-track and normal components respectively;

[0116] Secondly, based on the position formula of the slave satellite relative to the master satellite, the observation matrix H is obtained k ;

[0117]

[0118] In the formula, H k Indicates t k The observation matrix at time;

[0119] Finally, based on the observation matrix, the observation equation of the relative position δr is obtained;

[0120] δr k =H k δα k +V k

[0121] Where V k represents the observation noise vector; δα k is the state vector of the system.

[0122] The first state module is used to predict the intersatellite state of the satellite formation based on the observation equation;

[0123] The third calculation module is used to calculate the prediction error variance and the filter gain based on the state equation and the observation equation;

[0124] Its working process includes:

[0125] Calculate the forecast error variance using the following formula:

[0126]

[0127] Where P k,k-1 represents the prediction error variance matrix; P k-1 Indicates t k-1 The forecast error variance matrix at time Φ k,k-1 Indicates t k-1 Time to t k The n×n dimensional state transfer matrix at the moment, Γ k,k-1 represents the n×p dimensional noise input matrix, Q k-1 Indicates t k-1 The process noise W at time k The p×p dimensional symmetric non-negative definite variance matrix of ;

[0128] Calculate the filter gain using the following formula:

[0129]

[0130] In the formula, K k is the gain matrix, represents the m×n dimensional observation matrix, R k represents the observation noise V k The m-by-m symmetric positive definite variance matrix of .

[0131] The second state module is used to obtain an estimate of the inter-satellite state of the satellite formation based on the inter-satellite prediction state of the satellite formation, the filter gain and the error variance;

[0132] Its working process includes:

[0133] Calculate the filter equation, the formula is as follows:

[0134]

[0135] Based on the above filtering equation, the estimation of the intersatellite state of the satellite formation is obtained

[0136] Wherein, the third calculation module includes an error variance calculation unit and a filter gain calculation unit;

[0137] The error variance calculation unit is used to calculate the prediction error variance;

[0138] The filter gain calculation unit is used to calculate the filter gain.

[0139] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A method for continuous estimation of the flat-root state of a satellite formation based on Kalman filtering, characterized in that the steps include: S1. Based on the relative orbit elements of the satellite formation, a state equation of the relative orbit elements of the satellite formation is obtained; S2. Based on the state equation, an observation equation of relative position is obtained; S3. Predicting the intersatellite state of the satellite formation based on the observation equation; S4. Calculate the prediction error variance and the filter gain based on the state equation and the observation equation; S5. Obtaining an estimate of the inter-satellite state of the satellite formation based on the inter-satellite prediction state of the satellite formation, the filter gain and the error variance; The method for obtaining the state equation in S1 includes: S11. Based on the perturbation theory, the expression of the relative orbital elements of the satellite formation changing with time under the first-order long-term perturbation is obtained; S12. Based on the expression of the relative orbital elements of the satellite formation changing with time under the first-order long-term perturbation, a state equation of the relative orbital elements of the satellite formation is obtained; The relative orbital elements of the satellite formation in S11 are defined as δα=[δa, δu, δe x ,δe y ,δi x ,δi y ] T According to the perturbation theory, we can get the relative orbital element number of the satellite formation δα=[δa,δu,δe x ,δe y ,δi x ,δi y ] T The expression of the time variation under the first-order long-term perturbation of J2 is: Where a is the semi-major axis, u is the satellite latitude argument, e is the eccentricity, ω is the perigee argument, and e is the x =ecosω,e y =esinω, i represents the orbital inclination, i x =icosω,i y =isinω, u0 represents the initial value of latitude argument, γ represents a constant that depends on the J2 perturbation, represents the relative inclination vector phase angle, represents the differential of the phase angle of the E vector with respect to the latitude argument u of the formation master star, The method for obtaining the state equation of the satellite formation relative orbital elements in S12 includes: S121, calculating the state transfer matrix; In the formula, k represents the time, Φ k,k-1 Indicates t k-1 Time to t k Time state transfer matrix; u k Indicates t k The latitude argument at the time, u k-1 Indicates t k-1 The latitude argument at the time; S122, calculating the state equation of the relative orbital elements of the satellite formation; yes k =Φ k,k-1 yes k-1 +W k-1 In the formula, δα k is the state vector of the system, W k-1 is the process noise vector; The calculation method of the observation equation in S2 includes: S21. Based on the kinematic equations described by the orbital elements, the position formula of the slave satellite relative to the master satellite is obtained: In the formula, r represents the orbital vector diameter of the slave star relative to the master star, r r , r t , r n Respectively represent the radial, along-track and normal components; S22. Based on the position formula of the slave star relative to the master star, obtain the observation matrix H k ; In the formula, H k Indicates t k The observation matrix at time; S23, based on the observation matrix, obtaining an observation equation of the relative position δr; δr k =H k yes k +V k Where V k represents the observation noise vector; δα k is the state vector of the system; The prediction error variance and filter gain calculation method in S4 includes: The prediction error variance calculation formula is as follows: Where P k,k-1 represents the prediction error variance matrix; P k-1 Indicates t k-1 The forecast error variance matrix at time Φ k,k-1 Indicates t k-1 Time to t k The n×n dimensional state transfer matrix at the time; Γ k,k-1 represents the n×p dimensional noise input matrix, Q k-1 Indicates t k-1 The process noise W at time k The p×p dimensional symmetric non-negative definite variance matrix of ; The filter gain calculation formula is as follows: In the formula, K k is the gain matrix, represents the m×n dimensional observation matrix, R k represents the observation noise V k The m-by-m symmetric positive definite variance matrix of .

2. A satellite formation flat root state continuous estimation system based on Kalman filtering, the system is used to implement the method of claim 1, characterized in that: include: A first calculation module, a second calculation module, a first state module, a third calculation module and a second state module; The first calculation module is used to obtain a state equation of the relative orbit elements of the satellite formation based on the relative orbit elements of the satellite formation; The second calculation module is used to obtain an observation equation of the relative position based on the state equation; The first state module is used to predict the inter-satellite state of the satellite formation based on the observation equation; The third calculation module is used to calculate the prediction error variance and the filter gain based on the state equation and the observation equation; The second state module is used to obtain an estimate of the inter-satellite state of the satellite formation based on the inter-satellite prediction state of the satellite formation, the filter gain and the error variance.

3. The satellite formation flat root state continuous estimation system based on Kalman filtering according to claim 2, characterized in that: The third calculation module includes an error variance calculation unit and a filter gain calculation unit; The error variance calculation unit is used to calculate the prediction error variance; The filter gain calculation unit is used to calculate the filter gain.

Images