A Decentralized Processing Schema for Efficient and Robust Real-time Multi-GNSS Satellite Clock Estimation

Real-time multi-GNSS precise point positioning (PPP) requires the support of high-rate satellite clock corrections. Due to the large number of ambiguity parameters, it is difficult to update clocks at high frequency in real-time for a large reference network. With the increasing number of satellites of multi-GNSS constellations and the number of stations, real-time high-rate clock estimation becomes a big challenge. In this contribution, we propose a decentralized clock estimation (DECE) strategy, in which both undifferenced (UD) and epoch-differenced (ED) mode are implemented but run separately in different computers, and their output clocks are combined in another process to generate a unique product. While redundant UD and/or ED processing lines can be run in offsite computers to improve the robustness, processing lines for different networks can also be included to improve the clock quality. The new strategy is realized based on the Position and Navigation Data Analyst (PANDA) software package and is experimentally validated with about 110 real-time stations for clock estimation by comparison of the estimated clocks and the PPP performance applying estimated clocks. The results of the real-time PPP experiment using 12 global stations show that with the greatly improved computational efficiency, 3.14 cm in horizontal and 5.51 cm in vertical can be achieved using the estimated DECE clock.

combination of analysis center (AC) products at the IGS real-time combination center [24]. As each of the processing lines can be run in different computers, we named it decentralized clock estimation (DECE) in this study.
The processing strategy is realized based on the Position and Navigation Data Analyst (PANDA) software for validation [25,26]. A network of about 110 real-time stations with the tracking capacity of GPS, GLONASS, Galileo, and BeiDou four systems is employed and real-time estimated clocks, as well as the computation time using the UD processing mode and the DECE mode, are compared and the PPP performance is also evaluated.
First, we introduce the current research status for multi-GNSS and real-time clock estimation, and then the multi-GNSS real-time clock estimation model is discussed. In addition, we present the strategy of data processing and the real-time clock estimation system based on PANDA software. The real-time clock product estimated by the DECE strategy is also validated by clock comparison and their application to multi-GNSS PPP. Finally, a summary of the results and corresponding conclusions are given.

Multi-GNSS (Global Navigation Satellite Systems) Clock Estimation
As UD and ED observations are the basic observations used in most of the current clock estimation as well as in this study, the corresponding observation equations for multi-GNSS clock estimation will be introduced and discussed in this section followed by the process for clock combination.

Undifferenced Observation Model
In the clock estimation using UD observations, both carrier phase and pseudorange observations are involved, and the ionospheric-free linear combination is formulated to eliminate the first-order ionospheric delay. Thus, the UD observation model for a given epoch is expressed as where s, r are subscript for the satellite and receiver, and the subscript G, R, E, C represent GPS, GLONASS, Galileo, and BeiDou satellite, respectively; PC and LC denote ionospheric-free pseudorange and phase observations, respectively; dt r and dT s are receiver and satellite clock corrections, respectively; b r and b s are the signal delay bias in receiver and satellite, respectively; d trop represents the zenith tropospheric delay while m is the corresponding mapping coefficient; N denotes the ambiguity of the ionosphere-free phase observations. v and l are the post-fit residuals pre-fit residuals of the ionospheric-free observations, respectively. In addition, the antenna phase center corrections, the phase windup, as well as the code bias between different types of codes, e.g., differential code biases (DCB) of P1 and C1, should be considered in Equation (1). Concerning the satellite DCB, which is stable enough as demonstrated by Dach et al. [27], it can normally be absorbed by clock corrections in the data processing. In other words, the final clock estimation is a combination of the "true" clock and satellite DCB. In addition, since the initial information in normal equation is mostly provided by pseudorange observation and GPS time is usually chosen as the time reference, an additional Inter-System Bias (ISB) parameter can be introduced for R, E, and C to compensate its offset from GPS. The Inter-Frequency phase Bias (IFB) of GLONASS Remote Sens. 2019, 11, 2595 4 of 20 signals are estimated in the orbit determination and corrected as a constant in clock estimation. Thus, Equation (1) can be expressed as where ISB R , ISB E , and ISB C represent the Inter-System Bias of GLONASS, Galileo, and BeiDou, relative to GPS, and d T s is the satellite clock absorbing the satellite-related code bias. Due to the fact that the receiver clock terms are linearly dependent with the satellite clock, any bias in the receiver clock estimate can be compensated by the satellite clock estimates. In addition, since an individual receiver clock parameter for each satellite system is estimated, the datum deficiency exists for each system. Thus, the following constraint is introduced for each satellite system [14]: As presented in Equation (3), the sum of the satellite clock error for each system is set as zero to separate the receiver clock and the corresponding satellite clock.
In the solution of the satellite clock based on Equations (2) and (3), all parameters listed in Table 1 are going to be estimated simultaneously. Suppose that n sys , n site , and n sat are the numbers of navigation systems, stations, and total multi-GNSS satellites, respectively, while the averaged satellite number that is tracked at each station each epoch for a single system is n p . Thus, the parameters to be estimated with Equation (2) are listed in Table 1. Table 1. Parameters need to be estimated at each epoch in the undifferenced (UD) model.

Item Parameters Number
Receiver clock and Inter-System Bias (ISB) n sys × n site Satellite clock corrections n sat Zenith troposphere delay n site Undifferenced ambiguity n sys × n p × n site Total n tot (n sys + 1 + n sys × n p ) × n site + n sat Assuming the numbers of satellites for GPS, GLONASS, Galileo, and BeiDou are 32, 27, 30, and 35, respectively. Without loss of generality, it is assumed that n p = 8, while for the GPS-only solution, n sys = 1, n sat = 32, and thus the total number of parameters at each epoch is n tot = 10 × n site + 32, in the case of a four system multi-GNSS solution, n sys = 4, n sat = 32 + 27 + 30 + 35 = 124, and n tot = 37 × n site + n sat .

Epoch-Differenced Model
As is well known, the majority of parameters in the UD model are ambiguities which are increasing linearly with the number of satellites and stations. It is clearly the main reason for the heavy computational burden in the UD model. The motivation of ED model is to eliminate the ambiguity parameters, and consequently reduce the computation load. The ED observation equations at each epoch can be expressed as follows, where ∆ is the difference operator between two adjacent epochs and the meanings of the other parameters are the same as in the UD model. Since the ambiguity parameters are removed, there are only three types of parameters left in the equation and they have been shown in Table 2. Since the code bias between GPS and other systems are stable over a day (Zhang et al., 2016), the ∆ISB can be neglected in Equation (4). The equations can be expressed as, All parameters left in the equation are shown in Table 2. It should be noted that the ED zenith troposphere delays cannot be ignored although they are very small, because they have a systematic trend along with time and will bias the estimates if ignored [15,17].

Item Parameters Number
Receiver clocks n site Satellite clocks n sat Zenith troposphere delay n site Total n tot 2 × n site + n sat Under the same assumption of the number of satellites for each system above, the total number of unknowns is 2 × n site + 32 and 2 × n site + 124 for GPS-only and multi-GNSS solutions, respectively. However, it is suggested by Equation (5) that only the variation of the clock parameter between epochs, i.e., ∆t, can be estimated in the ED model.
To get an intuitive impression on the comparison of the number of unknowns, Figure 1 illustrates the varying total number of parameters along with the number of stations n sta . Let us take n sta = 100 Remote Sens. 2019, 11, 2595 6 of 20 as an example, the total number of unknowns is 232, 324, 1032, 3824 for ED of GPS-only, ED of multi-GNSS, UD of GPS-only, and UD of multi-GNSS, respectively. It is concluded that the inclusion of the new satellite systems implies a significant increase in computation load for the UD model, whereas this dilemma can be mitigated by using the ED model. To get an intuitive impression on the comparison of the number of unknowns, Figure 1 illustrates the varying total number of parameters along with the number of stations sta n . Let us take 100 sta n = as an example, the total number of unknowns is 232, 324, 1032, 3824 for ED of GPS-only, ED of multi-GNSS, UD of GPS-only, and UD of multi-GNSS, respectively. It is concluded that the inclusion of the new satellite systems implies a significant increase in computation load for the UD model, whereas this dilemma can be mitigated by using the ED model.

Clock Combination
Although the ED model is computationally efficient enough to support a high-rate update, it can only provide epoch-differenced clocks. Therefore, the UD model must be involved in order to obtain the absolute clock. Thus, a combination process is essential to generate the final clock product that takes advantage of both ED and UD clocks. Theoretically, we can also consider the combination of several ED and/or UD clock products, such as UD clocks from different software packages or different networks, etc. However, for clarification in the following study, we concentrated on the combination of a single UD and ED products. A simple but efficient approach was implemented by using the latest available UD clocks to align the high-rate clock changes estimated by the ED model to generate the final high-rate absolute clocks.
Assuming that the update interval of the UD model is k times of the ED model, for example, if the update rate of the UD is 20 sec and that of ED 5 sec, then k = 4, and the combination for post-processing can be expressed as follows [16], where UD dT , ED dT represent the UD and ED clock at each epoch, dT represents the combined clock, t is the epoch time, and the symbol mod means remainder operator. From Equation (6), it can be seen that the essence of this approach is using the last UD estimate to calibrate the ED estimates at the following epochs.

Clock Combination
Although the ED model is computationally efficient enough to support a high-rate update, it can only provide epoch-differenced clocks. Therefore, the UD model must be involved in order to obtain the absolute clock. Thus, a combination process is essential to generate the final clock product that takes advantage of both ED and UD clocks. Theoretically, we can also consider the combination of several ED and/or UD clock products, such as UD clocks from different software packages or different networks, etc. However, for clarification in the following study, we concentrated on the combination of a single UD and ED products. A simple but efficient approach was implemented by using the latest available UD clocks to align the high-rate clock changes estimated by the ED model to generate the final high-rate absolute clocks.
Assuming that the update interval of the UD model is k times of the ED model, for example, if the update rate of the UD is 20 s and that of ED 5 s, then k = 4, and the combination for post-processing can be expressed as follows [16], where dT UD , dT ED represent the UD and ED clock at each epoch, dT represents the combined clock, t is the epoch time, and the symbol mod means remainder operator. From Equation (6), it can be seen that the essence of this approach is using the last UD estimate to calibrate the ED estimates at the following epochs. However, in the case of real-time clock estimation, the approach must be adapted, as the computation of the UD model usually takes a much longer time than the ED model. Thus, only the latest available UD clocks could be used to calibrate the ED clocks to absolute ones. For real-time processing, the update rate of the UD and ED model should be properly selected, i.e., the processing/update interval must be larger than the processing time needed and it must be k times the ED interval. The processing schema is illustrated in Figure 2. In the Figure, j is the epoch index of UD processing, the related epoch time is t j×k , and within two adjacent UD epochs, there are k ED epochs. However, in the case of real-time clock estimation, the approach must be adapted, as the computation of the UD model usually takes a much longer time than the ED model. Thus, only the latest available UD clocks could be used to calibrate the ED clocks to absolute ones.
For real-time processing, the update rate of the UD and ED model should be properly selected, i.e., the processing/update interval must be larger than the processing time needed and it must be k times the ED interval. The processing schema is illustrated in Figure 2. In the Figure  According to the above criteria for selecting k, we can assume that at each ED epoch within  (7) if With the following definition, Equation (7) can be rewritten as In the combination, accumulated ED clocks over the UD interval are computed, then the combined clocks are calculated using Equation (9), according to which the UD clock is available for the working ED epoch. It should be pointed out that the accumulated ED clocks over the latest two UD epoch intervals must be saved for continuous propagation of the clock reference.
Besides this, it should be noted that the matching of the UD and ED processes is simplified to an addition operation and the noise of them are not considered in this study. We take the results of the UD process as a constraint which generates the absolute clock offset, whereas the ED process generates a set of relative clocks.

Parameter Estimation Methodology
In this contribution, we use Square Root Information Filter (SRIF) to estimate real-time clock parameters both in UD and ED models. SRIF is an enhanced Kalman filter, which is suggested in GNSS data processing as it is numerically stable. The core principle of SRIF is to update the filter state using the Householder orthogonal transformation [28]. The a priori information of the state can be expressed as According to the above criteria for selecting k, we can assume that at each ED epoch within t j×k and t ( j+1)×k , let us say epoch t j×k+m , the UD clock at epoch t ( j−1)×k is always available, and that at epoch t j×k it may be available depending on the computation efficiency of the UD process. Thus, the combined clock reads as the first expression of Equation (7) if dT UD (t j·k ) is unavailable at epoch t j×k while reads as the second expression of Equation (7) With the following definition, Equation (7) can be rewritten as In the combination, accumulated ED clocks over the UD interval are computed, then the combined clocks are calculated using Equation (9), according to which the UD clock is available for the working ED epoch. It should be pointed out that the accumulated ED clocks over the latest two UD epoch intervals must be saved for continuous propagation of the clock reference.
Besides this, it should be noted that the matching of the UD and ED processes is simplified to an addition operation and the noise of them are not considered in this study. We take the results of the UD process as a constraint which generates the absolute clock offset, whereas the ED process generates a set of relative clocks.

Parameter Estimation Methodology
In this contribution, we use Square Root Information Filter (SRIF) to estimate real-time clock parameters both in UD and ED models. SRIF is an enhanced Kalman filter, which is suggested in GNSS data processing as it is numerically stable. The core principle of SRIF is to update the filter state using the Householder orthogonal transformation [28]. The a priori information of the state can be expressed as where A i is the coefficient matrix, z i is the observation vector, v i is the corresponding observation error vector, and e i is the post-fit residual. By using the Household transformation matrix T, the a priori information matrix R i is updated toR i . Then the updated state vector can be expressed aŝ Specific to the application of real-time clock estimation, we express state vectorx i as two parts which are the constant parameter vector p and the dynamic parameter vector y [29]. Considering the state transition equation, we can get the following state equation of p i+1 and p i , where w is the process noise and R w is the square root of the covariance of w that is R −1 w R −T w = D(w). Based on the Household transformation and similar to Equation (11), we can express the parameter update of real-time clock estimation as Then, we can predict the state vector as Based on Equation (15) and the new observations, the real-time clock parameters are estimated in a recursive way epoch by epoch using SRIF.

Decentralized Clock Estimation (DECE)
Instead of integrating ED and UD and combination processing in a single process in terms of threads, we propose a decentralized clock estimation mode in which UD, ED, and the combination can be run in separated computers with a network connection. Figure 3 shows a flow chart of the simplest DECE configuration with one UD, and one ED and a combination process each running in a server. Ultra-rapid orbit and real-time observations are needed for each server group and we fix site coordinates to Solution INdependent EXchange (SINEX) format solution [30,31]. Both clock products from the ED and UD processes are transferred to the server where the combination is running to generate the final clock products simultaneously according to Equation (9). It should be noted that the orbit switching on each server is synchronous.
One of the advantages of the DECE is that the data exchange among the processes is very small and mainly for transmitting the estimated clocks, and thus can be easily realized through a TCP/UDP port. A much more important advantage is that DECE can significantly improve the robustness of the real-time precise positioning service by running redundant processing processes.
Although only two clock estimation processes are deployed in Figure 3 and used in this study, it can be extended to include more with the same strategy and stations as redundant processing.  One of the advantages of the DECE is that the data exchange among the processes is very small and mainly for transmitting the estimated clocks, and thus can be easily realized through a TCP/UDP port. A much more important advantage is that DECE can significantly improve the robustness of the real-time precise positioning service by running redundant processing processes.
Although only two clock estimation processes are deployed in Figure 3 and used in this study, it can be extended to include more with the same strategy and stations as redundant processing.
Processing processes with different strategies or different station sub-networks can also be involved for using more stations. All the processing processes can be deployed in servers either of a local network or in an offsite network. The combination processes automatically collect the necessary clock estimates according to specified weights to generate a combined product. It is obvious that the robustness of the processing or the availability of the clock products can be dramatically improved by such a distributed system, as the system can always provide high-rate clock product if one of the redundant ED solutions and one of the redundant UD solutions are working well. The crash of the clock combination process is negligible, since it can provide well-qualified combined clocks instantaneously after restarting.

Real-time Experimental Evaluation
In order to validate the DECE strategy, we use the Position and Navigation Data Analyst (PANDA) software package [25,26] as the base, and the real-time estimator is adapted to a distributed environment and a tool for real-time clock combinations developed for the DECE processing. First, we describe the whole data processing of the real-time precise positioning service based on the PANDA software package, and then the network and finally the processing parameters and processing schema used in the validation. Processing processes with different strategies or different station sub-networks can also be involved for using more stations. All the processing processes can be deployed in servers either of a local network or in an offsite network. The combination processes automatically collect the necessary clock estimates according to specified weights to generate a combined product. It is obvious that the robustness of the processing or the availability of the clock products can be dramatically improved by such a distributed system, as the system can always provide high-rate clock product if one of the redundant ED solutions and one of the redundant UD solutions are working well. The crash of the clock combination process is negligible, since it can provide well-qualified combined clocks instantaneously after restarting.

Real-Time Experimental Evaluation
In order to validate the DECE strategy, we use the Position and Navigation Data Analyst (PANDA) software package [25,26] as the base, and the real-time estimator is adapted to a distributed environment and a tool for real-time clock combinations developed for the DECE processing. First, we describe the whole data processing of the real-time precise positioning service based on the PANDA software package, and then the network and finally the processing parameters and processing schema used in the validation.

Data Processing System
In our system, the real-time Precise Orbit Determination (POD) is carried out in batch-processing mode using the observations from MGEX and IGS networks. All parameters in PDO are estimated by least squares adjustment and previous 24-h observations are used to generate the rapid orbit. Then the real-time orbit is predicted based on the rapid orbits in a batch-processing mode by using orbit integrator for at least six hours in a similar way to that of the IGS ultra-rapid data processing. Then the real-time clock estimation processes will do parameter estimation epoch by epoch based on the fixed real-time orbit using the methodology introduced in Section 2.4.
The BKG Ntrip Client (BNC) software (version 2.12) [32] is used to receive real-time observations as well as ephemeris data to feed the filter for clock estimation, and is also used to generate state-space representation corrections from the estimated orbits and clocks. Real-time data are transported using the NTRIP protocol and encoded/decoded following the RTCM 3.3 standard [11]. The structure of all systems is shown in Figure 4.
batch-processing mode using the observations from MGEX and IGS networks. All parameters in PDO are estimated by least squares adjustment and previous 24-hour observations are used to generate the rapid orbit. Then the real-time orbit is predicted based on the rapid orbits in a batch-processing mode by using orbit integrator for at least six hours in a similar way to that of the IGS ultra-rapid data processing. Then the real-time clock estimation processes will do parameter estimation epoch by epoch based on the fixed real-time orbit using the methodology introduced in section 2.4.
The BKG Ntrip Client (BNC) software (version 2.12) [32] is used to receive real-time observations as well as ephemeris data to feed the filter for clock estimation, and is also used to generate state-space representation corrections from the estimated orbits and clocks. Real-time data are transported using the NTRIP protocol and encoded/decoded following the RTCM 3.3 standard [11]. The structure of all systems is shown in Figure 4. In order to evaluate the performance of the DECE processing strategy, the "Real-time Clock" function in Figure 4 can be replaced by a single UD processing or by the DECE in Figure 3 to carry out UD processing and DECE processing for comparison.  In order to evaluate the performance of the DECE processing strategy, the "Real-time Clock" function in Figure 4 can be replaced by a single UD processing or by the DECE in Figure 3 to carry out UD processing and DECE processing for comparison.

Processing Parameters
The data processing parameters are shown in Table 3 and are kept the same for both UD and ED mode if it is not specified. It should be mentioned that there is no official Phase Center Offset (PCO) and Phase Center Variation (PCV) information for Beidou and Galileo satellite and receiver antennas at present, so we use GPS antenna parameters instead.

Processing Parameters
The data processing parameters are shown in Table 3 and are kept the same for both UD and ED mode if it is not specified. It should be mentioned that there is no official Phase Center Offset (PCO) and Phase Center Variation (PCV) information for Beidou and Galileo satellite and receiver antennas at present, so we use GPS antenna parameters instead.

Results
The experimental real-time precise positioning service has been run operationally and we take the orbit and DECE clock products from Day of Year (DOY) 274 to 280, 2018 for the validation. In the experiment, we use a 5-s update rate for ED clock estimator and 20 s for UD clock estimator to generate 5 s final DECE clock products. The CPU we used is Intel(R) Xeon(R) CPU E5-2637 v4 @ 3.50GHz.

Decentralized Clock Estimation (DECE) Compared to Undifferenced (UD)
Taking the results of DOY 251 in 2018 as an example, Figure 6 shows the time consumed using UD and DECE at each epoch. The computing time of the DECE is 0.07 s on average, which is much less than that of the UD mode, which is about 5.66 s. It should be noted that it is difficult to update the UD clock within five seconds at each epoch, especially when quality control is involved. What is more, with the completion of Galileo and Beidou in the near future, 124 satellites will be used and the burden of calculation will further increase. On the contrary, the DECE model can solve this computational problem easily.

Troposphere
Piecewise constant zenith delay for each station every 2 hours with a constraint of 2 cm/√h Inter-System Bias (ISB) Estimated as constant/None for ED Inter-Frequency Phase Bias (IFB) Estimated in the orbit determination and corrected as a constant

Results
The experimental real-time precise positioning service has been run operationally and we take the orbit and DECE clock products from Day of Year (DOY) 274 to 280, 2018 for the validation. In the experiment, we use a 5-sec update rate for ED clock estimator and 20 sec for UD clock estimator to generate 5 s final DECE clock products. The CPU we used is Intel(R) Xeon(R) CPU E5-2637 v4 @ 3.50GHz.

Decentralized Clock Estimation (DECE) Compared to Undifferenced (UD)
Taking the results of DOY 251 in 2018 as an example, Figure 6 shows the time consumed using UD and DECE at each epoch. The computing time of the DECE is 0.07 seconds on average, which is much less than that of the UD mode, which is about 5.66 seconds. It should be noted that it is difficult to update the UD clock within five seconds at each epoch, especially when quality control is involved. What is more, with the completion of Galileo and Beidou in the near future, 124 satellites will be used and the burden of calculation will further increase. On the contrary, the DECE model can solve this computational problem easily. According to the combination algorithm described by (9), the difference between DECE and UD clocks is mainly caused by the non-synchronization of the clock datum. Taking the results of DOY 251 in 2018 as an example, Figure 7 shows the time series of differences between the DECE and UD products and each subplot shows the differences of all satellites of one system. According to the combination algorithm described by (9), the difference between DECE and UD clocks is mainly caused by the non-synchronization of the clock datum. Taking the results of DOY 251 in 2018 as an example, Figure 7 shows the time series of differences between the DECE and UD products and each subplot shows the differences of all satellites of one system.
From the results, it is obvious that most of the differences between DECE and UD are very subtle, normally smaller than 0.02 nsec (about 6 mm). For multi-GNSS real-time PPP, the influence of such differences can be ignored. Figure 8 shows the detailed behavior of the difference between the UD and the DECE results, which takes a 5-min time series as an example. Since the interval of UD results is 20 s, which cannot fit with the 5-s DECE clock at each epoch, the linear interpolation has implemented to the UD clock to ensure there is always a corresponding point fit to DECE. In addition, it will not be clear enough if the results of many satellites are plotted together in one figure, so here we plot the single satellite's series and take G02 and G06 as an example for analysis. From the results, it is obvious that most of the differences between DECE and UD are very subtle, normally smaller than 0.02 nsec (about 6 mm). For multi-GNSS real-time PPP, the influence of such differences can be ignored. Figure 8 shows the detailed behavior of the difference between the UD and the DECE results, which takes a 5-minute time series as an example. Since the interval of UD results is 20 seconds, which cannot fit with the 5-seconds DECE clock at each epoch, the linear interpolation has implemented to the UD clock to ensure there is always a corresponding point fit to DECE. In addition, it will not be clear enough if the results of many satellites are plotted together in one figure, so here we plot the single satellite's series and take G02 and G06 as an example for analysis.  From the results, it is obvious that most of the differences between DECE and UD are very subtle, normally smaller than 0.02 nsec (about 6 mm). For multi-GNSS real-time PPP, the influence of such differences can be ignored. Figure 8 shows the detailed behavior of the difference between the UD and the DECE results, which takes a 5-minute time series as an example. Since the interval of UD results is 20 seconds, which cannot fit with the 5-seconds DECE clock at each epoch, the linear interpolation has implemented to the UD clock to ensure there is always a corresponding point fit to DECE. In addition, it will not be clear enough if the results of many satellites are plotted together in one figure, so here we plot the single satellite's series and take G02 and G06 as an example for analysis.  Figure 8 is interesting because there appears to be a periodical slip every 4 epochs (20 seconds) in the series, but not always obviously. This is probably because normally the DECE process needs to adjust to the UD process as a datum every 4 epochs. According to the combination methodology described in section 2.3, the difference between DECE and UD clocks present in the series is mainly caused by the non-synchronization of the clock datum. However, more detailed property of these slips still needs further research.

Decentralized Clock Estimation (DECE) Compared to Post-Processing Products
As for further evaluation of the DECE products, they are also compared with the MGEX final products provided by GFZ (GBM) which is supposed to have the best quality among the existing multi-GNSS products [37]. The clock difference and the Signal-In-Space Ranging Error (SISRE) which is a measure of the joint effect of orbits and clocks will be analyzed below.  Figure 8 is interesting because there appears to be a periodical slip every 4 epochs (20 s) in the series, but not always obviously. This is probably because normally the DECE process needs to adjust to the UD process as a datum every 4 epochs. According to the combination methodology described in Section 2.3, the difference between DECE and UD clocks present in the series is mainly caused by the non-synchronization of the clock datum. However, more detailed property of these slips still needs further research.

Decentralized Clock Estimation (DECE) Compared to Post-Processing Products
As for further evaluation of the DECE products, they are also compared with the MGEX final products provided by GFZ (GBM) which is supposed to have the best quality among the existing multi-GNSS products [37]. The clock difference and the Signal-In-Space Ranging Error (SISRE) which is a measure of the joint effect of orbits and clocks will be analyzed below.

Clock Difference
Since receiver and satellite clocks are estimated epoch-wise as white noise, the clock datum could change from epoch to epoch and the inter-system range bias at receivers could also induce inter-system clock biases. Therefore, in the multi-GNSS clock comparison, a satellite clock from each system is selected as a reference, and single-differenced clocks with respect to the reference satellites are used, in which the datum differences between different products are removed. Here G01, R01, E01, and C11 are selected as the reference satellite for each system, respectively.
As the bias in the time series of the between-satellite clock differences can be absorbed by ambiguities in PPP, the Standard Deviation (STD) of the clock is a significant indicator to reflect the impact of the clock on phase-based positioning. STD values of the differenced clocks are calculated for each satellite and shown in Figure 9. Each subplot in Figure 9 is for a system indicated inside the ith system averaged STD value. For BDS to the left and right of the slash means the value for MEO/IGSO and GEO, respectively. From the statistical results above, it can be seen that the average STD of GPS is 0.08 nsec, while GLONASS, Galileo, BDS MEO/ISO, and BDS GEO are 0.24 nsec, 0.10 nsec, 0.20 nsec, and 0.52 nsec, respectively. The main reason for the larger difference of BDS GEO satellites is due to the insufficient tracking stations and the poor tracking geometry.

Signal-in-Space Ranging Error (SISRE)
In the clock estimation, most of the orbit biases in the radial direction can be absorbed by clock parameters. This means that satellite clocks are biased by orbit biases in the radial direction and the biases are complementary. The SISRE is a statistical measure for the impact of orbit and clock errors on the modeled pseudorange. The SISRE takes the orbit differences into consideration while comparing clocks of two products by projecting their satellite position difference on the line-of-sight direction from satellite to the user position. It is introduced as a more reliable indicator of the comprehensive influence of orbits and clocks. The formula of SISRE is given by Montenbruck et al. [38] as From the statistical results above, it can be seen that the average STD of GPS is 0.08 nsec, while GLONASS, Galileo, BDS MEO/ISO, and BDS GEO are 0.24 nsec, 0.10 nsec, 0.20 nsec, and 0.52 nsec, respectively. The main reason for the larger difference of BDS GEO satellites is due to the insufficient tracking stations and the poor tracking geometry.

Signal-in-Space Ranging Error (SISRE)
In the clock estimation, most of the orbit biases in the radial direction can be absorbed by clock parameters. This means that satellite clocks are biased by orbit biases in the radial direction and the biases are complementary. The SISRE is a statistical measure for the impact of orbit and clock errors on the modeled pseudorange. The SISRE takes the orbit differences into consideration while comparing clocks of two products by projecting their satellite position difference on the line-of-sight direction from satellite to the user position. It is introduced as a more reliable indicator of the comprehensive influence of orbits and clocks. The formula of SISRE is given by Montenbruck et al. [38] as where dR, dA, and dC represent the orbit differences in radial direction, along direction, and cross directions, while dt denotes the real-time clock error compared to final products, C is the speed of light in vacuum, a and b represent the SISRE coefficient according to satellite altitude. The coefficients for multi-GNSS SISRE are shown in Table 4 [38]. Obviously, we can find that the SISRE is mostly affected by the radial direction orbit difference and clock difference. The results of SISRE are shown in Figure 10. From the statistical results above, it can be seen that the Root Mean Square (RMS) of SISRE for GPS is 0.12 m, while GLONASS, Galileo, BDS MEO/ISO, and GEO 1.51 m, 0.11 m, 0.33 m, and 0.52 m, respectively. The SISRE of GLONASS is significantly larger than for other systems probably due to frequency division multiple access (FMDA) which needs further research. It should be noted that the SISRE mainly reflects the influence of the products on pseudorange positioning.

Precise Point Positioning (PPP) Validation
PPP is a convincing way to directly verify the orbit and clock products together. In the experiment, real-time kinematic PPP is also carried out for 12 global multi-GNSS stations from MGEX using the DECE products and the UD products, respectively. The data processing strategy and parameter model are the same as in Table 3, except that the satellite clock is fixed in PPP and coordinates are estimated in the kinematic mode as white noise. The real-time streams come from mgex.igs-ip.net:2101 and the sample rate is set to five seconds. PPP for all stations is processed using multi-GNSS observations.

Precise Point Positioning (PPP) Validation
PPP is a convincing way to directly verify the orbit and clock products together. In the experiment, real-time kinematic PPP is also carried out for 12 global multi-GNSS stations from MGEX using the DECE products and the UD products, respectively. The data processing strategy and parameter model are the same as in Table 3, except that the satellite clock is fixed in PPP and coordinates are estimated in the kinematic mode as white noise. The real-time streams come from mgex.igs-ip.net:2101 and the sample rate is set to five seconds. PPP for all stations is processed using multi-GNSS observations. Figure 11 illustrates the positioning error of the station ASCG in the north, east, and vertical components. The figure shows that the positioning accuracy of PPP with the DECE clock in a real-time mode is better than 10 cm in the horizontal and 20 cm in the vertical after initialization, which can meet the need most positioning applications. Similar accuracies were obtained for the other stations. The RMS is an important indicator for evaluating the positioning performance which can be calculated as In the experiment, the RMS of the position differences using DECE and UD products are shown in Figure 12, in which the legend H and V mean horizontal and vertical direction, respectively. All results are counted after PPP convergence. The statistics show that the average RMS is about 3.10 cm for using UD and 3.14 cm for DECE in horizontal, and 5.47 cm for using UD and 5.51 cm for DECE in vertical. In general, for PPP, centimeter-level positioning results can be achieved using the clock products. Normally the result of UD is a little better than DECE, but there is almost no difference in the accuracy between using DECE and UD products. The RMS is an important indicator for evaluating the positioning performance which can be calculated as where X ppp and X re f represent any of the calculated and reference coordinates and n is the epoch number used for PPP. In the experiment, the RMS of the position differences using DECE and UD products are shown in Figure 12, in which the legend H and V mean horizontal and vertical direction, respectively. All results are counted after PPP convergence. The statistics show that the average RMS is about 3.10 cm for using UD and 3.14 cm for DECE in horizontal, and 5.47 cm for using UD and 5.51 cm for DECE in vertical. In general, for PPP, centimeter-level positioning results can be achieved using the clock products. Normally the result of UD is a little better than DECE, but there is almost no difference in the accuracy between using DECE and UD products.

Discussion
With the great progress of multi-GNSS real-time PPP, the performance of real-time clock products becomes more and more important. Different from post-processing and near-RT modes, the computational efficiency of real-time clock estimation is critical, because the delay or missing of the products will directly affects the positioning accuracy of the PPP client. Moreover, since the recovery of real-time service normally takes a lot of time, especially in case of software crashes, higher requirements for the stability and continuously of real-time products are also necessary. In section 3, to solve the problems above, we proposed a DECE strategy to improve the computational efficiency and enhance the robustness of the real-time system. First, both low-rate UD and high-rate ED processes are implemented to guarantee the timely update of real-time clock products. From Figure 1 and Figure 7, it is noted that although with the launch of more new satellites and setup of more ground GNSS stations in near future, the processing pressure will still not increase significantly. In addition, when using the DECE strategy, more than one line can be the processed in different processing centers or processing using different strategies (ED/UD). This means the robustness of the processing or the availability of the clock products can be significantly improved by such a distributed system, as the system can always provide high-rate clock products if one of the redundant ED solutions and one of the redundant UD solutions are both working well. As the PPP experimental results show in section 5.3, centimeter accuracy position can be achieved using the DECE real-time products, which is a great improvement of computational efficiency and robustness. It also improves the experience of real-time positioning users.

Conclusions
We proposed the decentralized clock estimation approach to improve the computational efficiency and robustness of the real-time precise positioning service. In the new approach, both the UD and ED modes are implemented but run separately in different computers. The UD mode estimates clock offsets with a lower update rate because a great number of ambiguities parameters are included, while the ED mode determines clock variations with a higher update rate. The products of the two modes are combined to generate the final products. Redundant UD and/or ED processing lines can be scheduled even in offsite computers with an internet connection to improve

Discussion
With the great progress of multi-GNSS real-time PPP, the performance of real-time clock products becomes more and more important. Different from post-processing and near-RT modes, the computational efficiency of real-time clock estimation is critical, because the delay or missing of the products will directly affects the positioning accuracy of the PPP client. Moreover, since the recovery of real-time service normally takes a lot of time, especially in case of software crashes, higher requirements for the stability and continuously of real-time products are also necessary. In Section 3, to solve the problems above, we proposed a DECE strategy to improve the computational efficiency and enhance the robustness of the real-time system. First, both low-rate UD and high-rate ED processes are implemented to guarantee the timely update of real-time clock products. From Figures 1 and 7, it is noted that although with the launch of more new satellites and setup of more ground GNSS stations in near future, the processing pressure will still not increase significantly. In addition, when using the DECE strategy, more than one line can be the processed in different processing centers or processing using different strategies (ED/UD). This means the robustness of the processing or the availability of the clock products can be significantly improved by such a distributed system, as the system can always provide high-rate clock products if one of the redundant ED solutions and one of the redundant UD solutions are both working well. As the PPP experimental results show in Section 5.3, centimeter accuracy position can be achieved using the DECE real-time products, which is a great improvement of computational efficiency and robustness. It also improves the experience of real-time positioning users.

Conclusions
We proposed the decentralized clock estimation approach to improve the computational efficiency and robustness of the real-time precise positioning service. In the new approach, both the UD and ED modes are implemented but run separately in different computers. The UD mode estimates clock offsets with a lower update rate because a great number of ambiguities parameters are included, while the ED mode determines clock variations with a higher update rate. The products of the two modes are combined to generate the final products. Redundant UD and/or ED processing lines can be scheduled even in offsite computers with an internet connection to improve the robustness of the processing system. More processing lines for different networks can also be included to improve the clock quality.
The new approach is realized for the experimental evaluation based on the PANDA software package for GNSS data processing and BNC software for data communication. The experiment was carried out in real time with about 110 stations for clock estimation and 12 stations for PPP. The clock comparison of the new approach with GFZ MGEX product shows that the STD of GPS is 0.08 nsec, while GLONASS, Galileo, BDS MEO/ISO, and BDS GEO are 0.24 nsec, 0.10 nsec, 0.20 nsec, and 1.52 nsec, respectively. The STD of signal-in-space range error (SISRE) of the clock product of GPS is 0.02 m, while GLONASS, Galileo, BDS MEO/ISO, and BDS GEO are 1.51 m, 0.11 m, 0.33 m, and 0.52 m, respectively. Using the estimated clocks and corresponding orbit products, real-time kinematic multi-GNSS precise point positioning can be realized with an averaged RMS of about 3.14 cm in horizontal and 5.51 cm in vertical.
All these results confirm that the decentralized clock estimation can provide comparable clock products as the most used undifferenced method but with a much higher computational efficiency and robustness. It is a suitable approach for multi-GNSS real-time clock estimation with the increasing number of satellites and stations and signals of different frequency bands.
Author Contributions: All authors contributed to the work. X.J. conceived and designed the experiments; X.J. and S.G. performed the experiments; X.J. and P.L. analyzed the data; X.J. and P.L. wrote the paper; S.G., P.L., M.G. and H.S. helped modify the paper.