Low earth orbit mega constellation orbit determination method based on vernal equinox point number

Abstract

The application discloses a low-orbit giant constellation orbit determination method based on vernal equinox point roots, and the method is used for Walker-delta constellation composed of thousands of satellites, under the perturbation model of non-spherical J2-J4 terms and MSIS00 atmospheric resistance, combined with inter-satellite / earth-space two types of ranging data, the vernal equinox point orbit roots, resistance coefficient b* and system error coefficient of all satellites are estimated at one time in the weighted least square framework; through two-stage iteration of 'orbit improvement-precise orbit determination', the orbit state is extrapolated by using RKF4(5) / RKF5(6) numerical integral, the non-zero elements of sparse Jacobian matrix are calculated based on the vernal equinox point root analytical formula, and then the conjugate gradient method is used to quickly solve the equation A T WAΔx=A T Wb, and the convergence of 8000-level constellation is completed in the general CPU platform 1300s, and the end error is better than 0.15km. The application can realize fast and precise orbit determination of giant constellation under limited ground control resources without additional hardware, and provides data support for space situation awareness, collision warning and accurate pointing of inter-satellite link.

Description

Technical Field

[0001] This invention relates to the field of low-Earth orbit giant constellation orbit determination technology based on vernal equinox roots, and particularly to a method for low-Earth orbit giant constellation orbit determination based on vernal equinox roots. Background Technology

[0002] As a crucial component of the future integrated land-sea-air-space network, emerging mega-constellations in low Earth orbit (LEO) have seen significant development in recent years due to their advantages such as wide coverage, low latency, and integrated space-ground capabilities. High-precision orbit determination for these mega-constellations has become a hot research topic. Among traditional orbit determination algorithms, dynamic orbit determination algorithms face a significant increase in computational load and complexity when dealing with large-scale LEO constellations, limiting their real-time performance. Geometric methods do not require precise initial satellite position information and have advantages for rapid orbit determination of LEO satellites, but they are susceptible to errors from the ionosphere and troposphere. Some new artificial intelligence algorithms learn and predict satellite orbit data based on neural networks, mining implicit information and patterns in the data to quickly and accurately estimate satellite orbital states and handle complex nonlinear relationships, but they require large amounts of training data and computational resources.

[0003] Inter-satellite links, typically established via microwave or laser, are used for inter-satellite ranging. Addressing the large-scale intelligent orbit calculation requirements of mega-constellations, this study comprehensively analyzes the perturbation mechanisms and technical characteristics of inter-satellite links. Orbit determination of mega-constellations can be achieved based on ranging data from inter-satellite links and satellite-to-ground links. Orbit determination using link ranging data is characterized by its massive scale, high accuracy requirements, and scarcity of ground-based measurement resources. Unlike traditional orbit determination based on radar and other tracking and control equipment, the satellite's orbit cannot be directly obtained from link ranging data; mathematically, it represents a typical inverse problem. Summary of the Invention

[0004] The purpose of this invention is to provide a method for determining the orbit of a low-Earth orbit mega-constellation based on the vernal equinox roots, thereby solving the aforementioned problems existing in the prior art.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] A method for determining the orbit of a large low-Earth orbit constellation based on the vernal equinox roots includes:

[0007] Step 1) Select the Walker-δ constellation configuration, with configuration code N / P / F, and obtain the initial vernal equinox orbital elements of all satellites in the constellation. The vernal equinox orbital elements are composed of the semi-major axis a, eccentric components h and k, inclination components p and q, and mean longitude λ.

[0008] Step 2) Establish an orbital dynamics model. The model is in the form of a set of differential equations, which only considers the non-spherical perturbation terms J2, J3, and J4 and the MSIS00 atmospheric drag term, and is fully decoupled from other calculation modules.

[0009] Step 3) Construct two types of observation models:

[0010] The ground-based ranging observation model acquires ranging data from the core node satellite and introduces system error coefficients a and b that drift linearly with time.

[0011] The inter-satellite ranging observation model acquires two-way ranging data from adjacent satellites in front and behind and left and right on the same orbit, and introduces linear drift over time and systematic error coefficients a, b, and c related to the satellite's z-coordinate;

[0012] Step 4) Embed the models from Step 2) and Step 3) together into a weighted least squares estimation framework. Use the vernal equinox orbital elements, drag coefficient b, and system error coefficients from Step 3) as the state vector to be estimated, and establish the observation equation y = H(X) + ε. The weight matrix W adds higher weights to the ground-based measurement data based on the principle that the inter-satellite ranging accuracy is higher than or equivalent to the satellite-to-ground ranging accuracy.

[0013] Step 5) Use a "one-time" orbit determination strategy to solve the problem in two iterative stages:

[0014] In the first stage, the drag coefficient b and the system error coefficient in step 3) are locked, and only the satellite orbital elements are estimated to obtain the improved orbit solution.

[0015] The second stage uses the improved orbit solution as the initial value, and jointly estimates the orbital elements, b* and the systematic error coefficient to obtain the precise orbit determination solution.

[0016] Step 6) In each iteration, the orbital state at the observation time is extrapolated using RKF4(5) or RKF5(6) numerical integration, and the non-zero elements corresponding to the inter-satellite / satellite-ground geometric correlation are calculated only based on the vernal equinox root analytical formula to obtain the sparse Jacobian matrix A of the orbital parameters of the observation pair.

[0017] Step 7) Transform the equation WAΔx=Wb into the normal equation A T WAΔx=A T Wb uses the conjugate gradient method to solve for the parameter correction Δx by multiplying only the non-zero elements, and updates the state vector X←X+Δx until rms satisfies the convergence condition;

[0018] Step 8) Output the converged orbital elements of the vernal equinox to complete the overall orbit determination of the giant constellation. Based on the final Jacobian matrix and residual vector, calculate the time consumption, initial error and final error to achieve a comprehensive evaluation of the orbit determination accuracy.

[0019] Preferably, the vernal equinox orbital elements in step 1) are used to eliminate the singularities of small inclination and small eccentricity through the following algebraic transformation:

[0020] .

[0021] Preferred,

[0022] The ground-based ranging observation model for step 3) is as follows: ;

[0023] The inter-satellite ranging observation model is as follows:

[0024] ;

[0025] In the formula, , c is the error coefficient. Zero-mean Gaussian white noise, To find the Euclidean distance function, i.e. , These are the coordinates of the survey station.

[0026] Preferably, in step 6), the sparse Jacobian matrix A only computes the following non-zero positions:

[0027] Inter-satellite ranging only takes partial derivatives with respect to the orbital parameters of adjacent satellites in the top, bottom, left, and right directions;

[0028] The satellite-to-ground ranging line only takes partial derivatives with respect to the orbital parameters of the corresponding node star;

[0029] The systematic error parameters correspond one-to-one with their respective observation data, and all other positions are set to zero;

[0030] In calculating A T WA and A T When Wb is used, only non-zero elements are traversed, and A is used. T The property that WA is a symmetric positive definite matrix is ​​that only the non-zero elements in the upper triangular region are calculated.

[0031] Preferably, in the conjugate gradient method solution process of step 7), sparse matrix multiplication skips zero-element operations through a pre-generated non-zero element index table and is implemented purely in software using a CPU general-purpose computing platform.

[0032] Preferably, in a scenario with 8,000 satellites at an orbital altitude of 550 km, the method achieves an error of ≤26.9 km at the end of the orbit improvement stage, an error of ≤0.15 km at the end of the precise orbit determination stage, and a total process time of ≤1336 s.

[0033] In some specific embodiments, a computer-readable storage medium stores a computer program that, when executed by a processor, implements the above-described method for determining the orbit of a low-Earth orbit mega-constellation based on the vernal equinox roots.

[0034] In some specific embodiments, a low-Earth orbit mega-constellation orbit determination system includes:

[0035] The orbital dynamics module is used to perform the integration of the differential equations in step 2) above;

[0036] The observation modeling module is used to perform the observation model construction in step 3) above;

[0037] The sparse Jacobian derivative module is used to perform the nonzero element calculations described above.

[0038] A conjugate gradient solver is used to solve the linear equations described above.

[0039] The accuracy assessment module is used to perform the orbit determination accuracy assessment in step 8).

[0040] The system couples the modules through a data bus to achieve a one-time orbit determination for the entire constellation.

[0041] The beneficial effects of this invention are:

[0042] This invention provides a method for determining the orbit of a large low-Earth orbit constellation based on the vernal equinox roots. Under limited ground measurement resources, it makes full use of inter-satellite measurement data to achieve rapid and accurate orbit determination of the large constellation. Attached Figure Description

[0043] Figure 1 This is a schematic diagram of the Walker-δ constellation of the present invention;

[0044] Figure 2 This is a flowchart of the single-stage differential improvement calculation of the present invention;

[0045] Figure 3 This is a diagram of the Jacobian matrix structure of the present invention;

[0046] Figure 4 This is a graph showing the relationship between the spacecraft's position and the vernal equinox roots according to the present invention;

[0047] Figure 5 This is a flowchart of the conjugate gradient method algorithm of the present invention;

[0048] Figure 6 This is a graph showing the changes in orbit improvement and precise orbit determination error of more than 3,000 satellites as a function of the number of iterations. Detailed Implementation

[0049] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0050] Reference Figures 1 to 6 The method shown is for determining the orbit of a low-Earth orbit giant constellation based on the vernal equinox roots, including:

[0051] Step 1) Select the Walker-δ constellation configuration, with configuration code N / P / F, and obtain the initial vernal equinox orbital elements of all satellites in the constellation. The vernal equinox orbital elements are composed of the semi-major axis a, eccentric components h and k, inclination components p and q, and mean longitude λ.

[0052] Step 2) Establish an orbital dynamics model. The model is in the form of a set of differential equations, which only considers the non-spherical perturbation terms J2, J3, and J4 and the MSIS00 atmospheric drag term, and is fully decoupled from other calculation modules.

[0053] Step 3) Construct two types of observation models:

[0054] The ground-based ranging observation model acquires ranging data from the core node satellite and introduces system error coefficients a and b that drift linearly with time.

[0055] The inter-satellite ranging observation model acquires two-way ranging data from adjacent satellites in front and behind and left and right on the same orbit, and introduces linear drift over time and systematic error coefficients a, b, and c related to the satellite's z-coordinate;

[0056] Step 4) Embed the models from Step 2) and Step 3) together into a weighted least squares estimation framework. Use the vernal equinox orbital elements, drag coefficient b, and system error coefficients from Step 3) as the state vector to be estimated, and establish the observation equation y = H(X) + ε. The weight matrix W adds higher weights to the ground-based measurement data based on the principle that the inter-satellite ranging accuracy is higher than or equivalent to the satellite-to-ground ranging accuracy.

[0057] Step 5) Use a "one-time" orbit determination strategy to solve the problem in two iterative stages:

[0058] In the first stage, the drag coefficient b and the system error coefficient in step 3) are locked, and only the satellite orbital elements are estimated to obtain the improved orbit solution.

[0059] The second stage uses the improved orbit solution as the initial value, and jointly estimates the orbital elements, b* and the systematic error coefficient to obtain the precise orbit determination solution.

[0060] Step 6) In each iteration, the orbital state at the observation time is extrapolated using RKF4(5) or RKF5(6) numerical integration, and the non-zero elements corresponding to the inter-satellite / satellite-ground geometric correlation are calculated only based on the vernal equinox root analytical formula to obtain the sparse Jacobian matrix A of the orbital parameters of the observation pair.

[0061] Step 7) Transform the equation WAΔx=Wb into the normal equation A T WAΔx=A T Wb uses the conjugate gradient method to solve for the parameter correction Δx by multiplying only the non-zero elements, and updates the state vector X←X+Δx until rms satisfies the convergence condition;

[0062] Step 8) Output the converged orbital elements of the vernal equinox to complete the overall orbit determination of the giant constellation. Based on the final Jacobian matrix and residual vector, calculate the time consumption, initial error and final error to achieve a comprehensive evaluation of the orbit determination accuracy.

[0063] This invention provides a method for orbit determination of a large low-Earth orbit constellation based on vernal equinox elements. Designed for the future Walker-δ constellation consisting of thousands of satellites, this method achieves rapid and precise orbit determination of the entire constellation using a purely software-based approach, even in environments with scarce ground-based telemetry and control resources but abundant inter-satellite link ranging data. The method first selects the constellation configuration code N / P / F according to mission requirements and obtains the initial vernal equinox orbit elements for all satellites. These elements are described by six parameters: semi-major axis a, eccentricity components h and k, inclination components p and q, and mean longitude λ. Algebraic transformations are used to completely eliminate the numerical singularities caused by small inclinations and eccentricities in traditional Kepler elements, providing a numerically stable basis for subsequent large-scale matrix operations.

[0064] Subsequently, an orbital dynamics model was established, containing only non-spherical perturbation terms J2, J3, and J4 and the MSIS00 atmospheric drag term, implemented as a system of differential equations. This model is fully decoupled from other computational modules, facilitating subsequent replacement or upgrading to higher-order perturbation models as needed. At the observation end, this invention simultaneously utilizes ground-based ranging data from the core node satellite and inter-satellite ranging data from adjacent satellites within the constellation to construct two types of observation models: the ground-based ranging model introduces system error coefficients a and b that drift linearly with time, while the inter-satellite ranging model further incorporates system error coefficients a, b, and c related to the satellite's z-coordinate. Each coefficient is treated as a constant within a single sampling arc, thus estimating the system error along with the orbital parameters, avoiding residual biases caused by prior corrections in traditional methods.

[0065] The aforementioned dynamics and observation models are jointly embedded into a weighted least squares estimation framework. The vernal equinox orbital elements, atmospheric drag coefficient b*, and two types of system error coefficients of all satellites in the constellation are used as the state vectors to be estimated. The observation equation y = H(X) + ε is established. The weight matrix W is based on the principle that "inter-satellite ranging accuracy is higher than or equivalent to satellite-to-ground ranging accuracy". It assigns higher weights to ground-based measurement data to ensure that the absolute positioning accuracy of the constellation as a whole relative to the Earth reference system can still be maintained when the geometric strength of inter-satellite data is insufficient.

[0066] To address the challenge of estimating ultra-high-dimensional nonlinearities caused by the coupling of thousands of satellites, this invention employs a "one-time" orbit determination strategy and designs a two-stage iterative process: The first stage locks the drag coefficient b* and all system error coefficients, estimating only the orbital elements to quickly obtain an improved orbital solution; the second stage uses this improved solution as initial values ​​to jointly estimate the orbital elements, b*, and system error coefficients, entering the precise orbit determination stage. This two-stage initial value transfer mechanism significantly reduces the degree of nonlinearity and the number of iterations.

[0067] Within each iteration, the orbital state at the extrapolated observation time is first propagated with high precision using RKF4(5) or RKF5(6) high-order adaptive numerical integration to obtain the satellite's position and velocity in the inertial frame. Then, based on the vernal equinox root analysis formula, only the non-zero partial derivative elements corresponding to the inter-satellite / satellite-to-ground geometric association are calculated to generate a sparse Jacobian matrix A in real time. The positions of its non-zero elements are predetermined by the Walker-δ constellation topology—the inter-satellite ranging row is only associated with the orbital parameters of the adjacent satellites above, below, left, and right, and the satellite-to-ground ranging row is only associated with the parameters of the corresponding node satellite. The system error parameters correspond one-to-one with their respective observation data, and all other positions are set to zero. This strategy reduces the Jacobian dimension of the 8000-satellite scale from tens of millions to millions of non-zero elements, with a memory footprint of less than 2KB per satellite.

[0068] After obtaining the sparse matrix A, the original equation WAΔx=Wb is transformed into the symmetric positive definite equation AᵀWAΔx=AᵀWb, which is solved using the conjugate gradient method. Multiplication operations only traverse the pre-generated non-zero element index table, avoiding invalid calculations of zero elements. Simultaneously, leveraging the symmetry of AᵀWA, only the upper triangular non-zero blocks are stored and computed, further reducing computational load. After iterative convergence, the vernal equinox orbital elements of all satellites in the constellation are output. Based on the final Jacobian matrix and residual vector, the initial error, termination error, and total time consumption are automatically calculated, enabling real-time assessment of orbit determination accuracy.

[0069] This invention can complete the orbit determination of a constellation of 8,000 satellites at an altitude of 550 km on a pure CPU general-purpose server. The error at the end of the orbit improvement phase is better than 27 km, and the error at the end of the precise orbit determination phase is better than 0.15 km, with a total time of approximately 1,300 seconds. Furthermore, the computation time increases approximately linearly with the size of the constellation, demonstrating excellent scalability. Since no additional hardware acceleration or constellation splitting / dimensional reduction is required, this invention, under limited ground-based telemetry and control resources, fully utilizes inter-satellite measurement data to achieve rapid and precise overall orbit determination of a giant constellation. This provides a high-precision orbital foundation for space situational awareness resource scheduling, on-orbit spacecraft collision early warning, millimeter-level pointing correction for inter-satellite links, and high-speed, high-capacity data transmission, significantly improving the operational safety and communication assurance capabilities of a large low-Earth orbit constellation.

[0070] Preferably, the vernal equinox orbital elements in step 1) are used to eliminate the singularities of small inclination and small eccentricity through the following algebraic transformation:

[0071] .

[0072] Preferred,

[0073] The ground-based ranging observation model for step 3) is as follows: ;

[0074] The inter-satellite ranging observation model is as follows:

[0075] ;

[0076] In the formula, , c is the error coefficient. Zero-mean Gaussian white noise, To find the Euclidean distance function, i.e. , These are the coordinates of the survey station.

[0077] Preferably, in step 6), the sparse Jacobian matrix A only computes the following non-zero positions:

[0078] Inter-satellite ranging only takes partial derivatives with respect to the orbital parameters of adjacent satellites in the top, bottom, left, and right directions;

[0079] The satellite-to-ground ranging line only takes partial derivatives with respect to the orbital parameters of the corresponding node star;

[0080] The systematic error parameters correspond one-to-one with their respective observation data, and all other positions are set to zero;

[0081] In calculating A T WA and A T When Wb is used, only non-zero elements are traversed, and A is used. T The property that WA is a symmetric positive definite matrix is ​​that only the non-zero elements in the upper triangular region are calculated.

[0082] Preferably, in the conjugate gradient method solution process of step 7), sparse matrix multiplication skips zero-element operations through a pre-generated non-zero element index table and is implemented purely in software using a CPU general-purpose computing platform.

[0083] Preferably, in a scenario with 8,000 satellites at an orbital altitude of 550 km, the method achieves an error of ≤26.9 km at the end of the orbit improvement stage, an error of ≤0.15 km at the end of the precise orbit determination stage, and a total process time of ≤1336 s.

[0084] In some specific embodiments, a computer-readable storage medium stores a computer program that, when executed by a processor, implements the above-described method for determining the orbit of a low-Earth orbit mega-constellation based on the vernal equinox roots.

[0085] In some specific embodiments, a low-Earth orbit mega-constellation orbit determination system includes:

[0086] The orbital dynamics module is used to perform the integration of the differential equations in step 2) above;

[0087] The observation modeling module is used to perform the observation model construction in step 3) above;

[0088] The sparse Jacobian derivative module is used to perform the nonzero element calculations described above.

[0089] A conjugate gradient solver is used to solve the linear equations described above.

[0090] The accuracy assessment module is used to perform the orbit determination accuracy assessment in step 8).

[0091] The system couples the modules through a data bus to achieve a one-time orbit determination for the entire constellation.

[0092] In some specific embodiments, the present invention provides a computer-readable storage medium, which can be any of the following forms: solid-state drive, hard disk drive, flash memory card, optical disk, or cloud-based distributed file system; or it can be a persistent storage unit embedded in a ground control server, onboard computer, or edge computing node. A compiled executable program is embedded on the storage medium. After the program code is called by the processor, it first reads the raw ranging data transmitted from the ground control network or inter-satellite link, and then sequentially completes the entire process, including Walker-δ constellation configuration analysis, vernal equinox orbital element initialization, numerical integration of non-spherical and atmospheric drag perturbations, construction of inter-satellite / ground-satellite observation models, non-zero element differentiation of sparse Jacobian matrices, conjugate gradient iteration of normal equations, and self-evaluation of orbit determination accuracy. During program execution, no manual intervention is required, enabling 24-hour near real-time orbit determination of up to 8000 satellites on a single general-purpose CPU server or containerized cluster. It also supports horizontal expansion by adding nodes to meet future constellation expansion needs.

[0093] In other specific embodiments, the present invention also provides a low-Earth orbit giant constellation orbit determination system. This system is designed based on the principles of modularity, low coupling, and easy expansion, and organically integrates various functional units through a high-speed data bus. The orbit dynamics module encapsulates the differential equations of the J2, J3, and J4 terms of non-spherical perturbation and the MSIS00 atmospheric drag model. It supports adaptive switching of the RKF4(5) or RKF5(6) higher-order integral algorithm according to the orbital altitude, and reserves interfaces for higher-order models such as solar radiation pressure, three-body gravity, and tidal perturbation, enabling hot-swappable models without interrupting business operations. The observation modeling module has dual-channel data access capabilities for both ground-based ranging and inter-satellite ranging. The ground-based channel receives raw pseudoranges from global tracking and control stations via TCP / UDP sockets, while the inter-satellite channel parses the inter-satellite ranging packets using the data frame format recommended by CCSDS. The module automatically completes time synchronization, coordinate system transformation, ionospheric / tropospheric delay correction, and initial values ​​of system error coefficients a, b, and c, providing a consistent observation benchmark for subsequent joint estimation.

[0094] The sparse Jacobian derivative module is based on the analytical formula for the vernal equinox roots. At runtime, it generates only pre-labeled non-zero elements: inter-satellite ranging rows are associated only with the six-root partial derivatives of their adjacent satellites (up, down, left, and right); satellite-to-ground ranging rows are associated only with the corresponding node's satellite parameters; and system error coefficients correspond one-to-one with their respective observations. The remaining positions are not allocated storage space in memory, thus compressing a theoretically tens of millions-dimensional dense matrix into a million-level non-zero element scale. Internally, the module employs the SIMD instruction set and OpenMP multi-threaded parallel technology, completing partial derivative calculations for 8000 satellites within seconds. It also supports dynamically loading new satellite nodes or new ranging links at runtime without requiring a system restart.

[0095] The conjugate gradient solver module implements the symmetric positive definite normal equation A T WAΔx=A T The iterative solution of Wb uses a core computational unit that only traverses a pre-generated non-zero element index table. Multiplication operations are performed using a cache-friendly block compressed sparse row (BCSR) format. It also integrates diagonal block preconditioners and residual soft thresholding, enabling the relative residual to be reduced to 10⁻ within 50 iterations. 4 The solver provides multi-language bindings for C / C++, Python, and Fortran, facilitating seamless integration with existing measurement and control software stacks.

[0096] After each orbit determination convergence, the accuracy assessment module automatically reads the final Jacobian matrix and residual vector, calculates the time series RMS of the difference between the predicted and actual measurements, the maximum absolute deviation, and the constellation's average position uncertainty, and outputs performance metrics such as CPU time, peak memory usage, and number of iterations. The assessment results are pushed to the monitoring center in JSON format, supporting real-time display on visualization platforms such as Grafana. They can also trigger downstream services such as collision warnings, orbit maintenance, and inter-satellite link pointing corrections via message queues.

[0097] The system's modules are interconnected via Gigabit Ethernet or InfiniBand high-speed data buses, and use ZeroMQ or gRPC for low-latency message passing. It supports both single-machine containerized deployment and horizontal scaling with Kubernetes clusters. The overall architecture has no single point of failure; any module failure can be resolved within seconds via hot standby nodes, ensuring uninterrupted 24 / 7 orbit determination service. Through this hardware and software co-design, the system can perform precise orbit determination for 8000 satellites at an orbital altitude of 550km within a 24-hour arc on a 28-core commercial server with only 128GB of memory. The entire process takes approximately 1300 seconds, with a completion error better than 0.15km. Furthermore, the computation time increases approximately linearly with the constellation size, providing an efficient, reliable, and low-cost orbit determination solution for the operation and control management of future mega-constellations with tens of thousands of satellites.

[0098] In another embodiment, please refer to Figure 1-6 This invention provides a method for orbit determination of giant constellations based on the vernal equinox roots. Step S1: Comprehensively analyze the perturbation dynamics of the giant constellation and construct an appropriate orbital model according to application requirements. The Walker-δ constellation, due to its unique configuration and good performance, is universally applicable, flexibly meeting diverse needs, and is stable and reliable. It is widely used in communication, navigation, Earth observation, and other fields. Therefore, this project assumes the giant constellation configuration is Walker-δ, with configuration codes N / P / F, such as... Figure 1 As shown;

[0099] The orbital dynamics model primarily considers factors such as non-spherical perturbations, solar radiation pressure, and atmospheric drag. This study considers non-spherical perturbations (including J2, J3, and J4 terms) and atmospheric drag perturbations (MSIS00 model). The dynamics model is presented as a system of differential equations and is fully decoupled from other computational modules to ensure its scalability. In practical applications, the required dynamics model can be flexibly replaced for upgrades. The orbital dynamics equations considering the J2, J3, and J4 non-spherical perturbations are as follows:

[0100] (1)

[0101] In the formula, GM: 3.96004415 × 10 14 ;

[0102] J2: -1.082626690598×10 -3 ;

[0103] J3: -2.532435345754×10 -6 ;

[0104] J4: -1.619331205072×10 -6 .

[0105] When a satellite flies in a thin atmosphere, the aerodynamic force it experiences is mainly drag, the magnitude of which is related to the drag coefficient. Cross-sectional area that withstands resistance and atmospheric density It is directly proportional to the square of the satellite's velocity relative to the atmosphere, that is:

[0106] (2)

[0107] (3)

[0108] in, For spacecraft mass; This represents the velocity vector of the spacecraft in the inertial frame. This represents the spacecraft's position vector in the inertial frame. This is the Earth's rotational angular velocity vector, with a magnitude of 7.292 × 10⁻⁶. -5 Atmospheric density (rad / s) The calculation was performed using the MSIS00 model.

[0109] Step S2: Construct a link ranging observation model and a ground-based measurement data observation model that reflect the actual situation.

[0110] The constellation link establishment method in this study is as follows: the core node satellite establishes a satellite-to-ground link, and all constellation satellites establish inter-satellite links with their neighboring satellites in the same direction of motion. Therefore, the measurement data consists of two types: one is the ranging data of the node satellites from the ground control equipment, and the other is the ranging data between adjacent satellites.

[0111] Analysis revealed errors in the foundation measurement data. For time A linear function, with measurement error being time. and two satellites A linear function of coordinates. Let the state of a single satellite be... ,in Let be the drag coefficient. For a single satellite, the ground-based measurement data observation model is:

[0112] (4)

[0113] In the formula, and The error coefficient, Zero-mean Gaussian white noise, To find the Euclidean distance function, i.e. , Here are the coordinates of the station. For inter-satellite measurement data, we have...

[0114] (5)

[0115] Taking into account the adaptability and computational complexity of the model, this project assumes that the coefficients of the systematic error model for satellite-to-ground measurement data and inter-satellite processing data do not change over time.

[0116] If a giant constellation is a composite constellation of multiple Walker constellations, then each constellation will be treated separately. Additionally, if the constellation is a non-standard Walker constellation, the satellites within each orbital plane are not uniformly distributed (but the semi-major axes of all satellites remain the same, and their flight periods are the same), but the satellites between different orbital planes are still roughly aligned, which does not affect the overall orbit determination.

[0117] Step S3: Construct a constellation orbit determination method framework based on least squares, transforming the link ranging data orbit determination problem into the optimal estimation problem of satellite state parameters;

[0118] When using the least squares method for orbit improvement, dimensional vector express There are 100 observations. The observation equation is:

[0119] (6)

[0120] in This is called a regression function. For state parameters, The noise is zero-mean Gaussian white noise. This project adopts a constellation-based orbit determination scheme, assuming the constellation includes... 1 satellite, including node satellite Ground measuring equipment The inter-satellite measurement data consists of distance data between adjacent satellites. At each sampling time, there are L inter-satellite distance data points. The ground state parameters are...

[0121] (7)

[0122] (8)

[0123] (9)

[0124] The orbit improvement problem is then equivalent to finding the parameters. The optimal estimate is:

[0125] (10)

[0126] Taking into account the different precision levels of the observation data and applying weighting, the above formula becomes:

[0127] (11)

[0128] in For state vectors, The initial value of the state vector; Let W be the initial time, and W be the observation precision. This represents the residual vector between the observed and extrapolated values. The weight matrix W significantly impacts the algorithm's convergence speed and accuracy; generally, a suitable value is selected based on the accuracy estimate of the observed data. In orbit determination for giant constellations, inter-satellite ranging is generally considered to have higher accuracy or be comparable to ground-based ranging. However, in early explorations, the inter-satellite ranging information for the Walker constellation was incomplete in orbit determination, requiring ground-based measurement data to provide an "anchor point" relative to Earth's reference. Therefore, a higher weight should be assigned to the ground-based measurement data.

[0129] Step S4: Giant Constellation Orbit Determination Steps

[0130] Compared to traditional single-satellite orbit determination problems, the interconnected and coupled satellites within a constellation present significant challenges to problem decomposition and dimensionality reduction. To avoid the accuracy loss caused by forced decomposition, a constellation-wide, one-time orbit determination scheme is adopted to solve the problem of optimal estimation of model parameters for a large-scale nonlinear system. Let there be n observation data points, and X represent m-dimensional orbital parameters and model parameters. Here, A represents the parameter correction, and A is the Jacobian matrix (the matrix of partial derivatives of the calculated state values ​​with respect to the initial orbital state values) in the reference orbit. The value at that location, i.e.

[0131] (12)

[0132] This is the derivative of the current measurement with respect to the orbital elements at the vernal equinox at the current moment. For state transition, the overall orbit determination of the constellation can be transformed into solving equations.

[0133] (13)

[0134] Orbit determination is divided into two stages. The first stage is to lock in the drag coefficient. The first stage involves estimating only the satellite orbital elements, along with the systematic error coefficients of each measurement value. This stage can be termed the "orbit improvement stage." The second stage involves jointly estimating the satellite orbital elements, drag coefficient, and systematic error coefficients of each measurement value. This stage can be termed the "precision orbit determination stage."

[0135] The improved differential solution process is as follows: Figure 2 As shown, Step 1: For each observation data point, based on the initial state... Extrapolate the orbital state at the observed time; calculate the residual matrix b between the observed state and the extrapolated calculated state; calculate the Jacobian matrix A. Step 2: Solve for parameter corrections. Step 1: Determine the descent direction and step size; Step 2: Check if the RMS meets the convergence condition; Step 3: Update the state. Step 5: Repeat the loop until convergence to obtain the orbit determination result. Among them, the extrapolation of the orbit state at the observation time in Step 1, the calculation of the Jacobian matrix, and the solution of the parameter correction in Step 2 are the key.

[0136] Step S5: Extrapolation of orbital state at the observation time

[0137] Orbit extrapolation can be performed using analytical formulas or numerical calculations. However, as is well known, except for the simplest two-body problem, the differential equations describing the motion of space targets have not yielded rigorous solutions, and the formulas for high-precision small-parameter power series solutions are cumbersome and complex. Therefore, this study uses numerical methods for orbit extrapolation. The most commonly used numerical integration method for differential equations is the Runge-Kutta-Fehlberg method (RKF). Specifically, it selects either the higher-order adaptive RKF4 (5) or RKF5 (6) algorithm based on the orbital altitude and data. The RKF4 (5) method provides a fourth-order precision reference solution while ensuring fifth-order precision, supporting error estimation and adaptive step size control; the RKF5 (6) method provides higher sixth-order precision and achieves more accurate error control through a seventh-order reference solution, meeting the requirements for high-precision orbit numerical integration.

[0138] Differential equations with known initial values:

[0139] (14)

[0140] In the formula, For state vectors, The initial value of the state vector, Let's take the initial time as an example. Using RKF5(6), the integration step size is... Then there is

[0141] (15)

[0142] (16)

[0143] No. The local truncation error of the step is

[0144] (17)

[0145] Step S6: Calculation of the Jacobian matrix

[0146] The Jacobian matrix is ​​the partial derivative matrix of the observations with respect to orbital parameters (such as orbital elements, perturbation parameters, and systematic errors). It reflects the degree of influence of small changes in orbital parameters on the observations and is the core of constructing the state transition equations and error equations. Its accuracy directly affects the accuracy of parameter estimation and the convergence speed of the algorithm. Within the overall orbit determination framework of this project, the structure of the matrix is ​​as follows: Figure 3 As shown;

[0147] As mentioned earlier, orbit determination is divided into two stages: "orbit improvement" and "precise orbit determination". To improve the speed and convergence of iterative calculations, the initial value of the derivative calculation in the second stage is set to the result of the calculation in the first stage.

[0148] Observations are the direct basis for orbit determination, and errors in their calculations directly lead to residual distortion. The Jacobian matrix, as a linearization tool, can be iteratively corrected as long as it approximates the influence of parameters on observations. Analytical calculation of the Jacobian matrix requires differentiating complex dynamic models (e.g., calculating partial velocity derivatives of atmospheric drag models) or using numerical differencing methods (e.g., recalculating observations after small parameter perturbations), and the computational load increases dramatically with the dimensionality of orbital parameters. Therefore, considering the role of the Jacobian matrix in orbit determination calculations... The computational accuracy requirement is relatively low. To improve the speed of numerical differentiation, the derivative is calculated using the state transition calculation and the approximate analytical formula of the observation model, thereby increasing the speed. As mentioned above, the Jacobian matrix can be calculated according to formulas (18) to (23) and (37) to (39). Relevant elements. Due to the locality of satellite correlation, the Jacobian matrix is ​​a large-scale sparse matrix. To avoid unnecessary calculations, it is necessary to fully analyze the distribution of non-zero elements and solve for the derivative using numerical methods.

[0149] Step S7: Fast sparse matrix multiplication based on nonzero element analysis

[0150] Research has revealed that in the improved differential algorithm, the massive scale of giant constellations leads to an excessive Jacobian matrix. The dimension is very high, the computational cost is large, and the normal equation matrix is... and It is also a high-dimensional matrix, involving a large amount of computation. Due to the locality of the dependence of measured values ​​on different initial orbital values ​​and system parameters, matrix, Matrix and The matrix is ​​a sparse matrix with a large number of non-zero elements. When calculating sparse matrix multiplication, directly calling the sparse matrix multiplication algorithm will result in trial calculations of all terms due to the lack of prior information. This project aims to fully analyze the problem structure, establish a network model of the relationships between various parameters, deduce and understand the positions of non-zero elements in the sparse matrix, and eliminate a large number of invalid calculations.

[0151] Calculate the Jacobian matrix The derivatives in the equation are used to solve only the non-zero elements, as shown in Table 1.

[0152] Table 1. Acceleration of Jacobian matrix differentiation

[0153]

[0154] As mentioned above, It is also a sparse matrix, and it is easy to further conclude that this matrix is ​​a symmetric positive definite matrix. Therefore, the calculation of this matrix can be performed by only multiplying the non-zero elements of the upper triangular matrix, as shown in Table 2.

[0155] Table 2 Sparse matrix A T A Multiplication Acceleration

[0156]

[0157] The calculation is similar, because The sparsity can be calculated by simply computing the sparsity of the data. Non-zero elements and The product of these factors significantly reduces redundant calculations.

[0158] Step S8: Calculate the derivative of the position quantity with respect to the orbital elements at the vernal equinox.

[0159] To avoid numerical singularities in traditional Kepler elements for small-inclination and small-eccentricity orbits when describing spacecraft motion, and considering both numerical solution efficiency and stability, this project uses the vernal equinox orbital elements to characterize the orbital state. The vernal equinox orbital elements include parameters such as the semi-major axis, two eccentric components, two inclination components, and the mean anomaly angle. The relationship between spacecraft position quantities and the vernal equinox elements is as follows: Figure 4 As shown, the relationship between the vernal equinox orbital elements and the Keplerian orbital elements is as follows:

[0160] (18)

[0161] in Longitude is average. Longitude deviation is defined as follows:

[0162] (19)

[0163] The following steps can be used to determine the position and velocity of a spacecraft based on the vernal equinox elements:

[0164] (20)

[0165] but Location Roots of the vernal equinox The following relationships exist (subscripts omitted). ):

[0166] (twenty one)

[0167] ① right The derivative is

[0168] (twenty two)

[0169] ② right The derivative is

[0170] (twenty three)

[0171] ③ right The derivative is

[0172] (twenty four)

[0173] ④ right The derivative is

[0174] (25)

[0175] ⑤ right The derivative is

[0176] (26)

[0177] ⑥ right The derivative is

[0178] (27)

[0179] Step S9: Vernal Equinox Root State Transition Calculation

[0180] The calculation of the state transition of the roots at the vernal equinox needs to be calculated. The derivative of the vernal equinox roots with respect to the initial vernal equinox roots. Considering the J2 perturbation, we have

[0181] (28)

[0182] but

[0183] (29) (30) (31) (32)

[0184] (33)

[0185] for

[0186] (34)

[0187] Therefore, there is

[0188] (35)

[0189] because ,therefore

[0190] (36)

[0191] (37)

[0192] (38)

[0193] (39)

[0194] (40)

[0195] Therefore, for each satellite, the partial derivatives of the roots at the vernal equinox are arranged into a matrix as follows: (41)

[0196] in .

[0197] Step S10: The derivative of the observed quantity with respect to the system error coefficients

[0198] Inter-satellite measurement data depends only on the position of the corresponding satellite; for each measurement data... ,have:

[0199] (42)

[0200] in Satellite-to-ground measurement data depends only on the positions of the satellite and the corresponding measurement station; for each measurement data point... have:

[0201] (43)

[0202] in .

[0203] As described above, the Jacobian matrix can be calculated using formulas (22) to (27) and formulas (41) to (43). Related elements.

[0204] Step S11: Parameter correction amount Calculation

[0205] Equation First convert to normal equations Then, the conjugate gradient method is used to solve the equation. The algorithm flow of the conjugate gradient method is as follows: Figure 5 As shown. The conjugate gradient method is an iterative method for solving symmetric positive definite linear equation systems. For the normal equations... Weight Since the value is fixed, the equation can be transformed into because Since the matrix is ​​symmetric positive definite, the conjugate gradient method can be used to solve it. The basic idea of ​​the conjugate gradient method is to construct a set of conjugate vectors to progressively approximate the solution to the system of equations. Each iteration searches along the conjugate direction, ensuring that the projection of the residual vector onto that direction is zero. Compared to traditional iterative methods, the conjugate gradient method has the advantages of fast convergence and low storage requirements. However, the conjugate gradient method is sensitive to initial values ​​and may have a slower convergence speed for ill-conditioned problems.

[0206] Step S12: Determine the final Jacobian matrix and residuals based on the orbit, and conduct a comprehensive evaluation of the orbit determination accuracy of the giant constellation.

[0207] This project adopts a process of first improving the orbit and then determining the precise orbit, with both processes based on the conjugate gradient method. Traditionally, the initial solution uses guessed values. To improve algorithm stability and convergence speed, the initial orbit data is used as the initial value in the orbit improvement step, and the orbit data after the orbit improvement calculation is used as the initial value in the precise orbit determination step. Therefore, the algorithm consists of several epochs (diff_correction_epoch), each epoch containing several steps. Based on measurement data, orbit determination calculations were performed using pure CPU on constellations of 1152, 3600, 5760, and 8100 stars (the algorithm verification computer processor was: Intel® Core™ Ultra7-155H 1.40GHz, 22 cores, 32GB memory; the competition computer processor was: Intel® Core™ i7-14700 2.10GHz, 28 cores, 128GB memory). The algorithm iteration convergence results are shown in Table 3 and... Figure 6 As shown in the table. Due to a lack of standard test data, the data for 8100 satellites and the 300km range were generated by a self-developed simulation program. In the table, the initial error and final error refer to the difference between the predicted and actual measurements in the differential correction algorithm. The initial error is the error at the end of the first iteration, and the final error is the error at the end of the last iteration. The algorithm itself has no special requirements on the constellation size; based on testing, the orbit determination algorithm can be scaled up to 8000+ satellites.

[0208] Table 3. Algorithm Iteration Convergence Status

[0209]

[0210] By adopting the above-disclosed technical solution of this invention, the following beneficial effects are obtained:

[0211] Eliminating singularities: The orbital state is represented by the orbital elements at the vernal equinox. Through algebraic transformations, numerical singularities caused by small inclination angles and small eccentricities are completely avoided, thereby improving the overall iterative stability of large-scale constellations.

[0212] Computational complexity is controllable: by utilizing the locality of inter-satellite / satellite-ground geometric correlations, only non-zero elements of the Jacobian matrix are calculated, and only upper triangular non-zero elements are stored through the symmetric positive definite property, resulting in a memory footprint of <2KB / satellite, which is 90% less memory usage for the same scale.

[0213] The project is easy to implement: the dynamic model, observation model, and sparse matrix module are fully decoupled and can be upgraded and replaced independently, making it convenient to connect to higher-order perturbation models such as solar radiation pressure and attitude coupling.

[0214] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for determining the orbit of a low-Earth orbit giant constellation based on the vernal equinox roots, characterized in that, include: Step 1) Select the Walker-δ constellation configuration, with configuration code N / P / F, and obtain the initial vernal equinox orbital elements of all satellites in the constellation. The vernal equinox orbital elements are composed of the semi-major axis a, eccentric components h and k, inclination components p and q, and mean longitude λ. Step 2) Establish an orbital dynamics model, which is a set of differential equations that only considers the non-spherical perturbation terms J2, J3, and J4 and the MSIS00 atmospheric drag term, and is fully decoupled from other calculation modules; Step 3) Construct two types of observation models: A ground-based ranging observation model is used to acquire ranging data from the core node satellite and introduces system error coefficients a and b that drift linearly with time. An inter-satellite ranging observation model is used to acquire two-way ranging data from adjacent satellites in front and behind and left and right on the same orbit, and introduces linear drift over time and systematic error coefficients a, b, and c related to the satellite's z-coordinate. Step 4) Embed the models from Step 2) and Step 3) together into a weighted least squares estimation framework. Use the vernal equinox orbital elements, drag coefficient b, and the system error coefficients described in Step 3) as the state vector to be estimated, and establish the observation equation y = H(X) + ε. The weight matrix W adds higher weights to the ground-based measurement data based on the principle that the inter-satellite ranging accuracy is higher than or equivalent to the satellite-to-ground ranging accuracy. Step 5) Use a "one-time" orbit determination strategy to solve the problem in two iterative stages: In the first stage, the drag coefficient b and the system error coefficient mentioned in step 3) are locked, and only the satellite orbital elements are estimated to obtain the improved orbit solution. The second stage uses the improved orbit solution as the initial value, and jointly estimates the orbital elements, b* and the systematic error coefficient to obtain the precise orbit determination solution. Step 6) In each iteration, the orbital state at the observation time is extrapolated using RKF4(5) or RKF5(6) numerical integration, and the non-zero elements corresponding to the inter-satellite / satellite-ground geometric correlation are calculated only based on the vernal equinox root analytical formula to obtain the sparse Jacobian matrix A of the orbital parameters of the observation pair. Step 7) Transform the equation WAΔx=Wb into the normal equation A T WAΔx=A T Wb uses the conjugate gradient method to solve for the parameter correction Δx by multiplying only the non-zero elements, and updates the state vector X←X+Δx until rms satisfies the convergence condition; Step 8) Output the converged orbital elements of the vernal equinox to complete the overall orbit determination of the giant constellation. Based on the final Jacobian matrix and residual vector, calculate the time consumption, initial error and final error to achieve a comprehensive evaluation of the orbit determination accuracy.

2. The method according to claim 1, characterized in that, The vernal equinox orbital elements obtained in step 1) are used to eliminate the singularities of small inclination and small eccentricity through the following algebraic transformation: 。 3. The method according to claim 1, characterized in that, The ground-based ranging observation model for step 3) is as follows: ; The inter-satellite ranging observation model is as follows: ； In the formula, , c is the error coefficient. Zero-mean Gaussian white noise, To find the Euclidean distance function, i.e. , These are the coordinates of the survey station.

4. The method according to claim 1, characterized in that, The sparse Jacobian matrix A in step 6) only computes the following non-zero positions: Inter-satellite ranging only takes partial derivatives with respect to the orbital parameters of adjacent satellites in the top, bottom, left, and right directions; The satellite-to-ground ranging line only takes partial derivatives with respect to the orbital parameters of the corresponding node star; The systematic error parameters correspond one-to-one with their respective observation data, and all other positions are set to zero; In calculating A T WA and A T When Wb is used, only the non-zero elements are traversed, and A is used. T The property that WA is a symmetric positive definite matrix is ​​that only the non-zero elements in the upper triangular region are calculated.

5. The method according to claim 1, characterized in that, In the conjugate gradient method solution process of step 7), sparse matrix multiplication skips zero-element operations through a pre-generated non-zero element index table and is implemented purely in software using a CPU general-purpose computing platform.

6. The method according to claim 5, characterized in that, In a scenario with 8,000 satellites at an orbital altitude of 550 km, the method achieves an error of ≤26.9 km at the end of the orbit improvement stage, ≤0.15 km at the end of the precise orbit determination stage, and a total process time of ≤1336 seconds.

7. A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for determining the orbit of a low-Earth orbit mega-constellation based on the vernal equinox roots as described in any one of claims 1-6.

8. A low-Earth orbit mega-constellation orbit determination system, comprising: The orbital dynamics module is used to perform the integration of the differential equation in step 2 of claim 1; An observation modeling module is used to perform the observation model construction in step 3) of claim 1. A sparse Jacobian derivative module is used to perform the nonzero element computation of claim 4; A conjugate gradient solver is used to solve the linear equation system of claim 5; The accuracy assessment module is used to perform the orbit determination accuracy assessment in step 8) of claim 1. The system couples the modules through a data bus to achieve a one-time orbit determination for the entire constellation.

Images