A Parallel Approach for Multi-GNSS Ultra-Rapid Orbit Determination

: Global Navigation Satellite System (GNSS) ultra-rapid orbit is critical for geoscience and real-time engineering applications. To improve the computational efﬁciency and the accuracy of predicted orbit, a parallel approach for multi-GNSS ultra-rapid orbit determination is proposed based on Message Passing Interface (MPI)/Open Multi Processing (OpenMP). This approach, compared with earlier efﬁcient methods, can improve the efﬁciency of multi-GNSS ultra-rapid orbit solution without changing the original observation data and retaining the continuity and consistency of the original parameters to be estimated. To obtain high efﬁciency, three steps are involved in the approach. First and foremost, the normal equation construction is optimized in parallel based on MPI. Second, equivalent reduction of the estimated parameters is optimized using OpenMP parallel method. Third, multithreading is used for parallel orbit extrapolation. Thus, GNSS ultra-rapid orbit determination is comprehensively optimized in parallel, and the computation efﬁciency is greatly improved. Based on the data from MGEX and IGS stations, experiments are carried out to analyze the performance of the proposed approach in computational efﬁciency, accuracy and stability. The results show that the approach greatly improves the efﬁciency of satellite orbit determination. It can realize 1-h update frequency for the multi-GNSS ultra-rapid orbit determination using 88 stations with four-system observations. The accuracy of the GPS, GLONASS, Galileo and BDS ultra-rapid orbit with 1-h update frequency using the parallel approach is approximately 33.4%, 31.4%, 40.1% and 32.8% higher than that of the original orbit, respectively. The root mean squares (RMS) of GPS, GLONASS, Galileo and BDS predicted orbit are about 3.2 cm, 5.1 cm, 5.6 cm and 11.8 cm. Moreover, the orbit provided by the proposed method has a better stability. The precision loss of all parallel optimization can be negligible and the original correlation between the parameters is fully preserved. are selected for the processing. The approach is implemented by Fortran and intel MPI/OpenMP parallel library. The experiments are carried out on a Linux computer cluster.


Introduction
Global Navigation Satellite System (GNSS) ultra-rapid orbit is the basis of geodetic research and real-time applications, such as real-time atmosphere monitoring [1,2], realtime PPP [3][4][5], real-time satellite clock estimation [6], and real-time orbit determination of LEO satellite [7]. With the development of Galileo and implementation of BDS global service plan, the number of satellites in the constellation has increased dramatically to more than 130 [8]. Currently, the International GNSS Service (IGS), the Multi-GNSS Experiment and Pilot Project (MGEX) and the International GNSS Monitoring and Assessment System (iGMAS) have provided more than 500 global distributed stations. A combined multi-GNSS solution can effectively reduce systematic errors and improve the accuracy and reliability of satellite orbit. However, on the other hand, the increase in the number of satellite constellations and stations makes the processing model more complicated and results in larger computational burden [9], which decrease the computational efficiency of The algorithm of precise orbit determination (POD) mainly uses dynamics and reduced dynamics methods [14][15][16][17][18]. The methods above are first applied to GPS satellite orbit determination and are widely used in GLONASS, Galileo, BDS and other satellite systems, as well as the combination of them with the development of GNSS [19][20][21][22][23][24]. Based on this, the predicted orbit of ultra-rapid products is generally obtained by integration of initial orbit and force model parameters [25]. The accuracy of orbit for real-time applications is directly related to the accuracy of satellite orbit initial values and the length of extrapolation time. In order to obtain accurate solutions and achieve higher computational efficiency, many scholars have optimized the theories and methods of orbit determination by dividing subnets, grouping combined systems, dividing arc length of observation, and reducing observations. , for the first time, proposed the method of merging the past three consecutive single-day arcs to orbit determination in 1999 [26]. Yao (2008) further used a dynamic smoothing algorithm to smooth the orbit which can improve the accuracy and reliability of the comprehensive algorithm solution. In order to improve the efficiency of ultra-rapid orbit determination [27], Lou (2008) deeply studied the problem of ambiguity processing in real-time orbit determination and proposed a sliding window method [28]. Through the superposition of short-arc normal equations, the orbit can be calculated with higher efficiency, and the accuracy is improved when the ambiguities are fixed. In addition, TUM used a two-step approach for generating Galileo and QZSS products [29]. Furthermore, a three-step approach was applied by WHU to orbit product calculation: GPS and GLONASS orbit and ERPs are solved first, then station coordinates, zenith tropospheric delay (ZTD) and receiver clock are estimated, and the orbit of other satellites are solved in the last step. Then, a block recursive least squares method was proposed to improve the accuracy of GNSS ultra-rapid orbits. The update frequency of ultra-rapid orbits was increased from 6-h to 1-h, and the accuracy was significantly improved [30]. Li (2019) realized 1-h update multi-GNSS ultra-rapid orbit determination by multi-threading and changed observation intervals [31]. The accuracy of predicted orbit was improved and the kinematic PPP solutions could achieve a shorter time to first fix and higher positioning accuracy.
To promote the accuracy of the predicted multi-GNSS orbit, we proposed a new approach for multi-GNSS hourly ultra-rapid orbit determination. This approach improves the efficiency of multi-GNSS ultra-rapid orbit solution without changing the original observation data or processing interval, and retains the continuity and consistency of the original parameters to be estimated using the multi-process and multi-thread method based on Message Passing Interface (MPI)/Open Multi Processing (OpenMP). MPI is a message passing library specification which provides a powerful and portable way for expressing parallel programs. OpenMP is a set of compiler directives as well as an API for programs that provides support for parallel programming in shared-memory environments. Therefore, MPI/OpenMP can achieve higher accuracy for multi-GNSS ultra-rapid orbit as the initial orbit is precise, but the length of extrapolation time is shortened significantly. To realize that, first and foremost, the construction of the normal equation for initial value of satellite orbit is divided into multi-arc parallel processing by the multi-process method based on MPI, and each process achieves consistency through processing the arc data and interacting matrix. In addition, the multi-thread method is used in parameter elimination and orbit extrapolation to further improve the efficiency. In this approach, the GPS+GLONASS+Galileo+BDS observations are processed together in one common procedure with the double-differenced network model. The promotion of computational efficiency, accuracy and stability is validated based on the data from MGEX and IGS stations.
In the following sections, the principle and main steps of the modified parallel multi-GNSS obit determination approach will be introduced. Then, some processing strategies are presented. Afterward, using MGEX/IGS data and related products, the improvement of computational efficiency of the new approach is analyzed. At last, the accuracy of multi-GNSS obit using the proposed approach is evaluated, and the results indicate that the proposed parallel approach is viable.

Methodology
The double-differenced network solution is used in the proposed parallel multi-GNSS ultra-rapid orbit determination. After the data preprocessing, the basic procedure of multi-GNSS ultra-rapid orbit determination mainly includes three steps. First, the double-differenced normal equations are constructed based on un-differenced observation equations. Second, the normal equations are solved after eliminating the redundant parameters. Finally, the predicted orbits are obtained by extrapolation on the basis of initial orbit and model parameters using orbit integration. We optimize the above core processing to achieve higher computational efficiency using a parallel computing method based on MPI/OpenMP. Since the original observational model and stochastic model have not been changed, the approach retains the original correlation of the estimated parameters and the consistency of the ambiguity parameters.

Multi-Arc Parallel Normal Equations Construction Based on MPI
In the multi-GNSS ultra-rapid orbit determination, the sequential batch processing method is adopted. The double-differenced carrier observation equation can be expressed as: ambiguity, λ is the wavelength of ionospheric-free observations, and ε Φ is the observation noise. The observation equation (1) is linearized and expressed as: where B = ∂ρ ij kl (t)/∂X r is the partial derivative of the measured value with respect to the state parameters, Ψ(t, t 0 ) is the transition matrix which can be obtained by solving the variational equation, X r0 , δX r are the approximate coordinates of the station and their corrections, respectively, X s (t) is the state of satellite at time t, which is obtained by solving the differential equation of satellite motion, ρ ij kl0 (t, X r0 , X s (t)) is the approximate value of the double-differenced satellite-ground geometric distance at time t, and δX s 0 is the correction of the satellite state at the reference time. Other symbols have the same meaning as in (1). Based on the observation (2), the GNSS satellite reduce-dynamic precise orbit determination (RDPOD) can be performed. Table 2 shows the main parameter settings of the RDPOD strategy. In this paper, we have no a priori box-wing model used for the orbit determination, which may generate orbit errors, especially for Galileo satellites due to the elongated shape of the satellite buses. The a priori box-wing model may be a solution for twice-per-revolution empiri-cal parameters from the CODE's new solar radiation pressure model are used for GNSS orbit determination [32]. In the paper, we adopt the one-step orbit determination of the GPS, GLONASS, Galileo and BDS in order to study that the parallel algorithm improves the efficiency of satellite orbit prediction, which does not affect the computational efficiency of the parallel algorithm. For BDS, no official metadata for a priori box-wing models are available. However, estimated parameters can be still used [33]. Equation (2) is the basic observation equation for ultra-rapid orbit determination. In order to determine the ultra-rapid orbit, the traditional sequential batch processing method traverses the observation data in all arcs by epoch to construct the normal equations. In this study, a parallel method based on MPI is used in this process. The main idea is to establish the multi-arc normal equation in the orbit determination by adopting the multiprocessing method. Assuming that the single-day data is divided into three arcs, the length of each arc is 8 h. Each sub-process processes the measured data of the current arc and then calculates the observation equation coefficients. Based on the current time period, sub-normal equations including initial satellite state, force model parameters, earth rotation parameters, atmospheric parameters, station coordinates and ambiguity parameters are generated. Each process sends the sub-normal equation matrices and auxiliary parameters to the main process to obtain the total normal equations.
The parameter sequence of each period is consistent with that of the entire original period. The mathematical model and the parameters have not changed, so there is no additional conversion and connection. Supposing that the entire arc is divided into k sub-arcs, the least squares problem can be solved by accumulating the normal equations. The sub-normal equations and the entire normal equation can be expressed as: where N k and W k are the related matrix of sub-normal equations and δx is unknown parameter vector. Based on this entire normal equation, the elimination of parameters and least squares (LS) estimation can be then carried out. It is worth noting that the number of time periods are designed according to the computer hardware and time requirements.

Parallel Parameter Elimination Based on OpenMP
The parameters to be estimated in the GNSS observation equation of POD include station coordinate parameters, satellite orbit parameters, tropospheric parameters, ambiguity parameters, and earth rotation parameters. Usually, only some of these parameters are of concern to us, and the other parameters can be eliminated to improve the computational efficiency [41]. Equivalent elimination is a classic algorithm that can reduce the dimensionality of normal equations without a precision loss (Yao 2004). The error equation can be expressed as: where L is the observation vector of dimension m, A and B are the design matrices of dimension m × (n − r) and m × r,x 1 andx 2 are unknown parameters of dimension n − r and r, v is the residual error matrix of dimension m, n is the number of unknowns, m is the number of observations, and P is the weight matrix of dimension m × m.
The parameters in the equation are decomposed into x 1 x 2 as needed, taking into account the correlation between the observations. The normal equations are: where If only x 1 is concerned, x 2 can be eliminated in order to reduce the dimension of the normal equation; 1 22 , and N −1 22 is the inverse matrix of N 22 . The left and right sides of (8) are multiplied by the matrix as follows: Then the normal equation after elimination is: The solution of the normal equation is: (13) Equation (12) is substituted into (5): When the number of observation equations and parameter dimensions is very large, the parameter elimination is time-consuming. In order to further improve computational efficiency, parallel optimization using OpenMP technology for parameter elimination is adopted.
When the initial orbit parameters and the force model parameters are determined by using LS estimation, the orbit of the entire arc will be updated in the order of epoch and satellite by using orbit integration. As with the the above-mentioned, we use the parallel orbit integration method to update the final orbit sequence.

Parallel Orbit Determination Process
The parallel approach proposed for the ultra-rapid orbit determination of the doubledifferenced network solution mainly adopts the MPI/OpenMP hybrid parallel programming model. The optimized orbit determination process is shown in Figure 1. After reading observation files, navigation ephemeris files, and some auxiliary model correction files, data preprocessing includes gross error detection, basic model correction, forming double-differenced baseline, and floating solution. Then, the entire time period to be calculated is divided into n segments, and multiple processes are developed based on MPI to read and construct sub-normal equations in parallel. Each process read from the first epoch to the last epoch of the current period so that the consistency of the current arc parameter sequence can be maintained. The above multiple sub-normal equations are accumulated to obtain the total normal equation. Next, in order to further improve the efficiency of orbit determination, the OpenMP multithreading model is used for the parameter elimination process. Finally, after obtaining the precise dynamic parameters by the LS estimation, the same multi-thread model is used to realize the parallel orbit integration to obtain the final satellite POD results. The entire parallel optimization does not destroy the correlation between the parameters without the ambiguity segmentation problem. The input and output of orbit determination do not need to be adjusted greatly, which is easy to implement. It can be applied to single satellite POD as well as multi-GNSS POD.
ciency of orbit determination, the OpenMP multithreading model is used for the parameter elimination process. Finally, after obtaining the precise dynamic parameters by the LS estimation, the same multi-thread model is used to realize the parallel orbit integration to obtain the final satellite POD results. The entire parallel optimization does not destroy the correlation between the parameters without the ambiguity segmentation problem. The input and output of orbit determination do not need to be adjusted greatly, which is easy to implement. It can be applied to single satellite POD as well as multi-GNSS POD.

Computation and Comparison
In order to verify the feasibility and effectiveness of the proposed GNSS parallel ultra-rapid orbit determination approach, multiple sets of experiments are carried out. The precision loss of the parallel POD, the improvement of computational efficiency, and the orbit accuracy after using the parallel approach are analyzed in detail.
Two sets of tracking stations from IGS/MGEX networks around the world are selected to analyze the accuracy and efficiency of multi-GNSS parallel ultra-rapid POD. The networks are shown in Figure 2. Relevant data from day of year (DOY) 182 to 189, 2018 are selected for the processing. The approach is implemented by Fortran and intel MPI/OpenMP parallel library. The experiments are carried out on a Linux computer cluster.
precision loss of the parallel POD, the improvement of computational efficiency, and the orbit accuracy after using the parallel approach are analyzed in detail.
Two sets of tracking stations from IGS/MGEX networks around the world are selected to analyze the accuracy and efficiency of multi-GNSS parallel ultra-rapid POD. The networks are shown in Figure 2. Relevant data from day of year (DOY) 182 to 189, 2018 are selected for the processing. The approach is implemented by Fortran and intel MPI/OpenMP parallel library. The experiments are carried out on a Linux computer cluster.

Impact of Parallel Approach on the Orbit Solution
Before further analyzing the efficiency and accuracy of the parallel approach, it is necessary to analyze the accuracy of the non-parallel GNSS POD method. The observation data from DOY 182 to 188, 2018 of the above 131 IGS/MGEX stations are used to carry out the multi-GNSS POD. Analysis of the errors in the radial, along and cross directions of the multi-GNSS orbits using the non-parallel method is performed. As a result, the along track direction accuracy is the lowest, followed by the cross track, and the radial is the highest. Except for slightly larger fluctuations in the BDS in the sense of the orbit RMS, the other systems have little change and the accuracy is relatively stable. The accuracy of each system is further analyzed, which is shown in Table 3. It can be seen from the table that the 1D mean RMS (1 = ( + + ) 3 ⁄ ) of GPS, GLONASS, Galileo and BDS are 1.3 cm, 2.7 cm, 2.7 cm, and 5.4 cm, respectively, which are at the same accuracy level as reported by international institutions [42]. The results show the correctness of the non-parallel POD processing, which lays the foundation for the proposed parallel POD.

Impact of Parallel Approach on the Orbit Solution
Before further analyzing the efficiency and accuracy of the parallel approach, it is necessary to analyze the accuracy of the non-parallel GNSS POD method. The observation data from DOY 182 to 188, 2018 of the above 131 IGS/MGEX stations are used to carry out the multi-GNSS POD. Analysis of the errors in the radial, along and cross directions of the multi-GNSS orbits using the non-parallel method is performed. As a result, the along track direction accuracy is the lowest, followed by the cross track, and the radial is the highest. Except for slightly larger fluctuations in the BDS in the sense of the orbit RMS, the other systems have little change and the accuracy is relatively stable. The accuracy of each system is further analyzed, which is shown in Table 3. It can be seen from the table that the 1D mean RMS (1DRMS = σ 2 R + σ 2 A + σ 2 C /3) of GPS, GLONASS, Galileo and BDS are 1.3 cm, 2.7 cm, 2.7 cm, and 5.4 cm, respectively, which are at the same accuracy level as reported by international institutions [42]. The results show the correctness of the non-parallel POD processing, which lays the foundation for the proposed parallel POD. Table 3. Accuracy statistics of multi-GNSS orbit using the non-parallel method (cm). It is necessary to further analyze the impact of the parallel approach on the covariance matrix and solution of the normal equations, as well as the precision loss of orbit. Based on MPI, the entire 480 epochs are divided into 2 periods. Statistics of the difference between the matrix elements of the parallel method and those of the non-parallel method is shown in Table 4. It can be seen from these figures that the error impact of the parallel approach on the covariance matrix and unknown parameters are approximately 10 −13 and 10 −6 , respectively.

Satellite
Using the same set of data, the normal equation coefficients of the above-mentioned parallel approach and the normal equation coefficients of the non-parallel approach are saved. After the procedures of parameter estimation, orbit update, and result output, the two sets of orbits are obtained. The comparison of the two sets of orbits show that their difference is below the millimeter level, which indicates that the precision loss caused by parallel method can be negligible.

Improvement of Computational Efficiency of Parallel Approach
In order to verify the characteristics of the parallel approach in terms of efficiency, the parallel approach is compared with the non-parallel method. The efficiency improvement of the approach in many aspects is further analyzed, such as multi-system, number of threads, and number of processes. In addition, the computational time changing with the number of stations is analyzed, which reflects the advantage of the parallel approach in the improvement of efficiency and provides a reference for the selection of the number of stations in the 1-h update frequency ultra-rapid POD.

Improvement of Efficiency with Different Number of Satellite Systems
The computational efficiency of GPS, GPS+GLONASS and GPS+GLONASS+Galileo +BDS parallel POD are compared to analyze the influence of the number of satellites. Two processes and ten threads are used in the parallel scheme. The statistical results of computational time are shown in Table 5, and the comparison between non-parallel and parallel approach is shown in Figure 3. Table 5. Time-consuming statistical results of GNSS parallel and non-parallel POD (unit: min).

Satellite System
Non-Parallel POD

Parallel Initial Orbit Estimation
Parallel POD It can be seen from Table 5 that the more satellite systems that are involved, the longer it takes for POD. After using MPI parallelism, the calculation time of this process is significantly reduced, indicating that it is reasonable to optimize this process in parallel. It can be seen from Figure 3 that the efficiency of GPS, GPS+GLONASS and GPS+GLONASS+Galileo+BDS POD have been obviously improved by the proposed parallel approach. As the number of satellites increases, the efficiency of the parallel POD has It can be seen from Figure 3 that the efficiency of GPS, GPS+GLONASS and GPS+ GLONASS+Galileo+BDS POD have been obviously improved by the proposed parallel approach. As the number of satellites increases, the efficiency of the parallel POD has been improved more significantly, and the efficiency of the four-system combined POD has been improved the most. The parameters to be estimated and the dimension of the equation both increase nonlinearly, so the advantages of the parallel POD are more obvious.

Improvement of Efficiency with Different Number of Threads
The computational efficiency of the parallel parameter elimination process with a different number of threads based on OpenMP is further analyzed. According to the hardware configuration of the platform, the number of threads is sequentially increased from 2 to 16 at an interval of 2, and the computational time for the parameter elimination process of each parallel scheme is compared, which is shown in Figure 4. It can be seen from Figure 3 that the efficiency of GPS, GPS+GLONASS and GPS+GLONASS+Galileo+BDS POD have been obviously improved by the proposed parallel approach. As the number of satellites increases, the efficiency of the parallel POD has been improved more significantly, and the efficiency of the four-system combined POD has been improved the most. The parameters to be estimated and the dimension of the equation both increase nonlinearly, so the advantages of the parallel POD are more obvious.

Improvement of Efficiency with Different Number of Threads
The computational efficiency of the parallel parameter elimination process with a different number of threads based on OpenMP is further analyzed. According to the hardware configuration of the platform, the number of threads is sequentially increased from 2 to 16 at an interval of 2, and the computational time for the parameter elimination process of each parallel scheme is compared, which is shown in Figure 4. It can be seen from Figure 4 that as the number of threads increases, the computational time of parameter elimination in GNSS POD gradually decreases and the scheme with 16 threads has the shortest computational time. However, the speed of efficiency improvement gradually decreases, most likely because when the number of threads is close to or greater than the number of cores, queuing round calculations occur, and the opening and ending of threads take more time. Therefore, further increasing the number It can be seen from Figure 4 that as the number of threads increases, the computational time of parameter elimination in GNSS POD gradually decreases and the scheme with 16 threads has the shortest computational time. However, the speed of efficiency improvement gradually decreases, most likely because when the number of threads is close to or greater than the number of cores, queuing round calculations occur, and the opening and ending of threads take more time. Therefore, further increasing the number of threads beyond the hardware core is of little significance to the improvement of parallel efficiency.

Improvement of Efficiency with Different Number of Processes
In order to further reflect the improvement of efficiency by the parallel POD approach, the influence of the number of processes on efficiency is analyzed. Multi-process is mainly used to construct normal equations in parallel. All schemes use the same number of threads, the same processing strategy, and select the number of processes as the only control variable to analyze the effect of increasing the number of processes on the computational time. The time consumption is shown in Table 6, and the speedup is shown in Figure 5. As can be seen from Figure 5 and Table 6, compared to a single-process, the speedup ratios of two-process, three-process, four-process, and five-process are 1.788, 2.525, 3.102, and 3.406, respectively. As the number of processes increases, the speedup ratio increases. When five processes are used, the speedup in POD reaches the maximum. The speedup shows a nonlinear trend probably because the cost of opening and closing processes and inter-process communication increases as the number of processes increases, and depends on the proportion of parallel programs to the total program. In general, if the hardware can support more processes, the computational efficiency of POD can be further improved, which is similar to multithreading.
In order to further reflect the improvement of efficiency by the parallel POD approach, the influence of the number of processes on efficiency is analyzed. Multi-process is mainly used to construct normal equations in parallel. All schemes use the same number of threads, the same processing strategy, and select the number of processes as the only control variable to analyze the effect of increasing the number of processes on the computational time. The time consumption is shown in Table 6, and the speedup is shown in Figure 5.  As can be seen from Figure 5 and Table 6, compared to a single-process, the speedup ratios of two-process, three-process, four-process, and five-process are 1.788, 2.525, 3.102, and 3.406, respectively. As the number of processes increases, the speedup ratio increases. When five processes are used, the speedup in POD reaches the maximum. The speedup shows a nonlinear trend probably because the cost of opening and closing processes and inter-process communication increases as the number of processes increases, and depends on the proportion of parallel programs to the total program. In general, if the hardware can support more processes, the computational efficiency of POD can be further improved, which is similar to multithreading.

Improvement of Efficiency with Different Number of Stations
Observation data from 40 to 110 IGS and MGEX stations are selected to analyze the improvement of efficiency of parallel multi-GNSS POD. The results are shown in Figure  6. It can be seen from Figure 6 that the parallel POD is obviously more efficient than the non-parallel POD. As the number of IGS/MGEX stations increases, the serial calculation time of non-parallel and parallel POD both gradually increases. The speedup ratio of parallel POD increases with the increase in the number of stations. The efficiency of POD using 80 to 90 stations has been significantly improved. Moreover, the parallel POD of four satellite system can be controlled within one hour, which provides the possibility to increase the update frequency to 1-h for multi-GNSS parallel POD.

Analysis of the Accuracy of Parallel Ultra-Rapid POD
Aiming at the problem that the predicted orbit accuracy of GNSS ultra-rapid products becomes worse as the extrapolation time is longer, we apply the parallel approach to the multi-GNSS ultra-rapid POD. Based on the results of the improvement of efficiency with different number of stations in the previous section, 88 IGS/MGEX stations are selected for multi-GNSS ultra-rapid POD. The distribution of stations is shown in Figure 2. 3.3.1. Improvement of Accuracy with Increased Update Frequency by the Parallel POD The update frequency of the original ultra-rapid orbit is 6 h, and the update fre- It can be seen from Figure 6 that the parallel POD is obviously more efficient than the non-parallel POD. As the number of IGS/MGEX stations increases, the serial calculation time of non-parallel and parallel POD both gradually increases. The speedup ratio of parallel POD increases with the increase in the number of stations. The efficiency of POD using 80 to 90 stations has been significantly improved. Moreover, the parallel POD of four satellite system can be controlled within one hour, which provides the possibility to increase the update frequency to 1-h for multi-GNSS parallel POD.

Analysis of the Accuracy of Parallel Ultra-Rapid POD
Aiming at the problem that the predicted orbit accuracy of GNSS ultra-rapid products becomes worse as the extrapolation time is longer, we apply the parallel approach to the multi-GNSS ultra-rapid POD. Based on the results of the improvement of efficiency with different number of stations in the previous section, 88 IGS/MGEX stations are selected for multi-GNSS ultra-rapid POD. The distribution of stations is shown in Figure 2.

Improvement of Accuracy with Increased Update Frequency by the Parallel POD
The update frequency of the original ultra-rapid orbit is 6 h, and the update frequency can be increased to 1 h after applying the parallel POD. The accuracy of orbit with 6-h update frequency and 1-h update frequency are compared in this section. Observation data from DOY 182 to 189, 2018 are selected for multi-GNSS ultra-rapid POD. The accuracy of GPS, GLONASS, Galileo, and BDS predicted orbit with parallel and non-parallel POD are analyzed, which are shown in Figures 7-10.

Analysis of the Accuracy of Parallel Ultra-Rapid POD
Aiming at the problem that the predicted orbit accuracy of GNSS ultra-rapid products becomes worse as the extrapolation time is longer, we apply the parallel approach to the multi-GNSS ultra-rapid POD. Based on the results of the improvement of efficiency with different number of stations in the previous section, 88 IGS/MGEX stations are selected for multi-GNSS ultra-rapid POD. The distribution of stations is shown in Figure 2.

Improvement of Accuracy with Increased Update Frequency by the Parallel POD
The update frequency of the original ultra-rapid orbit is 6 h, and the update frequency can be increased to 1 h after applying the parallel POD. The accuracy of orbit with 6-h update frequency and 1-h update frequency are compared in this section. Observation data from DOY 182 to 189, 2018 are selected for multi-GNSS ultra-rapid POD. The accuracy of GPS, GLONASS, Galileo, and BDS predicted orbit with parallel and non-parallel POD are analyzed, which are shown in Figures 7-10.        After the update frequency is increased to 1-h by parallel approach, the accuracy of predicted orbit is higher and more stable at different extrapolation time. The statistical results of the accuracy of GPS/GLONASS/Galileo/BDS orbit are shown in Table 7. Selecting the precise orbit of IGS and CODE analysis center as reference, the along, normal, radial and 1D mean RMS values of predicted orbit of 1-h update frequency and 6-h update frequency are compared. It can be seen from the table that the accuracy of GPS, GLONASS, Galileo, and BDS 1-h updated orbits are improved in normal, radial and along direction, and the radial direction is improved the most compared to the 6-h updated orbit. The RMS of the predicted orbit of GPS, GLONASS, Galileo and BDS are 33.4%, 31.4%, 40.1% and 32.8% higher than those of the original 6-h results, respectively.

Evaluation of Predicted Orbit Accuracy of Ultra-Rapid POD by Parallel Approach
In order to further analyze the performance of the parallel approach, the accuracy of the ultra-rapid predicted orbit is evaluated. The GPS/GLONASS precise orbit of IGS and the Galileo/BDS precise orbit of CODE are selected as reference values. The experiment includes three groups. The first one uses the proposed parallel POD approach with a 1-h updated result (denoted as SDU-1h). The second one uses the non-parallel POD approach with a 6-h updated result (denoted as SDU-6h), and the third one is the 1-h updated result of Wuhan University analysis center (denoted as WHU-1h). The statistical results of accuracy are shown in Figure 11. From Figure 11, it can be seen that as the extrapolation time increases, the orbit accuracy of SDU-6h decreases in the orbit sequence. The phenomenon that the orbit accuracy fluctuates with extrapolation time does not exist in WUH-1h and SDU-1h. The accuracy of WUH-1h is better than that of SDU-6h result, and it is relatively stable except for occasionally large outliers. The predicted orbits of GPS, GLONASS, Galileo, and BDS of SDU-1h have the highest accuracy and the most stable performance. Combined with the statistical results in Table 8, the 1D mean RMS of the ultra-rapid predicted orbit with the proposed parallel approach is 3.2 cm, 5.1 cm, 5.6 cm and 11.8 cm, respectively. The accuracy of the predicted orbit of each satellite from GPS/GLONASS/Galileo/BDS are plotted in Figure 12, and BDS-2 MEO and BDS-2 IGSO of the BDS system are used in this processing. It can be seen from these figures that the predicted orbits of GPS, GLONASS, Galileo and BDS determined by the parallel approach with 1-h update frequency have the highest accuracy. The accuracies of GPS, GLONASS, Galileo and BDS predicted orbit are mostly better than 3 cm, 5 cm, 6 cm, and 12 cm, respectively. One of Figure 11. RMS of GPS, GLONASS, Galileo and BDS (from top to bottom) orbit of SDU-1h, SDU-6h and WHU-1h compared with CODE precise orbit.
From Figure 11, it can be seen that as the extrapolation time increases, the orbit accuracy of SDU-6h decreases in the orbit sequence. The phenomenon that the orbit accuracy fluctuates with extrapolation time does not exist in WUH-1h and SDU-1h. The accuracy of WUH-1h is better than that of SDU-6h result, and it is relatively stable except for occasionally large outliers. The predicted orbits of GPS, GLONASS, Galileo, and BDS of SDU-1h have the highest accuracy and the most stable performance. Combined with the statistical results in Table 8, the 1D mean RMS of the ultra-rapid predicted orbit with the proposed parallel approach is 3.2 cm, 5.1 cm, 5.6 cm and 11.8 cm, respectively. The accuracy of the predicted orbit of each satellite from GPS/GLONASS/Galileo/BDS are plotted in Figure 12, and BDS-2 MEO and BDS-2 IGSO of the BDS system are used in this processing. It can be seen from these figures that the predicted orbits of GPS, GLONASS, Galileo and BDS determined by the parallel approach with 1-h update frequency have the highest accuracy. The accuracies of GPS, GLONASS, Galileo and BDS predicted orbit are mostly better than 3 cm, 5 cm, 6 cm, and 12 cm, respectively. One of the GPS satellites, G18, has poor accuracy, which may be due to its large orbit error during the shadow period. In general, the accuracy of the ultra-rapid predicted orbit of each GNSS satellite with the parallel approach and 1-h updated frequency have greatly improved. the GPS satellites, G18, has poor accuracy, which may be due to its large orbit error during the shadow period. In general, the accuracy of the ultra-rapid predicted orbit of each GNSS satellite with the parallel approach and 1-h updated frequency have greatly improved.

Conclusions
A parallel ultra-rapid POD approach has been proposed in this study. The parallel approach properly makes use of parallel computing to improve the efficiency of multi-GNSS orbit determination. The construction of normal equations, parameter elimination and orbit extrapolation in the GNSS POD algorithm are optimized based on MPI/OpenMP. The following conclusions can be drawn: (1) The precision loss and feasibility of parallel POD are verified. The results show that 1D mean RMS of GPS, GLONASS, Galileo and BDS orbit in non-parallel mode are 1.3 cm, 2.7 cm, 2.7 cm, and 5.4 cm, respectively. Moreover, the precision loss of using the parallel approach is at the sub-millimeter level, which can be ignored, thus verifying the feasibility of the parallel method.
(2) The improvement of efficiency by the parallel approach is analyzed in detail, including multiple factors such as the number of processes, the number of threads, the number of stations and the number of satellite systems. The results show that with the increase of the number of threads and processes, the efficiency of the parallel approach increases more obviously. The approach is more applicable to large reference networks and multiple satellite systems.
(3) The superiority of the parallel approach used in multi-GNSS ultra-rapid POD is verified. The improvement of accuracy of predicted orbit with an increased update frequency of 1-h is analyzed in detail. The 1-h updated orbit determined by the parallel approach is compared with the 6-h updated orbit determined by the non-parallel approach, and the results show that the accuracy of GPS, GLONASS, Galileo and BDS are 33.4%, 31.4%, 40.1% and 32.8% higher than the 6-h results of the non-parallel approach, respectively. Additionally, the accuracy of the 1-h updated orbit is evaluated by comparing with IGS and CODE precise orbit products. The 1-h orbit product of Wuhan University is selected for comparison. The results show that the orbit with 1-h update frequency by the parallel approach has higher accuracy and stability. Further, statistical results show that the 1D mean RMS of the ultra-rapid predicted orbit of GPS, GLONASS, Galileo, and BDS with the proposed parallel algorithm is 3.2 cm, 5.1 cm, 5.6 cm and 11.8 cm, respectively.
In summary, the proposed parallel approach can greatly improve the efficiency of multi-GNSS ultra-rapid orbit determination and enhance the update frequency from 6-h to 1-h, which can help for the improvement of the ultra-rapid orbit accuracy and stability. The proposed parallel approach makes it possible to use a denser network of reference stations for multi-GNSS ultra-rapid POD with high update frequency for real-time users.
It should be noted that we adopted the ECOM1 SRP model for the one-step orbit determination of the GPS, GLONASS, Galileo and BDS in this study. For the specific elongated shape of Galileo satellites, orbit accuracy will be improved by using a prior model. In addition, clock determination is not related in this study, making it difficult to get the Signal In Space Range Error (SISRE). In the next studies, we will use an a priori model to reduce the orbit perturbations for Galileo, and we will combine the orbit and clocks through the processing in order to calculate clocks and orbital SISRE for the analysis of the algorithm.