High-orbit satellite autonomous orbit determination method based on combination of inter-satellite relative position vector measurement and GNSS

Abstract

The invention provides a method for determining an autonomous orbit of a high-orbit satellite based on combination of inter-satellite relative position vector measurement and a GNSS (Global Navigation Satellite System), and belongs to the field of autonomous navigation. The method comprises the following specific steps: constructing a GEO orbit (geostationary orbit) environment of a satellite, and establishing a gravity model by using an EGM2008 model; establishing a relative position vector measurement model between the spacecrafts and a relative distance measurement model between the spacecrafts and the GNSS satellites; analyzing and fusing the system and measurement errors; performing linearization processing on the measurement equation; and estimating the state of the high-orbit spacecraft by adopting an extended Kalman filter, and determining the absolute position and speed of the spacecraft. On the basis of determining the orbit of the spacecraft by using a GNSS satellite traditionally, direction finding measurement is added between high-orbit spacecrafts, so that the orbit determination precision of the high-orbit spacecrafts is improved, and high-precision orbit determination is realized.

Description

Technical Field

[0001] The invention relates to the field of satellite autonomous navigation and provides a fully autonomous orbit determination method for a high-orbit satellite based on inter-satellite relative position vector measurement. Background Art

[0002] High-orbit satellites play an important role in the fields of communication, navigation, earth observation and meteorology. With the increasing demand for high-orbit satellite applications in various countries, the number of high-orbit satellites has gradually increased, and the requirements for the orbit determination accuracy of high-orbit satellites have become increasingly higher. So far, the autonomous navigation and positioning of high-orbit satellites mainly rely on the sidelobe signals of the Global Navigation Satellite System (GNSS). Due to the constraints of extremely weak signals and susceptibility to interference, as well as the low visibility of high-orbit satellites, high-precision and high-reliability navigation of high-orbit satellites has also been challenged. Therefore, it is considered to introduce new measurement methods to ensure the complete autonomous navigation of high-orbit satellites while improving the accuracy of their absolute orbit determination; and the satellite autonomous navigation method based on inter-satellite relative position vector measurement has become one of the most noteworthy satellite autonomous navigation measurement methods due to its simplicity and the maturity of existing technologies.

[0003] Relative position measurement navigation aims to accurately measure the relative position, distance, angle and relative speed of one object relative to another object through various technical means. With the advancement of sensor technology and computer vision technology, the accuracy and reliability of relative position measurement navigation are constantly improving, and the scope of application is also expanding. The basic idea of ​​satellite inter-satellite relative position vector measurement navigation is to determine the absolute position of two or more satellites by measuring the relative position vector from one satellite to another. This method not only measures the distance between satellites, but also includes azimuth and altitude angles, thereby obtaining complete three-dimensional relative position information. Relative position measurement can be achieved through inter-satellite links and star sensors. Inter-satellite links are used to measure relative distances, and star sensors are used to measure relative azimuths. In order to ensure the feasibility of absolute navigation, the relative position vector needs to refer to an inertial coordinate system, and a time series of relative position measurements is also required.

[0004] Based on the inter-satellite relative position vector measurement technology, the accuracy and reliability of autonomous navigation can be improved. Practical engineering applications include: The GRACE mission uses two satellites to infer the changes in the Earth's gravity field by accurately measuring the changes in the relative distance between them. This technology is not only used for scientific research, but also provides high-precision data support for satellite autonomous navigation. The GRAIL mission draws a high-resolution lunar gravity field image through the relative position measurement of two satellites, realizing the application of inter-satellite relative position measurement in autonomous navigation.

[0005] This patent proposes a method for high-precision autonomous navigation of high-orbit satellites by combining GNSS and relative position vector measurement. This method can further improve the navigation accuracy of high-orbit satellites. Summary of the invention

[0006] The technical problem solved by this application is to make up for the shortcomings of the existing high-orbit satellite autonomous navigation technology and provide a method for determining the absolute orbit of a high-orbit satellite based on the inter-satellite relative position vector measurement. Here, taking the relative position measurement of two high-orbit satellites and the relative distance measurement between the two and a GNSS satellite as an example, the extended Kalman filter is used as a sequential estimation method to explain the technical solution of the invention as follows:

[0007] 1. Establishing the state equation model based on orbital dynamics

[0008] The positions and velocities of the two high-orbit satellites 1 and 2 in the Earth Inertial Coordinate System (ECI) are:

[0009]

[0010] The subscript i indicates the corresponding high-orbit satellite.

[0011] The position and velocity of all high-orbit satellites constitute the state variable x, that is:

[0012]

[0013] For high-orbit satellites, their motion is mainly determined by the Earth's gravity, but is also affected by other perturbations including the gravity of the Sun and the Moon. Therefore, from the basic kinematics and orbital dynamics theory, it can be seen that the state equation satisfied by the state variable can be modeled as the following continuous nonlinear time-varying state equation with random input:

[0014]

[0015] Among them, a im is the modeled acceleration, which may include the acceleration caused by the earth's gravity, three-body gravity, atmospheric drag, and solar radiation pressure. The specific modeling accuracy can be selected according to actual needs and constraints; it is random noise to compensate for the unmodeled acceleration and can generally be modeled as zero-mean Gaussian white noise.

[0016] By discretizing the continuous state equation, we can get the following discrete state equation:

[0017] x k =f(x k-1 ,t k ,t k-1 )+w k-1 (4)

[0018] Among them, w k-1 is the process noise in discrete form, x k Indicates t k The state variables of the spatial target at a certain moment. The function f(·) represents the deterministic part of the state variable transfer relationship between two discrete moments. Generally, there is no explicit expression, but the function value can be obtained by numerical integration based on the input state variables at the previous moment.

[0019] In the extended Kalman filter, the state transfer matrix is ​​used to propagate the covariance of the posterior estimate of the state variable at the previous moment to obtain the covariance of the prior estimate of the state variable at the current moment. k-1 The state transition matrix from time to time t is denoted as Φ(t,t k-1 ), then it satisfies the following matrix differential equation:

[0020] Where A(t) is the coefficient matrix obtained after linearization of the state equation.

[0021]

[0022] Then A(t) can be expressed as

[0023]

[0024] At the same time, the state transfer matrix has the following properties:

[0025] Φ(t k-1 ,t k-1 )=I 12×12 (8)

[0026] By numerically integrating the matrix differential equation satisfied by the above state transfer matrix, we can obtain t k-1 Time to t k The state transfer matrix Φ at time k =Φ(t k ,t k-1 ).

[0027] In addition, the process noise w k-1 The covariance matrix of is also needed in the Kalman filter algorithm. k-1 The process noise matrix from time to time t is recorded as Q(t,t k-1 ), then it satisfies the following matrix differential equation:

[0028]

[0029] Where Q0 is a constant matrix. Let the continuous process noise wit The covariances of (i=1,2) are And order

[0030]

[0031] Then Q0 can be expressed as

[0032]

[0033] At the same time, the discrete form of the process noise covariance matrix satisfies the following initial value conditions:

[0034] Q(t k-1 ,t k-1 )=0 12×12 (12)

[0035] By numerically integrating the matrix differential equation satisfied by the above process noise covariance matrix, we can obtain t k-1 Time to t k The process noise covariance matrix Q at time k =Q(t k ,t k-1 ).

[0036] 2. Establish a measurement model based on relative orientation and distance information

[0037] The measurement information includes the relative orientation observations between two high-orbit satellites, the relative distance between the high-orbit satellites, and the relative distance observations between the high-orbit satellites and the GNSS satellite.

[0038] First, a relative direction finding model between high-orbit satellites is established. The relative direction finding information is the two-dimensional observation quantity z u , can have different specific representations. Relative azimuth observation z u It can be expressed as the right ascension α and declination δ of the relative line of sight vector between two high-orbit satellites in the ECI coordinate system. However, since the relative direction finding information is usually obtained through an optical camera, a more refined modeling method can also be used to express it as the two-dimensional coordinates (Y, Z) of the pixels of the observed satellite in the imaging plane of the optical camera. Next, the relative direction finding model is explained using the first method as an example.

[0039] When the relative azimuth observation is expressed by the right ascension α and declination δ of the relative line of sight vector between two high-orbit satellites in the ECI coordinate system, we have

[0040]

[0041] At the same time, the relative line of sight vector between two high-orbit satellites can be expressed as

[0042]

[0043] Then the relationship between the observed quantity and the relative sight vector is:

[0044]

[0045] Therefore, the measurement equation can be written as

[0046] z u =h u (x)+υ u (16)

[0047] where h u (x) is a nonlinear function of the state variable x, υ u is the measurement noise, and its covariance matrix is ​​R u .

[0048] At the same time, the partial derivative matrix of the relative azimuth observation to the position and velocity of each high-orbit satellite is

[0049]

[0050] The partial derivatives can be specifically expressed as follows

[0051]

[0052] Therefore, the partial derivative matrix of the relative orientation observation with respect to the state variable x is

[0053] H u =[H u1 H u2 ](19)

[0054] Secondly, a relative distance measurement model between high-orbit satellites is established. The relative distance observation between two high-orbit satellites is z r It can be expressed as

[0055] z r =h r (x)+υ r =r 21 +υ r =||r2-r1||+υ r (20)

[0056] Among them, h r (x) is a nonlinear function of the state variable x, v r is the measurement noise, whose covariance is R r .

[0057] At the same time, the partial derivative matrices of the relative distance observations between high-orbit satellites with respect to the position and velocity of each high-orbit satellite are respectively

[0058]

[0059] The partial derivatives can be specifically expressed as follows

[0060]

[0061] Therefore, the partial derivative matrix of the relative distance observation between high-orbit satellites with respect to the state variable x is

[0062] H r =[H r1 H r2 ](twenty three)

[0063] Finally, a relative distance measurement model between high-orbit satellites and GNSS satellites is established. ρ It can be expressed as

[0064]

[0065] Among them, h ρ (x) is a nonlinear function of the state variable x, υ ρ is the measurement noise, and its covariance matrix is ​​R ρ , r GNSS It is the position of the GNSS satellite in the ECI coordinate system that measures the relative distance between the two high-orbit satellites.

[0066] At the same time, the partial derivative matrices of the relative distance observation between each high-orbit satellite and the GNSS satellite with respect to the position and velocity of each corresponding high-orbit satellite are respectively

[0067]

[0068] The partial derivatives can be specifically expressed as follows

[0069]

[0070] Therefore, the partial derivative matrix of the relative distance observation between the high-orbit satellite and the GNSS satellite with respect to the state variable x is

[0071]

[0072] Now, we can build a complete measurement model by combining the above various types of observations. Let the complete observation z, measurement function h(x) and measurement noise υ be as follows:

[0073]

[0074] Then we have the complete measurement equation

[0075] z = h(x) + υ (29) where the measurement noise has the following covariance:

[0076]

[0077] And the partial derivative matrix of the observation with respect to the state variable is

[0078]

[0079] 3. Establish a sequential state estimation algorithm based on extended Kalman filter

[0080] Suppose the state posterior estimate after the k-1 time measurement update is The state estimation error covariance matrix is According to equations (5) to (12), by numerical integration, the state transfer matrix Φ can be obtained. k =Φ(t k ,t k-1 ) and the process noise matrix Q k =Q(t k ,t k-1 ), so the one-step prediction equation at time k is

[0081]

[0082] When obtaining the current moment observation z k Then, according to formula (31), the partial derivative matrix of the measurement equation at this moment with respect to the state variable is calculated as follows:

[0083]

[0084] Finally, the measurement update can be performed according to the extended Kalman filter formula:

[0085]

[0086] This is the a posteriori estimate of the position and velocity of each high-orbit satellite at the current time k and its error covariance.

[0087] The present invention has the following advantages: a method for combining GNSS with relative position vector measurement to perform high-precision autonomous navigation on a high-orbit satellite is proposed, and the method can further improve the navigation accuracy of the high-orbit satellite. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 The method flow chart of the present invention is DETAILED DESCRIPTION

[0089] The present invention proposes a method for autonomous orbit determination of high-orbit satellites based on inter-satellite relative position vector measurement and GNSS combination. Figure 1Here, taking two high-orbit satellites performing relative position measurement and two high-orbit satellites performing relative distance measurement with a GNSS satellite while using extended Kalman filtering as a sequential estimation method as an example, the specific implementation of the present invention is described as follows:

[0090] Step 1: Establish system state equation and measurement equation

[0091] The positions of the two high-orbit spacecraft in the ECI coordinate system are i With speed v i Describe their motion states respectively, and let the state variables be

[0092]

[0093] For each high-orbit satellite, select the appropriate force model to obtain the modeled acceleration a of the spacecraft im , using zero-mean Gaussian white noise w it Compensate for unmodeled acceleration and select appropriate acceleration noise covariance based on the acceleration modeling accuracy The following state equation with random input is obtained:

[0094]

[0095] The measurement equation is established based on the actual relative position information between the two high-orbit satellites and the relative distance information between the high-orbit satellite and the GNSS satellite. Here, the right ascension α and declination δ of the relative line of sight vector between the two high-orbit satellites in the ECI coordinate system are used as the relative azimuth observation z u , and select the relative azimuth observation z u , the relative distance observation between two high-orbit satellites z r And the relative distance observation z between the two high-orbit satellites and the GNSS satellite ρ The measurement noise covariance of u , R r With R ρ , so the measurement equation corresponding to each observation is:

[0096]

[0097] z r =h r (x)+υ r =r 21 +υ r =||r2-r1||+υ r ,,υ r ~N(0,R r ) (38)

[0098]

[0099] Let R = diag([R u R r R ρ ]), thus establishing the following complete measurement equation:

[0100]

[0101] This step corresponds to box ① in the figure, and then proceed to step 2.

[0102] Step 2: Initialize the filter

[0103] Initialize the state variable estimate and its error covariance at the initial moment, that is, select the initial estimate of the state variable and its error covariance This step corresponds to box ② in the figure, and then proceed to step 3.

[0104] Step 3: Time update of state estimate

[0105] According to the state estimate given in the previous step Numerically integrate the state equation to obtain the time update result of the current state estimate:

[0106]

[0107] make

[0108]

[0109] Numerically integrate the following:

[0110]

[0111] Get the state transfer matrix Φ from the previous moment to the current moment k =Φ(t k ,t k-1 ).make

[0112]

[0113] Then perform numerical integration on the following formula:

[0114]

[0115] Get the covariance matrix Q of the discrete process noise from the previous moment to the current moment k =Q(t k ,t k-1 ).

[0116] Finally, the propagation result of the state estimation covariance is calculated as follows:

[0117]

[0118] This step corresponds to box ③ in the figure, and then proceed to step 4.

[0119] Step 4: Calculate the partial derivative matrix of the measurement equation

[0120] First, the partial derivative matrix of each observation component to the position and velocity of the high-orbit satellite is established. The partial derivative matrix of the relative azimuth observation component to the position and velocity of the high-orbit satellite is:

[0121]

[0122] The partial derivative matrix of the relative distance observation between high-orbit satellites with respect to the position and velocity of each high-orbit satellite is:

[0123]

[0124] The partial derivative matrix of the relative distance observation between each high-orbit satellite and the GNSS satellite with respect to the position and velocity of each corresponding high-orbit satellite is:

[0125]

[0126] Secondly, calculate the partial derivative matrix of each measurement equation with respect to the state variable at the current moment. Calculate the partial derivative matrix H of the relative direction finding equation with respect to the state variable as follows: uk :

[0127]

[0128] The partial derivative matrix of the relative distance observation between high-orbit satellites with respect to the state variable x is calculated as follows:

[0129]

[0130] The partial derivative matrix of the relative distance observation between the high-orbit satellite and the GNSS satellite with respect to the state variable x is calculated as follows:

[0131]

[0132] The above partial derivative matrix calculations are all based on Conducted at.

[0133] Finally, the above partial derivative matrices are grouped to obtain the following complete measurement matrix

[0134]

[0135] This step corresponds to box ④ in the figure, and then proceed to step 5.

[0136] Step 5: Measure and update the state estimate

[0137] Based on the current observation z k , the current state prior estimate and its error covariance The partial derivative matrix H of the measurement equation k , measurement noise covariance matrix R k , the state estimate and its error covariance are measured and updated according to the EKF formula:

[0138]

[0139] This gives us the posterior estimate of the state variable at the current moment: and its error covariance

[0140] This step corresponds to box ⑤ in the figure. If the orbit determination task is completed, the algorithm stops; otherwise, it returns to step 3 to update the state estimation time.

Claims

1. A method for autonomous orbit determination of a high-orbit satellite based on inter-satellite relative position vector measurement combined with GNSS, characterized in that: The method comprises the following steps: Step 1: Establish system state equation and measurement equation Firstly, the positions and velocities of all high-orbit satellites in the ECI coordinate system are taken as the state variable x. According to the required accuracy of the dynamic model, appropriate force models such as the earth's gravity, three-body gravity, atmospheric drag and solar radiation pressure are selected to obtain the modeled acceleration of each spacecraft, and random noise is introduced to compensate for the unmodeled acceleration. Based on the basic kinematics and orbital dynamics theory, the state equation with random input is established. Secondly, the measurement information mainly includes the inter-satellite relative position observation z u , the relative distance between stars observed z r And the relative distance observation z between the high-orbit satellite and GNSS ρ When establishing the measurement equation, the observed quantity can be modeled according to the actual acquisition method of the observation information. For example, the relative position observation quantity z between satellites u It can be modeled as the right ascension α and declination δ of the relative line of sight vector between two high-orbit satellites in the ECI coordinate system. However, since the relative direction-finding information is usually obtained through an optical camera, a more refined modeling method can also be used to express it as the two-dimensional coordinates (Y, Z) of the pixels of the observed satellite in the imaging plane of the optical camera. By grouping various types of observations, the following complete measurement equation can be established: z=h(x)+υ Among them, z is the u 、z r With z ρ The complete observable is constructed, h(x) is the mapping function from the state variable to the theoretical measurement value, and υ is the measurement noise. Step 2: Initialize the filter Initialize the state variable estimate and its error covariance at the initial moment, that is, select the initial estimate of the state variable and its error covariance Then proceed to step three. Step 3: Time update of state estimate According to the state estimate given in the previous step The state equation is numerically integrated, and the matrix differential equation satisfied by the state transfer matrix and the discrete form of the process noise covariance matrix is ​​numerically integrated to obtain the time update result of the current state estimate. Then proceed to step 4. Step 4: Calculate the partial derivative matrix of the measurement equation Record the relative position between satellites z u , the relative distance between stars observed z r And the relative distance observation z between the high-orbit satellite and GNSS ρ The partial derivative matrix with respect to the state variables is The values ​​of H uk , H rk With H ρk , then the partial derivative matrix of the measurement equation to the state variable at the current moment is Step 5: Calculate measurement updates According to the current observation z k , the current state prior estimate and its error covariance The partial derivative matrix H of the measurement equation k , measurement noise covariance matrix R k , the state estimate and its error covariance are measured and updated according to the EKF formula: This gives us the posterior estimate of the state variable at the current moment: and its error covariance If the orbit determination task is completed, the algorithm stops; otherwise, it returns to step 3 to update the state estimation time.

Images