Optical formation satellite orbit determination method and system

Abstract

The invention provides an optical formation satellite orbit determination method and system, and the method comprises the steps: S1, building a kinetic model of a formation satellite system according to the stress condition of a satellite; the formation satellite system is used for formation satellites; s2, acquiring an inter-satellite distance and an inter-satellite angle based on the kinetic model; s3, establishing a state equation and an observation equation according to the inter-satellite distance and the inter-satellite angle; s4, carrying out nonlinear transformation on the random variable to obtain a nonlinear transformation result; carrying out statistical averaging on the nonlinear transformation result, and then filtering to obtain an unbiased result; and S5, constructing an iterative forecasting equation based on the unbiased result, the state equation and the observation equation, and further determining the orbit of the formation satellite. According to the method, distance measurement and angle measurement information is utilized, the orbit determination precision is improved, the principle is simple, programming calculation is convenient, the method is suitable for the field of deep space formation satellite orbit determination, and compared with traditional VLBI measurement, the orbit determination precision is improved.

Description

Technical Field

[0001] The present invention belongs to the field of aerospace optical technology, and specifically, relates to an optical formation satellite orbit determination method and system, and more specifically, relates to a high-precision orbit determination method for space gravitational wave detection formation satellites. Background Art

[0003] Regarding the problem of formation satellite orbit determination, previous researchers have conducted the following research.

[0004] Reference 1. Shao Kai [Shao Kai, Zhang Houzhe, Qin Xianping, Huang Zhiyong, Yi Bin, Gu Defeng "Precise Absolute and Relative Orbit Determination of Distributed InSAR Formation Satellites", Acta Geodaetica et Cartographica Sinica, 2021, 50(5)] used the measured data of GPS onboard the formation satellite and improved the absolute and relative orbit determination accuracy by modeling the orbital maneuvers with constant acceleration.

[0005] Reference 2. Song Yezhi [Song Yezhi, "DRO orbit and its navigation in Earth-Moon space flight with LEO formation", the 12th National Multibody Dynamics and Control and the 7th National Space Dynamics and Control and the 15th National Joint Academic Conference on Analytical Mechanics, 2020] developed the LEO-DRO formation flight orbit determination software, which improved the orbit determination accuracy of DRO orbit and its formation mission with LEO.

[0006] Reference 3. Gu Defeng [Gu Defeng, An Zicong, Zhao Yue, Yang Chengyun, Li Kangkang, Zhu Jubo, Liu Yuan, Song Jianing, "Tianqin Satellite Orbit Determination Technology", Journal of Sun Yat-sen University, 2021, 60(1-2)] studied the problem of high-precision satellite orbit determination of the Tianqin satellite, and analyzed and prospected the technology from the perspectives of modeling, compensation, micro-perturbation and multi-measurement fusion.

[0007] Reference 4. Li Kangkang [Li Kangkang, Gu Defeng, An Zicong, Song Jianing, Zhu Jubo, "Orbit Determination and Simulation Analysis of Space Gravitational Wave Detectors", Journal of Sun Yat-sen University, 2021, 8(2)] proposed a method to improve the orbit determination accuracy of the LISA detector using nonlinear weighted least squares and Monte Carlo methods.

[0008] Reference 5, Zhuo Li [Zhuo Li, Jianhua Zheng. Orbit determination for a space-based gravitational wave observatory. 2021, 185: 170-178] proposed an orbit determination method based on extended Kalman filtering for the Taiji space gravitational wave detection mission using inter-satellite ranging and angle measurement.

[0009] Among them, References 1 and 2 provide orbit determination methods for single-star problems, but they cannot be used directly for formation satellites with more than two stars; References 3, 4, and 5 solve the orbit determination problem of space gravitational wave detection formation missions; among them, Reference 3 mainly proposes a solution from the perspective of modeling and multi-measurement fusion, and Reference 4 uses the least squares and Monte Carlo method. References 3 and 4 study the orbit determination problem of geocentric orbit space gravitational wave detection satellites and use GNSS. However, for the orbit determination of heliocentric orbit space gravitational wave detection satellites, GNSS cannot play a role because they are far away from the earth. Reference 5 is based on the extended Kalman filter, and its limitation is that it works better for linear systems. Summary of the invention

[0010] In view of the defects in the prior art, an object of the present invention is to provide a method and system for determining the orbits of optical formation satellites.

[0011] A method for determining an orbit of an optical formation satellite provided by the present invention comprises:

[0012] Step S1: establishing a dynamic model of a formation satellite system according to the force conditions of the satellites; the formation satellite system is used for formation satellites;

[0013] Step S2: Based on the dynamic model, obtaining the inter-satellite distance and inter-satellite angle;

[0014] Step S3: establishing a state equation and an observation equation according to the inter-satellite distance and the inter-satellite angle;

[0015] Step S4: nonlinearly transform the random variable to obtain a nonlinear transformation result; statistically average the nonlinear transformation result, and then filter it to obtain an unbiased result;

[0016] Step S5: Based on the unbiased result and the state equation and observation equation, an iterative prediction equation is constructed to determine the orbits of the formation satellites.

[0017] Preferably, in step S1, the dynamic model uses the heliocentric J2000 coordinate system as the reference coordinate system, and the mathematical expression is:

[0018]

[0019] in, is the satellite velocity vector; v k represents the velocity vector of the formation satellite, k = 1, 2, 3;

[0020]

[0021] in, is the satellite’s acceleration vector; μ 0represents the solar gravitational constant; r k represents the position vector of the formation satellite, k = 1, 2, 3; r k is the satellite position vector; R i represents the position vectors of the eight planets, namely Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and the Moon, i = 1, 2, ..., 8, 9; μ i Represent the gravitational constants of the eight planets and the moon, i=1,2,...,8,9; R i A scalar representing the positions of the eight planets.

[0022] Preferably, in step S2, the inter-satellite angle includes an inter-satellite azimuth angle and an inter-satellite altitude angle;

[0023] The mathematical expression of the inter-satellite azimuth is:

[0024]

[0025] Where A represents the inter-satellite azimuth; ρ x , y , z is the position component of the slave star in the master star's own system, determined based on the dynamic model;

[0026] The mathematical expression of the intersatellite elevation angle is:

[0027]

[0028] Where E represents the inter-satellite elevation angle; ρ represents the inter-satellite distance;

[0029] The mathematical expression of the inter-satellite distance is:

[0030]

[0031] in, represents partial derivative; ρ represents the intersatellite distance vector; r represents the position vector; B I represents the conversion matrix from J2000 to the local system; T represents the transposition symbol;

[0032] The conversion matrix from J2000 to the local system is expressed as:

[0033] B I =M z [u c ]⊙M x [i c ]⊙M z [Ω c ] (7)

[0034] Among them, uc 、i c With Ω c Respectively represent the latitude argument, orbital inclination and right ascension of the ascending node of the reference star; M z Indicates rotation around the Z axis; M x represents rotation around the X-axis, and ⊙ represents the rotation order.

[0035] Preferably, in step S3, the state equation is the system state equation; the measurement equation is the system measurement equation;

[0036] The system state equation, mathematical expression is:

[0037] X k+1 =f(X k ,w k ,k) (1)

[0038] Among them, X k+1 is the state variable at time k+1; f is the dynamic equation; k is the state variable at time k; Z k is the position velocity of the three formation satellites, which is an 18-dimensional vector; k is the model error of the system;

[0039] The system measurement equation, mathematical expression is:

[0040] Y k+1 =h(X k+1 ,e k+1 ,k+1) (2)

[0041] Among them, Y k+1 The position, velocity, inter-satellite distance and angle measurement information of the three formation satellites; k+1 is the measurement error, w k ~N(0,Q k ),e k+1 ~N(0,R k+1 );Q k is the error variance, R k+1 is the measurement error variance; h is the observation equation;

[0042] The formation satellites include: satellite 1, satellite 2 and satellite 3;

[0043] According to the formation satellites and the observation matrix, the mathematical expression of the observation equation is:

[0044] h=(X 1 ,Y 1 ,Z 1 ,V x1 ,V y1 ,V z1 ,X2 ,Y 2 ,Z 2 ,V x2 ,V y2 ,V z2 ,X 3 ,Y 3 ,Z 3 ,V x3 ,V y3 ,V z3 ,

[0045] ρ 12 ,A 12 ,E 12 ,ρ 23 ,A 23 ,E 23 ,ρ 31 ,A 31 ,E 31 ) T (12)

[0046] Among them, ρ 12 is the distance between satellite 1 and satellite 2; ρ 23 is the distance between satellite 2 and satellite 3; ρ 31 is the distance between satellite 3 and satellite 1; A 12 is the azimuth between satellite 1 and satellite 2; A 23 is the azimuth between satellite 2 and satellite 3; A 31 is the azimuth between satellite 3 and satellite 1; E 12 is the elevation angle between satellite 1 and satellite 2; E 23 is the elevation angle between satellite 2 and satellite 3; E 31 is the elevation angle between satellite 3 and satellite 1; X 1 , Y 1 , Z 1 is the position vector of satellite 1; X 2 , Y 2 , Z 2 is the position vector of satellite 2; X 3 , Y 3 , Z 3 is the position vector of satellite 3; V x1 、V y1 、V z1 is the velocity vector of satellite 1; V x2 、V y2 、V z2 is the velocity vector of satellite 2; V x3 、V y3 、V z3 is the velocity vector of satellite 3;

[0047] The mathematical expression of the observation matrix is:

[0048]

[0049] Among them, H b represents the observation matrix; x, y, z represent the satellite position scalar; V x 、V y 、V z A scalar representing the satellite's velocity.

[0050] In step S4, the nonlinear transformation random variable is expressed as follows:

[0051] y=f(x)

[0052] Among them, y represents the result of nonlinear transformation; x represents random variable; f represents nonlinear transformation;

[0053] In step S4, the nonlinear transformation result is statistically averaged by Unscented transformation, that is, 2L+1 extension vectors are constructed according to the random variable, and the expression is:

[0054]

[0055] Among them, X i represents the extended vector, i=0,1,…2L; represents the mean of the random variable; P x represents the variance of the random variable; L represents the dimension; k represents the scalar;

[0056] An Unscented transformation point set is formed according to the 2L+1 extension vectors; and an unbiased result is obtained after Kalman filtering is performed on the Unscented transformation point set.

[0057] Preferably, in the step S5, an iterative prediction equation is constructed based on the unbiased result, the system state equation and the system measurement equation; the iterative prediction equation includes: an extended state prediction equation, a system state prediction equation, an error covariance prediction equation, an extended state prediction equation and a measurement prediction equation;

[0058] The extended state prediction equation has the following mathematical expression:

[0059] χ k+1 / k,i =f(χ k+1,i ,k),i=0,1,2,…2L

[0060] Among them, χ k+1 / k,i is the result of the extended state forecast; k+1,i Indicates the current state;

[0061] The system state prediction equation has the mathematical expression:

[0062]

[0063] in, is the result of system state prediction; W i represents weight;

[0064] The error covariance prediction equation has the mathematical expression:

[0065]

[0066] Among them, P x,k+1 / k is the result of error covariance prediction; Q k is the model error variance, Predict the results of system status;

[0067] The extended state prediction equation has the following mathematical expression:

[0068]

[0069] in, Forecast results for extended states;

[0070] The mathematical expression of the measurement prediction equation is:

[0071]

[0072] Among them, Y k+1 / k To measure the results of the prediction equation.

[0073] An optical formation satellite orbit determination system provided by the present invention comprises:

[0074] Module M1: Establishing a dynamic model of a formation satellite system according to the force conditions of the satellites; the formation satellite system is used for formation satellites;

[0075] Module M2: based on the dynamic model, obtaining the inter-satellite distance and inter-satellite angle;

[0076] Module M3: establishing a state equation and an observation equation according to the inter-satellite distance and the inter-satellite angle;

[0077] Module M4: nonlinearly transforming random variables to obtain nonlinear transformation results; statistically averaging the nonlinear transformation results, and then filtering to obtain unbiased results;

[0078] Module M5: Based on the unbiased results and the state equation and observation equation, an iterative prediction equation is constructed to determine the orbits of the formation satellites.

[0079] Preferably, in the module M1, the dynamic model uses the heliocentric J2000 coordinate system as the reference coordinate system, and the mathematical expression is:

[0080]

[0081] in, is the satellite velocity vector; v k represents the velocity vector of the formation satellite, k = 1, 2, 3;

[0082]

[0083] in, is the satellite’s acceleration vector; μ 0 represents the solar gravitational constant; r k represents the position vector of the formation satellite, k = 1, 2, 3; r k is the satellite position vector; R i represents the position vectors of the eight planets, namely Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and the Moon, i = 1, 2, ..., 8, 9; μ i Respectively represent the gravitational constants of the eight planets and the moon, i = 1, 2, ..., 8, 9. i A scalar representing the positions of the eight planets.

[0084] Preferably, in the module M2, the inter-satellite angle includes the inter-satellite azimuth angle and the inter-satellite altitude angle;

[0085] The mathematical expression of the inter-satellite azimuth is:

[0086]

[0087] Where A represents the inter-satellite azimuth; ρ x , y , z is the position component of the slave star in the master star's own system, determined based on the dynamic model;

[0088] The mathematical expression of the intersatellite elevation angle is:

[0089]

[0090] Where E represents the inter-satellite elevation angle; ρ represents the inter-satellite distance;

[0091] The mathematical expression of the inter-satellite distance is:

[0092]

[0093] in, represents partial derivative; ρ represents the intersatellite distance vector; r represents the position vector; B I represents the conversion matrix from J2000 to the local system; T represents the transposition symbol;

[0094] The conversion matrix from J2000 to the local system is expressed as:

[0095] B I =M z [u c ]⊙M x [i c ]⊙M z [Ω c ] (7)

[0096] Among them, u c 、i c With Ω c Respectively represent the latitude argument, orbital inclination and right ascension of the ascending node of the reference star; M z Indicates rotation around the Z axis; M x represents rotation around the X-axis, and ⊙ represents the rotation order.

[0097] Preferably, in the module M3, the state equation is the system state equation; the measurement equation is the system measurement equation;

[0098] The system state equation, mathematical expression is:

[0099] X k+1 =f(X k ,w k ,k) (1)

[0100] Among them, X k+1 is the state variable at time k+1; f is the dynamic equation; k is the state variable at time k; X k is the position velocity of the three formation satellites, which is an 18-dimensional vector; k is the model error of the system;

[0101] The system measurement equation, mathematical expression is:

[0102] Y k+1 =h(X k+1 ,e k+1 ,k+1) (2)

[0103] Among them, Y k+1 The position, velocity, inter-satellite distance and angle measurement information of the three formation satellites; k+1 is the measurement error, w k ~N(0,Q k ),e k+1 ~N(0,Rk+1 );Q k is the error variance, R k+1 is the measurement error variance; h is the observation equation;

[0104] The formation satellites include: satellite 1, satellite 2 and satellite 3;

[0105] According to the formation satellites and the observation matrix, the mathematical expression of the observation equation is:

[0106] h=(X 1 ,Y 1 ,Z 1 ,V x1 ,V y1 ,V z1 ,X 2 ,Y 2 ,Z 2 ,V x2 ,V y2 ,V z2 ,X 3 ,Y 3 ,Z 3 ,V x3 ,V y3 ,V z3 ,

[0107] ρ 12 ,A 12 ,E 12 ,ρ 23 ,A 23 ,E 23 ,ρ 31 ,A 31 ,E 31 ) T (12)

[0108] Among them, ρ 12 is the distance between satellite 1 and satellite 2; ρ 23 is the distance between satellite 2 and satellite 3; ρ 31 is the distance between satellite 3 and satellite 1; A 12 is the azimuth between satellite 1 and satellite 2; A 23 is the azimuth between satellite 2 and satellite 3; A 31 is the azimuth between satellite 3 and satellite 1; E 12 is the elevation angle between satellite 1 and satellite 2; E 23 is the elevation angle between satellite 2 and satellite 3; E 31 is the elevation angle between satellite 3 and satellite 1; X 1 , Y 1 , Z 1 is the position vector of satellite 1; X 2 , Y2 , Z 2 is the position vector of satellite 2; X 3 , Y 3 , Z 3 is the position vector of satellite 3; V x1 、V y1 、V z1 is the velocity vector of satellite 1; V x2 、V y2 、V z2 is the velocity vector of satellite 2; V x3 、V y3 、V z3 is the velocity vector of satellite 3;

[0109] The mathematical expression of the observation matrix is:

[0110]

[0111] Among them, H b represents the observation matrix; x, y, z represent the satellite position scalar; V c 、V y 、V z represents the satellite's velocity scalar;

[0112] In the module M4, the nonlinear transformation random variable is expressed as follows:

[0113] y=f(x)

[0114] Among them, y represents the result of nonlinear transformation; x represents random variable; f represents nonlinear transformation;

[0115] In the module M4, the nonlinear transformation results are statistically averaged by Unscented transformation, that is, 2L+1 extension vectors are constructed according to the random variables, and the expression is:

[0116]

[0117] Among them, X i represents the extended vector, i=0,1,…2L; represents the mean of the random variable; P x represents the variance of the random variable; L represents the dimension; k represents the scalar;

[0118] An Unscented transformation point set is formed according to the 2L+1 extension vectors; and an unbiased result is obtained after Kalman filtering is performed on the Unscented transformation point set.

[0119] Preferably, in the module M5, an iterative prediction equation is constructed based on the unbiased result, the system state equation and the system measurement equation; the iterative prediction equation includes: an extended state prediction equation, a system state prediction equation, an error covariance prediction equation, an extended state prediction equation and a measurement prediction equation;

[0120] The extended state prediction equation has the following mathematical expression:

[0121] χ k+1 / k,i =f(χ k+1,i ,k),i=0,1,2,…2L

[0122] Among them, χ k+1 / k,i is the result of the extended state forecast; k+1,i Indicates the current state;

[0123] The system state prediction equation has the mathematical expression:

[0124]

[0125] in, is the result of system state prediction; W i represents weight;

[0126] The error covariance prediction equation has the mathematical expression:

[0127]

[0128] Among them, P x,k+1 / k is the result of error covariance prediction; Q k is the model error variance, Predict the results of system status;

[0129] The extended state prediction equation has the following mathematical expression:

[0130]

[0131] in, Forecast results for extended states;

[0132] The mathematical expression of the measurement prediction equation is:

[0133]

[0134] Among them, Y k+1 / k To measure the results of the prediction equation.

[0135] Compared with the prior art, the present invention has the following beneficial effects:

[0136] 1. The present invention utilizes distance measurement and angle measurement information to improve the accuracy of orbit determination. The principle is simple and easy to program and calculate. It is applicable to the field of deep space formation satellite orbit determination. Compared with traditional VLBI measurement, the orbit determination accuracy is improved.

[0137] 2. The present invention utilizes the high-precision distance and angle measurement information of the space gravitational wave detection mission itself, without introducing other variables.

[0138] 3. Compared with the linearization method, the present invention does not need to solve the Jacobian matrix of the dynamic equation, which simplifies the difficulty of solving this task. BRIEF DESCRIPTION OF THE DRAWINGS

[0139] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings:

[0140] Figure 1 The present invention provides a flow chart of the method. DETAILED DESCRIPTION

[0141] The present invention is described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those of ordinary skill in the art, several changes and improvements can also be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0142] This scheme establishes high-precision orbit equations for space gravitational wave detection formation satellites, including system state equations and system measurement equations. In the system state equations, the forces on the formation satellites need to consider the gravity of the sun, the perturbations of the eight planets in the solar system and the moon. The system measurement equations include the position velocities of the three satellites and the inter-satellite distances and inter-satellite angles; secondly, the unscented Kalman filter is used for sampling transformation to calculate the sampling weights. Finally, recursion is performed according to the formula.

[0143] A method for determining an orbit of an optical formation satellite provided by the present invention comprises:

[0144] Step 1: Establish a dynamic model of the formation satellite system based on the satellite force conditions;

[0145] Step 2: Obtaining inter-satellite distance and inter-satellite angle information through the dynamic model;

[0146] Step 3: Establish the state equation and observation equation of the formation satellite system;

[0147] Step 4: Based on step 3, the random variables are subjected to nonlinear transformation and then statistically averaged, and then Kalman filtering is performed to eliminate the deviation caused by linearization.

[0148] Step 5: Construct the iterative prediction equation based on steps 3 and 4.

[0149] Step 1: Establish a dynamic model of the formation satellite system based on the satellite force conditions;

[0150] Specifically, taking the heliocentric J2000 coordinate system as the reference coordinate system, the position velocity vectors of satellite 1, satellite 2, and satellite 3 are expressed as (X 1 ,Y 1 ,Z 1 ,V x1 ,V y1 ,V z1 ) T 、(X 2 ,Y 2 ,Z 2 ,V x2 ,V y2 ,V z2 ) T and (X 3 ,Y 3 ,Z 3 ,V x3 ,V y3 ,V z3 ) T .

[0151] The mathematical expression of the dynamic model in the heliocentric J2000 coordinate system is:

[0152]

[0153] in, is the satellite velocity vector; v k Represents the velocity vectors of the three satellites, k = 1, 2, 3;

[0154]

[0155] in, is the satellite’s acceleration vector; μ 0 represents the solar gravitational constant; r k Represents the position vectors of three satellites, k = 1, 2, 3; r k is the satellite position vector; R i represents the position vectors of the eight planets, namely Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and the Moon, i = 1, 2, ..., 8, 9; μ i Respectively represent the gravitational constants of the eight planets and the moon, i = 1, 2, ..., 8, 9. i A scalar representing the positions of the eight planets.

[0156] Step 2: further obtaining inter-satellite distance and inter-satellite angle information through the dynamic model;

[0157] Specifically, the inter-satellite angles include azimuth and altitude, which are converted from the reference star coordinate system to the inertial coordinate system.

[0158] In other words, the inter-satellite angle information is the information of A and E, where A and E are the inter-satellite azimuth and altitude angles respectively. The mathematical expression is:

[0159]

[0160]

[0161] Among them, ρ x , y , z is the position component of the slave star in the master star's own system; ρ represents the inter-satellite distance;

[0162] Let u c 、i c ,Ω c Represent the latitude argument, orbital inclination and right ascension of the reference star respectively. Then the transformation matrix from the J2000 coordinate system to the reference star coordinate system is expressed as:

[0163] B I =M z [u c ]⊙M x [i c ]⊙M z [Ω c ] (7)

[0164] Among them, B I is the conversion matrix from J2000 to the local system; M z Indicates rotation around the Z axis; M x represents rotation around the X-axis, and ⊙ represents the rotation order.

[0165] The following process variables are then obtained for calculation;

[0166]

[0167] Among them, ρ is the inter-satellite distance vector, and the mathematical expression is:

[0168]

[0169] in, represents partial derivative; ρ represents the intersatellite distance vector; r represents the position vector; B I represents the conversion matrix from J2000 to the local system; T represents the transposition symbol;

[0170] Step 3: Establish the state equation and observation equation of the formation satellite system;

[0171] Specifically, the state equation and the measurement equation correspond to the system state equation and the system measurement equation, respectively;

[0172] The system state equation, mathematical expression is:

[0173] X k+1 =f(X k ,w k ,k) (1)

[0174] Among them, X k+1 is the state variable at time k+1; f is the dynamic equation; k is the state variable at time k;

[0175] X k is the position velocity of the three formation satellites, which is an 18-dimensional vector; k is the model error of the system;

[0176] The system measurement equation, mathematical expression is:

[0177] Y k+1 =h(X k+1 ,e k+1 ,k+1) (2)

[0178] Among them, Y k+1 Including the position and speed of the three formation satellites, as well as inter-satellite ranging and angle measurement information; k+1 is the measurement error, w k ~N(0,Q k ),e k+1 ~N(0,R k+1 );Q k is the error variance of the system model, R k+1 is the measurement error variance of the system; h represents the observation equation of the system;

[0179] This proposal attempts to use the ranging and angle measurement information in space gravitational wave detection to improve the accuracy of orbit determination. Therefore, the system measurement equation includes the position velocity of the three satellites and the distance and angle between them.

[0180] In other words, the measurement contents of the system measurement equation include: the positions, speeds, inter-satellite distances and inter-satellite angles of the three satellites.

[0181] h=(X 1 ,Y 1 ,Z 1 ,V x1 ,V y1 ,V z1 ,X2 ,Y 2 ,Z 2 ,V x2 ,V y2 ,V z2 ,X 3 ,Y 3 ,Z 3 ,V x3 ,V y3 ,V z3 ,

[0182] ρ 12 ,A 12 ,E 12 ,ρ 23 ,A 23 ,E 23 ,ρ 31 ,A 31 ,E 31 ) T (12)

[0183] Among them, h is the observation equation, that is, the measurement vector; ρ 12 is the distance between satellite 1 and satellite 2; 2 23 is the distance between satellite 2 and satellite 3; ρ 31 is the distance between satellite 3 and satellite 1; A 12 is the azimuth between satellite 1 and satellite 2; A 23 is the azimuth between satellite 2 and satellite 3; ρ 31 is the azimuth between satellite 3 and satellite 1; E 12 is the elevation angle between satellite 1 and satellite 2; E 23 is the elevation angle between satellite 2 and satellite 3; E 31 is the elevation angle between satellite 3 and satellite 1; X, Y, Z, V x ,V y ,V z are the components of the three satellite position velocities;

[0184] Specifically, X 1 , Y 1 , Z 1 is the position vector of satellite 1; X 2 , Y 2 , Z 2 is the position vector of satellite 2;

[0185] X 3 , Y 3 , Z 3 is the position vector of satellite 3; V x1 、V y1 、V z1 is the velocity vector of satellite 1; V x2、V y2 、V z2 is the velocity vector of satellite 2; V x3 、V y3 、V z3 is the velocity vector of satellite 3;

[0186] The observation matrix H about the inter-satellite distance ρ and the inter-satellite angles A and E n , which is part of h, and its mathematical expression is:

[0187]

[0188] Among them, h n represents the observation equation related to the inter-satellite distance and inter-satellite angle; x represents the position vector; x, y, z represent the position scalar; V x 、V y 、V z represents the velocity scalar;

[0189] Step 4: Based on step 3, the random variables are subjected to nonlinear transformation and then statistically averaged, and then Kalman filtering is performed to eliminate the deviation caused by linearization.

[0190] Suppose that the L-dimensional variable x is a mean Variance is P x Gaussian random vector. y is a nonlinear transformation function of the random variable x:

[0191] y=f(x)

[0192] Among them, y represents the result of nonlinear transformation; x represents random variable; f represents nonlinear transformation;

[0193] The statistical characteristics of y can be obtained by the following Unscented transformation. For random variables x, use the mean and variance P x Expand and get 2L+1 new vectors X i ,i=0,1,…2L, where:

[0194]

[0195] Where k is a scalar that controls the distance from each point to the mean. Let the corresponding weight W of the weighted sum be i is a set of scalars and is the following expression

[0196]

[0197] and

[0198] Let the expansion vector X iThe nonlinear transformation value is The expression is

[0199]

[0200] Then the mean of the nonlinear transformation function y of the random variable x is and variance P y It can be calculated by the following formula

[0201]

[0202] The 2L+1 new vector X formed by the yellow formula i The set of is the Unscented transformation point set.

[0203] Assume that the state equation and measurement equation of the discrete nonlinear system in the step from k to k+1 are:

[0204] X k+1 =f(X k ,w k ,k)

[0205] Y k+1 =h(X k+1 ,e k+1 ,k+1)

[0206] Among them, X k+1 is the state variable at time k+1; f is the dynamic equation; k is the state variable at time k;

[0207] X k is the position velocity of the three formation satellites, which is an 18-dimensional vector; k is the model error of the system;

[0208] Y k+1 It includes the position and velocity of the three formation satellites, as well as the inter-satellite ranging and angle measurement information. k+1 is the measurement error, w k ~N(0,Q k ),e k+1 ~N(0,R k+1 );Q k is the system model error variance, R k+1 is the system measurement error variance; h represents the measurement equation.

[0209] The system iteration step is t, and the state estimation value of the previous step is calculated at time (k+1)T. Perform Unscented transformation to obtain the extended vector X i , when filtering, let the state estimate of the previous step That is, the current state mean, with the state estimate variance P x,k+1That is, the expanded equation is as follows; the system is the formation satellite system.

[0210] The mathematical expression of the extended equation is:

[0211]

[0212] Provide system status forecast.

[0213] Step 5: Construct the iterative prediction equation based on steps 3 and 4.

[0214] The mathematical expression of the extended state prediction equation is:

[0215] χ k+1 / k,i =f(χ k+1,i ,k),i=0,1,2,…2L

[0216] Among them, χ k+1 / k,i is the result of the extended state forecast; k+1,i Indicates the current state;

[0217] The system state prediction equation, mathematical expression is:

[0218]

[0219] in, is the result of system state prediction; i represents weight;

[0220] Error covariance prediction equation, mathematical expression is:

[0221]

[0222] Among them, P x,k+1 / k is the result of error covariance prediction; Q k is the model error variance, Predict the results of system status.

[0223] The mathematical expression of the extended state prediction equation is:

[0224]

[0225] Predict results for extended states.

[0226] The measurement prediction equation, mathematical expression is:

[0227]

[0228] Among them, Y k+1 / k To measure the results of the prediction equation.

[0229] The prediction process ends, according to the measured value Y k+1 and the measurement error variance R k+1 Perform status updates and corrections.

[0230] The measurement error variance update equation is expressed as follows:

[0231]

[0232] P y,k+1 is the result of the measurement error variance update equation, P k+1 is the measurement error variance;

[0233] The state and measurement error covariance update equation is expressed as:

[0234]

[0235] Among them, P xy,k+1 is the result of updating the covariance of state and measurement error.

[0236] The Kalman gain matrix, mathematical expression is:

[0237]

[0238] Among them, K k+1 is the Kalman gain matrix; The inverse of the measurement error variance update value.

[0239] The state correction equation, mathematical expression is:

[0240]

[0241] in, is the result of state correction;

[0242] The state estimation error variance correction equation is expressed as follows:

[0243]

[0244] Among them, P x,k+1 is the result of state estimation error variance correction.

[0245] The present invention also provides an optical formation satellite orbit determination system, which can be implemented by executing the process steps of the optical formation satellite orbit determination method, that is, those skilled in the art can understand the optical formation satellite orbit determination method as a preferred implementation of the optical formation satellite orbit determination system.

[0246] An optical formation satellite orbit determination system provided by the present invention comprises:

[0247] Module M1: Establishing a dynamic model of a formation satellite system according to the force conditions of the satellites; the formation satellite system is used for formation satellites;

[0248] Module M2: based on the dynamic model, obtaining the inter-satellite distance and inter-satellite angle;

[0249] Module M3: establishing a state equation and an observation equation according to the inter-satellite distance and the inter-satellite angle;

[0250] Module M4: nonlinearly transforming random variables to obtain nonlinear transformation results; statistically averaging the nonlinear transformation results, and then filtering to obtain unbiased results;

[0251] Module M5: Based on the unbiased results and the state equation and observation equation, an iterative prediction equation is constructed to determine the orbits of the formation satellites.

[0252] Those skilled in the art know that, in addition to realizing the system and its various devices, modules, and units provided by the present invention in a purely computer-readable program code, the system and its various devices, modules, and units provided by the present invention can be completely realized in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps to realize the same functions. Therefore, the system and its various devices, modules, and units provided by the present invention can be considered as a hardware component, and the devices, modules, and units included therein for realizing various functions can also be regarded as structures within the hardware component; the devices, modules, and units for realizing various functions can also be regarded as both software modules for realizing the method and structures within the hardware component.

[0253] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which does not affect the essence of the present invention. In the absence of conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. A method for determining the orbit of optical formation satellites, characterized in that: include: Step S1: establishing a dynamic model of a formation satellite system according to the force conditions of the satellites; the formation satellite system is used for formation satellites; Step S2: Based on the dynamic model, obtaining the inter-satellite distance and inter-satellite angle; Step S3: establishing a state equation and an observation equation according to the inter-satellite distance and the inter-satellite angle; Step S4: nonlinearly transform the random variable to obtain a nonlinear transformation result; Performing statistical averaging on the nonlinear transformation results and then filtering to obtain unbiased results; Step S5: Based on the unbiased result and the state equation and observation equation, an iterative prediction equation is constructed to determine the orbits of the formation satellites.

2. The optical formation satellite orbit determination method according to claim 1, characterized in that: In step S1, the dynamic model uses the heliocentric J2000 coordinate system as the reference coordinate system, and the mathematical expression is: in, is the satellite velocity vector; v k represents the velocity vector of the formation satellite, k = 1, 2, 3; in, is the satellite's acceleration vector; μ0 represents the solar gravitational constant; r k represents the position vector of the formation satellite, k = 1, 2, 3; r k is the satellite position vector; R i represents the position vectors of the eight planets, namely Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and the Moon, i = 1, 2, ..., 8, 9; μ i Represent the gravitational constants of the eight planets and the moon, i=1,2,...,8,9; R i A scalar representing the positions of the eight planets.

3. The optical formation satellite orbit determination method according to claim 2, characterized in that: In step S2, the inter-satellite angle includes the inter-satellite azimuth angle and the inter-satellite altitude angle; The mathematical expression of the inter-satellite azimuth is: Where A represents the inter-satellite azimuth; ρ x , y , z is the position component of the slave star in the master star's own system, determined based on the dynamic model; The mathematical expression of the intersatellite elevation angle is: Where E represents the inter-satellite elevation angle; ρ represents the inter-satellite distance; The mathematical expression of the inter-satellite distance is: in, represents partial derivative; ρ represents the intersatellite distance vector; r represents the position vector; B I represents the conversion matrix from J2000 to the local system; T represents the transposition symbol; The conversion matrix from J2000 to the local system is expressed as: B I =M z [u c ]⊙M x [i c ]⊙M z [Ω c ] (7) Among them, u c 、i c With Ω c Respectively represent the latitude argument, orbital inclination and right ascension of the ascending node of the reference star; M z Indicates rotation around the Z axis; M x represents rotation around the X-axis, and ⊙ represents the rotation order.

4. The optical formation satellite orbit determination method according to claim 3, characterized in that: In the step S3, the state equation is the system state equation; the measurement equation is the system measurement equation; The system state equation, mathematical expression is: X k+1 =f(X k ,w k ,k) (1) Among them, X k+1 is the state variable at time k+1; f is the dynamic equation; k is the state variable at time k; X k is the position velocity of the three formation satellites, which is an 18-dimensional vector; k is the model error of the system; The system measurement equation, mathematical expression is: Y k+1 =h(X k+1 ,e k+1 ,k+1) (2) Among them, Y k+1 The position, velocity, inter-satellite distance and angle measurement information of the three formation satellites; k+1 is the measurement error, w k ~N(0,Q k ),e k+1 ~N(0,R k+1 );Q k is the error variance, R k+1 is the measurement error variance; h is the observation equation; The formation satellites include: satellite 1, satellite 2 and satellite 3; According to the formation satellites and the observation matrix, the mathematical expression of the observation equation is: h=(X1,Y1,Z1,V x1 ,V y1 ,V z1 ,X2,Y2,Z2,V x2 ,V y2 ,V z2 ,X3,Y3,Z3,V x3 ,V y3 ,V z3 , r 12 ,A 12 ,E 12 ,r 23 ,A 23 ,E 23 ,r 31 ,A 31 ,E 31 ) T (12) Among them, ρ 12 is the distance between satellite 1 and satellite 2; ρ 23 is the distance between satellite 2 and satellite 3; ρ 31 is the distance between satellite 3 and satellite 1; A 12 is the azimuth between satellite 1 and satellite 2; A 23 is the azimuth between satellite 2 and satellite 3; A 31 is the azimuth between satellite 3 and satellite 1; E 12 is the elevation angle between satellite 1 and satellite 2; E 23 is the elevation angle between satellite 2 and satellite 3; E 31 is the elevation angle between satellite 3 and satellite 1; X1, Y1, Z1 are the position vector of satellite 1; X2, Y2, Z2 are the position vector of satellite 2; X3, Y3, Z3 are the position vector of satellite 3; V x1 、V y1 、V z1 is the velocity vector of satellite 1; V x2 、V y2 、V z2 is the velocity vector of satellite 2; V x3 、V y3 、V z3 is the velocity vector of satellite 3; The mathematical expression of the observation matrix is: Among them, H b represents the observation matrix; x, y, z represent the satellite position scalar; V x 、V y 、V z represents the satellite's velocity scalar; In step S4, the nonlinear transformation random variable is expressed as follows: y=f(x) Among them, y represents the result of nonlinear transformation; x represents random variable; f represents nonlinear transformation; In step S4, the nonlinear transformation result is statistically averaged by Unscented transformation, that is, 2L+1 extension vectors are constructed according to the random variable, and the expression is: Among them, X i represents the extended vector, i=0,1,…2L; represents the mean of the random variable; P x represents the variance of the random variable; L represents the dimension; k represents the scalar; An Unscented transformation point set is formed according to the 2L+1 extension vectors; and an unbiased result is obtained after Kalman filtering is performed on the Unscented transformation point set.

5. The optical formation satellite orbit determination method according to claim 4, characterized in that: In the step S5, an iterative prediction equation is constructed based on the unbiased result, the system state equation and the system measurement equation; The iterative prediction equations include: an extended state prediction equation, a system state prediction equation, an error covariance prediction equation, an extended state prediction equation and a measurement prediction equation; The extended state prediction equation has the following mathematical expression: χ k+1 / k,i =f(χ k+1,i ,k),i=0,1,2,…2L Among them, χ k+1 / k,i is the result of the extended state forecast; k+1,i Indicates the current state; The system state prediction equation has the mathematical expression: in, is the result of system state prediction; W i represents weight; The error covariance prediction equation has the mathematical expression: Among them, P x,k+1 / k is the result of error covariance prediction; Q k is the model error variance, Predict the results of system status; The extended state prediction equation has the following mathematical expression: in, Forecast results for extended states; The mathematical expression of the measurement prediction equation is: Among them, Y k+1 / k To measure the results of the prediction equation.

6. An optical formation satellite orbit determination system, characterized in that: include: Module M1: Establishing a dynamic model of a formation satellite system according to the force conditions of the satellites; the formation satellite system is used for formation satellites; Module M2: based on the dynamic model, obtaining the inter-satellite distance and inter-satellite angle; Module M3: establishing a state equation and an observation equation according to the inter-satellite distance and the inter-satellite angle; Module M4: nonlinearly transform random variables to obtain nonlinear transformation results; Performing statistical averaging on the nonlinear transformation results and then filtering to obtain unbiased results; Module M5: Based on the unbiased results and the state equation and observation equation, an iterative prediction equation is constructed to determine the orbits of the formation satellites.

7. The optical formation satellite orbit determination system according to claim 6, characterized in that: In the module M1, the dynamic model uses the heliocentric J2000 coordinate system as the reference coordinate system, and the mathematical expression is: in, is the satellite velocity vector; v k represents the velocity vector of the formation satellite, k = 1, 2, 3; in, is the satellite's acceleration vector; μ0 represents the solar gravitational constant; r k represents the position vector of the formation satellite, k = 1, 2, 3; r k is the satellite position vector; R i represents the position vectors of the eight planets, namely Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and the Moon, i = 1, 2, ..., 8, 9; μ i Represent the gravitational constants of the eight planets and the moon, i=1,2,...,8,9; R i A scalar representing the positions of the eight planets.

8. The optical formation satellite orbit determination system according to claim 7, characterized in that: In the module M2, the inter-satellite angle includes the inter-satellite azimuth angle and the inter-satellite altitude angle; The mathematical expression of the inter-satellite azimuth is: Where A represents the inter-satellite azimuth; ρ x , y , z is the position component of the slave star in the master star's own system, determined based on the dynamic model; The mathematical expression of the intersatellite elevation angle is: Where E represents the inter-satellite elevation angle; ρ represents the inter-satellite distance; The mathematical expression of the inter-satellite distance is: in, represents partial derivative; ρ represents the intersatellite distance vector; r represents the position vector; B I represents the conversion matrix from J2000 to the local system; T represents the transposition symbol; The conversion matrix from J2000 to the local system is expressed as: B I =M z [u c ]⊙M x [i c ]⊙M z [Ω c ] (7) Among them, u c 、i c With Ω c Respectively represent the latitude argument, orbital inclination and right ascension of the ascending node of the reference star; M z Indicates rotation around the Z axis; M x represents rotation around the X-axis, and ⊙ represents the rotation order.

9. The optical formation satellite orbit determination system according to claim 8, characterized in that: In the module M3, the state equation is the system state equation; the measurement equation is the system measurement equation; The system state equation, mathematical expression is: X k+1 =f(X k ,w k ,k) (1) Among them, X k+1 is the state variable at time k+1; f is the dynamic equation; k is the state variable at time k; X k is the position velocity of the three formation satellites, which is an 18-dimensional vector; k is the model error of the system; The system measurement equation, mathematical expression is: Y k+1 =h(X k+1 ,e k+1 ,k+1) (2) Among them, Y k+1 The position, velocity, inter-satellite distance and angle measurement information of the three formation satellites; k+1 is the measurement error, w k ~N(0,Q k ),e k+1 ~N(0,R k+1 );Q k is the error variance, R k+1 is the measurement error variance; h is the observation equation; The formation satellites include: satellite 1, satellite 2 and satellite 3; According to the formation satellites and the observation matrix, the mathematical expression of the observation equation is: h=(X1,Y1,Z1,V x1 ,V y1 ,V z1 ,X2,Y2,Z2,V x2 ,V y2 ,V z2 ,X3,Y3,Z3,V x3 ,V y3 ,V z3 , r 12 ,A 12 ,r 12 ,r 23 ,A 23 ,E 23 ,r 31 ,A 31 ,E 31 ) T (12) Among them, ρ 12 is the distance between satellite 1 and satellite 2; ρ 23 is the distance between satellite 2 and satellite 3; ρ 31 is the distance between satellite 3 and satellite 1; A 12 is the azimuth between satellite 1 and satellite 2; A 23 is the azimuth between satellite 2 and satellite 3; A 31 is the azimuth between satellite 3 and satellite 1; E 12 is the elevation angle between satellite 1 and satellite 2; E 23 is the elevation angle between satellite 2 and satellite 3; E 31 is the elevation angle between satellite 3 and satellite 1; X1, Y1, Z1 are the position vector of satellite 1; X2, Y2, Z2 are the position vector of satellite 2; X3, Y + , Z3 is the position vector of satellite 3; V x1 、V y1 、V z1 is the velocity vector of satellite 1; V x2 、V y2 、V z2 is the velocity vector of satellite 2; V x3 、V y3 、V z3 is the velocity vector of satellite 3; The mathematical expression of the observation matrix is: Among them, H b represents the observation matrix; x, y, z represent the satellite position scalar; V x 、V y 、V z represents the satellite's velocity scalar; In the module M4, the nonlinear transformation random variable is expressed as follows: y=f(x) Among them, y represents the result of nonlinear transformation; x represents random variable; f represents nonlinear transformation; In the module M4, the nonlinear transformation results are statistically averaged by Unscented transformation, that is, 2L+1 extension vectors are constructed according to the random variables, and the expression is: Among them, X i represents the extended vector, i=0,1,…2L; represents the mean of the random variable; P x represents the variance of the random variable; L represents the dimension; k represents the scalar; An Unscented transformation point set is formed according to the 2L+1 extension vectors; and an unbiased result is obtained after Kalman filtering is performed on the Unscented transformation point set.

10. The optical formation satellite orbit determination system according to claim 9, characterized in that: In the module M5, an iterative prediction equation is constructed based on the unbiased result, the system state equation and the system measurement equation; The iterative prediction equations include: an extended state prediction equation, a system state prediction equation, an error covariance prediction equation, an extended state prediction equation and a measurement prediction equation; The extended state prediction equation has the following mathematical expression: χ k+1 / k,i =f(χ k+1,i ,k),i=0,1,2,…2L Among them, χ k+1 / k,i is the result of the extended state forecast; k+1,i Indicates the current state; The system state prediction equation has the mathematical expression: in, is the result of system state prediction; W i represents weight; The error covariance prediction equation has the mathematical expression: Among them, P x,k+1 / k is the result of error covariance prediction; Q k is the model error variance, Predict the results of system status; The extended state prediction equation has the following mathematical expression: in, Forecast results for extended states; The mathematical expression of the measurement prediction equation is: Among them, Y k+1 / k To measure the results of the prediction equation.

Images