Parallel least squares recursive parameter estimation method for GNSS station network observation data processing

Abstract

The present invention relates to a parallel least squares recursive parameter estimation method for processing GNSS station network observation data, belonging to the technical field of solution processing of large-scale GNSS station network observation data, and solves the problems of low hardware resource utilization efficiency, low matrix operation efficiency, and long data processing time in the prior art. The present invention derives an estimator based on the mathematical formula of least squares in the navigation algorithm, realizes multi-threaded parallel operation, and optimizes the CPU parallel processing of each epoch data matrix, thereby improving the utilization efficiency of hardware resources and matrix operation efficiency, with significant optimization effect; shortening data processing time, and maintaining complete consistency in the optimized calculation results, with a very stable optimization effect; combining the technical characteristics of GSTAR with multi-threaded optimization, the invention has many expansion directions, is universal, and is convenient for transplanting to different systems.

Description

Technical Field

[0001] The present invention relates to the technical field of solution processing of large-scale GNSS station network observation data, and in particular to a parallel least squares recursive parameter estimation method for GNSS station network observation data processing. Background Art

[0002] With the accelerated progress of my country's Beidou system's global coverage and the increasing density and upgrades of its ground-based reference station network, large-scale GNSS networks are experiencing new development opportunities. As a critical national ground infrastructure, many countries and regions have established GNSS (Global Navigation Satellite System) reference station networks comprised of hundreds or even thousands of continuously operating reference stations. While observation data from large-scale global GNSS networks provides abundant computing resources, it also presents significant data processing challenges. As the volume of GNSS network observation data increases, the numerous large-scale linear matrix operations involved increase exponentially in data processing time. This leads to computational inefficiency and even impossibility in processing global GNSS network observation data. Consequently, high-efficiency data processing models and methods are urgently needed to fully leverage the advantages of the GSTAR network and enhance the accuracy, reliability, and real-time performance of its services. The pursuit of efficient and rapid processing of large-scale GNSS network data has become a hot topic and is receiving increasing attention and attention.

[0003] Distributed computing technology, with its efficient resource utilization and computational efficiency, is increasingly being applied to geodesy. Since the essence of high-performance computing is parallel computing, decomposing traditional GNSS network observation data processing techniques and designing appropriate parallel algorithms is paramount for achieving rapid GNSS network observation data processing. Most existing optimization methods focus on matrix operations and high-performance computing, or utilize matrix partitioning to improve computational efficiency. However, these improvements alone offer limited efficiency gains, with little consideration given to optimizing program execution and thus failing to effectively utilize computer resources. Faced with the high computational cost of GNSS network observation data processing, one traditional strategy is the subnet solution, which breaks down the heavy workload of large network processing into several smaller tasks. However, theoretically, subnet solutions cannot guarantee the most comprehensive results compared to integrated solutions. Traditional GNSS network observation data processing programs are often written in a serial mode for single-processor architectures, inefficiently utilizing the performance of current data processing hardware platforms. For GNSS estimators, mainstream programs currently use the traditional least-squares algorithm. These programs typically run serially. However, the mathematical derivation of the least-squares recursive algorithm strongly depends on the matrix data before and after, making it difficult to parallelize the program on a large scale. Consequently, due to the nature of the estimator algorithm, few researchers have studied the program's operational structure. Starting from the mathematical formulas for least-squares in navigation algorithms, they identify the matrix and parameter expressions required for each step and then combine them with the GSTAR program for parallelization. Since least-squares involves iterative data computation, it is difficult to directly parallelize the program on a large scale. Instead, it is necessary to identify subdivided parts of each step to share data in parallel, thereby amortizing the cost of moving data between the CPU and memory and improving efficiency. As one of the most fundamental linear algebra kernels, the Basic Linear Algebra Subroutines (BLAS) provide a set of matrix-matrix operations with a portable interface. This is designed to amortize the cost of moving data between the processor and memory on modern microprocessors. However, from a computational perspective, there is still room for optimization.

[0004] In summary, the existing technology has the problems of low utilization efficiency of hardware resources, low matrix operation efficiency, and long data processing time. Summary of the Invention

[0005] In view of the above problems, the present invention provides a parallel least squares recursive parameter estimation method for GNSS station network observation data processing, which solves the problems of low hardware resource utilization efficiency, low matrix operation efficiency and long data processing time in the existing technology.

[0006] The present invention provides a parallel least squares recursive parameter estimation method for GNSS station network observation data processing, characterized in that a parallel least squares recursive algorithm is used to perform time update and measurement update; wherein the time update includes the following steps:

[0007] Step S1. The GSTAR program starts running and reads the GNSS observation data of the nth epoch and preprocesses it according to the .cfg configuration file;

[0008] Step S2. Processing the pre-processed GNSS observation data of the nth epoch using a parallel least squares recursive algorithm to calculate the calculation parameters;

[0009] According to the dimension of the calculated parameters, the observation coefficient matrix N is obtained (n) ; For the observation coefficient matrix N (n) Initialize the blocks and then assign values ​​to them in parallel through a for loop;

[0010] Step S3. Perform the first write and read operations in parallel; wherein the first write includes storing the equation of the recursive expression and storing the equation coefficient information generated by the nth epoch according to the recursive expression through time update; the recursive expression is used to represent the calculation parameters and the observation value coefficient matrix N (n) The functional relationship between the initialization blocks;

[0011] The first reading includes calculating the equation coefficient information generated by the nth epoch by calling the calculation function to obtain the observation coefficient matrix after the time update Including the equation coefficient information generated by the nth epoch after the time update and the intermediate value parameters of the recursive expression after the time update;

[0012] Step S4. Performing a second parallel write and read operation; wherein the second write operation includes storing equation coefficient information generated at the nth epoch after the time update and intermediate value parameters of the recursive expression after the time update to an intermediate file; the second write operation also includes determining the size of the intermediate file and dividing the intermediate file;

[0013] The second read includes the updated observation coefficient matrix according to the intermediate value parameters of the recursive expression after time update. Divide into blocks, and then reassign Make updates;

[0014] Step S5. Continue to use the parallel least squares recursive algorithm to process the GNSS observation data of the nth epoch after preprocessing until the time of the n+1th epoch is updated, accumulate the value of n, and return to step S1 until the processing of the GNSS observation data is completed.

[0015] Furthermore, the time update is performed using a parallel least squares recursive algorithm and is implemented by parallel processing of the RemoveXData() function.

[0016] Furthermore, in step S2, the calculation parameters specifically include constant parameter Y and random parameter X (n) Among them, the constant parameter Y is the parameter that does not change with the epoch update, including orbital elements and force model parameters, station coordinates and inter-system deviation parameters; random parameter X (n) are the random parameters of the nth epoch, including clock error parameters, ionospheric parameters, tropospheric parameters and ambiguity parameters.

[0017] Furthermore, in step S2, the observation coefficient matrix N (n) Is a symmetric matrix; for the observation coefficient matrix N (n) Initialize the blocks and obtain the upper triangular matrix N 11(n) 、N 22(n) and matrix N 12(n) , the situation is as follows:

[0018] The initial normal equation is constructed based on the initial error equation; where,

[0019] The initial error equation is expressed as:

[0020] V=AY+BX (n) -l;

[0021] Where A is the coefficient matrix of Y, B is the coefficient matrix of X (n) The coefficient matrix of , l is the observation vector;

[0022] The initial normal equation is expressed as:

[0023]

[0024] in, P is the precision of the true observation;

[0025] For the observation coefficient matrix N (n) Perform initialization block and get:

[0026]

[0027] Among them, W1 (n) A T Pl;W 2(n) For B T Pl.

[0028] Furthermore, in step S2, the observation coefficient matrix N is obtained by a for loop. (n) The initialization block assignment includes: (n) Read N 11(n) 、N 22(n) and N 12(n), the reading process creates three child processes Thread, and encapsulates the assignment code and three for loops into functions, N 11(n) 、N 22(n) 、N 12(n) 、W 1(n) and W 2(n) Use the structure to pass the child thread and run the child thread at the same time.

[0029] Furthermore, in step S3, a recursive expression is constructed based on the least squares algorithm, specifically:

[0030]

[0031] Among them, Φ (n) is the state transfer matrix of the parameters, that is, the parameters X between epochs (n) and X (n-1) The recursive relationship between P w is the accuracy of the virtual observation value, that is, X (n) and X (n-1) process noise between parameters;

[0032] The equation coefficient information generated by the nth epoch is:

[0033]

[0034] The intermediate value parameters of the recursive expression are:

[0035]

[0036] Furthermore, in step S3, the calculation function is read and called from the BLAS library for the first time, specifically including matrix inversion and matrix multiplication.

[0037] Furthermore, in step S3, when performing the first parallel operation of writing and reading, conflicting data is found, a temporary variable is created, the pointer address and value of the temporary variable are copied using memcpy, and then the temporary variable is passed to the child thread.

[0038] Furthermore, in step S4, when performing the second parallel operation of writing and reading, for the writing phase, a new pointer variable is created to point to the data memory space, and the data is modified directly to ensure data consistency.

[0039] Furthermore, step S5, the time update of the nth epoch is completed, specifically including the block matrix N in the RemoveXData() function. 11(n) 、N 22(n) and N 12(n) The allocated space is reclaimed so that the parameter information of the nth epoch is correctly used when allocating memory space for the new block matrix during the time update phase of the n+1th epoch.

[0040] Compared with the prior art, the present invention has at least the following beneficial effects:

[0041] (1) The present invention optimizes the parallelization of the estimator, starting from the mathematical formula derivation of least squares in the navigation algorithm. Multi-threaded parallel work currently completes the CPU parallel optimization of the data matrix processing of each epoch, improving the utilization efficiency of hardware resources and the efficiency of matrix operations. The optimization effect is significant, and each epoch has an efficiency increase of about 5%-30% compared with the source program. The optimization effect varies depending on the hard disk. The longer the solution time, that is, the more epochs there are, the better the optimization effect.

[0042] (2) The present invention shortens the data processing time. Taking the experimental results of one day at an interval of 30 seconds as an example, the total time is about 5%-15% more efficient than the source program, and the calculation results remain completely consistent after optimization, and the optimization effect is very stable.

[0043] (3) This invention improves the efficiency of the computational program from another perspective of the traditional navigation algorithm, namely, from the perspective of data structure and code logic. Starting from the mathematical formula derivation of least squares in the navigation algorithm, the matrix and parameter expressions required for each step are found, and the multi-threaded parallelization of the GNSS estimator is realized. Combining the technical features of GSTAR with multi-threaded optimization, this method has many expansion directions, is universal, and is easy to port to different systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] The drawings are only for purposes of illustrating particular embodiments and are not to be considered limiting of the invention.

[0045] Figure 1 Flowchart of the least squares recursive algorithm of the present invention in GNSS measurement;

[0046] Figure 2 This is a flow chart of a traditional serial algorithm corresponding to the least squares method in the estimator of the present invention;

[0047] Figure 3 This is a flow chart after the parallelization analysis of the least squares method in the estimator of the present invention;

[0048] Figure 4 A schematic diagram of outputting the running time results of each epoch for the present invention;

[0049] Figure 5 The following is a comparison chart of serial and parallel results between the No-BLAS version and the BLAS version. DETAILED DESCRIPTION

[0050] In order to more clearly understand the above-mentioned objects, features and advantages of the present invention, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments can be combined with each other. In addition, the present invention can also be implemented in other ways different from those described herein. Therefore, the scope of protection of the present invention is not limited by the specific embodiments disclosed below.

[0051] GSTAR's estimation layer achieves highly efficient parameter estimation. Currently, the estimator layer implements two estimation methods: the square root information filter (SRIF) method for real-time applications and the least squares method (LSQ) for post-processing applications. The estimator layer consists of a single library called Engine, which includes both estimation methods. Because parameter estimation is generally the most time-consuming part, the feasibility of parallelization is analyzed by deriving the least squares method formula.

[0052] According to Figure 3 The least squares process shown in the figure first calculates the dimensions of the random parameter vector and the constant parameter vector, and then calculates the initial error equation and the normal equation:

[0053] V=AY+BX (n) -l;

[0054] Where A is the coefficient matrix of Y, B is the coefficient matrix of X (n) The coefficient matrix of , l is the observation vector;

[0055]

[0056] in, P is the precision of the true observation value, which is expressed as the inverse matrix of the observation variance.

[0057] Then the observation coefficient matrix N (n) Perform block initialization and simplify the normal equation as follows:

[0058]

[0059] Then, according to the GNSS algorithm, the calculation of each epoch requires two updates:

[0060] First time update and new error formula:

[0061] X (2) =Φ (2) X (1) +ω;

[0062] V x =Φ (2) X (1) -X (2).

[0063] The virtual error equation and virtual normal equation are obtained as follows:

[0064]

[0065]

[0066] By adding equation (1) and equation (2) together, we can get the recursive expression of the least squares recursion:

[0067]

[0068] It is not difficult to see that formula (3) only contains (2) and X (2) related.

[0069] After substitution, we can get the matrix after time update

[0070]

[0071] in:

[0072]

[0073]

[0074] Store the parameters into an intermediate file for subsequent reverse smoothing.

[0075] The second is a measurement update, which introduces data from the next epoch:

[0076] V (2) =A (2) Y+B (2) X (2) -l (2) ; Get the new error equation and normal equation:

[0077]

[0078]

[0079] Adding equations (4) and (5) yields:

[0080]

[0081] Simplified to:

[0082]

[0083] Equation (6) can be similarly viewed as the new recursive term of Equation (1).

[0084] And so on, we can finally get:

[0085]

[0086] Then we get the general formula, which is the recursive expression based on the least squares algorithm:

[0087]

[0088] Traditional least squares recursive algorithms, because they involve a large number of matrix calculations, are typically written in C language to perform large matrix inversions and matrix multiplications, resulting in extremely low efficiency. Therefore, the present invention combines the above mathematical formula derivation and utilizes the linear matrix calculation tool libraries BLAS and LAPACK to accelerate matrix operations. As can be seen from the above mathematical derivation, during each time and measurement update phase, large matrix inversions and a large number of matrix multiplications are required. By utilizing the third-level functions in the BLAS library, namely matrix-matrix operations, computer performance can be fully utilized to calculate intermediate results more quickly, thereby achieving more efficient matrix operations.

[0089] Currently, almost all programs on the market that use the least squares recursive algorithm adopt a serial operation structure, such as Figure 2 As shown, the serial process is as follows:

[0090] At each epoch, the matrix parameters must first be calculated and the matrix blocks initialized according to the parameters. After allocating space, values ​​must also be assigned to prepare for subsequent block matrix calculations.

[0091] Store the parameter values ​​and addresses of the block matrix of this epoch into the intermediate file.

[0092] Perform matrix calculations according to the least squares recursive algorithm steps and update the block matrix.

[0093] Write the updated matrix to a temporary file.

[0094] The updated block matrix is ​​written back into the original matrix to complete the update of the original matrix for this epoch.

[0095] Combining the serial program with the derivation and analysis of the above-mentioned least squares algorithm, it can be concluded that when solving large-scale GNSS data, the corresponding parameters need to be stored at the same time for the subsequent smooth reverse solution. Therefore, it is necessary to find the sharing relationship between the data to prevent data conflicts. The Thread multi-thread library can be used to achieve parallel reading and writing, thereby improving the operating efficiency. Therefore, the second technical feature of the present invention is to achieve parallel reading and writing in the traditional least squares recursive algorithm. Based on the results derived from the above formula: First, after the time is updated, the operations of matrix multiplication, inversion, and writing the new block matrix to the intermediate file are parallelized to improve the CPU processing efficiency; then, after the measurement is updated, the two operations of updating the values ​​of the block and the large matrix and writing the new matrix after the measurement update to the intermediate file are parallelized, which can also achieve efficiency improvement. In order to solve the data conflict problem, we design a structure and pass it into the child thread for write operation to ensure the stability of the data results.

[0096] The process of least squares recursive algorithm in GNSS measurement is as follows Figure 1 As shown, the least squares method is divided into two parts: time update and measurement update. The time update part involves a large number of complex operations such as matrix multiplication and matrix inversion and file storage operations, while the measurement update is more of a matrix addition and does not take much time. Therefore, the present invention starts from the RemoveXData() function, which is the most time-consuming time update, and performs parallel analysis. A multi-system non-differential and non-combined operation program with a 30s interval is selected, and a 200-station multi-satellite system is selected to ensure stable operation time and facilitate statistics.

[0097] The present invention discloses a parallel least squares recursive parameter estimation method for processing GNSS station network observation data, which uses a parallel least squares recursive algorithm to perform time update and measurement update. The measurement update is part of the least squares algorithm and is a prior art. The present invention mainly focuses on the parallel processing of time update. The measurement update involves a large number of matrix additions but does not involve a large number of matrix multiplications and inversions, so it is less time-consuming and does not require parallel processing. The time update includes the following steps:

[0098] Step S1. The GSTAR program begins running. According to the .cfg configuration file, it reads and preprocesses the GNSS observation data for the nth epoch in preparation for subsequent calculations. A parallel least squares recursive algorithm is designed based on this. The time update for each epoch is implemented in parallel using the RemoveXData() function to update the parametric equations.

[0099] Step S2: Use the parallel least squares recursive algorithm to process the GNSS observation data of the nth epoch after preprocessing. Figure 3 , the function can be divided into several stages, first calculating the calculation parameters.

[0100] Specifically, the parallel least squares recursive algorithm is implemented through the RemoveXData() function.

[0101] According to the dimension of the calculated parameters, the observation coefficient matrix N is obtained (n) ; For the observation coefficient matrix N (n) Initialize the blocks and then assign values ​​to them in parallel through a for loop.

[0102] The calculation parameters specifically include constant parameter Y and random parameter X (n) Among them, the constant parameter Y is the parameter that does not change with the epoch update, including orbital elements and force model parameters, station coordinates and inter-system deviation parameters; random parameter X (n) are the random parameters of the nth epoch, including clock error parameters, ionospheric parameters, tropospheric parameters and ambiguity parameters.

[0103] Observation coefficient matrix N (n) is a symmetric matrix; initialize the observation coefficient matrix N into blocks, specifically obtaining the upper triangular matrix N 11(n) 、N 22(n) and matrix N 12(n) , the situation is as follows:

[0104] The initial normal equation is constructed based on the initial error equation; where,

[0105] The initial error equation is expressed as:

[0106] V=AY+BX (n) -l;

[0107] Where A is the coefficient matrix of Y, B is the coefficient matrix of X (n) The coefficient matrix of , l is the observation vector.

[0108] The initial normal equation is expressed as:

[0109]

[0110] in, P is the precision of the true observations, expressed as the inverse matrix of the observation variance.

[0111] For the observation coefficient matrix N (n) Perform initialization block and get:

[0112]

[0113] Among them, W 1(n) A T Pl;W 2(n) For B T Pl.

[0114] Parallel observation coefficient matrix N through for loop (n) The initialization block assignment includes: (n) Read N 11(n) 、N 22(n) and N 12(n) ,because It is not used alone, so it is treated as a whole and stored in a variable. Since only the data at different positions of the same array are read, there is no conflict between readings. Since the matrix dimension is large, the serial reading for loop takes a long time, so the reading process creates three sub-processes Thread, and the assignment code is encapsulated into a function to establish three for loops. 11(n) 、N 22(n) 、N 12(n) 、W 1(n) and W 2(n) Parameters such as matrix parameters are passed to child threads using structures, and child threads are run simultaneously to improve efficiency. The more matrix parameters there are and the larger the for loop is, the better the parallel effect.

[0115] Step S3. Perform the first write and read operations in parallel; wherein the first write includes storing the equation of the recursive expression and storing the equation coefficient information generated by the nth epoch according to the recursive expression through time update; the recursive expression is used to represent the calculation parameters and the observation value coefficient matrix N (n) The functional relationship between the initialization blocks.

[0116] The recursive expression is constructed based on the least squares algorithm, specifically:

[0117]

[0118] Among them, Φ (n) is the state transfer matrix of the parameters, that is, the parameters X between epochs (n) and X (n-1) The recursive relationship between P w is the accuracy of the virtual observation value, that is, X (n) and X (n-1) Process noise between parameters.

[0119] The equation coefficient information generated by the nth epoch is:

[0120]

[0121] This information will be calculated during the time update and then stored in the intermediate file, so a file write operation will be generated in this function.

[0122] The first reading includes calculating the equation coefficient information generated by the nth epoch by calling the calculation function to obtain the observation coefficient matrix after the time update This includes the equation coefficient information generated at the nth epoch after the time update and the intermediate value parameters of the recursive expression after the time update.

[0123] The intermediate value parameters of the recursive expression are:

[0124]

[0125] The first read calls computational functions from the BLAS library, specifically matrix inversion and matrix multiplication.

[0126] After feasibility analysis, the first write and read can be processed in parallel.

[0127] During the first parallel write and read operation, if direct parallel processing is used, that is, a shallow copy of the parameters is performed, the program will encounter a sharing conflict in accessing some data, resulting in a write error and abnormal final results. The solution is to identify the conflicting data, create a temporary variable, copy the pointer address and value of the temporary variable using memcpy, and then pass the temporary variable to the child thread to ensure that the data is not changed during the child thread's write process, ensuring the consistency of the results.

[0128] Step S4. Perform a second parallel operation of writing and reading; wherein, the second writing includes storing the equation coefficient information generated by the nth epoch after the time update and the intermediate value parameters of the recursive expression after the time update.

[0129] The second read includes the updated observation coefficient matrix according to the intermediate value parameters of the recursive expression after time update. Divide into blocks, and then reassign Update to facilitate calculation of the next epoch.

[0130] When performing the second parallel write and read operation, the file writing and computation operations can also be parallelized. During the write phase, due to the limited size of intermediate files, file partitioning is necessary when processing large amounts of data. However, when analyzing data consistency, the file size parameter is a global parameter that changes continuously with epoch updates. Unaddressed data update issues encountered in child threads can lead to global data initialization, resulting in incorrect file partitioning and affecting the correctness of the results. To synchronously modify global variables in child threads, the necessary measures are to create a pointer variable pointing to the data memory space and then make modifications to ensure data consistency.

[0131] Step S5: Finalize the time update of the nth epoch.

[0132] Continue to use the parallel least squares recursive algorithm to process the GNSS observation data of the nth epoch after preprocessing until the time of the n+1th epoch is updated, accumulate the value of n, and return to step S1 until the processing of the GNSS observation data is completed.

[0133] Finally, the running time results of each epoch are output as follows Figure 4 As shown; the average time comparison between the parallel operation of the present invention and the traditional serial operation is shown in Table 1.

[0134] Table 1 Comparison of average time between the parallel operation of the present invention and the traditional serial operation

[0135] Average serial running time Average serial running time Percentage increase 6.583369s 4.589524s 30.29％

[0136] The serial and parallel results of the No-BLAS version and the BLAS version are compared. Figure 5 shown.

[0137] Compared with the existing technology, the present invention optimizes the parallelization of the estimator, starting from the mathematical formula derivation of least squares in the navigation algorithm. Multi-threaded parallel work currently completes the CPU parallel optimization of each epoch data matrix processing, improves the utilization efficiency of hardware resources and the efficiency of matrix operations, and has a significant optimization effect. Each epoch has an efficiency increase of about 5%-30% compared with the source program. The optimization effect varies depending on the hard disk. The longer the solution time, that is, the more epochs, the better the optimization effect. The present invention shortens the data processing time. Taking 30-second intervals and one-day experimental results as an example, the total time has an efficiency increase of about 5%-15% compared with the source program, and the calculation results after optimization remain completely consistent, and the optimization effect is very stable. The technical characteristics of GSTAR are combined with multi-threaded optimization, with more expansion directions, versatility, and convenient porting to different systems.

[0138] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by any technician familiar with this technical field within the technical scope disclosed by the present invention should be covered by the scope of protection of the present invention.

Claims

1. A parallel least squares recursive parameter estimation method for GNSS station network observation data processing, characterized in that: A parallel least squares recursive algorithm is used to perform time update and measurement update; wherein the time update includes the following steps: Step S1. The GSTAR program starts running and reads the GNSS observation data of the nth epoch and preprocesses it according to the .cfg configuration file; Step S2. Processing the pre-processed GNSS observation data of the nth epoch using a parallel least squares recursive algorithm to calculate the calculation parameters; According to the dimension of the calculated parameters, the observation coefficient matrix N is obtained (n) ; For the observation coefficient matrix N (n) Initialize the blocks and then assign values ​​to them in parallel through a for loop; Step S3. Perform the first write and read operations in parallel; wherein the first write includes storing the equation of the recursive expression and storing the equation coefficient information generated by the nth epoch according to the recursive expression through time update; the recursive expression is used to represent the calculation parameters and the observation value coefficient matrix N (n) The functional relationship between the initialization blocks; The first reading includes calculating the equation coefficient information generated by the nth epoch by calling the calculation function to obtain the observation coefficient matrix after the time update Including the equation coefficient information generated by the nth epoch after the time update and the intermediate value parameters of the recursive expression after the time update; Step S4. Performing a second parallel write and read operation; wherein the second write operation includes storing equation coefficient information generated at the nth epoch after the time update and intermediate value parameters of the recursive expression after the time update to an intermediate file; the second write operation also includes determining the size of the intermediate file and dividing the intermediate file; The second read includes the updated observation coefficient matrix according to the intermediate value parameters of the recursive expression after time update. Divide into blocks, and then reassign Make updates; Step S5. Continue to use the parallel least squares recursive algorithm to process the GNSS observation data of the nth epoch after preprocessing until the time of the n+1th epoch is updated, accumulate the value of n, and return to step S1 until the processing of the GNSS observation data is completed.

2. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 1, characterized in that: Time update is performed using a parallel least squares recursive algorithm and is implemented through parallel processing of the RemoveXData() function.

3. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 2, characterized in that: Step S2: Calculate the parameters including constant parameter Y and random parameter X. (n) Among them, the constant parameter Y is the parameter that does not change with the epoch update, including orbital elements and force model parameters, station coordinates and inter-system deviation parameters; random parameter X (n) are the random parameters of the nth epoch, including clock error parameters, ionospheric parameters, tropospheric parameters and ambiguity parameters.

4. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 3, characterized in that: Step S2, observation coefficient matrix N (n) Is a symmetric matrix; for the observation coefficient matrix N (n) Initialize the blocks and obtain the upper triangular matrix N 11(n) 、N 22(n) and matrix N 12(n) ,as follows: The initial normal equation is constructed based on the initial error equation; where, The initial error equation is expressed as: V=AY+BX (n) -l; Where B is the coefficient matrix of Y, B is the coefficient matrix of X (n) The coefficient matrix of , l is the observation vector; The initial normal equation is expressed as: in, P is the precision of the true observation; For the observation coefficient matrix N (n) Perform initialization block and get: Among them, W 1(n) A T ol;W 2(n) For B T Pl.

5. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 4, characterized in that: Step S2, through the for loop and the observation coefficient matrix N (n) The initialization block assignment includes: (n) Read N 11(n) 、N 22(n) and N 12(n) , the reading process creates three child processes Thread, and encapsulates the assignment code and three for loops into functions, N 11(n) 、N 22(n) 、N 12(n) 、W 1(n) and W 2(n) Use the structure to pass the child thread and run the child thread at the same time.

6. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 5, characterized in that: In step S3, a recursive expression is constructed based on the least squares algorithm, specifically: Among them, Φ (n) is the state transfer matrix of the parameters, that is, the parameters X between epochs (n) and X (n-1) The recursive relationship between P w is the accuracy of the virtual observation value, that is, X (n) and X (n-1) process noise between parameters; The equation coefficient information generated by the nth epoch is: The intermediate value parameters of the recursive expression are:

7. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 6, characterized in that: Step S3: First read and call calculation functions from the BLAS library, specifically including matrix inversion and matrix multiplication.

8. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 7, characterized in that: Step S3, when performing the first parallel operation of writing and reading, find the conflicting data, create a temporary variable, use memcpy to copy the pointer address and value of the temporary variable, and then use the temporary variable to pass it to the child thread.

9. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 8, characterized in that: Step S4, when performing the second parallel operation of writing and reading, for the writing phase, a new pointer variable is created to point to the data memory space, and the modification is performed directly to ensure data consistency.

10. The parallel least squares recursive parameter estimation method for GNSS station network observation data processing according to claim 9, characterized in that: Step S5, finishing the time update of the nth epoch, specifically including the block matrix N in the RemoveXData() function 11(n) 、N 22(n) and N 12(n) The allocated space is reclaimed so that the parameter information of the nth epoch is correctly used when allocating memory space for the new block matrix during the time update phase of the n+1th epoch.

Images