Low-orbit giant constellation autonomous orbit determination method based on star sensor and inter-satellite angle measurement

Abstract

The invention relates to the technical field of satellite orbit determination, in particular to a low-orbit giant constellation autonomous orbit determination method based on a star sensor and inter-satellite angle measurement. The method comprises the following steps: step 1, determining a constellation configuration by using inter-satellite quantity measurement: expressing each satellite state quantity in a constellation by using six elements, constructing a partial derivative matrix of inter-satellite quantity measurement relative to state information, constructing a linear equation set according to the state quantity, the quantity measurement and the partial derivative matrix, and determining a constellation configuration; and solving the satellite state quantity describing the constellation configuration through an improved least square method. Step 2, correcting constellation drift by using a star sensor and a camera quantity measurement and combining the constellation configuration in the step 1: acquiring a true value of an included angle between an adjacent satellite direction vector and a certain fixed star direction vector by using the star sensor and the camera, then calculating an estimated value of the included angle between the adjacent satellite direction vector and the certain fixed star direction vector, and correcting the constellation drift; and the difference between the real value and the estimated value of the included angle is minimized to realize constellation correction. By adopting the scheme, orbit determination can be realized without relying on a ground observation station, and data is easy to manage.

Description

Technical Field

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

[0002] With the formation of low-orbit giant constellation systems, the demand for constellation measurement and control has exploded. At present, most low-orbit giant constellations rely on global distributed ground stations for continuous observation. Through the communication between the ground stations and satellites, accurate orbit data can be obtained. However, due to limited ground observation resources, it is impossible to continuously track each satellite, resulting in a significant decrease in orbit determination accuracy. The scale of giant constellations has exacerbated the severity of the scarcity of ground stations.

[0003] In order to solve the defects of relying on ground stations, many researchers have tried to use intersatellite links to obtain mutual measurements between satellites to obtain relative position and velocity information, thereby reducing dependence on ground stations. However, the traditional method of using intersatellite measurements faces problems such as constellation drift, and has little effect on improving the orbit determination accuracy and timeliness of giant constellations; on the other hand, the complexity of intersatellite communication and data processing also puts higher requirements on constellation management.

[0004] The defects of the existing technology can be attributed to the following two points:

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

[0006] Although intersatellite links are gradually developing, current orbit determination still relies heavily on observations from ground stations. This dependence is particularly evident in some areas due to the uneven distribution and coverage of ground observation resources. There are also external factors such as weak signal interference, which may cause orbit observation interruptions and affect accuracy.

[0007] 2. The complexity and scalability of intersatellite link management are poor

[0008] For giant constellations, the number of inter-satellite links is very large. As the number of satellites increases, the maintenance and management of the links becomes more complicated. In particular, when the relative motion between satellites is constantly changing, existing technical solutions are difficult to handle giant constellations. Existing technical solutions are also targeted at small constellations or single satellites, making it difficult to ensure the stability of giant constellation links and perform effective data management. Summary of the invention

[0009] In order to solve the above technical problems, the present invention proposes an autonomous orbit determination method for low-orbit giant constellations based on star sensors and inter-satellite angle measurement. It is an autonomous orbit determination method for low-orbit giant constellations that does not require any ground measurement stations at all, can achieve orbit determination without relying on ground measurement stations and is easy to manage data.

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

[0011] A method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement is carried out in the following steps:

[0012] Step 1: Determine the constellation configuration using inter-satellite measurement;

[0013] Step 2: Use star sensors and cameras to measure and correct constellation drift in combination with the constellation configuration in step 1.

[0014] Furthermore, the determination of the constellation configuration by using inter-satellite quantity measurement in step 1 comprises the following steps:

[0015] Step 1.1, use six numbers to express the state of each satellite in the constellation;

[0016] Step 1.2, construct the partial derivative matrix of the intersatellite quantity measurement relative to the state information, and construct a linear equation system according to the state quantity, quantity measurement and partial derivative matrix;

[0017] Step 1.3: Solve the satellite state variables describing the constellation configuration using the improved least squares method.

[0018] Furthermore, the six numbers used in step 1.1 to express the state of each satellite in the constellation are specifically:

[0019] m represents the total number of intersatellite measurements, n represents the total number of satellites in the constellation, D = [D1D2…D m ] T It represents the quantity measurement generated when satellites measure each other; D is related to the position of the satellite in the constellation; the state quantities of all satellites in the constellation are represented by six orbital numbers, namely a, e, i, Ω, ω, M, where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the longitude of the ascending node, ω is the pericenter angle, and M is the mean anomaly angle;

[0020] The status of each satellite in the constellation is represented by a vector form:

[0021]

[0022] In the formula, X represents a state quantity to be estimated, and X represents the vector composed of all state quantities.

[0023] Furthermore, in step 1.2, the partial derivative matrix of the intersatellite quantity measurement relative to the state information is constructed, and a linear equation system is constructed according to the state quantity, quantity measurement, and partial derivative matrix, specifically:

[0024]

[0025] Where R represents the partial derivative matrix of inter-satellite quantity measurement with respect to state information, and D represents the quantity measurement generated when satellites measure each other.

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

[0027] Using the LM method, we get the following solution formula:

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

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

[0030] Furthermore, the constellation drift is corrected by using the star sensor and camera to measure in step 2, combined with the constellation configuration in step 1, including the following steps:

[0031] Step 2.1, using a star sensor and camera, obtain the true value of the angle between the direction vector of an adjacent satellite and the direction vector of a certain star;

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

[0033] Step 2.3: Use an iterative method to minimize the difference between the true angle value and the estimated angle value to achieve constellation correction.

[0034] Furthermore, the star sensor and camera described in step 2.1 are used to obtain the true value of the angle between the adjacent satellite direction vector and a certain star direction vector. The specific contents are as follows:

[0035] Let B be the transformation matrix from the inertial system to the satellite body fixed coordinate system, and let A(α,β,γ)=C x (α)C y (β)C z (γ) is a coordinate rotation matrix to be solved. The real coordinates of a satellite in the constellation are column vector 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 vector AQ. The satellite coordinates obtained by intersatellite orbit determination are Q. When observing distant stars from an inertial system, it is assumed that the unit line of sight vector from any satellite to the distant star is a constant vector

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

[0037]

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

[0039]

[0040] 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 distant star in the inertial system, is the unit line of sight vector from satellite p to satellite q, B p represents the coordinate rotation matrix of the selected satellite from the local system to the inertial system, QP represents the difference between the coordinate vectors of the selected satellite and its adjacent satellites, θ true For L 0 and The angle is the true value of the angle between the adjacent satellite direction vector and a certain star direction vector.

[0041] Furthermore, in step 2.2, the estimated value of the angle between the adjacent satellite direction vector and a certain star direction vector is calculated according to the satellite state quantity obtained in step 1, specifically:

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

[0043]

[0044] Where cosθ esti It represents the estimated value of the angle between the direction vector of an adjacent satellite and the direction vector of a certain star. is the unit sight vector of the satellite pointing to the distant star in the inertial system, B p It represents the coordinate rotation matrix from the local system of the selected satellite to the inertial system, and QP represents the difference between the coordinate vectors of the selected satellite and its adjacent satellites.

[0045] Furthermore, the iterative method described in step 2.3 is used to minimize the difference between the true angle value and the estimated angle value to achieve constellation correction. The specific process is as follows:

[0046] Define the residual r i is cosθ esti,i and cosθ true,i The absolute value of the difference is:

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

[0048] The subscript i represents the i-th group of data, and the value range is 1, 2, ... N

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

[0050]

[0051] We want to minimize F(A), and this problem is solved using the Gauss-Newton method;

[0052] Define the initial θ0 and convert the residual r i Expanding at θ = θ0 and ignoring the second-order small terms, we get

[0053]

[0054] use represents the partial derivative of the residual with respect to each parameter,

[0055]

[0056] Taking the partial derivative of the above formula with respect to δθ and setting the result to 0, we get:

[0057]

[0058] in The calculation formula is:

[0059]

[0060] To consider multiple sets of data, construct the Jacobian matrix and stack all the residuals:

[0061]

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

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

[0064]

[0065] Using the iterative formula:

[0066]

[0067] Repeat multiple times until the preset convergence condition or accuracy requirement is reached, and α0, β0, γ0 are iterated to the actual rotation angle to obtain more accurate constellation orbit parameters.

[0068] The present invention also provides a low-orbit giant constellation autonomous orbit determination system based on star sensors and inter-satellite angle measurement, which is used to implement the above-mentioned low-orbit giant constellation autonomous orbit determination method based on star sensors and inter-satellite angle measurement.

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

[0070] 1. The present invention makes the most of inter-satellite measurement data and gets rid of the requirement for ground stations. Furthermore, the present invention uses star sensors and on-board measurement equipment instead of ground measurement data to eliminate the overall drift caused by using inter-satellite measurement data, and can achieve autonomous orbit determination that is completely free of external data.

[0071] 2. The present invention first uses the ranging information between satellites to perform relative orbit determination, and then uses the star sensor and intersatellite measurement equipment to rotate as a whole, making the data in the information fusion process easier to manage. At the same time, the accuracy is higher than the prior art when the error form is complex; when the measurement error form is normal noise, the performance of this scheme is the same as that of the prior art; and when the measurement error is in a more complex form, the orbit determination accuracy of this scheme is better than that of the prior art. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0074] The present invention innovatively proposes a method for autonomous orbit determination of giant constellations without the need for ground stations. It uses the ranging information between giant constellation satellites to perform relative orbit determination and determine the constellation configuration; then the constellation is regarded as a whole, and the star position information provided by the star sensor and the neighboring satellite position information provided by the camera are used to rotate the constellation as a whole, offsetting the rotation error caused by using only inter-satellite ranging, so as to achieve efficient and accurate giant constellation orbit calculation problems. This solution is particularly suitable for the scenario of giant constellation (6000+ satellites) orbit calculation (the scale is similar to the low-orbit Walker constellation or hybrid constellation).

[0075] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention are described in detail below.

[0076] A method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement is carried out in the following steps:

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

[0078] Six numbers are used to express the state quantities of each satellite in the constellation, and the partial derivative matrix of inter-satellite quantity measurement relative to the state information is constructed. A set of linear equations is constructed based on the state quantities, quantity measurements, and partial derivative matrices. The satellite state quantities describing the constellation configuration are solved using the improved least squares method.

[0079] Step 1.1: Use six numbers to express the state of each satellite in the constellation, specifically:

[0080] m represents the total number of intersatellite measurements, n represents the total number of satellites in the constellation, D = [D1D2…D m ] T It represents the quantity measurement generated when satellites measure each other; D is related to the position of the satellite in the constellation; the state quantities of all satellites in the constellation are represented by six orbital numbers, namely a, e, i, Ω, ω, and M, where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the longitude of the ascending node, ω is the periapsis angle, and M is the mean anomaly angle.

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

[0082]

[0083] In the formula, X represents a state quantity to be estimated, and X represents the vector composed of all state quantities.

[0084] Step 1.2: Construct the partial derivative matrix of intersatellite quantity measurement relative to state information, and construct a linear equation system based on state quantity, quantity measurement, and partial derivative matrix, specifically:

[0085]

[0086] Where R represents the partial derivative matrix of inter-satellite quantity measurement with respect to state information, and D represents the quantity measurement generated when satellites measure each other.

[0087] Step 1.3: Solve the satellite state quantity describing the constellation configuration by using the improved least squares method, specifically:

[0088] Using the LM method, we get the following solution formula:

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

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

[0091] If λE is not added, then R T R is not rank sufficient and cannot be inverted. The fact that the R matrix is ​​not rank sufficient means that the system of equations has infinite solutions, and adding a small amount allows one of the solutions to be obtained with higher accuracy.

[0092] It can be proved that the satellite positions obtained in sequence differ from the actual satellite positions by only one rotation around the center of the earth, and this rotation is time-invariant, thus obtaining the constellation configuration.

[0093] Step 2: Use star sensors and cameras to measure and correct constellation drift in combination with the constellation configuration in step 1.

[0094] Using star sensors and cameras, the true value of the angle between the direction vector of an adjacent satellite and the direction vector of a certain star is obtained. The estimated value of the angle between the direction vector of an adjacent satellite and the direction vector of a certain star is calculated based on the satellite state quantity obtained in step 1. The difference between the true value of the angle and the estimated value of the angle is minimized using an iterative method to achieve constellation correction.

[0095] Step 2.1: Use the star sensor and camera to obtain the true value of the angle between the adjacent satellite direction vector and a certain star direction vector. The specific contents are as follows:

[0096] Let B be the transformation matrix from the inertial system to the satellite body fixed coordinate system, and let A(α,β,γ)=C x (α)C y (β)C z (γ) is a coordinate rotation matrix to be solved. The real coordinates of a satellite in the constellation are column vector 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 vector AQ. The satellite coordinates obtained by intersatellite orbit determination are Q. When observing distant stars from an inertial system, it is assumed that the unit line of sight vector from any satellite to the distant star is a constant vector

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

[0098]

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

[0100]

[0101] 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 distant star in the inertial system, is the unit line of sight vector from satellite p to satellite q, B p represents the coordinate rotation matrix of the selected satellite from the local system to the inertial system, QP represents the difference between the coordinate vectors of the selected satellite and its adjacent satellites, θ true For L 0 and The angle is the true value of the angle between the adjacent satellite direction vector and a certain star direction vector.

[0102] Step 2.2: Calculate the estimated value of the angle between the adjacent satellite direction vector and a certain star direction vector based on the satellite state quantity obtained in step 1, specifically:

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

[0104]

[0105] Where cosθ esti It represents the estimated value of the angle between the direction vector of an adjacent satellite and the direction vector of a certain star. is the unit sight vector of the satellite pointing to the distant star in the inertial system, B p It represents the coordinate rotation matrix from the local system of the selected satellite to the inertial system, and QP represents the difference between the coordinate vectors of the selected satellite and its adjacent satellites.

[0106] Step 2.3: Use an iterative method to minimize the difference between the true angle value and the estimated angle value to achieve constellation correction. The specific process is as follows:

[0107] Define the residual r i is cosθ esti,i and cosθ true,i The absolute value of the difference is:

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

[0109] The subscript i represents the i-th group of data, and the value range is 1, 2, ... N

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

[0111]

[0112] We want to minimize F(A). This problem can be solved using the Gauss-Newton method.

[0113] Define the initial θ0 and convert the residual r i Expanding at θ = θ0 and ignoring the second-order small terms, we get

[0114]

[0115] use represents the partial derivative of the residual with respect to each parameter,

[0116]

[0117] Taking the partial derivative of the above formula with respect to δθ and setting the result to 0, we get:

[0118]

[0119] in The calculation formula is:

[0120]

[0121] To consider multiple sets of data, construct the Jacobian matrix and stack all the residuals:

[0122]

[0123] Write (17) and (18) in vector form:

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

[0125]

[0126] Using the iterative formula:

[0127]

[0128] Repeat multiple times until the preset convergence condition or accuracy requirement is reached, and α0, β0, γ0 are iterated to the actual rotation angle to obtain more accurate constellation orbit parameters.

[0129] Based on the same inventive concept, the present invention also provides a low-orbit giant constellation autonomous orbit determination system based on star sensors and inter-satellite angle measurement, which is used to implement the above-mentioned low-orbit giant constellation autonomous orbit determination method based on star sensors and inter-satellite angle measurement.

[0130] Simulation verification:

[0131] The following low-orbit giant constellation is used for simulation verification, see Table 1:

[0132] Table 1 Parameters of low-orbit giant constellation

[0133]

[0134] There are 1,296 satellites in the constellation, 16 of which are equipped with star sensors and the required inter-satellite angle measurement equipment.

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

[0136] Table 2 Geophysical constants

[0137]

[0138] The vector X composed of satellite states contains a total of 7776 unknowns. Each satellite has ranging data with the four satellites in front, behind, left and right. Each measurement time point generates 2592 sets of inter-satellite measurement data. Taking the partial derivative of the inter-satellite ranging vector at N time points with respect to the vector X, we get Using formula (3), the satellite state quantity describing the constellation configuration is calculated. The whole network position error (expressed as root mean square average rms) is 155m, and the whole network speed error (rms) is 0.17m / s. The 16 node satellites in the constellation measure the stars and the satellites in front of them through sensors, and obtain the cosθ described in formula (5) p,q ; Calculate the above situation with a position error of 155m. These 16 node satellites measure the stars and the satellites in front of them through sensors, and get the cosθ described in formula (6): esti . Minimize cosθ by iteration p,q and cosθ esti The difference between the two values ​​is obtained, and the changing curves of position error and velocity error are obtained as the number of iterations increases, as shown in Figure 1 and Figure 2 shown.

[0139] Figure 1 , Figure 2 They are respectively obtained by applying the method of the present invention, and are trend diagrams showing that the average position deviation of the entire network and the average speed deviation of the entire network decrease with the number of iterations. Figure 1 In the figure, the horizontal axis is the number of iterations, and the vertical axis is the square root average of the entire network position; Figure 2 In the figure, the horizontal axis is the number of iterations, and the vertical axis is the square root average of the speed of the entire network. Both indicators represent the orbit determination accuracy. The smaller the value, the higher the accuracy. As can be seen from the figure, during the 8 iterations, the rms of both the position and the speed of the entire network continued to decrease, and finally reached a very low value, thus proving the effectiveness and stability of the present invention in the precise orbit determination of giant constellations.

[0140] It should be noted that the embodiments of the present invention are preferred implementations rather than limiting the present invention. It should be noted that, for those skilled in the art, without departing from the principles of the present invention, the specific implementations may be modified or some technical features may be replaced by equivalents, which shall all be deemed to be within the scope of the present invention.

Claims

1. A method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement, characterized in that: Follow these steps: Step 1: Determine the constellation configuration using inter-satellite measurement; Step 2: Use star sensors and cameras to measure and correct constellation drift in combination with the constellation configuration in step 1.

2. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 1, characterized in that: Determining the constellation configuration by using inter-satellite quantity measurement in step 1 includes the following steps: Step 1.1, use six numbers to express the state of each satellite in the constellation; Step 1.2, construct the partial derivative matrix of the intersatellite quantity measurement relative to the state information, and construct a linear equation system according to the state quantity, quantity measurement and partial derivative matrix; Step 1.3: Solve the satellite state variables describing the constellation configuration using the improved least squares method.

3. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 2, characterized in that: The six numbers described in step 1.1 are used to express the state of each satellite in the constellation, specifically: m represents the total number of intersatellite measurements, n represents the total number of satellites in the constellation, D = [D1D2…D m ] T It represents the quantity measurement generated when satellites measure each other; D is related to the position of the satellite in the constellation; the state quantities of all satellites in the constellation are represented by six orbital numbers, namely a, e, i, Ω, ω, M, where a is the semi-major axis, e is the eccentricity, i is the orbital inclination, Ω is the longitude of the ascending node, ω is the pericenter angle, and M is the mean anomaly angle; The status of each satellite in the constellation is represented by a vector form: In the formula, X represents a state quantity to be estimated, and X represents the vector composed of all state quantities.

4. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 2, characterized in that: In step 1.2, the partial derivative matrix of intersatellite quantity measurement relative to state information is constructed, and a linear equation system is constructed according to the state quantity, quantity measurement, and partial derivative matrix, specifically: Where R represents the partial derivative matrix of inter-satellite quantity measurement with respect to state information, and D represents the quantity measurement generated when satellites measure each other.

5. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 2, characterized in that: The satellite state quantity describing the constellation configuration is solved by the improved least squares method described in step 1.3, specifically: Using the LM method, we get the following solution formula: X=(R T R+λE) -1 R T D(3) Wherein, the subscript T represents the matrix transpose, λ represents the damping factor, and E is the unit matrix.

6. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 1, characterized in that: The star sensor and camera measurements described in step 2 are combined with the constellation configuration in step 1 to correct the constellation drift, including the following steps: Step 2.1, using a star sensor and camera, obtain the true value of the angle between the direction vector of an adjacent satellite and the direction vector of a certain star; Step 2.2, calculating the estimated value of the angle between the adjacent satellite direction vector and a certain star direction vector according to the satellite state quantity obtained in step 1; Step 2.3: Use an iterative method to minimize the difference between the true angle value and the estimated angle value to achieve constellation correction.

7. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 6, characterized in that: As described in step 2.1, the star sensor and camera are used to obtain the true value of the angle between the adjacent satellite direction vector and a certain star direction vector. The specific contents are as follows: Let B be the transformation matrix from the inertial system to the satellite body fixed coordinate system, and let A(α,β,γ)=C x (α)C y (β)C z (γ) is a coordinate rotation matrix to be solved. The real coordinates of a satellite in the constellation are column vector 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 vector AQ. The satellite coordinates obtained by intersatellite orbit determination are Q. When observing distant stars from an inertial system, it is assumed that the unit line of sight vector from any satellite to the distant star is a constant vector Then the unit line of sight vector L of the star in the satellite body fixed coordinate system is 0 for: The calculation method for the true value of the angle between the adjacent satellite direction vector and a certain star direction vector is: 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 sight vector of the satellite pointing to the distant star in the inertial system, is the unit line of sight vector from satellite p to satellite q, B p represents the coordinate rotation matrix of the selected satellite from the local system to the inertial system, QP represents the difference between the coordinate vectors of the selected satellite and its adjacent satellites, θ true For L 0 and The angle is the true value of the angle between the adjacent satellite direction vector and a certain star direction vector.

8. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 6, characterized in that: In step 2.2, the estimated value of the angle between the adjacent satellite direction vector and a certain star direction vector is calculated based on the satellite state quantity obtained in step 1, specifically: Given the attitude of a satellite at time t and its coordinates (B, P, Q) with its neighboring satellites, the estimated angle between the direction vector of the neighboring satellite and the direction vector of a star is cosθ esti for: Where cosθ esti It represents the estimated value of the angle between the adjacent satellite direction vector and a certain star direction vector. is the unit sight vector of the satellite pointing to the distant star in the inertial system, B p It represents the coordinate rotation matrix from the local system of the selected satellite to the inertial system, and QP represents the difference between the coordinate vectors of the selected satellite and its adjacent satellites.

9. The method for autonomous orbit determination of a low-orbit giant constellation based on star sensors and inter-satellite angle measurement according to claim 6, characterized in that: In step 2.3, the iterative method is used to minimize the difference between the true angle value and the estimated angle value to achieve constellation correction. The specific process is as follows: Define the residual r i is cosθ esti,i and cosθ true,i The absolute value of the difference is: r i =|cosθ esti,i -cosθ true,i |,i=1,2,……N(7) The subscript i represents the i-th group of data, and the value range is 1, 2, ... N Define the objective function F(A) as: We want to minimize F(A), and this problem is solved using the Gauss-Newton method; Define the initial θ0 and convert the residual r i Expanding at θ = θ0 and ignoring the second-order small terms, we get use represents the partial derivative of the residual with respect to each parameter, Taking the partial derivative of the above formula with respect to δθ and setting the result to 0, we get: in The calculation formula is: To consider multiple sets of data, construct the Jacobian matrix and stack all the residuals: Write (17) and (18) in vector form: K T R=-(K T K)sth(19) Using the iterative formula: Repeat multiple times until the preset convergence condition or accuracy requirement is reached, and α0, β0, γ0 are iterated to the actual rotation angle to obtain more accurate constellation orbit parameters.

10. A low-orbit giant constellation autonomous orbit determination system based on star sensors and inter-satellite angle measurement, characterized in that: Used to implement the autonomous orbit determination method for low-orbit giant constellation based on star sensors and inter-satellite angle measurement as described in claims 1 to 9.

Images