Fast orbit determination system based on global large-scale GNSS station network data

Abstract

The application relates to the technical field of satellite orbit calculation, and discloses a fast orbit determination system based on global large-scale GNSS station network data, which comprises a precise orbit determination execution module and three-layer acceleration modules, the three-layer acceleration modules are embedded into multi-round iterations of the precise orbit determination execution module; the precise orbit determination execution module is used for acquiring observation data of a global tracking station network and finally outputs multi-system precise orbit and satellite clock difference products; the three-layer acceleration modules comprise a process-level parallel acceleration module, a thread-level parallel acceleration module and an algorithm-level block elimination optimization module. The fast orbit determination system based on global large-scale GNSS station network data shortens the calculation time consumption while ensuring the centimeter-level orbit determination precision through the collaborative design of the process-level, thread-level and algorithm-level three-layer parallel acceleration architecture.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation, positioning, and orbit determination data processing technology, and specifically to a rapid orbit determination system based on global large-scale GNSS station network data. Background Technology

[0002] Precise orbit determination for Global Navigation Satellite System (GNSS) is fundamental to satellite navigation applications. The accuracy and timeliness of satellite orbit and clock bias products directly determine the quality of user-end positioning, navigation, and timing services. With the comprehensive deployment of multiple systems such as GPS, GLONASS, Galileo, and BeiDou, the number of navigation satellites in orbit has now exceeded 120, and the global GNSS tracking station network has expanded from dozens to hundreds, with the volume and types of observation data continuously increasing. Against this backdrop, how to efficiently process massive amounts of observation data while ensuring solution accuracy and stability, and achieve high-precision and rapid orbit determination for multiple navigation satellite systems, has become a significant challenge in the field of GNSS data processing.

[0003] Existing GNSS precise orbit determination methods are mainly based on batch least squares estimation techniques. A typical processing flow includes observation data preprocessing, orbit numerical integration, orbit fitting, parameter estimation, residual editing, and ambiguity fixing. When processing data from over 200 global tracking stations, more than 100 satellites across multiple systems, and 24-hour arc data with a sampling interval of 30 seconds, the number of parameters to be estimated can reach hundreds of thousands, and the dimension of the normal equations can rise to hundreds of thousands of orders. Constrained by the rapid expansion of station network size and constellation numbers, traditional batch orbit determination processes typically suffer from the following three levels of efficiency bottlenecks: (1) Insufficient process-level parallelism: There are a large number of independent tasks in the orbit determination process that can be split at the satellite or station level, but traditional implementations often use serial or low-parallel processing methods, which make it difficult to make full use of distributed computing resources, and the overall processing efficiency is easily limited by the computing and storage capabilities of a single node. (2) Insufficient use of thread-level parallelism: The orbit fitting and parameter estimation stages require the construction and solution of large-scale matrices, which is the most computationally intensive core link. Traditional software often does not give full play to the thread parallelism capabilities of multi-core processors and high-performance linear algebra libraries, resulting in low efficiency in operations such as normal equation accumulation, decomposition and solution and residual calculation. (3) The algorithm-level parameter elimination mechanism is imperfect: the orbit determination solution contains a large number of local parameters (such as ambiguity) that increase with the advancement of the arc segment. Traditional parameter elimination or uniform retention solution strategies are prone to increase in computational complexity and matrix dimension expansion, thus limiting the space for overall efficiency improvement.

[0004] For example, the International GNSS Service (IGS) Analysis Center typically takes several hours or even longer to process a day's worth of global multi-system data, which severely restricts the timeliness of precision orbit products.

[0005] To address the aforementioned issues, existing research has attempted to introduce thread parallelism or distributed processing strategies in localized stages. However, these approaches often only cover a single level or accelerate only some modules, making it difficult to simultaneously achieve distributed task parallelism, operator-level thread acceleration, and parameter elimination algorithm optimization at the system level. Consequently, the overall acceleration effect remains limited.

[0006] In summary, existing GNSS precise orbit determination methods have not yet systematically integrated acceleration technologies at three levels—process-level parallelism, thread-level parallelism, and algorithm-level optimization—when dealing with large-scale global GNSS station network data processing. Therefore, there is an urgent need for a multi-level collaborative parallel rapid orbit determination system to meet the practical application requirements of efficient processing of large-scale global GNSS station network data. Summary of the Invention

[0007] The purpose of this invention is to address the shortcomings of the aforementioned technologies by providing a rapid orbit determination system based on large-scale global GNSS station network data. Through the collaborative design of a three-layer parallel acceleration architecture at the process, thread, and algorithm levels, the system significantly reduces computation time while maintaining centimeter-level orbit determination accuracy.

[0008] To achieve the above objectives, the present invention relates to a rapid orbit determination system based on global large-scale GNSS station network data, comprising a precise orbit determination execution module and a three-layer acceleration module, wherein the three-layer acceleration module is embedded in the multiple iterations of the precise orbit determination execution module; The precision orbit determination execution module is used to acquire observation data from the global tracking station network, complete the preprocessing and quality control of the observation data, orbit numerical integration and orbit fitting, use the least squares parameter estimation method to solve the equations, and combine residual editing and iterative convergence control to iteratively update the least squares parameter estimation, orbit parameter update and residual check. After the iteration is completed, the ambiguity is fixed and re-estimated, and finally the multi-system precision orbit and satellite clock error products are output. The three-layer acceleration module includes: The process-level parallel acceleration module is used to execute computational tasks that can be split into satellite-dimensional or station-dimensional components serving the precise orbit determination execution module using the MPI distributed parallel mechanism. The computational tasks include orbital numerical integration, observation data preprocessing and quality control, as well as frequency code deviation calculation and station load displacement calculation that provide correction information to the precise orbit determination execution module. The thread-level parallel acceleration module is used to implement operator-level parallel acceleration for the construction of normal equations, Cholesky decomposition, residual calculation and matrix operations involved in trajectory fitting, least squares parameter estimation and block elimination. The algorithm-level block elimination optimization module is used to dynamically eliminate inactive local parameters during the least squares parameter estimation process using a block matrix elimination strategy based on Schur complement, and saves the back-substitution recovery information to control the dimension of the normal equation.

[0009] Preferably, the process-level parallel acceleration module includes: MPI Process Initialization and Task Allocation Unit: Used to initialize the MPI computing environment and start... Each MPI process sets up the overall computing task, including... Each subtask is assigned to an MPI process in a round-robin fashion. The round-robin strategy ensures that the number of subtasks assigned to each MPI process differs by at most one, thus guaranteeing an approximately balanced load distribution. Parallel computation unit for satellite orbit numerical integration: Used by each MPI process to independently perform orbit numerical integration on the satellite subtask assigned to it, based on the task allocation results; Parallel computing unit for preprocessing station observation data: Based on the task allocation results, each MPI process independently performs GNSS observation data preprocessing and quality control on the station sub-tasks assigned to it. The observation data preprocessing process of each station only involves the observation files of that station. The stations are independent of each other and there is no data exchange between processes. Intra-frequency code deviation parallel calculation unit: used to independently execute the intra-frequency code deviation differential observation construction for each station subtask assigned to it based on the task allocation results; Parallel calculation unit for station load displacement: Based on the task allocation results, each MPI process independently performs load deformation displacement correction calculations for the station subtasks assigned to it; Process synchronization and result aggregation unit: After the satellite orbit numerical integration parallel calculation unit, station observation data preprocessing parallel calculation unit, intra-frequency code deviation parallel calculation unit and station load displacement parallel calculation unit have been executed, process synchronization is performed, and an MPI process is selected as the main process. The local result files generated by each MPI process are aggregated to form a complete result file for subsequent orbit determination calculation.

[0010] Preferably, the MPI process initialization and task allocation unit: Startup Each MPI process obtains the total number of MPI processes in the global communication domain. With its own MPI process number r Set up the overall computing task, including Each subtask is assigned sequentially to an MPI process using a round-robin strategy. r The set of task indexes under your responsibility Represented as: , in, r Indicates the MPI process number. , This represents the total number of subtasks in the current overall computing task, when broken down by satellite dimension. When splitting by station dimension , This indicates the total number of satellites participating in orbit determination. This indicates the total number of stations involved in orbit determination.

[0011] Preferably, the satellite orbit numerical integration parallel computing unit is used by the MPI process according to the task allocation results. r For task index set Satellites assigned to their own sub-missions independently perform orbital numerical integration; for those assigned to the MPI process... satellite sub-mission Its six-dimensional state vector The equation of motion that is satisfied is expressed as: , in, Indicates satellite s The three-dimensional position vector in the geocentric inertial coordinate system. This represents the corresponding three-dimensional velocity vector. t For time, F This represents a nonlinear dynamic function composed of several perturbation forces, including Earth's gravitational field, the gravitational pull of the Sun and Moon, solar radiation pressure, Earth's tidal perturbations, and relativistic effects. Since the equations of motion for each satellite depend only on its own state vector... Different satellites are independent of each other, and there is no need for inter-process data exchange during the integration process of each MPI process.

[0012] Preferably, the intra-frequency code deviation parallel calculation unit is used by the MPI process according to the task allocation result. r For task index set Each station subtask assigned to it independently reads pseudorange observation data and constructs pseudorange differential observations between different signal attributes on the same carrier frequency; for those assigned to the MPI process... r Station subtask The observed satellite s in epoch The differential observation equation for its intra-frequency code bias is expressed as: , in, and Representing the stations j In the calendar Received satellite s In frequency f The above are signal attributes and Encoded pseudorange observations, Indicates the satellite end-frequency code deviation. Indicates the intra-frequency code deviation at the receiver end. This represents observation noise. Each MPI process independently performs data reading and differential observation construction for the station subtasks assigned to it, and different stations are independent of each other.

[0013] Preferably, the parallel calculation unit for station load displacement is used by the MPI process based on the task allocation results. r For task index set The subtasks assigned to each station within the system independently perform load deformation displacement correction calculations, while those assigned to the MPI process... r Station subtask In the calendar The total load displacement correction vector is expressed as: , in, This indicates the displacement caused by ocean tidal loads. This indicates the displacement caused by the polar tide load. This indicates the displacement caused by atmospheric tidal loads. This represents displacement caused by non-tidal atmospheric and ocean loads. The load displacement correction of each station depends only on its own coordinates and external geophysical model parameters. The stations are independent of each other and are executed in complete parallel.

[0014] Preferably, the thread-level parallel acceleration module includes: Multi-threaded environment initialization unit: used to set the number of parallel threads and enable the multi-threaded linear algebra subroutine during the orbit fitting and least squares parameter estimation stages; A multi-threaded normal equation construction unit is used to call the multi-threaded linear algebra subroutine to perform thread-level parallel acceleration of matrix operations in the normal equation construction process under the number of parallel threads, so as to construct the normal equation coefficient matrix. Multi-threaded Cholesky decomposition solution unit: used to call the multi-threaded linear algebra subroutine to perform Cholesky decomposition of the coefficient matrix of the normal equation and subsequent previous and back substitution solutions under the number of parallel threads; Multi-threaded residual calculation unit: used to call the multi-threaded linear algebra subroutine to perform thread-level parallel acceleration of matrix-vector operations in the residual recovery process under the number of parallel threads, so as to recover the observation residuals.

[0015] Preferably, the multi-threaded environment initialization unit is used to set the number of parallel threads during the orbit fitting and least squares parameter estimation stages. This allows subsequent linear algebra subroutines involving matrix operations to be executed within the process. The linear algebra subroutines are executed concurrently by multiple threads. The linear algebra subroutines are provided by the multi-threaded basic linear algebra subroutine library, covering subroutines for three-level matrix-matrix operations, two-level matrix-vector operations, and decomposition and solution of linear equation systems. Multi-threaded method equation building unit: In the orbit fitting and least squares parameter estimation stages, let the cumulative observation equation of all epochs be: ,in for The design matrix for The vector of differences between observed and calculated values. The observation weight matrix is ​​a symmetric rank matrix. The update operation constructs the normal equation matrix in one step, and the normal equation matrix is ​​specifically represented as follows: , in, Represents the prior constraint matrix. express The coefficient matrix of the normal equations, express The vector on the right-hand side of the normal equation, Represents the total number of observations. Represents the total number of parameters to be estimated, symmetric rank The update operation internally will observe the equations Assigned to Multiple threads concurrently compute and accumulate local outer products, achieving thread-level parallel acceleration of normal equation construction; Multi-threaded Cholesky decomposition solution unit: used for solving the coefficient matrix of the normal equations Perform Cholesky decomposition, which is represented as: , in, express The upper triangular matrix, based on the Cholesky decomposition results, yields parameter corrections through prior and back-substitution solutions. The previous generation solution and the back generation solution are expressed as: , in, Representing the intermediate vector, the Cholesky decomposition internally employs a block-based algorithm, dividing the upper triangular matrix... Assigned to Multiple threads concurrently perform decomposition and update operations to achieve thread-level parallel acceleration of parameter solving; Multi-threaded residual calculation unit: used to calculate parameter corrections. Recovery of observation residual vectors using matrix-vector multiplication Observation residual vector Represented as: , in, express Calculate the weighted sum of squared residuals from the observed residual vectors. Used for unit weight variance estimation Unit weighted variance estimation Represented as: , Matrix-vector multiplication internally designs the matrix Assigned to Multiple threads perform concurrent calculations and accumulations, achieving thread-level parallel acceleration for residual recovery.

[0016] Preferably, the algorithm-level block elimination optimization module includes: Inactive parameter identification unit: used to dynamically identify inactive local parameters during the epoch-by-epoch processing of the least squares parameter estimation stage and add them to the parameter index set to be eliminated; Normal equation block unit: used to perform row and column permutations on the normal equation matrix, dividing it into a matrix of retained parameter coefficients, a matrix of parameter coefficients to be eliminated, and a cross-coupling coefficient matrix; Schur complement elimination unit: used to perform Cholesky decomposition on the coefficient matrix of the parameter to be eliminated, and calculate Schur complement based on the decomposition result to update the reduced normal equation coefficient matrix and right-hand vector with retained parameters, thereby reducing the dimension of the normal equation; The back-substitution information storage unit is used to calculate and store the back-substitution recovery matrix for subsequent recovery of the estimated value of eliminated parameters; Normal equation compression unit: used to reclaim the storage space freed up by eliminated parameters; The parameter back-substitution recovery unit is used to recover the estimated values ​​of the eliminated parameters in each elimination process in batches, based on the saved back-substitution recovery information, after the solution of the global reduction method equation is completed in the least squares parameter estimation stage.

[0017] Preferably, the algorithm-level block elimination optimization module specifically includes: Inactive Parameter Identification Unit: During the parameter estimation stage, inactive local parameters are dynamically identified in the epoch-by-epoch processing and included in the parameter index set to be eliminated. The inactive local parameters include carrier phase ambiguity parameters where the tracking arc has ended, tropospheric delay parameters that have left the observation window, and satellite and station clock bias parameters for the corresponding epochs. The number of parameters to be eliminated is [number missing]. ; Normal equation block unit: to eliminate the The rows and columns corresponding to the parameters are moved to the bottom right corner of the matrix, giving the normal equation a block structure. The permuted normal equation is specifically represented as follows: , in, The dimension is The retained parameter vector includes orbital parameters, station coordinate parameters, zenith tropospheric delay parameters, and ambiguity parameters that are still active. The dimension is The parameter vector to be eliminated, for The retained parameter coefficient matrix, for The parameter coefficient matrix to be eliminated, and for and The cross-coupling matrix, and This is the corresponding right-hand vector; Schur complemented elimination element: the coefficient matrix of the parameters to be eliminated Perform Cholesky decomposition and calculate Schur complement based on the decomposition results to complete block elimination. Cholesky decomposition is represented as: , in, for Construct the cross matrix from the upper triangular matrix. Its dimensions are Then, perform the previous generation solution and the symmetric rank k update operation in sequence: First, the intermediate matrix is ​​obtained through the previous solution. The previous solution can be expressed as: , Then, the Schur complement matrix is ​​obtained through a symmetric rank k update operation. Schur complement matrix Represented as: , Schur complement matrix Divide into blocks, The block form is represented as: ; Using Schur complement matrix A one-time update of the reduced normal equations with preserved parameters, resulting in the reduced normal equation coefficient matrix. and the right-hand vector Represented as: , Synchronously update the weighted sum of squared residuals The cumulative value, and the update formula is specifically expressed as: , in, This represents the cumulative sum of squared weighted residuals updated after eliminating the parameters to be eliminated. After elimination, the dimension of the normal equation is changed from... Reduced to ; Back-substitution information storage unit: Used to calculate and store the back-substitution recovery matrix Z, which is used for subsequent estimation of eliminated parameters. The back-substitution recovery matrix Z is obtained through back-substitution solution. ; Normal equation compression unit: Used to perform compression processing on the reduced normal equations, reclaiming the storage space freed up by the eliminated parameters, so that the storage dimension of the normal equations is consistent with the effective parameter dimension. Assume a total of [number] operations are performed. K Next block elimination, the first k elimination times With parameters, the effective dimension of the normal equation remains constant: , in, This represents the total number of parameters introduced during the parameter estimation process. This represents the cumulative number of parameters eliminated in each elimination round. Through the dynamic block elimination strategy, the dimension of the normal equation is constrained to be much smaller than... This level of computation avoids memory overflow and significantly reduces the computational complexity of subsequent Cholesky decomposition; Parameter back-substitution recovery unit: After solving the global reduction method equations in the least squares parameter estimation stage, it uses the saved back-substitution recovery information to recover the estimated values ​​of the eliminated parameters in each elimination round in batches, for the th... k The parameter vector eliminated in the second elimination. Its back-substitution recovery formula is specifically expressed as: , in, Indicates the first k The back-substitution recovery matrix preserved by elimination express The last column is , express The former column , Indicates the first k The estimated values ​​of the retained parameters obtained in the global solution during the elimination process are strictly followed in the reverse order of elimination (the 1st elimination). k The process is repeated from the first to the second time to ensure that the estimated values ​​of the retained parameters required for each batch of recovery have been obtained.

[0018] Compared with the prior art, the present invention has the following advantages: 1. By adopting a message passing interface process-level parallel method, independent subtasks such as orbit integration, observation data preprocessing, intra-frequency code deviation calculation and station load displacement calculation are evenly distributed to multiple processes for concurrent execution according to the satellite or station dimension, thereby achieving data-level parallel acceleration. This makes the execution time of each calculation module decrease nearly linearly with the number of processes, effectively eliminating the computational bottleneck caused by hundreds of global stations and multi-system satellite constellations. 2. For computationally intensive matrix operations such as normal equation construction, Cholesky decomposition and residual recovery, the multi-threaded linear algebra subroutine library is used to achieve operator-level parallel acceleration. The matrix is ​​distributed to multiple threads for concurrent computation by row blocks or column blocks. The multi-core processor resources are fully utilized within the process, which significantly shortens the execution time of matrix operations. 3. A dynamic block elimination strategy based on Schur complement is adopted to identify and eliminate inactive local parameters in real time along the epoch progression direction. At the same time, the back-substitution recovery matrix is ​​saved so that the estimated values ​​of all eliminated parameters can be restored in reverse order after global solution. This constrains the effective dimension of the normal equation to a level much smaller than the total number of parameters, fundamentally reducing the computational complexity and memory consumption of Cholesky decomposition, enabling the system to handle large-scale orbit determination problems with hundreds of thousands of parameters. 4. The system integrates a three-layer parallel architecture that combines process-level data parallelism, thread-level operator parallelism, and algorithm-level block elimination optimization. These three layers form a nested and complementary relationship, which significantly improves computational efficiency and system scalability while maintaining centimeter-level precision orbit determination accuracy. This provides an efficient and reliable technical solution for the near real-time generation and release of precision orbit products for global navigation satellite systems. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the architecture of the rapid orbit determination system based on large-scale global GNSS station network data of the present invention; Figure 2 Comparison of calculation time for precise orbit determination of BDS3 single system; Figure 3 Comparison of calculation time for precise orbit determination of four systems. Detailed Implementation

[0020] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings and embodiments. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] like Figure 1 As shown, a rapid orbit determination system based on global large-scale GNSS station network data includes a precise orbit determination execution module and a three-layer acceleration module, with the three-layer acceleration module embedded in the multiple iterations of the precise orbit determination execution module; The precision orbit determination execution module is used to acquire observation data from the global tracking station network, complete the preprocessing and quality control of the observation data, orbit numerical integration and orbit fitting, use least squares parameter estimation to solve the equations, and combine residual editing and iterative convergence control to iteratively update the least squares parameter estimation, orbit parameter update and residual check. After the iteration is completed, ambiguity is fixed and re-estimated, and finally the multi-system precision orbit and satellite clock error products are output. The three-layer acceleration module includes: The process-level parallel acceleration module is used to execute computational tasks that can be split into satellite-level or station-level tasks serving the precision orbit determination execution module using the MPI distributed parallel mechanism. The computational tasks include orbit numerical integration, observation data preprocessing and quality control, as well as frequency code deviation calculation and station load displacement calculation to provide correction information for the precision orbit determination execution module. The thread-level parallel acceleration module is used to implement operator-level parallel acceleration for the construction of normal equations, Cholesky decomposition, residual calculation and matrix operations involved in trajectory fitting, least squares parameter estimation and block elimination. The algorithm-level block elimination optimization module is used to dynamically eliminate inactive local parameters during the least squares parameter estimation process using a block matrix elimination strategy based on Schur complement, and saves the back-substitution recovery information to control the dimension of the normal equation.

[0022] The orbital numerical integration includes the following: Satellite state vector: , in, Represents the satellite state vector (usually 6-dimensional: position + velocity). Represents the satellite velocity vector. This represents the satellite's position vector.

[0023] Dynamic differential equation: in, Let be the derivative of the state vector with respect to time. For the dynamic model function, the output is... Includes acceleration model, This represents the parameter vector of the dynamic model.

[0024] Linearization (state transition): in, This is the state correction (the increment relative to the reference orbit). Here is the state transition matrix. This is the initial state correction. The parameter sensitivity matrix is ​​for dynamic parameters. This is the correction amount for the dynamic parameters.

[0025] Numerical integration: given initial values In the solar arc Integrating yields the reference orbit: , in Represents the numerical integration operator. The reference orbital state (the prior / current iteration orbital obtained by numerical integration).

[0026] The orbit fitting process involves fitting the orbit to an external orbit and updating the initial orbit state and dynamic parameters.

[0027] Least squares parameter estimation includes the following: Pseudorange observation equation: , Carrier phase observation equation: , in: In order to observe epoch time, j For receivers or measuring stations, s For satellites, speed of light in vacuum, For the station j To satellite s geometric distance, For the station j To satellite s At frequency f pseudorange observations, For the station j To satellite s At frequency fThe carrier observations (if in weeks, multiply by ) (conversion) For receiver / station clock bias, For satellite clock bias, For tropospheric delay, In frequency point f Ionospheric delay, This is the phase-related comprehensive deviation term. This is the comprehensive deviation term related to pseudorange. This is due to pseudorange observation noise and unmodeled errors. This includes carrier observation noise and unmodeled errors. For the station j With satellite s At frequency f Integer blur, For frequency point f The corresponding wavelength.

[0028] Linearization formation: ,in, For OC vector (by , (Observations minus theoretical values ​​calculated using a reference orbit and current parameters) These are corrections to the parameters to be estimated (initial orbital values / dynamic parameters / clock error / troposphere / ambiguity / deviation term, etc.). For designing the matrix (including partial derivatives with respect to orbit, clock error, troposphere, ambiguity, etc.).

[0029] Weighted least squares objective function: ,in, Let be the weight matrix. Normal equation: ,in, Let be the normal matrix. Let be the vector on the right-hand side of the normal equation. Residuals and statistics: ,in, This represents the residual vector. Unit weight variance factor: ,in, The standard deviation factor is the unit weight. The ratio of the number of observation equations to the total number of observations. This is the number of unknown parameters divided by the total dimension of the parameters to be estimated.

[0030] In this embodiment, numerical integration is performed based on the satellite dynamics model and the updated orbit parameters to obtain the reference orbit, state transition matrix, and parameter sensitivity matrix. Then, the residuals are calculated based on the floating-point solution and residual editing is performed to determine whether the convergence or quality control conditions are met. If not, the orbit numerical integration is performed again for iterative updates.

[0031] In the above process, a three-layer acceleration module is embedded. Specifically, the process-level parallel acceleration module includes: MPI Process Initialization and Task Allocation Unit: Used to initialize the MPI computing environment and start... Each MPI process sets up the overall computing task, including... Each subtask is assigned to a specific MPI process using a round-robin strategy. Specifically, the process starts... Each MPI process obtains the total number of MPI processes in the global communication domain. With its own MPI process number r Set up the overall computing task, including Each subtask is assigned sequentially to an MPI process using a round-robin strategy. r The set of task indexes under your responsibility Represented as: , in, r Indicates the MPI process number. , This represents the total number of subtasks in the current overall computing task, when broken down by satellite dimension. When splitting by station dimension , This indicates the total number of satellites participating in orbit determination. This indicates the total number of stations participating in orbit determination; The satellite orbit numerical integration parallel computing unit is used to enable each MPI process to independently perform orbit numerical integration on the satellite subtasks assigned to it, based on the task allocation results. Specifically, the satellite orbit numerical integration parallel computing unit is used to enable each MPI process to independently perform orbit numerical integration on the satellite subtasks assigned to it, based on the task allocation results. r For task index set Satellites assigned to their own sub-missions independently perform orbital numerical integration; for those assigned to the MPI process... satellite sub-mission Its six-dimensional state vector The equation of motion that is satisfied is expressed as: , in, Indicates satellite s The three-dimensional position vector in the geocentric inertial coordinate system. This represents the corresponding three-dimensional velocity vector. t For time, F It represents a nonlinear dynamic function composed of the Earth's gravitational field, the gravitational pull of the Sun and Moon, solar radiation pressure, Earth's tidal perturbations, and relativistic effects. Different satellites are independent of each other, and there is no need for inter-process data exchange during the integration process of each MPI process. Parallel computing unit for preprocessing station observation data: Based on the task allocation results, each MPI process independently performs GNSS observation data preprocessing and quality control on the station sub-tasks assigned to it; Intra-frequency code offset parallel calculation unit: This unit is used by each MPI process to independently construct the intra-frequency code offset differential observations for its assigned station subtasks based on the task allocation results. Specifically, the intra-frequency code offset parallel calculation unit is used by the MPI processes to construct the intra-frequency code offset differential observations based on the task allocation results. r For task index set Each station subtask assigned to it independently reads pseudorange observation data and constructs pseudorange differential observations between different signal attributes on the same carrier frequency; for those assigned to the MPI process... r Station subtask The observed satellite s in epoch The differential observation equation for its intra-frequency code bias is expressed as: , in, and Representing the stations j In the calendar Received satellite s In frequency f The above are signal attributes and Encoded pseudorange observations, Indicates the satellite end-frequency code deviation. Indicates the intra-frequency code deviation at the receiver end. To represent observation noise, each MPI process independently performs data reading and differential observation construction for the station subtasks assigned to it, and different stations are independent of each other; The station load-displacement parallel calculation unit is used to independently perform load deformation displacement correction calculations for the station subtasks assigned to it by each MPI process, based on the task allocation results. Specifically, the station load-displacement parallel calculation unit is used to perform load deformation displacement correction calculations for each MPI process based on the task allocation results. r For task index set The subtasks assigned to each station within the system independently perform load deformation displacement correction calculations, while those assigned to the MPI process... r Station subtask In the calendar The total load displacement correction vector is expressed as: , in, This indicates the displacement caused by ocean tidal loads. This indicates the displacement caused by the polar tide load. This indicates the displacement caused by atmospheric tidal loads. This represents displacement caused by non-tidal atmospheric and ocean loads. The load displacement correction for each station depends only on the station's own coordinates and external geophysical model parameters. The stations are independent of each other and are executed in complete parallel. Process Synchronization and Result Summarization Unit: This unit synchronizes the processes after the completion of the satellite orbit numerical integration parallel calculation unit, the station observation data preprocessing parallel calculation unit, the intra-frequency code deviation parallel calculation unit, and the station load displacement parallel calculation unit. It selects one MPI process as the master process and summarizes the local result files generated by each MPI process to form a complete result file for subsequent orbit determination calculations. The theoretical speedup of process-level parallelism is [not specified]. Specifically, it is expressed as follows: , Parallel efficiency of process-level parallelism Specifically, it is expressed as follows: in, Indicates all The total time taken for each subtask to be executed serially Indicates adoption The total time taken by executing MPI processes in parallel, when Can be When divisible, This achieves ideal load balancing.

[0032] In addition, the thread-level parallel acceleration module includes: Multithreaded environment initialization unit: Used to set the number of parallel threads in the orbit fitting and least squares parameter estimation stages, and to enable the multithreaded linear algebra subroutine. Specifically, it is used to set the number of parallel threads in the orbit fitting and least squares parameter estimation stages. This allows subsequent linear algebra subroutines involving matrix operations to be executed within the process. The threads execute concurrently, and the linear algebra subroutines are provided by the multithreaded basic linear algebra subroutine library, covering subroutines for three-level matrices and matrix operations, two-level matrices and vector operations, and decomposition and solution of linear equation systems. The multi-threaded normal equation construction unit is used to call multi-threaded linear algebra subroutines to perform thread-level parallel acceleration of matrix operations in the normal equation construction process, under the condition of parallel threads, in order to construct the normal equation coefficient matrix. Specifically, in the orbit fitting and least squares parameter estimation stage, let the cumulative observation equation of all epochs be... ,in for The design matrix for The vector of differences between observed and calculated values. The observation weight matrix is ​​a symmetric rank matrix. The update operation constructs the normal equation matrix in one step, and the normal equation matrix is ​​specifically represented as follows: , in, Represents the prior constraint matrix. express The coefficient matrix of the normal equations, express The vector on the right-hand side of the normal equation, Represents the total number of observations. Represents the total number of parameters to be estimated, symmetric rank The update operation internally will observe the equations Assigned to Multiple threads concurrently compute and accumulate local outer products, achieving thread-level parallel acceleration of normal equation construction; Multi-threaded Cholesky decomposition unit: This unit accelerates the Cholesky decomposition of the normal equation coefficient matrix and subsequent prior and back substitution problems by calling multi-threaded linear algebra subroutines under parallel thread conditions. Specifically, it is used for the Cholesky decomposition of the normal equation coefficient matrix... Perform Cholesky decomposition, which is represented as: , in, express The upper triangular matrix, based on the Cholesky decomposition results, yields parameter corrections through prior and back-substitution solutions. The previous generation solution and the back generation solution are expressed as: , in, Representing the intermediate vector, Cholesky decomposition internally employs a block-based algorithm, dividing the upper triangular matrix... Assigned to Multiple threads concurrently perform decomposition and update operations to achieve thread-level parallel acceleration of parameter solving; Multi-threaded residual calculation unit: Used to perform thread-level parallel acceleration of matrix-vector operations in the residual recovery process by calling multi-threaded linear algebra subroutines under parallel thread conditions, in order to recover the observation residuals. Specifically, it is used to calculate parameter corrections. Recovery of observation residual vectors using matrix-vector multiplication Observation residual vector Represented as: , in, express Calculate the weighted sum of squared residuals from the observed residual vectors. Used for unit weight variance estimation Unit weighted variance estimation Represented as: , Matrix-vector multiplication internally designs the matrix Assigned to Multiple threads perform concurrent calculations and accumulations, achieving thread-level parallel acceleration for residual recovery.

[0033] Finally, the algorithm-level block elimination optimization module includes: Inactive Parameter Identification Unit: This unit dynamically identifies inactive local parameters during the epoch-by-epoch processing of the least squares parameter estimation stage and adds them to the parameter index set to be eliminated. Specifically, during the parameter estimation stage, inactive local parameters are dynamically identified during the epoch-by-epoch processing and added to the parameter index set to be eliminated. Inactive local parameters include carrier phase ambiguity parameters where the tracking arc has ended, tropospheric delay parameters that have left the observation window, and satellite and station clock bias parameters for the corresponding epochs. The number of parameters to be eliminated is [number missing]. ; Normal equation block unit: Used to perform row and column permutations on the normal equation matrix, dividing it into a matrix of retained parameter coefficients, a matrix of parameter coefficients to be eliminated, and a cross-coupling coefficient matrix. Specifically, the matrix of parameter coefficients to be eliminated... The rows and columns corresponding to the parameters are moved to the bottom right corner of the matrix, giving the normal equation a block structure. The permuted normal equation is specifically represented as follows: , in, The dimension is The retained parameter vector includes orbital parameters, station coordinate parameters, zenith tropospheric delay parameters, and ambiguity parameters that are still active. The dimension is The parameter vector to be eliminated, for The retained parameter coefficient matrix, for The parameter coefficient matrix to be eliminated, and for and The cross-coupling matrix, and This is the corresponding right-hand vector; Schur complement elimination unit: Used to perform Cholesky decomposition on the coefficient matrix of the parameters to be eliminated, and calculate Schur complement based on the decomposition result to update the reduced normal equation coefficient matrix and right-hand side vector with preserved parameters, thereby reducing the dimension of the normal equation. Specifically, for the coefficient matrix of the parameters to be eliminated... Perform Cholesky decomposition and calculate Schur complement based on the decomposition results to complete block elimination. Cholesky decomposition is represented as: , in, for Construct the cross matrix from the upper triangular matrix. Its dimensions are Then, perform the previous generation solution and the symmetric rank k update operation in sequence: First, the intermediate matrix is ​​obtained through the previous solution. The previous solution can be expressed as: , Then, the Schur complement matrix is ​​obtained through a symmetric rank k update operation. Schur complement matrix Represented as: , Schur complement matrix Divide into blocks, The block form is represented as: ; Using Schur complement matrix A one-time update of the reduced normal equations with preserved parameters, resulting in the reduced normal equation coefficient matrix. and the right-hand vector Represented as: , Synchronously update the weighted sum of squared residuals The cumulative value, and the update formula is specifically expressed as: , in, This represents the cumulative sum of squared weighted residuals updated after eliminating the parameters to be eliminated. After elimination, the dimension of the normal equation is changed from... Reduced to ; The back-substitution information storage unit is used to calculate and store the back-substitution recovery matrix Z, which is used for subsequent recovery of the estimated value of the eliminated parameters. The back-substitution recovery matrix Z is obtained through back-substitution solution. ; Normal equation compression unit: Used to perform compression processing on the reduced normal equations, reclaiming the storage space freed up by the eliminated parameters, so that the storage dimension of the normal equations is consistent with the effective parameter dimension. Assume a total of [number] operations are performed. K Next block elimination, the first k elimination times With parameters, the effective dimension of the normal equation remains constant: , in, This represents the total number of parameters introduced during the parameter estimation process. This represents the cumulative number of parameters eliminated in each elimination round. Through the dynamic block elimination strategy, the dimension of the normal equation is constrained to be much smaller than... This level of computation avoids memory overflow and significantly reduces the computational complexity of subsequent Cholesky decomposition; Parameter back-substitution recovery unit: After solving the global reduction method equations in the least squares parameter estimation stage, it uses the saved back-substitution recovery information to recover the estimated values ​​of the eliminated parameters in each elimination round in batches, for the th... k The parameter vector eliminated in the second elimination. Its back-substitution recovery formula is specifically expressed as: , in, Indicates the first k The back-substitution recovery matrix preserved by elimination express The last column is , express The former column , Indicates the first k The estimated values ​​of the retained parameters obtained in the global solution during the elimination process are strictly followed in the reverse order of elimination (the 1st elimination). k The process is repeated from the first to the second time to ensure that the estimated values ​​of the retained parameters required for each batch of recovery have been obtained.

[0034] The above embodiments also relate to a computer-readable storage medium storing a computer program that, when executed by a processor, implements the function of the rapid orbit determination system based on global large-scale GNSS station network data of the present invention.

[0035] Based on observational data from 40 to 200 IGS monitoring stations worldwide, precise orbit determination calculations were performed for two typical scenarios: a single system (BDS3 series, 24 satellites) and a four-system system (GPS / GLONASS / Galileo / BDS, approximately 120 satellites). The orbit determination arc was set to 24 hours, and the data sampling interval was set to 300 seconds. The hardware platform employed a high-performance computing server equipped with dual Intel Xeon Gold 6258R processors. The Intel Xeon Gold 6258R processors have a base clock speed of 2.70GHz and a maximum turbo frequency of 4.00GHz. Each processor has 28 physical computing cores, for a total of 56 physical computing cores across both processors.

[0036] Under the same experimental conditions, the traditional single-process serial scheme uses only a single processing core to execute all calculation steps such as data preprocessing, orbit integration, and parameter estimation in sequence; the three-layer parallel scheme proposed in this embodiment simultaneously enables process-level parallel acceleration, thread-level parallel acceleration, and algorithm-level block elimination optimization, making full use of all the computing core resources of the dual-processor.

[0037] Combination Figure 2 As shown in Table 1, a comparison of the calculation time for precise orbit determination using the BDS3 single system is presented.

[0038] Table 1 Comparison of calculation time for precise orbit determination using BDS3 single system With a single system and 24 satellites, the computation time of the serial scheme increases from 30 minutes to approximately 71 minutes with the increase in the number of stations, an increase of about 2.4 times. After adopting the three-layer parallel scheme, the computation time slowly increases from 10.9 minutes to 22.27 minutes, an increase of only 2.0 times, showing a significantly slower growth trend. With 150 stations, the parallel scheme improves computational efficiency by approximately 68.6% compared to the serial scheme, achieving a speedup of 3.19 times. However, due to the relatively small number of satellites in a single system, the task allocation granularity at the process level is limited, thus restricting the speedup effect to some extent.

[0039] It should be noted that the speedup ratio in Table 1 does not increase monotonically with the number of stations; this is a normal phenomenon in parallel computing. The speedup ratio equals the serial execution time divided by the parallel execution time. The parallel execution time increases steadily with the number of stations, while the serial execution time is affected by factors such as differences in the amount of observation data from each station, fixed overhead, and computer cache and memory, resulting in irregular growth. This fluctuation is reflected in the speedup ratio, causing it to fluctuate. However, across the entire range of 40 to 150 stations, this system consistently achieved a stable speedup of 2.75 to 4.27 times, indicating that both the speedup effect and the experimental data are reliable.

[0040] Combination Figure 3 As shown in Table 2, the computation time for different station sizes under the four-system full constellation configuration is compared.

[0041] Table 2 Comparison of computation time for different station sizes under the four-system full constellation configuration With a configuration of approximately 120 satellites across four systems, the computation time of the serial scheme exhibits a significant superlinear growth characteristic with the increase in the number of stations: it increases dramatically from approximately 263 minutes with 40 stations to approximately 2541 minutes (approximately 42.3 hours) with 200 stations, an increase of 9.7 times. In contrast, the computation time using a three-layer parallel scheme only increases from 49.77 minutes to 140.53 minutes (approximately 2.3 hours), an increase of 2.8 times, with a near-linear growth trend. The difference in computation time between the two schemes widens dramatically with the expansion of the station scale: the speedup is 5.29 times with 40 stations, increases to 9.59 times with 100 stations, and further increases to 18.08 times with 200 stations, resulting in a 94.5% improvement in computational efficiency and a breakthrough in computing speed.

[0042] The experimental results of the two sets above show that the parallel acceleration effect of the method in this embodiment is positively correlated with both the number of satellites and the scale of the stations. In the small-scale scenario of a single system with 24 satellites, the speedup ratio remains between 2.75 and 4.27 times due to the limited number of parallelizable tasks. However, in the large-scale scenario of a four-system system with approximately 120 satellites, a large number of satellite orbit integration and station data preprocessing tasks can be fully allocated to multiple processes for concurrent execution. At the same time, the significantly increased dimension of the normal equations in the parameter estimation stage allows the parallelism of thread-level matrix operations and the optimization of algorithm-level block elimination to be fully utilized. The speedup ratio continues to climb with the increase in the scale of the stations, reaching a maximum of 18.08 times.

[0043] Therefore, the rapid orbit determination system based on global large-scale GNSS station network data can effectively solve the computational efficiency bottleneck under the conditions of large-scale station network and multi-system satellite constellation. It has good scalability and can significantly shorten the calculation time while maintaining centimeter-level orbit determination accuracy, thereby significantly improving timeliness and meeting the needs of near real-time generation and release of precise orbit products of global navigation satellite systems.

[0044] It should be noted that the above description of the technical solutions is exemplary, and this specification may be embodied in different forms and should not be construed as limiting it to the technical solutions set forth herein. Rather, providing these descriptions will ensure that the disclosure of this invention is thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art. Furthermore, the technical solutions of this invention are defined only by the scope of the claims.

[0045] The aspects disclosed in this specification and claims are merely examples, and therefore, this specification and claims are not limited to the details shown. In the foregoing description, detailed descriptions of relevant known functions or configurations have been omitted where it would unnecessarily obscure the focus of this specification and claims.

[0046] Finally, it should be noted that the above description, in conjunction with specific embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, any simple substitutions made without departing from the inventive concept should be considered within the scope of protection of the present invention. The above embodiments are merely representative examples of the present invention. Obviously, the present invention is not limited to the above embodiments and many variations are possible. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A rapid orbit determination system based on global large-scale GNSS station network data, characterized in that: It includes a precision orbit determination execution module and a three-layer acceleration module, wherein the three-layer acceleration module is embedded in the multiple iterations of the precision orbit determination execution module; The precision orbit determination execution module is used to acquire observation data from the global tracking station network, complete the preprocessing and quality control of the observation data, orbit numerical integration and orbit fitting, use the least squares parameter estimation method to solve the equations, and combine residual editing and iterative convergence control to iteratively update the least squares parameter estimation, orbit parameter update and residual check. After the iteration is completed, the ambiguity is fixed and re-estimated, and finally the multi-system precision orbit and satellite clock error products are output. The three-layer acceleration module includes: The process-level parallel acceleration module is used to execute computational tasks that can be split into satellite-dimensional or station-dimensional components serving the precise orbit determination execution module using the MPI distributed parallel mechanism. The computational tasks include orbital numerical integration, observation data preprocessing and quality control, as well as frequency code deviation calculation and station load displacement calculation that provide correction information to the precise orbit determination execution module. The thread-level parallel acceleration module is used to implement operator-level parallel acceleration for the construction of normal equations, Cholesky decomposition, residual calculation and matrix operations involved in trajectory fitting, least squares parameter estimation and block elimination. The algorithm-level block elimination optimization module is used to dynamically eliminate inactive local parameters during the least squares parameter estimation process using a block matrix elimination strategy based on Schur complement, and saves the back-substitution recovery information to control the dimension of the normal equation.

2. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 1, characterized in that: The process-level parallel acceleration module includes: MPI Process Initialization and Task Allocation Unit: Used to initialize the MPI computing environment and start... Each MPI process sets up the overall computing task, including Each subtask is assigned to an MPI process in a round-robin fashion. Parallel computation unit for satellite orbit numerical integration: Used by each MPI process to independently perform orbit numerical integration on the satellite subtask assigned to it, based on the task allocation results; Parallel computing unit for preprocessing station observation data: Based on the task allocation results, each MPI process independently performs GNSS observation data preprocessing and quality control on the station sub-tasks assigned to it; Intra-frequency code deviation parallel calculation unit: used to independently execute the intra-frequency code deviation differential observation construction for each station subtask assigned to it based on the task allocation results; Parallel calculation unit for station load displacement: Based on the task allocation results, each MPI process independently performs load deformation displacement correction calculations for the station subtasks assigned to it; Process synchronization and result aggregation unit: After the satellite orbit numerical integration parallel calculation unit, station observation data preprocessing parallel calculation unit, intra-frequency code deviation parallel calculation unit and station load displacement parallel calculation unit have been executed, process synchronization is performed, and an MPI process is selected as the main process. The local result files generated by each MPI process are aggregated to form a complete result file for subsequent orbit determination calculation.

3. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 2, characterized in that: MPI process initialization and task allocation unit: Startup Each MPI process obtains the total number of MPI processes in the global communication domain. With its own MPI process number r Set up the overall computing task, including Each subtask is assigned sequentially to an MPI process using a round-robin strategy. r The set of task indexes under your responsibility Represented as: ， in, r Indicates the MPI process number. , This represents the total number of subtasks in the current overall computing task, when broken down by satellite dimension. When splitting by station dimension , This indicates the total number of satellites participating in orbit determination. This indicates the total number of stations involved in orbit determination.

4. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 3, characterized in that: The satellite orbit numerical integration parallel computing unit is used by the MPI process based on the task allocation results. r For task index set Satellites assigned to their own sub-missions independently perform orbital numerical integration; for those assigned to the MPI process... satellite sub-mission Its six-dimensional state vector The equation of motion that is satisfied is expressed as: ， in, Indicates satellite s The three-dimensional position vector in the geocentric inertial coordinate system. This represents the corresponding three-dimensional velocity vector. t For time, F It represents a nonlinear dynamic function composed of the Earth's gravitational field, the gravitational pull of the Sun and Moon, solar radiation pressure, Earth's tidal perturbations, and relativistic effects. Different satellites are independent of each other, and there is no need for inter-process data exchange during the integration process of each MPI process.

5. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 3, characterized in that: The intra-frequency code deviation parallel calculation unit is used by the MPI process based on the task allocation results. r For task index set Each station subtask assigned to itself independently reads pseudorange observation data and constructs pseudorange differential observations between different signal attributes on the same carrier frequency; For allocation to MPI process r Station subtask The observed satellite s in epoch The differential observation equation for its intra-frequency code bias is expressed as: ， in, and Representing the stations j In the calendar Received satellite s In frequency f The above are signal attributes and Encoded pseudorange observations, Indicates the satellite end-frequency code deviation. Indicates the intra-frequency code deviation at the receiver end. This represents observation noise. Each MPI process independently performs data reading and differential observation construction for the station subtasks assigned to it, and different stations are independent of each other.

6. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 3, characterized in that: The parallel calculation unit for station load displacement is used by the MPI process based on the task allocation results. r For task index set The subtasks assigned to each station within the system independently perform load deformation displacement correction calculations, while those assigned to the MPI process... r Station subtask In the calendar The total load displacement correction vector is expressed as: ， in, This indicates the displacement caused by ocean tidal loads. This indicates the displacement caused by the polar tide load. This indicates the displacement caused by atmospheric tidal loads. This represents displacement caused by non-tidal atmospheric and ocean loads. The load displacement correction of each station depends only on its own coordinates and external geophysical model parameters. The stations are independent of each other and are executed in complete parallel.

7. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 1, characterized in that: The thread-level parallel acceleration module includes: Multi-threaded environment initialization unit: used to set the number of parallel threads and enable the multi-threaded linear algebra subroutine during the orbit fitting and least squares parameter estimation stages; A multi-threaded normal equation construction unit is used to call the multi-threaded linear algebra subroutine to perform thread-level parallel acceleration of matrix operations in the normal equation construction process under the number of parallel threads, so as to construct the normal equation coefficient matrix. Multi-threaded Cholesky decomposition solution unit: used to call the multi-threaded linear algebra subroutine to perform Cholesky decomposition of the coefficient matrix of the normal equation and subsequent previous and back substitution solutions under the number of parallel threads; Multi-threaded residual calculation unit: used to call the multi-threaded linear algebra subroutine to perform thread-level parallel acceleration of matrix-vector operations in the residual recovery process under the number of parallel threads, so as to recover the observation residuals.

8. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 7, characterized in that: Multi-threaded environment initialization unit: used to set the number of parallel threads during the orbit fitting and least squares parameter estimation stages. This allows subsequent linear algebra subroutines involving matrix operations to be executed within the process. The linear algebra subroutines are executed concurrently by multiple threads. The linear algebra subroutines are provided by the multi-threaded basic linear algebra subroutine library, covering subroutines for three-level matrix-matrix operations, two-level matrix-vector operations, and decomposition and solution of linear equation systems. Multi-threaded method equation building unit: In the orbit fitting and least squares parameter estimation stages, let the cumulative observation equation of all epochs be: ,in for The design matrix for The vector of differences between observed and calculated values. The observation weight matrix is ​​a symmetric rank matrix. The update operation constructs the normal equation matrix in one step, and the normal equation matrix is ​​specifically represented as follows: ， in, Represents the prior constraint matrix. express The coefficient matrix of the normal equations, express The vector on the right-hand side of the normal equation, Represents the total number of observations. Represents the total number of parameters to be estimated, symmetric rank The update operation internally will observe the equations Assigned to Multiple threads concurrently compute and accumulate local outer products, achieving thread-level parallel acceleration of normal equation construction; Multi-threaded Cholesky decomposition solution unit: used for solving the coefficient matrix of the normal equations Perform Cholesky decomposition, which is represented as: ， in, express The upper triangular matrix, based on the Cholesky decomposition results, yields parameter corrections through prior and back-substitution solutions. The previous generation solution and the back generation solution are expressed as: ， in, Representing the intermediate vector, the Cholesky decomposition internally employs a block-based algorithm, dividing the upper triangular matrix... Assigned to Multiple threads concurrently perform decomposition and update operations to achieve thread-level parallel acceleration of parameter solving; Multi-threaded residual calculation unit: used to calculate parameter corrections. Recovery of observation residual vectors using matrix-vector multiplication Observation residual vector Represented as: ， in, express Calculate the weighted sum of squared residuals from the observed residual vectors. Used for unit weight variance estimation Unit weighted variance estimation Represented as: ， Matrix-vector multiplication internally designs the matrix Assigned to Multiple threads perform concurrent calculations and accumulations, achieving thread-level parallel acceleration for residual recovery.

9. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 1, characterized in that: The algorithm-level block elimination optimization module includes: Inactive parameter identification unit: used to dynamically identify inactive local parameters during the epoch-by-epoch processing of the least squares parameter estimation stage and add them to the parameter index set to be eliminated; Normal equation block unit: used to perform row and column permutations on the normal equation matrix, dividing it into a matrix of retained parameter coefficients, a matrix of parameter coefficients to be eliminated, and a cross-coupling coefficient matrix; Schur complement elimination unit: used to perform Cholesky decomposition on the coefficient matrix of the parameter to be eliminated, and calculate Schur complement based on the decomposition result to update the reduced normal equation coefficient matrix and right-hand vector with retained parameters, thereby reducing the dimension of the normal equation; The back-substitution information storage unit is used to calculate and store the back-substitution recovery matrix for subsequent recovery of the estimated value of eliminated parameters; Normal equation compression unit: used to reclaim the storage space freed up by eliminated parameters; The parameter back-substitution recovery unit is used to recover the estimated values ​​of the eliminated parameters in each elimination process in batches, based on the saved back-substitution recovery information, after the solution of the global reduction method equation is completed in the least squares parameter estimation stage.

10. The rapid orbit determination system based on large-scale global GNSS station network data as described in claim 9, characterized in that: The algorithm-level block elimination optimization module specifically includes: Inactive Parameter Identification Unit: During the parameter estimation stage, inactive local parameters are dynamically identified in the epoch-by-epoch processing and included in the parameter index set to be eliminated. The inactive local parameters include carrier phase ambiguity parameters where the tracking arc has ended, tropospheric delay parameters that have left the observation window, and satellite and station clock bias parameters for the corresponding epochs. The number of parameters to be eliminated is [number missing]. ; Normal equation block unit: to eliminate the The rows and columns corresponding to the parameters are moved to the bottom right corner of the matrix, giving the normal equation a block structure. The permuted normal equation is specifically represented as follows: ， in, The dimension is The retained parameter vector includes orbital parameters, station coordinate parameters, zenith tropospheric delay parameters, and ambiguity parameters that are still active. The dimension is The parameter vector to be eliminated, for The retained parameter coefficient matrix, for The parameter coefficient matrix to be eliminated, and for and The cross-coupling matrix, and This is the corresponding right-hand vector; Schur complemented elimination element: the coefficient matrix of the parameters to be eliminated Perform Cholesky decomposition and calculate Schur complement based on the decomposition results to complete block elimination. Cholesky decomposition is represented as: ， in, for Construct the cross matrix from the upper triangular matrix. Its dimensions are Then, perform the previous generation solution and the symmetric rank k update operation in sequence: First, the intermediate matrix is ​​obtained by solving the previous solution. The previous solution can be expressed as: ， Then, the Schur complement matrix is ​​obtained through a symmetric rank k update operation. Schur complement matrix Represented as: ， Schur complement matrix Divide into blocks, The block form is represented as: ； Using Schur complement matrix A one-time update of the reduced normal equations with preserved parameters, resulting in the reduced normal equation coefficient matrix. and the right-hand vector Represented as: ， Synchronously update the weighted sum of squared residuals The cumulative value, and the update formula is specifically expressed as: ， in, This represents the cumulative sum of squared weighted residuals updated after eliminating the parameters to be eliminated. After elimination, the dimension of the normal equation is changed from... Reduced to ; Back-substitution information storage unit: Used to calculate and store the back-substitution recovery matrix Z, which is used for subsequent estimation of eliminated parameters. The back-substitution recovery matrix Z is obtained through back-substitution solution. ， in, The dimension is The back-substitution recovery matrix Z, along with the names of the eliminated parameters, their respective station and satellite numbers, effective time periods, and the row and column permutation index table, are written into a temporary file for storage. Normal equation compression unit: Used to perform compression processing on the reduced normal equations, reclaiming the storage space freed up by the eliminated parameters, so that the storage dimension of the normal equations is consistent with the effective parameter dimension. Assume a total of [number] operations are performed. K Next block elimination, the first k elimination times With parameters, the effective dimension of the normal equation remains constant: ， in, This represents the total number of parameters introduced during the parameter estimation process. This indicates the cumulative number of parameters eliminated in each elimination round; Parameter back-substitution recovery unit: After solving the global reduction method equations in the least squares parameter estimation stage, it uses the saved back-substitution recovery information to recover the estimated values ​​of the eliminated parameters in each elimination round in batches, for the th... k The parameter vector eliminated in the second elimination. Its back-substitution recovery formula is specifically expressed as: ， in, Indicates the first k The back-substitution recovery matrix preserved by elimination express The last column is , express The former column , Indicates the first k The estimated values ​​of the retained parameters obtained in the global solution during the elimination process are strictly executed in reverse order of elimination to ensure that the estimated values ​​of the retained parameters required for each batch of recovery are obtained.

Images