Giant constellation orbit determination method based on configuration rotation

Abstract

The invention provides a giant constellation orbit determination method based on configuration rotation, which comprises the following steps of: configuration solution under a deficit equation: constructing a partial derivative matrix of inter-satellite observed quantity to a state, constructing an equation set, and solving the equation set by using an L-M algorithm to obtain a constellation configuration in which a constellation wholly rotates around the earth center, therefore, the state of each satellite in the constellation is obtained; determining a rotation Euler angle: establishing a function relationship between satellite-ground ranging and a satellite position vector, a survey station position vector and offset, performing Taylor expansion to obtain a linear relationship between ranging and the rotation Euler angle, constructing a linear equation set, and obtaining the rotation Euler angle through least square calculation; and constellation configuration correction: constructing a rotation matrix by using the obtained rotation Euler angle, and correcting the constellation configuration to obtain accurate orbit information of the constellation. The method is low in degree of dependence on the number of satellite-earth observation arc segments; the degree of dependence on the length of a satellite-earth observation arc section is low; and observation resources are saved.

Description

Technical Field

[0001] The invention belongs to the field of aerospace, and in particular relates to a giant constellation orbit determination method based on configuration rotation. Background Art

[0002] At present, the orbit determination problem of giant constellations is still generally solved in the same way as that of small constellations, that is, the measurements between satellites and the measurements between satellites and ground stations are used to make the best estimate, and the orbit is determined by the least squares or Kalman filtering method (joint satellite-ground orbit determination). For small constellations with dozens of satellites such as GPS and Beidou constellations, limited ground stations can observe each satellite. However, for giant constellations consisting of thousands or tens of thousands of constellations in the future, it is impossible to observe each satellite given the limited ground observation resources.

[0003] In the satellite-ground joint orbit determination method, since intersatellite measurements are estimated together with satellite-ground measurements, the orbit determination accuracy often decreases sharply with the reduction of satellite-ground measurement arcs. When satellite-ground measurement data is lacking, the constellation will "drift" as a whole. At the same time, the duration of each satellite-ground measurement arc will also affect the orbit determination accuracy of the traditional method. Finally, although the accuracy of the traditional orbit determination method is related to the richness of the measurement data, it cannot explain how much measurement data is required to complete the orbit determination of the giant star constellation, and it is obviously impossible to guarantee the orbit determination accuracy when ground measurement data is scarce.

[0004] Therefore, for giant constellations, ensuring the minimum ground observation resources for orbit determination accuracy and ensuring orbit determination accuracy under the premise of scarce ground measurement data are urgent issues to be resolved. Summary of the invention

[0005] In order to solve the problems existing in the above-mentioned prior art, the present invention proposes a giant constellation orbit determination method based on configuration rotation, which specifically includes the following steps:

[0006] Step 1: Solve the configuration under the rank-deficient equation: Construct the partial derivative matrix of the satellite intersatellite observation quantity with respect to the state, use the partial derivative matrix and the deviation between the observation quantity and the actual value to construct the equation group, use the LM algorithm to solve the equation group, and obtain the constellation configuration of the entire constellation rotating around the center of the earth, thereby obtaining the state of each satellite in the constellation;

[0007] Step 2, determination of the rotation Euler angle: establish a functional relationship between the satellite-to-ground ranging and the satellite position vector, the station position vector, and the offset, perform Taylor expansion on the function, obtain the linear relationship between the ranging and the rotation Euler angle, construct a linear equation system using n sets of data, and calculate the rotation Euler angle through least squares;

[0008] Step 3: constellation configuration correction: construct a rotation matrix using the obtained rotation Euler angles, and use the rotation matrix to correct the constellation configuration obtained in step 1 to obtain the accurate orbit information of the constellation.

[0009] The partial derivative matrix of the constellation observations in step 1 is:

[0010]

[0011] Where D represents the distance between each satellite, r represents the number of inter-satellite measurements, c represents the number of satellites in the constellation, x represents the x-axis coordinate of the satellite in the inertial system, y represents the y-axis coordinate of the satellite in the inertial system, and v z , represents the speed of the satellite in the z-axis direction in the inertial system; the partial derivatives of D with respect to the z-direction coordinate, x-direction speed, and y-direction speed are omitted in the formula.

[0012] The steps for constructing and solving the orbit determination equations described in step 1 are as follows:

[0013] Autonomous orbit determination using the least squares method:

[0014] X=(B'B) -1 B'Y(2)

[0015] By adding the LM algorithm damping factor, equation (2) is improved, and the improved equation group is:

[0016] X=(B'B+εI) -1 B'Y(3)

[0017] Among them, I is the unit matrix, and ε is the damping factor, which ranges from 0.01 to 0.001.

[0018] The process of establishing the functional relationship between satellite-to-ground ranging and offset in step 2 is as follows:

[0019] In the inertial system, assume that the true coordinates of a satellite in the constellation are column vectors P, and its offset is represented by the rotation matrix A. The satellite coordinates obtained in step 1 are AP, and the station coordinates are column vectors R. The measured value of the distance between the satellite and the station is D. true , the distance calculated by AP is D init ;

[0020] D true , D init The following relationship exists between the station coordinates R, satellite coordinates P, and coordinate rotation matrix A:

[0021] D true 2 =(PR) T (PR)(4)

[0022] Dinit 2 =(AP-R) T (AP-R)(5)

[0023] Subtracting equation (4) from equation (5), we get the following relationship:

[0024]

[0025] The left side of equation (6) is denoted as L true .

[0026] The process of deriving the linear relationship between the measured quantity and the Euler angle described in step 2 is as follows:

[0027] Assume A T It is a coordinate rotation matrix consisting of three angles rotated in the order of zyx:

[0028] A T =Cx(α)Cy(β)Cz(γ)(7)

[0029] Assume that the initial value of the iteration is α0β0γ0, and A0 is calculated from α0β0γ0. According to formula (6), when the coordinate rotation matrix A0 is known, the corresponding L0 is calculated as follows:

[0030] L0=P T (A0 T -E)R(8)

[0031] Subtracting equation (6) from equation (8), we have:

[0032] L true -L0=P T (A T -A0 T )R(9)

[0033] A T In A0 T Expanding the three angles and ignoring the second-order small terms, we have:

[0034]

[0035] Substituting formula (10) into formula (9), we have:

[0037]

[0038] Formula (11) is written as a linear equation:

[0039]

[0040] The process of constructing the equation group in step 2 to solve the rotation Euler angle is as follows:

[0041] When n sets of satellite-earth observation data are used, the row vector in equation (12) is expanded to:

[0042]

[0043] The equations for solving the three Euler angles when using n sets of satellite-earth observation data are obtained:

[0044]

[0045] Equation (14) is a statically indeterminate system and can be solved by least squares as follows:

[0046]

[0047]

[0048] After k iterations, α0β0γ0 is iterated to the actual rotation angle.

[0049] The constellation configuration correction process described in step 3 is as follows:

[0050] For the state of each satellite in the constellation obtained in step 1, the following correction can be performed to correct the entire constellation:

[0051]

[0052] Among them, x,y,z,vx,vy,v z Represents the true position and velocity of each satellite in the constellation in the inertial system. The six state quantities on the right side of the equation represent the state of the satellite calculated in step 1. The A matrix is ​​a coordinate rotation matrix that rotates α, β, and γ in the three directions of zyx in sequence.

[0053] The beneficial effects of the present invention compared with the prior art are as follows:

[0054] 1. Low dependence on the number of satellite-to-earth observation arcs.

[0055] Traditional methods estimate inter-satellite measurements and satellite-to-ground measurements simultaneously. When there are few satellite-to-ground measurements, in order to ensure the optimal estimate, the orbit determination result of the overall constellation drift will appear, that is, after orbit determination, the difference between the inter-satellite measurement value and the actual value is small, while the difference between the satellite-to-ground measurement value and the actual value is large, because a small amount of satellite-to-ground observation data is difficult to suppress the overall drift of the constellation. In the present invention, satellite-to-ground measurements are used separately from inter-satellite measurements, which ensures that the satellite-to-ground observation data does not affect the convergence of the optimal estimate of the constellation configuration, and the inter-satellite observation data does not hinder the rotation matrix from converging to the real rotation matrix. Therefore, this method has the advantage of not needing to use a large amount of satellite-to-ground measurement data to suppress rotation.

[0056] 2. Low dependence on the length of the satellite-to-earth observation arc.

[0057] The traditional orbit determination method requires that each arc segment has a certain length, so that each satellite connected to the station can be accurately determined to obtain an accurate orbit, and then determine the status of the entire constellation. In this method, since the entire constellation is considered as a whole, no matter which arc segment is used, the overall orbit determination accuracy can be improved. Therefore, as long as different arc segments are selected, even if each arc segment is only a short length, this method can still ensure the orbit determination accuracy.

[0058] 3. Save observation resources.

[0059] Theoretically, the minimum number of measurements is 3. After the constellation configuration is determined by relative orbit determination, the entire giant constellation can be regarded as a rigid body. A rigid body in any space has six degrees of freedom, three for the position of the center of mass and three for the overall attitude. Taking advantage of the fact that the center of mass of the giant constellation is at the center of the earth, there are only three degrees of freedom, that is, three independent variables to be determined. For three unknowns, only three independent equations are needed to solve. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 The orbit determination accuracy curve of the present invention under different numbers of observation arc segments;

[0061] Figure 2 It is the orbit determination accuracy curve under different observation conditions when the observation arc length is reduced in the present invention;

[0062] Figure 3 This is the orbit determination accuracy curve for different numbers of observation arcs when the length of the observation arc is reduced using the traditional method. DETAILED DESCRIPTION

[0063] The present invention is described in detail below with reference to the accompanying drawings.

[0064] First, it is proved that the position of the constellation determined by autonomous orbit determination for a configuration always differs from the true position of the constellation by the same rotation, i.e., the rotation is time-invariant.

[0065] Take any pair of satellites m and n in the constellation with ranging, then the coordinates of these two satellites in the inertial system can be expressed by six numbers:

[0066]

[0067] a, e, i, ω, Ω, θ are the semi-major axis, eccentricity, orbital inclination, perigee angular distance, ascending node right ascension, and true anomaly of the six elements. Rotating from the orbital plane coordinate system to the inertial system is to rotate around the z axis to rotate the negative perigee angular distance, around the x axis to rotate the negative orbital inclination, and around the z axis to rotate the negative ascending node right ascension.

[0068] In the two-body model, after T time, since the satellite's semi-major axis, eccentricity, orbital inclination, right ascension of ascending node, and perigee angular distance remain unchanged, only the true anomaly angle will change with time, so the coordinates of the two in the inertial system are expressed as:

[0069]

[0070] Where Δθ represents the true anomaly angle of the satellite that changes during the time T.

[0071] If both satellites rotate around the center of the earth by the same angle at the initial moment, their coordinates at the initial moment become:

[0072]

[0073] Where A is the coordinate rotation matrix.

[0074] From equations ⑤ and ⑥, we can see that after the satellite rotates around the center of the earth, the semi-major axis, eccentricity, and initial true anomaly angle do not change. Therefore, after the two satellites rotate around the center of the earth for T time at the initial moment, the change in true anomaly angle remains unchanged, and the coordinates of the two become:

[0075]

[0076] Combining equations ⑤ to ⑧, we have:

[0077] P m,t ′-P n,t ′=A(P m,t -P n,t )⑨

[0078] Since A is the coordinate rotation matrix, the inter-satellite distance is:

[0079] |P m,t ′-P n,t ′|=|A||P m,t -P n,t |=|P m,t -P n,t |=D⑩

[0080] Therefore, after the constellation rotates as a whole at the initial moment, the distance observation at any moment remains unchanged, and at any moment it differs from the real constellation by the same coordinate rotation matrix. Since the A matrix is ​​determined by three independent angles, and any rotation satisfies the original equations, it can be concluded that the constellation obtained by orbit determination using only intersatellite measurements differs from the real constellation only by the rotation around the center of the earth.

[0081] Based on the above proof, the present invention proposes a giant constellation orbit determination method based on configuration rotation, which specifically includes the following steps:

[0082] Step 1: Solve the configuration under the rank-deficient equation: Construct the partial derivative matrix of the satellite intersatellite observation quantity with respect to the state, use the partial derivative matrix and the deviation between the observation quantity and the actual value to construct the equation group, use the LM algorithm to solve the equation group, and obtain the constellation configuration of the entire constellation rotating around the center of the earth, thereby obtaining the state of each satellite in the constellation;

[0083] The partial derivative matrix of the constellation observations in step 1 is:

[0084]

[0085] Where D represents the distance between each satellite, r represents the number of inter-satellite measurements, c represents the number of satellites in the constellation, x represents the x-axis coordinate of the satellite in the inertial system, y represents the y-axis coordinate of the satellite in the inertial system, and v z , represents the speed of the satellite in the z-axis direction in the inertial system; the partial derivatives of D with respect to the z-direction coordinate, x-direction speed, and y-direction speed are omitted in the formula.

[0086] The steps for constructing and solving the orbit determination equations described in step 1 are as follows:

[0087] Autonomous orbit determination using the least squares method:

[0088] X=(B'B) -1 B'Y(2)

[0089] By adding the LM algorithm damping factor, equation (2) is improved, and the improved equation group is:

[0090] X=(B'B+εI) -1 B'Y(3)

[0091] Among them, I is the unit matrix, and ε is the damping factor, which ranges from 0.01 to 0.001.

[0092] The process of establishing the functional relationship between satellite-to-ground ranging and offset in step 2 is as follows:

[0093] In the inertial system, assume that the true coordinates of a satellite in the constellation are column vectors P, and its offset is represented by the rotation matrix A. The satellite coordinates obtained in step 1 are AP, and the station coordinates are column vectors R. The measured value of the distance between the satellite and the station is D. true , the distance calculated by AP is D init ;

[0094] D true , D init The following relationship exists between the station coordinates R, satellite coordinates P, and coordinate rotation matrix A:

[0095] D true 2=(PR) T (PR)(4)

[0096] D init 2 =(AP-R) T (AP-R)(5)

[0097] Subtracting equation (4) from equation (5), we get the following relationship:

[0098]

[0099] The left side of equation (6) is denoted as L true .

[0100] Step 2, determination of the rotation Euler angle: establish a functional relationship between the satellite-to-ground ranging and the satellite position vector, the station position vector, and the offset, perform Taylor expansion on the function, obtain the linear relationship between the ranging and the rotation Euler angle, construct a linear equation system using n sets of data, and calculate the rotation Euler angle through least squares;

[0101] The process of deriving the linear relationship between the measured quantity and the Euler angle described in step 2 is as follows:

[0102] Assume A T It is a coordinate rotation matrix consisting of three angles rotated in the order of zyx:

[0103] A T =Cx(α)Cy(β)Cz(γ)(7)

[0104] Assume that the initial value of the iteration is α0β0γ0, and A0 is calculated from α0β0γ0. According to formula (6), when the coordinate rotation matrix A0 is known, the corresponding L0 is calculated as follows:

[0105] L0=P T (A0 T -E)R(8)

[0106] Subtracting equation (6) from equation (8), we have:

[0107] L true -L0=P T (A T -A0 T )R(9)

[0108] A T In A0 T Expanding the three angles and ignoring the second-order small terms, we have:

[0109]

[0110] Substituting formula (10) into formula (9), we have:

[0112]

[0113] Formula (11) is written as a linear equation:

[0114]

[0115] The process of constructing the equation group in step 2 to solve the rotation Euler angle is as follows:

[0116] When n sets of satellite-to-earth observation data are used, the row vector in equation (12) is expanded to:

[0117]

[0118] The equations for solving the three Euler angles when using n sets of satellite-earth observation data are obtained:

[0119]

[0120] Equation (14) is a statically indeterminate system and can be solved by least squares as follows:

[0121]

[0122] After k iterations, α0β0γ0 is iterated to the actual rotation angle.

[0123] Step 3: constellation configuration correction: construct a rotation matrix using the obtained rotation Euler angles, and use the rotation matrix to correct the constellation configuration obtained in step 1 to obtain the accurate orbit information of the constellation.

[0124] The constellation configuration correction process described in step 3 is as follows:

[0125] For the state of each satellite in the constellation obtained in step 1, the following correction can be performed to correct the entire constellation:

[0126]

[0127] Among them, x,y,z,vx,vy,v z Represents the true position and velocity of each satellite in the constellation in the inertial system. The six state quantities on the right side of the equation represent the state of the satellite calculated in step 1. The A matrix is ​​a coordinate rotation matrix that rotates α, β, and γ in the three directions of zyx in sequence.

[0128] Simulation test:

[0129] The following giant star constellations were used for simulation testing:

[0130] Table 1 Constellation parameter table

[0131]

[0132] In this constellation, there are 108 stars with satellite-ground measurements, and each star has measurements with four ground stations.

[0133] Each satellite measures the distance with the previous and next satellites in the same orbit every ten seconds, and measures the distance with two satellites in the same phase in adjacent orbital planes every ten seconds.

[0134] There are errors in intersatellite measurements in the following forms:

[0135]

[0136] Where ε is a normal noise. The constant term b is the dominant factor, so in the subsequent simulation process, only b is estimated.

[0137] The relevant geophysical constants are as follows:

[0138] Table 2 Geophysical constants

[0139]

[0140] Figure 1 The orbit determination accuracy curves of the present invention with different numbers of observation arcs are shown in the figure. In the figure, the horizontal axis represents the number of observation arcs used, and the vertical axis represents the root mean square error rms of the whole network. The lower the index, the higher the orbit determination accuracy. Figure 2 For adopting the present invention, when the length of the observation arc is reduced, the orbit determination accuracy curves under different observation angles are obtained; Figure 3 The orbit determination accuracy curves under different numbers of observation arcs are shown when the length of the observation arc is reduced using the traditional method. The simulation results show that the accuracy of the traditional orbit determination method is lost in the case of low number of arcs and short arcs, while the present invention can still maintain the orbit determination accuracy of the giant constellation.

Claims

1. A giant constellation orbit determination method based on configuration rotation, comprising the following steps: Step 1: Solve the configuration under the rank-deficient equation: Construct the partial derivative matrix of the satellite intersatellite observation quantity with respect to the state, use the partial derivative matrix and the deviation between the observation quantity and the actual value to construct the equation group, use the LM algorithm to solve the equation group, and obtain the constellation configuration of the entire constellation rotating around the center of the earth, thereby obtaining the state of each satellite in the constellation; Step 2, determination of the rotation Euler angle: establish a functional relationship between the satellite-to-ground ranging and the satellite position vector, the station position vector, and the offset, perform Taylor expansion on the function, obtain the linear relationship between the ranging and the rotation Euler angle, construct a linear equation system using n sets of data, and calculate the rotation Euler angle through least squares; Step 3: constellation configuration correction: construct a rotation matrix using the obtained rotation Euler angles, and use the rotation matrix to correct the constellation configuration obtained in step 1 to obtain the accurate orbit information of the constellation.

2. The method for determining orbits of a giant constellation based on configuration rotation according to claim 1, characterized in that: The partial derivative matrix of the constellation observations in step 1 is: Where D represents the distance between each satellite, r represents the number of inter-satellite measurements, c represents the number of satellites in the constellation, x represents the x-axis coordinate of the satellite in the inertial system, y represents the y-axis coordinate of the satellite in the inertial system, and v z, It represents the speed of the satellite in the z-axis direction in the inertial system. The partial derivative of D with respect to the z-direction coordinate, x-direction speed, and y-direction speed is omitted in the formula.

3. The method for giant constellation orbit determination based on configuration rotation according to claim 1, characterized in that: The steps for constructing and solving the orbit determination equations described in step 1 are as follows: Autonomous orbit determination using the least squares method: X=(B'B) -1 B'Y(2) By adding the LM algorithm damping factor, equation (2) is improved, and the improved equation group is: X=(B'B+εI) -1 B'Y(3) Among them, I is the unit matrix, and ε is the damping factor, which ranges from 0.01 to 0.

001.

4. The method for orbit determination of a giant constellation based on configuration rotation according to claim 1, characterized in that: The process of establishing the functional relationship between satellite-to-ground ranging and offset in step 2 is as follows: In the inertial system, assume that the true coordinates of a satellite in the constellation are column vectors P, and its offset is represented by the rotation matrix A. The satellite coordinates obtained in step 1 are AP, and the station coordinates are column vectors R. The measured value of the distance between the satellite and the station is D. true , the distance calculated by AP is D init ; D true , D init The following relationship exists between the station coordinates R, satellite coordinates P, and coordinate rotation matrix A: D true 2 =(P-R) T (P-R)(4) D init 2 =(AP-R) T (AP-R)(5) Subtracting equation (4) from equation (5), we get the following relationship: The left side of equation (6) is denoted as L true .

5. The method for giant constellation orbit determination based on configuration rotation according to claim 1, characterized in that: The process of deriving the linear relationship between the measured quantity and the Euler angle described in step 2 is as follows: Assume A T It is a coordinate rotation matrix consisting of three angles rotated in the order of zyx: A T =Cx(α)Cy(β)Cz(γ)(7) Assume that the initial value of the iteration is α0β0γ0, and A0 is calculated from α0β0γ0. According to formula (6), when the coordinate rotation matrix A0 is known, the corresponding L0 is calculated as follows: L0=P T (A0 T -E)R(8) Subtracting equation (6) from equation (8), we have: THE true -L0=P T (TO T -A0 T )R(9) A T In A0 T Expanding the three angles and ignoring the second-order small terms, we have: Substituting formula (10) into formula (9), we have: Formula (11) is written as a linear equation:

6. The method for orbit determination of a giant constellation based on configuration rotation according to claim 1, characterized in that: The process of constructing the equation group in step 2 to solve the rotation Euler angle is as follows: When n sets of satellite-to-earth observation data are used, the row vector in equation (12) is expanded to: The equations for solving the three Euler angles when using n sets of satellite-earth observation data are obtained: Equation (14) is a statically indeterminate system and can be solved by least squares as follows: After k iterations, α0β0γ0 is iterated to the actual rotation angle.

7. The method for giant constellation orbit determination based on configuration rotation according to claim 1, characterized in that: The constellation configuration correction process described in step 3 is as follows: For the state of each satellite in the constellation obtained in step 1, the following correction can be performed to correct the entire constellation: Among them, x,y,z,vx,vy,v z Represents the true position and velocity of each satellite in the constellation in the inertial system. The six state quantities on the right side of the equation represent the state of the satellite calculated in step 1. The A matrix is ​​a coordinate rotation matrix that rotates α, β, and γ in the three directions of zyx in sequence.

Images