Autonomous orbit determination method for low earth orbit mega constellation based on star sensor and inter-star angle measurement

Abstract

The present application relates to satellite orbit determination technical field, specifically to low-orbit giant constellation autonomous orbit determination method based on star sensor and inter-satellite angle measurement. The method comprises: step 1, determining constellation configuration by using inter-satellite measurement: using six number expression each satellite state quantity in constellation, constructing partial derivative matrix of inter-satellite measurement relative to state information, constructing linear equation set according to state quantity, measurement and partial derivative matrix, solving satellite state quantity describing constellation configuration by improved least square method. Step 2, using star sensor and camera measurement, combining constellation configuration in step 1, correcting constellation drift: using star sensor and camera, obtaining the real value of the angle between the direction vector of adjacent satellite and the direction vector of a certain star, then calculating the estimated value of the angle between the direction vector of adjacent satellite and the direction vector of a certain star, then minimizing the difference between the real value and the estimated value of the angle, realizing constellation correction. The above scheme can realize orbit determination without relying on ground station and is easy to manage data.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of satellite orbit determination, in particular to a low-orbit mega constellation autonomous orbit determination method based on star sensors and inter-satellite angle measurement. BACKGROUND

[0002] With the formation of low-orbit mega constellation system, the demand for constellation measurement and control has increased explosively. At present, most low-orbit mega constellations rely on global ground stations for continuous observation, and through the communication between ground stations and satellites, accurate orbit data can be obtained. However, due to the limited ground observation resources, it is impossible to continuously track each satellite, resulting in a significant decrease in orbit determination accuracy. The size of the mega constellation further exacerbates the severity of the ground station shortage.

[0003] To solve the defects of relying on ground stations, many researchers have tried to use inter-satellite links to obtain inter-satellite measurements to obtain relative position and velocity information, thereby reducing the dependence on ground stations. However, the traditional method of applying inter-satellite measurement faces problems such as constellation drift, and has little effect on improving the accuracy and timeliness of the orbit determination of mega constellations. On the other hand, the complexity of inter-satellite communication and data processing also puts higher requirements on the management of the constellation.

[0004] The defects of the prior art can be summarized as follows:

[0005] 1. Over-reliance on ground station measurement data

[0006] Although the inter-satellite link is gradually developing, the current orbit determination still relies heavily on ground station observations. Due to the uneven distribution and coverage of ground observation resources, this dependence is particularly evident in some areas; external factors such as weak signal interference can also cause orbit observation to be interrupted, affecting accuracy.

[0007] 2. Poor complexity and scalability of inter-satellite link management

[0008] For mega constellations, the number of inter-satellite links is very large, and as the number of satellites increases, the maintenance and management of the links also become more complex. Especially in the case of continuous changes in inter-satellite relative motion, existing technical solutions are difficult to handle mega constellations, and the objects targeted by existing technical solutions are small-scale constellations or individual satellites, making it difficult to ensure the stability of mega constellation links and effectively manage data. SUMMARY

[0009] To solve the above technical problems, the present application proposes a low-orbit mega constellation autonomous orbit determination method based on star sensors and inter-satellite angle measurement, which is a low-orbit mega constellation autonomous orbit determination method that does not require ground stations at all, and can realize orbit determination without relying on ground stations and easy data management.

[0010] To achieve the above object, the application adopts the following technical solutions:

[0011] A low-orbit giant constellation autonomous orbit determination method based on star sensors and inter-satellite angle measurement is performed according to the following steps:

[0012] Step 1, determining the constellation configuration by using inter-satellite measurement;

[0013] Step 2, correcting the constellation drift by using star sensors and camera measurement in combination with the constellation configuration in step 1.

[0014] Further, the step 1 of determining the constellation configuration by using inter-satellite measurement comprises the following steps:

[0015] Step 1.1, expressing the state quantity of each satellite in the constellation by using six elements;

[0016] Step 1.2, constructing the partial derivative matrix of inter-satellite measurement with respect to the state information, and constructing a linear equation set according to the state quantity, the measurement and the partial derivative matrix;

[0017] Step 1.3, solving the satellite state quantity describing the constellation configuration by using the improved least square method.

[0018] Further, the step 1.1 of expressing the state quantity of each satellite in the constellation by using six elements comprises the following steps:

[0019] Let m represent the total number of inter-satellite measurements, n represent the total number of satellites in the constellation, and D=[D1 D2 … D m ] T represent the measurement generated when measuring each other in the constellation; D is related to the position of the satellite in the constellation; the state quantity of all satellites in the constellation is expressed by the orbital six elements, i.e. a, e, i, Ω, ω and M, wherein a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the ascending node longitude, ω is the perihelion argument, and M is the mean anomaly;

[0020] The state expression vector of each satellite in the constellation is as follows:

[0021]

[0022] In the formula, X k represents a state quantity to be estimated, wherein k=1, 2, …, 6n, and X represents a vector composed of all state quantities.

[0023] Further, the step 1.2 of constructing the partial derivative matrix of inter-satellite measurement with respect to the state information, and constructing a linear equation set according to the state quantity, the measurement and the partial derivative matrix comprises the following steps:

[0024]

[0025] X = RD

[0026] Where R represents the partial derivative matrix of inter-satellite measurement with respect to state information, and D represents the measurement generated when measuring the inter-satellite measurement.

[0027] Further, the satellite state quantity describing the constellation configuration is solved by the improved least square method in step 1.3, specifically:

[0028] Using the L-M method, the following solving formula is obtained:

[0029] X = (R T R+λE) -1 R T D (3)

[0030] Where the subscript T represents the matrix transpose, λ represents the damping factor, and E is the unit matrix.

[0031] Further, the constellation drift is corrected by using the star sensor and camera measurement in step 2, combined with the constellation configuration in step 1, including the following steps:

[0032] Step 2.1, using the star sensor and camera to obtain the real value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star;

[0033] Step 2.2, calculating the estimated value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star according to the satellite state quantity obtained in step 1;

[0034] Step 2.3, using the iterative method to minimize the difference between the real value and the estimated value of the angle, to realize the constellation correction.

[0035] Further, the real value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star is obtained by using the star sensor and camera in step 2.1, and the specific content is as follows:

[0036] Let B be the conversion matrix from the inertial system to the satellite body-fixed coordinate system, and let A(α,β,γ) = C x (α)C y (β)C z (γ) be a coordinate rotation matrix to be solved, the real coordinates of a certain satellite in the constellation are column vectors AP, the coordinates of the corresponding satellite in the constellation obtained in step 1 are P, the real coordinates of another adjacent satellite are column vectors AQ, and the satellite coordinates obtained by inter-satellite orbit determination are Q; when observing a remote star from the inertial system, the unit line-of-sight vector pointing to the remote star from any satellite is considered as a constant vector

[0037] Then the unit line-of-sight vector L 0 of the star in the satellite body-fixed coordinate system is:

[0038]

[0039] The calculation method of the real value of the included angle between the adjacent satellite direction vector and a certain star direction vector is as follows:

[0040]

[0041] In formula (4) and formula (5), L 0 is a unit line-of-sight vector of a star in a satellite body-fixed coordinate system, is a unit line-of-sight vector of a star in an inertial system, is a unit line-of-sight vector from a p satellite to a q satellite, B p represents a coordinate rotation matrix of a selected satellite from a body system to an inertial system, Q-P represents a difference between a coordinate vector of the selected satellite and a coordinate vector of an adjacent satellite of the selected satellite, and θ true is an included angle between L 0 and , which is the real value of the included angle between the adjacent satellite direction vector and the certain star direction vector.

[0042] Further, the calculation of the estimated value of the included angle between the adjacent satellite direction vector and the certain star direction vector according to the satellite state quantity obtained in step 1 in step 2.2 is as follows:

[0043] Given the attitude of a certain satellite at a time t and the coordinates (B, P, Q) of the satellite and an adjacent satellite, the estimated value cosθ esti of the included angle between the adjacent satellite direction vector and the certain star direction vector is as follows:

[0044]

[0045] In the formula, cosθ esti represents the estimated value of the included angle between the adjacent satellite direction vector and the certain star direction vector, is a unit line-of-sight vector of a star in an inertial system, p represents a coordinate rotation matrix of a selected satellite from a body system to an inertial system, and Q-P represents a difference between a coordinate vector of the selected satellite and a coordinate vector of an adjacent satellite of the selected satellite.

[0046] Further, the use of an iterative method to minimize the difference between the real value and the estimated value of the included angle in step 2.3 to realize constellation correction is as follows:

[0047] The residual r i is defined as the absolute value of the difference between cosθ esti,i and cosθ true,i :

[0048] r i = |cosθ esti,i -cosθtrue,i |,i = 1,2, … N (7)

[0049] where the subscript i represents the i-th group of data, and the value range is 1, 2, … N

[0050] The objective function F(A) is defined as:

[0051]

[0052] It is desired to minimize F(A), which is solved by the Gauss-Newton method.

[0053] An initial θ0is defined, and the residual r i Expand at θ = θ0and ignore the second-order small terms, and we get

[0054]

[0055] Let denote the partial derivative of the residual with respect to each parameter,

[0056] Take the partial derivative of the above equation with respect to δθ, and set the result to 0, and we get:

[0057]

[0058] where The calculation formula of is:

[0059]

[0060] In order to consider multiple groups of data, construct the Jacobian matrix, and stack all the residuals:

[0061]

[0062] Write (17) and (18) into vector form:

[0063] K T R = - (K T K) δθ (19)

[0064]

[0065] Use the iterative formula:

[0066]

[0067] Iterate multiple times until the preset convergence condition or accuracy requirement is reached, and iterate α0, β0, γ0to the true rotation angle to obtain constellation orbit parameters with higher accuracy.

[0068] The application also provides a low-orbit mega constellation autonomous orbit determination system based on a star sensor and inter-satellite angle measurement, which is used for realizing the low-orbit mega constellation autonomous orbit determination method based on the star sensor and the inter-satellite angle measurement.

[0069] Compared with the prior art, the application has the following beneficial technical effects:

[0070] 1. The application makes full use of inter-satellite measurement data and gets rid of the requirement for ground stations. Further, the application uses a star sensor and on-satellite measurement equipment instead of ground measurement data to eliminate the overall drift caused by the full use of inter-satellite measurement data, and can realize autonomous orbit determination without external data.

[0071] 2. The application firstly uses inter-satellite ranging information to perform relative orbit determination, and then uses the star sensor and inter-satellite measurement equipment to rotate the whole constellation, so that the data in the information fusion process is more convenient to manage. Meanwhile, the accuracy is higher than that of the prior art in the case of complex error forms; the performance of the application is the same as that of the prior art when the measurement error form is normal noise; and the orbit determination accuracy of the application is better than that of the prior art when the measurement error form is more complex. BRIEF DESCRIPTION OF DRAWINGS

[0072] Figure 1 For the relationship between the whole network satellite position error and the number of iterations in the simulation verification;

[0073] Figure 2 For the relationship between the whole network satellite velocity error and the number of iterations in the simulation verification. DETAILED DESCRIPTION

[0074] The application innovatively proposes a mega constellation autonomous orbit determination method without ground stations. The method uses inter-satellite ranging information of the mega constellation to perform relative orbit determination and determine the constellation configuration; then, the constellation is regarded as a whole, the star sensor is used to provide star azimuth information, and the camera is used to provide adjacent satellite azimuth information to rotate the constellation as a whole, so as to offset the rotation error caused by the use of only inter-satellite ranging, and realize efficient and accurate mega constellation orbit calculation. The application is particularly suitable for the scenario of mega constellation (6000+ satellites) orbit calculation (the scale is similar to the low-orbit Walker constellation or the hybrid constellation).

[0075] In order to make the purpose, technical scheme and advantages of the application more clear, the embodiments of the application are described in detail below.

[0076] A low-orbit mega constellation autonomous orbit determination method based on a star sensor and inter-satellite angle measurement is performed according to the following steps:

[0077] Step 1: determining the constellation configuration by using inter-satellite measurement, including the following steps:

[0078] The state quantity of each satellite in the constellation is expressed by six roots, a partial derivative matrix of inter-satellite measurement with respect to state information is constructed, a linear equation set is constructed according to the state quantity, the measurement and the partial derivative matrix, and the state quantity of the satellite describing the constellation configuration is solved by the improved least square method.

[0079] Step 1.1, the state quantity of each satellite in the constellation is expressed by six roots, specifically:

[0080] Let m represent the total number of inter-satellite measurements, n represent the total number of satellites in the constellation, D = [D1 D2 … D m ] T represents the measurement generated when measuring each other in the constellation; D is related to the position of the satellite in the constellation; the state quantity of all satellites in the constellation is expressed by the orbital six roots, i.e. a, e, i, Ω, ω, m, wherein a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the longitude of the ascending node, ω is the argument of perigee, and M is the mean anomaly.

[0081] The state of each satellite in the constellation can be expressed in vector form:

[0082]

[0083] In the formula, X k represents a to-be-estimated state quantity, wherein k = 1, 2 … 6n, and X represents a vector composed of all state quantities.

[0084] Step 1.2, a partial derivative matrix of inter-satellite measurement with respect to state information is constructed, a linear equation set is constructed according to the state quantity, the measurement and the partial derivative matrix, and specifically:

[0085]

[0086] X = RD

[0087] In the formula, R represents a partial derivative matrix of inter-satellite measurement with respect to state information, and D represents the measurement generated when measuring each other in the constellation.

[0088] Step 1.3, the state quantity of the satellite describing the constellation configuration is solved by the improved least square method, and specifically:

[0089] The L-M method is used to obtain the following solving formula:

[0090] X = (R T R + λE) -1 R T D (3)

[0091] In the formula, the subscript T represents matrix transposition, λ represents a damping factor, and E is an identity matrix.

[0092] If λE is not added, R TR is not full rank, and cannot be inverted. R matrix not full rank means that the equation group has infinite solutions, and adding a small amount allows to get one of the solutions with higher accuracy.

[0093] It can be proved that the satellite positions obtained successively only differ from the real satellite positions by a rotation around the earth center, and the rotation is time-invariant, thus obtaining the constellation configuration.

[0094] Step 2, using star sensor and camera measurement, combining the constellation configuration in step 1, correcting the constellation drift.

[0095] Using the star sensor and the camera, the real value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star is obtained, and the estimated value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star is calculated according to the satellite state quantity obtained in step 1; the difference between the real value and the estimated value of the angle is minimized by using an iterative method, and the constellation correction is realized.

[0096] Step 2.1, using the star sensor and the camera, the real value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star is obtained, and the estimated value of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star is calculated according to the satellite state quantity obtained in step 1; the difference between the real value and the estimated value of the angle is minimized by using an iterative method, and the constellation correction is realized.

[0097] Let B be the conversion matrix from the inertial system to the satellite body-fixed coordinate system, and let A (α, β, γ) = C x (α)C y (β)C z (γ) be a coordinate rotation matrix to be solved, the real coordinates of a certain satellite in the constellation are column vectors AP, the coordinates of the corresponding satellite in the constellation obtained in step 1 are P, the real coordinates of another adjacent satellite are column vectors AQ, and the satellite coordinates obtained by inter-satellite orbit determination are Q; when observing a remote star from the inertial system, the unit line-of-sight vector pointing to the remote star from any satellite is considered as a constant vector

[0098] Then the unit line-of-sight vector L 0 of the star in the satellite body-fixed coordinate system is:

[0099]

[0100] The real value calculation method of the angle between the direction vector of the adjacent satellite and the direction vector of a certain star is:

[0101]

[0102] In formula (4) and formula (5), L 0 is the unit line-of-sight vector of the star in the satellite body-fixed coordinate system, is the unit line-of-sight vector of the satellite pointing to the remote star in the inertial system, is the unit line-of-sight vector from the p satellite to the q satellite, and B pQ-P is the difference between the selected satellite and its neighboring satellite coordinate vectors, and θ true is the angle between B 0 and P , which is the true value of the angle between the neighboring satellite direction vector and a certain star direction vector.

[0103] Step 2.2, calculate the estimated value of the angle between the neighboring satellite direction vector and a certain star direction vector according to the satellite state quantity obtained in step 1, specifically:

[0104] Given the attitude of a certain satellite at time t and its coordinates (B, P, Q) with neighboring satellites, the estimated value of the angle between the neighboring satellite direction vector and a certain star direction vector cosθ esti is:

[0105]

[0106] In the formula, cosθ esti represents the estimated value of the angle between the neighboring satellite direction vector and a certain star direction vector, is the unit line-of-sight vector of the satellite pointing to a distant star in the inertial system, B p represents the coordinate rotation matrix of the selected satellite from the body system to the inertial system, and Q-P represents the difference between the selected satellite and its neighboring satellite coordinate vectors.

[0107] Step 2.3, use the iterative method to minimize the difference between the true value of the angle and the estimated value of the angle, and realize constellation correction, the specific process is as follows:

[0108] Define the residual r i as the absolute value of the difference between cosθ esti,i and cosθ true,i :

[0109] r i = |cosθ esti,i -cosθ true,i |, i = 1, 2, ……N (7)

[0110] Where the subscript i represents the i-th group of data, and the value range is 1, 2, ……N

[0111] Define the objective function F(A) as:

[0112]

[0113] We hope to minimize F(A). This problem can be solved by Gauss-Newton method.

[0114] Define the initial θ0, and the residual r i is expanded at θ = θ0, and the second-order small term is ignored, which is

[0115]

[0116] with denotes the partial derivative of the residual with respect to each parameter,

[0117]

[0118] Take the partial derivative of the above formula with respect to δθ and set the result to 0, and obtain:

[0119]

[0120] where The calculation formula of is:

[0121]

[0122] In order to consider multiple groups of data, construct the Jacobian matrix, and stack all the residuals:

[0123]

[0124] Write (17) and (18) into vector form:

[0125] K T R = -(K T K) δθ (19)

[0126]

[0127] Use the iterative formula:

[0128]

[0129] Multiple iterations until the preset convergence condition or accuracy requirement is reached, and α0, β0, γ0 are iterated to the true rotation angle, and the constellation orbit parameters with higher accuracy are obtained.

[0130] Based on the same inventive concept, the application also provides a low-orbit giant constellation autonomous orbit determination system based on a star sensor and interstellar angle measurement, which is used to realize the low-orbit giant constellation autonomous orbit determination method based on the star sensor and interstellar angle measurement.

[0131] Simulation verification:

[0132] The following low-orbit giant constellation is used for simulation verification, and Table 1 is shown:

[0133] Table 1 Low-orbit giant constellation parameter table

[0134]

[0135] There are 1296 satellites in the constellation, of which 16 satellites are equipped with star sensors and the required inter-satellite angle measurement devices.

[0136] The relevant geophysical constants are shown in Table 2:

[0137] Table 2 Geophysical constants table

[0138]

[0139] The vector X composed of satellite states contains 7776 unknowns. Each satellite has distance measurement data with the four satellites in front, behind, left and right, and 2592 sets of inter-satellite measurement data are generated at each measurement time point. The partial derivative of the inter-satellite distance measurement vector at N time points with respect to the vector X is obtained as Using formula (3), the satellite state quantity describing the constellation configuration is calculated, and the whole network position error (expressed as root mean square rms) is 155 m, and the whole network velocity error (rms) is 0.17 m / s. The 16 node stars in the constellation measure the stars and the front satellites through sensors, and the cosθ p,q described in formula (5) is obtained. The 16 node stars in the constellation measure the stars and the front satellites through sensors under the condition that there is a 155 m position error, and the cosθ esti described in formula (6) is obtained. By iteratively minimizing the difference between cosθ p,q and cosθ esti , the curves of the change of position error and velocity error with the increase of iteration number are obtained, as shown in Figure 1 and Figure 2 .

[0140] Figure 1 , Figure 2 are respectively the trend graphs of the whole network average position deviation and the whole network average velocity deviation with the decrease of iteration number obtained by applying the method of the present application. Figure 1 In the figure, the horizontal axis is the iteration number, and the vertical axis is the root mean square of the whole network position; Figure 2 In the figure, the horizontal axis is the iteration number, and the vertical axis is the root mean square of the whole network velocity. Both indexes represent the orbit determination accuracy, and the smaller the value is, the higher the accuracy is. As can be seen from the figure, in the process of 8 iterations, the whole network rms of both position and velocity continues to decrease, and finally reaches a very low value, thereby proving the effectiveness and stability of the present application in the precise orbit determination of the giant constellation.

[0141] It should be noted that the embodiments of the present application are preferred embodiments, but not limitations. It should be noted that for ordinary skilled persons in the art, the specific embodiments can be modified or some technical features can be replaced equivalently without departing from the principles of the present application, which should be considered as the scope of the present application.

Claims

1. A low earth orbit mega constellation autonomous orbit determination method based on star sensor and interstellar angle measurement, characterized in that The steps are as follows: Step 1, determining the constellation configuration by using inter-satellite measurements; Step 1.1, expressing the state quantities of each satellite in the constellation by using six numbers; Step 1.2, constructing the partial derivative matrix of inter-satellite measurements with respect to state information, and constructing a linear equation set according to the state quantities, measurements, and the partial derivative matrix; Step 1.3, solving the state quantities of the satellites describing the constellation configuration by using an improved least square method; Step 2, correcting the constellation drift by using star sensor and camera measurements in combination with the constellation configuration in Step 1; Step 2.1, obtaining the real value of the included angle between the direction vector of an adjacent satellite and the direction vector of a certain star by using the star sensor and the camera; Step 2.2, calculating the estimated value of the included angle between the direction vector of an adjacent satellite and the direction vector of a certain star according to the state quantities of the satellites obtained in Step 1; Step 2.3, minimizing the difference between the real value and the estimated value of the included angle by using an iterative method to achieve constellation correction.

2. The method according to claim 1, wherein, In Step 1.1, the state quantities of each satellite in the constellation are expressed by using six numbers, which are specifically as follows: where m denotes the total number of inter-satellite measurements, n denotes the total number of satellites in the constellation, D = [D1 D2...Dm]Tdenotes the vector of inter-satellite measurements, and F = [F1 F2...Fn]Tdenotes the vector of satellite states. m ] T denotes the inter-satellite measurements; D is related to the positions of the satellites in the constellation; the states of all satellites in the constellation are represented by the orbital elements, i.e., a, e, i, Ω, ω, M, where a is the semi-major axis, e is the eccentricity, i is the inclination, Ω is the longitude of the ascending node, ω is the argument of the perigee, and M is the mean anomaly. The state representation vector of each satellite in the constellation is in the form of: where X k represents a certain state quantity to be estimated, where k = 1, 2, …, 6n, and X represents a vector composed of all state quantities.

3. The method according to claim 2, wherein, In Step 1.2, the partial derivative matrix of inter-satellite measurements with respect to state information is constructed, and a linear equation set is constructed according to the state quantities, measurements, and the partial derivative matrix, which are specifically as follows: In the formula, R represents the partial derivative matrix of inter-satellite measurements with respect to state information, X represents a vector composed of all state quantities, and D represents the measurements generated when measuring each other by inter-satellites.

4. The method according to claim 3, wherein, In Step 1.3, the state quantities of the satellites describing the constellation configuration are solved by using an improved least square method, which is specifically as follows: The L-M method is used to obtain the following solving formula: X = (R T R + λE) -1 R T D (3) In the formula, X represents a vector composed of all state quantities, the subscript T represents matrix transposition, λ represents a damping factor, E is a unit matrix, and D represents the measurements generated when measuring each other by inter-satellites.

5. The method of claim 1, wherein the low earth orbit mega constellation autonomous orbit determination based on star sensor and inter-star angle measurement, further comprises: In Step 2.1, the real value of the included angle between the direction vector of an adjacent satellite and the direction vector of a certain star is obtained by using the star sensor and the camera, and the specific content is as follows: Let B be the transformation matrix from the inertial frame to the satellite body-fixed coordinate frame, and let A(α, β, γ) = C x (α)C y (β)C z (γ) be a coordinate rotation matrix to be solved, the real coordinates of a satellite in the constellation are column vectors AP, and the coordinates of the corresponding satellite in the constellation obtained in step 1 are P, the real coordinates of another adjacent satellite are column vectors AQ, and the satellite coordinates obtained by inter-satellite orbit determination are Q; when observing a remote star from the inertial frame, it is considered that the unit line-of-sight vector from any satellite to the remote star is a constant vector Then the unit line-of-sight vector L 0 of the star in the satellite body-fixed coordinate system is: The real value of the included angle between the direction vector of an adjacent satellite and the direction vector of a certain star is calculated as follows: In formula (4) and formula (5), L 0 is a unit line-of-sight vector of a star in a satellite body-fixed coordinate system, is a unit line-of-sight vector of a star in an inertial coordinate system, is a unit line-of-sight vector from a p-th satellite to a q-th satellite, B p represents a coordinate rotation matrix from a body coordinate system of a selected satellite to an inertial coordinate system, Q-P represents a difference between coordinate vectors of the selected satellite and its adjacent satellites, θ true is an angle between L 0 and B , and is a true value of an angle between a direction vector of an adjacent satellite and a direction vector of a certain star.

6. The method of claim 5, wherein the low earth orbit mega constellation autonomous orbit determination based on star sensors and inter-star angular measurements is characterized by, In Step 2.2, the estimated value of the included angle between the direction vector of an adjacent satellite and the direction vector of a certain star is calculated according to the state quantities of the satellites obtained in Step 1, which is specifically as follows: where cos θ esti denotes the estimate of the angle between the direction vector of a neighboring satellite and a certain star direction vector, is the unit line-of-sight vector of the satellite pointing to a remote star in the inertial frame, B p denotes the coordinate rotation matrix from the body frame of the selected satellite to the inertial frame, Q-P denotes the difference between the coordinate vectors of the selected satellite and its neighboring satellites.

7. The method according to claim 6, wherein, In Step 2.3, the difference between the real value and the estimated value of the included angle is minimized by using an iterative method to achieve constellation correction, and the specific process is as follows: The definition of the residual r i is cos θ esti,i The absolute value of the difference between cos θ true,i and cos θ r i = |cos θ esti,i -cos θ true,i |, i = 1, 2,... N (7) Wherein, the subscript i represents the i-th group of data, and the value range is 1, 2, …, N The objective function F(A) is defined as: It is desired to minimize F(A), and the Gauss-Newton method is used to solve it; Define the initial θ0, and the residual r i Expanding at θ = θ0, and neglecting second order terms, we get with denotes the partial derivative of the residual with respect to the respective parameter, The partial derivative of the above formula with respect to δθ is taken, and the result is set to 0, which is: wherein The calculation formula is: In order to consider multiple sets of data, the Jacobian matrix is constructed, and all the residuals are stacked: The formula (17) and (18) are written in vector form as follows: K T R = -(K T K)δθ (19) The iterative formula is used as follows: Multiple iterations are performed until the preset convergence condition or accuracy requirement is reached, and α0, β0, γ0 are iterated to the true rotation angles to obtain higher-precision constellation orbit parameters.

8. A low earth orbit mega constellation autonomous orbit determination system based on star sensor and inter-star angle measurement, characterized in that, A method for implementing the low-orbit mega-constellation autonomous orbit determination based on star sensors and inter-satellite angle measurement according to any one of claims 1-7.

Images