A method for orbit determination of giant constellations based on configuration rotation

Abstract

This invention proposes a method for orbit determination of a giant constellation based on configuration rotation, comprising: Solving the configuration under rank-deficient equations: constructing the partial derivative matrix of inter-satellite observations with respect to the state, constructing a system of equations, and solving the system of equations using the L-M algorithm to obtain the constellation configuration of the entire constellation rotating around the Earth's center, thereby obtaining the state of each satellite in the constellation; Determining the rotation Euler angles: establishing the functional relationship between satellite-to-ground ranging and satellite position vectors, station position vectors, and offsets, performing Taylor expansion to obtain the linear relationship between ranging and rotation Euler angles, constructing a system of linear equations, and calculating the rotation Euler angles using least squares; Constellation configuration correction: constructing a rotation matrix using the obtained rotation Euler angles to correct the constellation configuration, obtaining accurate orbital information for the constellation. This invention has low dependence on the number of satellite-to-ground observation arc segments; low dependence on the length of satellite-to-ground observation arc segments; and saves observation resources.

Description

Technical Field

[0001] This invention belongs to the field of aerospace, specifically relating to a method for orbit determination of giant constellations based on configuration rotation. Background Technology

[0002] Currently, the common approach to determining the orbits of mega-constellations is still the same as that used for smaller constellations. This involves using measurements between satellites and between satellites and ground stations to perform optimal estimations, and then determining the orbit using least squares or Kalman filtering (joint satellite-ground orbit determination). For small constellations like GPS and BeiDou, which have dozens of satellites, a limited number of ground stations can observe each satellite. However, for future mega-constellations consisting of thousands or tens of thousands of satellites, given the limited ground observation resources, it will be impossible to observe every single satellite.

[0003] In satellite-ground joint orbit determination, because inter-satellite measurements are estimated together with satellite-ground measurements, the accuracy of orbit determination often decreases sharply as the number of satellite-ground measurement arcs decreases. When satellite-ground measurement data is lacking, the constellation will experience overall "drift." Furthermore, the duration of each satellite-ground measurement arc also affects the accuracy of traditional orbit determination methods. Finally, although the accuracy of traditional orbit determination methods is related to the abundance of measurement data, it cannot specify the minimum amount of measurement data required to determine the orbit of the constellation Gigantopithecus, and it is clearly unable to guarantee orbit determination accuracy when ground measurement data is scarce.

[0004] Therefore, for giant constellations, ensuring the minimum ground observation resources required to guarantee orbit determination accuracy, and ensuring orbit determination accuracy under the premise of scarce ground measurement data, are urgent problems to be solved. Summary of the Invention

[0005] To address the problems existing in the prior art, this invention proposes a method for orbit determination of giant constellations based on configuration rotation, specifically including the following steps:

[0006] Step 1, Solving the configuration under the rank-deficient equation: Construct the partial derivative matrix of the inter-satellite observations with respect to the state. Use this partial derivative matrix and the deviation between the observations and the actual values ​​to construct a system of equations. Use the LM algorithm to solve the system of equations to obtain the constellation configuration of the constellation rotating around the Earth's center, thereby obtaining the state of each satellite in the constellation.

[0007] Step 2, Determination of the rotating Euler angle: Establish the functional relationship between the satellite-to-ground ranging and the satellite position vector, the station position vector, and the offset. Perform a Taylor expansion on this function to obtain the linear relationship between the ranging and the rotating Euler angle. Construct a system of linear equations using n sets of data, and calculate the rotating Euler angle using least squares.

[0008] Step 3, Constellation Configuration Correction: Construct a rotation matrix using the obtained Euler angles, and use this rotation matrix to correct the constellation configuration obtained in Step 1 to obtain accurate orbital information of the constellation.

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

[0010]

[0011] Where D represents the inter-satellite distance, 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 frame, y represents the y-axis coordinate of the satellite in the inertial frame, and v z , represents the velocity of the satellite in the z-axis direction in the inertial frame; the partial derivatives of D with respect to the z-coordinate, x-velocity, and y-velocity 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 incorporating the damping factor of the LM algorithm, equation (2) is improved, and the improved set of equations is as follows:

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

[0017] Where I is the identity matrix and ε is the damping factor, with a value of 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 an inertial frame, assuming the true coordinates of a satellite in the constellation are column vector P, and its offset is represented by rotation matrix A, then the satellite coordinates obtained in step 1 are AP. Let the station coordinates be column vector R, and the measured distance between the satellite and the station be D. true The distance calculated by AP is D. init ;

[0020] D true D init The following relationships exist 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 obtain the following relationship:

[0024]

[0025] Let L be the left side of equation (6). true .

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

[0027] Let A T It is a coordinate rotation matrix consisting of three rotations in the order of zyx:

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

[0029] Let the initial value of the iteration be α0β0γ0. A0 is calculated from α0β0γ0. According to equation (6), when the coordinate rotation matrix A0 is known, the corresponding L0 is obtained at this time:

[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 on the three angles and ignoring second-order minterms, we have:

[0034]

[0035] Substituting equation (10) into equation (9), we get: ...

[0037] L true -L0=P T (CxCyCz*dα+CxCyCz*dβ+CxCyCz*dγ)R(11)

[0038] Equation (11) can be written as a linear equation:

[0039]

[0040] The process of constructing the system of equations and solving for the rotational Euler angles in step 2 is as follows:

[0041] When using n sets of satellite-to-ground observation data, the row vector in equation (12) is expanded as follows:

[0042]

[0043] The system of equations for solving the three Euler angles was obtained using n sets of satellite-to-ground observation data:

[0044]

[0045] The system of equations (14) is a statically indeterminate system of equations, which can be solved by least squares as follows:

[0046]

[0047]

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

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

[0050] By applying the following corrections to the state of each satellite in the constellation obtained in step 1, the overall constellation can be corrected:

[0051]

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

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

[0054] 1. Low dependence on the number of star-to-ground observation arcs.

[0055] Traditional methods estimate both inter-satellite and satellite-to-ground measurements simultaneously. When satellite-to-ground measurements are limited, to ensure optimal estimation, the overall constellation drifts in the orbit determination results. That is, after orbit determination, the inter-satellite measurements differ little from the actual values, while the inter-satellite measurements differ significantly, because the limited satellite-to-ground observation data is insufficient to suppress the overall constellation drift. In this invention, satellite-to-ground and inter-satellite measurements are used separately. This ensures that satellite-to-ground observation data does not affect the convergence of the optimal constellation configuration estimation, and that inter-satellite observation data does not prevent the rotation matrix from converging to the true rotation matrix. Therefore, this method has the advantage of not requiring a large amount of satellite-to-ground measurement data to suppress rotation.

[0056] 2. Low dependence on the length of the star-to-ground observation arc.

[0057] Traditional orbit determination methods require each arc segment to have a certain length to accurately determine the orbits of individual stars connected to the station, thus determining the overall constellation's status. However, this method treats the entire constellation as a whole, improving overall orbit determination accuracy regardless of the arc segment used. Therefore, by selecting different arc segments, even those with very short lengths, this method can still guarantee orbit determination accuracy.

[0058] 3. Saves observation resources.

[0059] Theoretically, the minimum number of measurements is 3. After determining the constellation configuration using relative orbit determination, the entire mega-constellation can be considered as a rigid body. A rigid body in any space has six degrees of freedom: three for the position of its center of mass and three for its overall orientation. Taking advantage of the fact that the mega-constellation's center of mass is located at the Earth's center, there are only three degrees of freedom (position), i.e., three independent variables to be determined. For these three unknowns, only three independent equations are needed to solve them. Attached Figure Description

[0060] Figure 1 The orbit determination accuracy curves of this invention under different numbers of observation arc segments;

[0061] Figure 2 This invention presents orbit determination accuracy curves under different observation shield lengths when the observation arc length is reduced.

[0062] Figure 3 This is a curve showing the orbit determination accuracy under different numbers of observation arcs when the length of the observation arc is reduced, using traditional methods. Detailed Implementation

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

[0064] First, it is proven that for a constellation position determined by autonomous orbital determination, the rotation is always the same as the actual position of the constellation, meaning that the rotation does not change over time.

[0065] Taking any pair of satellites m and n in a constellation that can measure distances, the coordinates of these two satellites in the inertial frame can be expressed using six roots:

[0066]

[0067] a, e, i, ω, Ω, θ represent the semi-major axis, eccentricity, orbital inclination, distance from perigee, right ascension of the ascending node, and true anomaly, respectively, of the six roots. Rotating from the orbital plane coordinate system to the inertial frame, the coordinates pass through the negative distance from perigee about the z-axis, the negative orbital inclination about the x-axis, and the negative right ascension of the ascending node about the z-axis.

[0068] In the two-body model, after time T, since the satellite's semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, and perigee distance remain unchanged, only the true perigee angle changes with time. Therefore, their coordinates in the inertial frame are represented as follows:

[0069]

[0070] Here, Δθ represents the true anomaly angle that the satellite changes over time T.

[0071] If both satellites have rotated by the same angle around the Earth's center at the initial moment, their coordinates at the initial moment would be:

[0072]

[0073] Where A is the coordinate rotation matrix.

[0074] As can be seen from equations ⑤ and ⑥, after the satellite rotates around the Earth's center, its semi-major axis, eccentricity, and initial true anomaly angle remain unchanged. Therefore, after time T, the change in the true anomaly angle remains constant for the two satellites rotating around the Earth's center at the initial moment, and their coordinates 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 a coordinate rotation matrix, for interstellar distances we have:

[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 the coordinate rotation matrix at any moment differs from the real constellation by the same amount. Furthermore, since matrix A is determined by three independent angles, and any rotation satisfies the original equations, it can be concluded that the constellation obtained using only interstellar measurements differs from the real constellation only in its rotation around the Earth's center.

[0081] Based on the above evidence, this invention proposes a method for orbit determination of giant constellations based on configuration rotation, specifically including the following steps:

[0082] Step 1, Solving the configuration under the rank-deficient equation: Construct the partial derivative matrix of the inter-satellite observations with respect to the state. Use this partial derivative matrix and the deviation between the observations and the actual values ​​to construct a system of equations. Use the LM algorithm to solve the system of equations to obtain the constellation configuration of the constellation rotating around the Earth's center, 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 inter-satellite distance, 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 frame, y represents the y-axis coordinate of the satellite in the inertial frame, and v z , represents the velocity of the satellite in the z-axis direction in the inertial frame; the partial derivatives of D with respect to the z-coordinate, x-velocity, and y-velocity 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 incorporating the damping factor of the LM algorithm, equation (2) is improved, and the improved set of equations is as follows:

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

[0091] Where I is the identity matrix and ε is the damping factor, with a value of 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 an inertial frame, assuming the true coordinates of a satellite in the constellation are column vector P, and its offset is represented by rotation matrix A, then the satellite coordinates obtained in step 1 are AP. Let the station coordinates be column vector R, and the measured distance between the satellite and the station be D. true The distance calculated by AP is D. init ;

[0094] D true D init The following relationships exist 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 obtain the following relationship:

[0098]

[0099] Let L be the left side of equation (6). true .

[0100] Step 2, Determination of the rotating Euler angle: Establish the functional relationship between the satellite-to-ground ranging and the satellite position vector, the station position vector, and the offset. Perform a Taylor expansion on this function to obtain the linear relationship between the ranging and the rotating Euler angle. Construct a system of linear equations using n sets of data, and calculate the rotating Euler angle using least squares.

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

[0102] Let A T It is a coordinate rotation matrix consisting of three rotations in the order of zyx:

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

[0104] Let the initial value of the iteration be α0β0γ0. A0 is calculated from α0β0γ0. According to equation (6), when the coordinate rotation matrix A0 is known, the corresponding L0 is obtained at this time:

[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 on the three angles and ignoring second-order minterms, we have:

[0109]

[0110] Substituting equation (10) into equation (9), we get: ...

[0112] L true -L0=P T (CxCyCz*dα+CxCyCz*dβ+CxCyCz*dγ)R(11)

[0113] Equation (11) can be written as a linear equation:

[0114]

[0115] The process of constructing the system of equations and solving for the rotational Euler angles in step 2 is as follows:

[0116] When using n sets of satellite-to-ground observation data, the row vector in equation (12) is expanded as follows:

[0117]

[0118] The system of equations for solving the three Euler angles was obtained using n sets of satellite-to-ground observation data:

[0119]

[0120] The system of equations (14) is a statically indeterminate system of equations, which can be solved by least squares as follows:

[0121]

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

[0123] Step 3, Constellation Configuration Correction: Construct a rotation matrix using the obtained Euler angles, and use this rotation matrix to correct the constellation configuration obtained in Step 1 to obtain accurate orbital information of the constellation.

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

[0125] By applying the following corrections to the state of each satellite in the constellation obtained in step 1, the overall constellation can be corrected:

[0126]

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

[0128] Simulation test:

[0129] The following simulation test was conducted using the constellation Giants:

[0130] Table 1. Zodiac Parameters

[0131]

[0132] In this constellation, 108 stars are used for star-to-ground measurements, with each star having measurements with four ground stations.

[0133] Each satellite measures distance to the preceding and following satellites in the same orbit every ten seconds, and measures distance to two satellites in the same phase in adjacent orbital planes every ten seconds.

[0134] Inter-satellite measurements are subject to the following forms of error:

[0135]

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

[0137] The relevant geophysical constants are as follows:

[0138] Table 2 Geophysical constants

[0139]

[0140] Figure 1 The figure shows the orbit determination accuracy curves under different numbers of observation arcs when using this invention. 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 entire network. The lower this index, the higher the orbit determination accuracy. Figure 2 To illustrate the use of this invention, orbit determination accuracy curves under different observation arc lengths were obtained when the observation arc length was reduced. Figure 3This paper presents orbit determination accuracy curves for different numbers of observation arcs when using traditional methods and reducing the length of the observation arcs. Simulation results demonstrate that traditional orbit determination methods suffer some accuracy loss in scenarios with low numbers and short arcs, while this invention can still maintain the orbit determination accuracy of giant constellations.

Claims

1. A method for orbit determination of giant constellations based on configuration rotation, comprising the following steps: Step 1, Solving the configuration under the rank-deficient equation: Construct the partial derivative matrix of the inter-satellite observations with respect to the state. Use this partial derivative matrix and the deviation between the observations and the actual values ​​to construct a system of equations. Use the LM algorithm to solve the system of equations to obtain the constellation configuration of the constellation rotating around the Earth's center, thereby obtaining the state of each satellite in the constellation. Step 2, Determination of the rotating Euler angle: Establish the functional relationship between the satellite-to-ground ranging and the satellite position vector, the station position vector, and the offset. Perform a Taylor expansion on this function to obtain the linear relationship between the ranging and the rotating Euler angle. Construct a system of linear equations using n sets of data, and calculate the rotating Euler angle using least squares. Step 3, Constellation Configuration Correction: Construct a rotation matrix using the obtained Euler angles, and use this rotation matrix to correct the constellation configuration obtained in Step 1 to obtain the accurate orbital information of the constellation; The partial derivative matrix of the constellation observations in step 1 is: in, D represents the inter-satellite distance, 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 frame, y represents the y-axis coordinate of the satellite in the inertial frame, and v z , represents the velocity of the satellite in the z-axis direction in the inertial frame; the partial derivatives of D with respect to the z-coordinate, x-velocity, and y-velocity are omitted in the formula; Step 2, establishing the functional relationship between satellite-to-ground ranging and offset, is as follows: In an inertial frame, assuming the true coordinates of a satellite in the constellation are column vector P, and its offset is represented by rotation matrix A, then the satellite coordinates obtained in step 1 are AP. Let the station coordinates be column vector R, and the measured distance between the satellite and the station be D. true The distance calculated by AP is D. init ; D true D init The following relationships exist between the station coordinates R, satellite coordinates P, and coordinate rotation matrix A: D true 2 =(P-R) T (P-R) (2) D init 2 =(AP-R) T (AP-R) (3) Subtracting equation (4) from equation (5), we obtain the following relationship: Let L be the left side of equation (6). true ; Step 2, deriving the linear relationship between the measured quantity and Euler angles, is as follows: Let A T It is a coordinate rotation matrix consisting of three rotations in the order of zyx: A T =Cx(α)Cy(β)Cz(γ) (5) Let the initial value of the iteration be α0β0γ0. A0 is calculated from α0β0γ0. According to equation (6), when the coordinate rotation matrix A0 is known, the corresponding L0 is obtained at this time: L0=P T (A0 T -E)R (6) Subtracting equation (6) from equation (8), we have: THE true -L0=P T (TO T -A0 T )R (7) A T In A0 T Expanding on the three angles and ignoring second-order minterms, we have: Substituting equation (10) into equation (9), we get: Equation (11) can be written as a linear equation:

2. The method for orbit determination of a giant constellation based on configuration rotation according to claim 1, characterized in that, The steps for constructing and solving the orbit determination equations in step 1 are as follows: Autonomous orbit determination using the least squares method: X=(B'B) -1 B'Y (11) By incorporating the damping factor of the LM algorithm, equation (2) is improved, and the improved set of equations is as follows: X=(B'B+εI) -1 B'Y (12) Where I is the identity matrix, and ε is the damping factor, with a value of 0.01 to 0.

001.

3. The method for orbit determination of a giant constellation based on configuration rotation according to claim 1, characterized in that, Step 2 involves constructing a system of equations to solve for the rotational Euler angles, as follows: When using n sets of satellite-to-ground observation data, the row vector in equation (12) is expanded as follows: The system of equations for solving the three Euler angles was obtained using n sets of satellite-to-ground observation data: The system of equations (14) is a statically indeterminate system of equations, which can be solved by least squares as follows: Through k iterations, α0β0γ0 is iterated to the true rotation angle.

4. The method for orbit determination of a giant constellation based on configuration rotation according to claim 1, characterized in that, Step 3, the constellation configuration correction process, is as follows: By applying the following corrections to the state of each satellite in the constellation obtained in step 1, the overall constellation can be corrected: Where x,y,z,v x ,v y ,v z The equation represents the true position and velocity of each satellite in the constellation in the inertial frame. The six state variables on the right side of the equation represent the state of the satellite calculated in step 1. Matrix A is a coordinate rotation matrix that rotates the satellite in the three directions zyx in sequence by α, β, and γ.

Images