A square root information filter for multi-GNSS real-time precise clock estimation

Real-time satellite orbit and clock estimations are the prerequisite for Global Navigation Satellite System (GNSS) real-time precise positioning services. To meet the high-rate update requirement of satellite clock corrections, the computational efficiency is a key factor and a challenge due to the rapid development of multi-GNSS constellations. The Square Root Information Filter (SRIF) is widely used in real-time GNSS data processing thanks to its high numerical stability and computational efficiency. In real-time clock estimation, the outlier detection and elimination are critical to guarantee the precision and stability of the product but could be time-consuming. In this study, we developed a new quality control procedure including the three standard steps: i.e., detection, identification, and adaption, for real-time data processing of huge GNSS networks. Effort is made to improve the computational efficiency by optimizing the algorithm to provide only the essential information required in the processing, so that it can be applied in real-time and high-rate estimation of satellite clocks. The processing procedure is implemented in the PANDA (Positioning and Navigation Data Analyst) software package and evaluated in the operational generation of real-time GNSS orbit and clock products. We demonstrated that the new algorithm can efficiently eliminate outliers, and a clock precision of 0.06 ns, 0.24 ns, 0.06 ns, and 0.11 ns can be achieved for the GPS, GLONASS, Galileo, and BDS-2 IGSO/MEO satellites, respectively. The computation time per epoch is about 2 to 3 s depending on the number of existing outliers. Overall, the algorithm can satisfy the IGS real-time clock estimation in terms of both the computational efficiency and product quality.


Introduction
Global Navigation Satellite System (GNSS) Global Real-Time Precise Positioning (RTPP) based on Precise Point Positioning (PPP) has been a hot topic for decades and it is recognized as the most promising system for future real-time precise positioning services. Currently several research institutions and commercial companies are providing their operational RTPP service (Caissy et al., 2011;Dixon, 2006;Leandro et al., 2011). It extends not only the geographical coverage of GNSS real-time precise positioning but also its application to a number of new areas, such as earthquake monitoring (Allen & Ziv, 2011), tsunami warning systems (Saito et al., 2016), and real-time atmospheric sounding (Li et al., 2014).
As a prerequisite of the RTPP service the provision of precise satellite clock corrections is challenging since satellite clocks must be estimated precisely using a globally distributed network and disseminated to users at an update interval of 5 s or higher. On the one hand, a dense network can provide better observation geometry and ensure a higher clock precision and stability. On the other hand, it is time-consuming to process such huge number of observations, especially for multi-GNSS with around Zuo et al. Satellite Navigation (2021) 2:28 120 operational satellites to be processed. Therefore, it is important to balance the station number and the computation efficiency, and to improve the computation efficiency is always a persistent pursue in real-time clock determination. There are two major aspects to speed up the clock estimation. One is reducing the number of the parameters to be estimated, especially the ambiguities, for example, using epoch-differenced observations (Ge et al., 2012;Jiang et al., 2019;Zhang et al., 2007Zhang et al., , 2011, or using the fixed integer ambiguities to convert phases to carrier-ranges (Chen et al., 2014a(Chen et al., , 2014bDai, Lou, et al., 2019;Laurichesse et al., 2009). The other one is utilizing the most advanced computation algorithms and techniques, such as fast matrix operation and parallel processing (Chen et al., 2014a(Chen et al., , 2014bGong et al., 2018;Fu et al., 2019).
Three estimators are commonly used in real-time clock estimation, i.e., the Kalman Filter (KF), the Least Squares adjustment (LSQ), and the Square Root Information Filter (SRIF, Bierman, 1975). Among them KF is the most widely used. For example, Hauschild & Montenbruck (2009) developed a real-time system using KF to generate GPS satellite clocks and applied them to the precise orbit determination of low-earth-orbiting satellites. The Centre National d'Études Spatiales (CNES) also adopted KF for multi-GNSS real-time satellite clock determination (Laurichesse et al., 2013). Usually, LSQ is more suitable for post-processing where the normal equation is accumulated and solved only once after the last epoch, as the inversion of the epoch-wise high-order normal equation to obtain the epoch-wise solution in real-time is time consuming. Nevertheless, using LSQ for real-time GNSS clock estimation was also investigated and the computation efficiency can also satisfy the real-time requirements (Fu et al., 2019). SRIF was developed for system control and spacecraft navigation because of its high numerical stability and computational efficiency (Bierman, 1975). It was implemented in the GIPSY software package for GPS data processing and later in the RTG (Real Time GIPSY) for real-time applications (Bertiger et al., 2020). It was also realized in the EPOS-RT (Ge et al. 2006) and PANDA software package (Liu & Ge, 2003;Gong et al., 2018) for real-time precise positioning services. Gong et al. (2018) further improved processing efficiency of SRIF by taking the advantage of the dense linear algebra algorithms and ambiguity elimination method.
Although these estimators were developed and realized for different purposes, they are based on the same principles, and thus suffer from the same problem of computational efficiency. It must be pointed out that the calculation of the state vector and its full covariance matrix is a must in both LSQ and KF for real-time clock estimation because of the KF state propagation and the quality control in LSQ. However, SRIF is more flexible in obtaining the estimates without loss of the opportunity of quality control. Besides the full solution as KF and LSQ, the other options could be no solution but only the updated square root information, or the estimated state vector only, or state vector with the standard deviations. These characteristics can be utilized to save computational time by providing a minimal solution, for example, we do not need the covariance matrix of the estimates for real-time clock estimation. More important is that choosing any of the above-mentioned solutions, the posterior observation residuals are always available, as they were generated in the measurement update before solving the equation, so that the quality control based on the posterior residuals can be carried out for all the solutions. Moreover, the possible numerical instability caused by indefinite matrix of KF and sequential LSQ can be overcome since SRIF is working on the form of triangularized Square Root Information (SRI) matrix instead of the covariance matrix or normal equation, which means it can always guarantee that the corresponding covariance matrix is positive and symmetric.
Online quality control is a critical and indispensable component of the real-time GNSS data processing in both the server and user end, as undetected outliers can degrade the product accuracy and reliability, even cause the estimator to crash and significantly deteriorate the stability and robustness of the online service, which will consequently affect the user's positioning. The quality control of GNSS real-time processing using LSQ and KF is quite mature and can be found in a number of references (Yang & Gao, 2006, Yang et al., 2010Fu et al., 2019), however it is not yet well demonstrated for SRIF. In the study by Lichten (1990), the predicted residual and its STD are employed to make the decision whether the observation contains an outlier or not according to the rule of three times of STD and the confirmed outliers are excluded in the processing. Shi et al. (2008) developed a new method based on a correlation analysis to scan outliers in the real-time positioning and its efficiency was demonstrated with both simulated and actual outliers in a kinematic experiment. Dai, Lou, et al. (2019), ) also applied SRIF for real-time precise orbit determination where an approach was developed for detecting irregular satellite movement, such as satellite maneuvers. However, the quality control within SRIF should be further optimized for satellite clock estimation in a large network with an update rate of every 5 s or even higher. Based on existing studies, we developed a complete set of quality control procedure for GNSS real-time clock estimation, including the Detection, Identification, and Adaption (DIA) steps, which is capable of detecting multiple outliers simultaneously. The computational efficiency was also improved by optimizing the computation of sensitivity vectors used in the quality control and properly applying the BLAS/LAPACK matrix operation algorithms. The processing procedure is implemented in the PANDA software package and evaluated in the operational real-time processing at the GFZ (German Geoscience Research Center) IGS (International GNSS Service) real-time analysis center.
This study is organized as follows. First a short description of the mathematical model of multi-GNSS clock estimation is given. Then the principle and the basic algorithm of SRIF including measurement and time updating and parameter eliminating are introduced. The procedure of the quality control for SRIF is presented and discussed in detail. Then the efficiency of this method is investigated with the numerical result of the operational running at the GFZ IGS real-time analysis center. Finally, the study is summarized, and conclusions are drawn.

Mathematical model for real-time clock estimation
As station coordinates and satellite orbits can be precisely derived in advance and fixed as known in the clock estimation, the linearized observation equations of ionospherefree phase and range can be expressed as: where s and r refer to the satellite and receiver, respectively; v lc and v pc denote "Observed Minus Computed" (OMC) residuals of phase and range observations; δt r and δt s are receiver and satellites clock parameters; b s r and d s r represent the receiver-dependent un-calibrated phase delay and code bias of satellite s; b s and d s are the satellite-dependent un-calibrated phase delay and code bias; B s r is the ambiguity of the ionosphere-free phase observation; T r and m s is the zenith wet tropospheric delay and its mapping function; ε lc and ε pc denote the sum of measurement noise and unmodeled error for phase and range observations. Note that the systematic error corrections such as phase windup effect, antenna phase center offset and variation, and the station displacements are already applied in the equation.
In order to obtain the satellite clock solution from Eq. (1), the bias parameters need to be reformed as some of them are strongly correlated. For example, the stable satellite code bias cannot be separated from the corresponding satellite clock, and also the un-calibrated phase delay will be absorbed by the corresponding ambiguity. Thus, Eq. (1) can be reformulated as: (2) v lc = δt r − δt s + B s r + m s · T r + ε lc v pc = δt r − δt s + m s · T r + ε pc where the estimated receiver and satellite clock are δt r = δt r + d s r and δt s = δt s + d s , respectively, and the ambiguity is expressed as B s r = B s r − b s + d s + b s r − d s r . For a multi-GNSS receiver the measurements could be biased systematically from satellite to satellite due to the signal delay caused by hardware. For the GNSS using Code Division Multiple Access (CDMA) signals, including GPS, BDS, and GALILEO, the differences between different satellites of the same system are very close to each other as they have the same frequency, and thus only the biases between different systems should be handled, that is, the Inter-System Biases (ISB). As for the GLONASS satellites using Frequency Division Multiple Access (FDMA) signals, different satellites have different frequencies, and thus the bias should be parameterized between different satellites, that is, the Inter-Frequency Bias (IFB). Using GPS time as time reference and parameterizing the ISB and IFB with respect to the GPS observations, the combined clock estimation models of GPS, GLONASS, BDS, and GALILEO can be expressed as: where G, C and E refer to the GPS, BDS, and GALILEO systems, respectively; R k denotes the GLONASS satellite with frequency factor k; The clock offset of a satellite and a receiver can be modeled as white noise or random walk process. The residual zenith tropospheric delay is usually parameterized as a random walk process on each station after a priori model correction, which removes the tropospheric dry component. The ISB and IFB parameters could be modeled as constant because of its stability or as random walk to consider possible unexpected variation. In order to solve the singularity caused by the full correlation between the satellite clocks from a single system and the corresponding ISB parameters of all stations, the constraint that the sum of all ISB parameters for each individual system equals zero, is imposed. For the same reason, the same constraint should be applied to the IFB parameters of GLONASS satellites in a network. Furthermore, in order to save computation time, IFB and ISB can be corrected using the values estimated in the previous days according to the stability of the estimates. Assuming that we have m observations and n parameters at epoch i , the observation equations can be expressed in the following simplified form as: where l i is the observation vector; H i is the design matrix; x i is the parameter vector; ε i is the measurement noise with a variance matrix D ε i .
In GNSS data processing, the parameters can be handled as either deterministic parameters y or the time varying stochastic parameters p , which are usually parameterized as white noise (for instance, satellite and receiver clocks) or random-walk process (for instance, tropospheric delays). The random-walk process can be described with the following state equation of the firstorder Gauss-Markov process: In SRIF, the noises in the observation equations and state equations should follow the standard normal distribution, that is N (0, I) . This can be achieved by adapting the corresponding equations. Taking the variance matrix of the observation as example, D ε i can be decomposed as D ε i = R −1 z i R −T z i where R z i is an upper triangular matrix and referred as square root information (SRI) of the observations. Then Eq. (4) can be normalized by multiplying R z i on the both sides as: The same normalization can be applied to the state equations and with D w i = R −1 w i R −T w i the normalized state equation of Eq. (5) reads: , the corresponding initial constraint can be expressed in form of pseudo-observations as: , the normalized observation equations should be: In summary, the Eqs. (4), (5) and (8) which are directly used in LSQ and KF, but must be normalized to Eqs. (6), (7) and (9) for SRIF.

Square root information filter
The basic operations of SRIF are the measurement update and time update using the Householder orthogonal transformation (Bierman, 1975). The principle of the SRIF and the Householder transformation, and the time and measurement update are summarized in this section, as they are essential for understanding the new quality control algorithm.

The principle of SRIF
In order to explain the principle of SRIF, we assume that the estimation problem includes only the initial information of parameter x in Eq. (9) and the observation equations Eq. (6), the state equations Eq. (7), that is, the time update will be presented later. Therefore, the observation equations of the estimation can be summarized as: By applying the Householder orthogonal transformation T which will be introduced in the next section to Eq. (10) instead of forming normal equation in the LSQ, we have: The least squares solution of Eq. (11) is to find the parameter vector x that satisfies: where R is also an upper triangular matrix and e is the a posteriori residual vector corresponding to the observations z . It should be noted that the a posteriori residual vector e also follows the standard normal distribution like the noise vector of the observations and initial information as T is an orthogonal transformation.
In order to minimize the norm of the residual vector, it is obvious that x must satisfy: Then, the estimates and their covariance matrix can be derived accordingly. Equation (13) can also be used as initial information to be combined with new observations or with state equations for sequential estimation. It is obvious that the orthogonal transformation T is the key of SRIF and will be used for both time and measurement update.

Householder transformation
The Householder orthogonal transformation is a computationally efficient and numerically stable approach among all the QR factorization methods (Demmel, 1997). It is carried out column by column from left to right of the matrix to be factorized and ended with an upperright triangular matrix. The main purpose to repeat this algorithm here is to explain how to calculate the sensitivity vectors needed for the detection of the possible outliers in the quality control.
Given is a unit vector u , its orthogonal space is U ⊥ , for any vector g , its reflection g ′ with respect to U ⊥ can be obtained by the Householder transformation T u (Bierman, 2006): with If g and g ′ are given, u can be calculated for getting the transformation matrix: The triangularization of an information matrix is carried out column by column from left to right sequentially. For the j-th column the transformation is applied to the low-right part of the matrix and the transformation matrix T j can be expressed by the corresponding vector u j , whereas the rest remains unchanged. The complete transformation matrix for the whole matrix is: In practice, when we apply T to vector a , instead of reconstructing T which is too time consuming we reform the matrix multiplication to vector addition as: Start from j = 1 with a 0 = a to j = n, the resulted a n is the transformed vector of a.

Time update
In real-time satellite clock estimation, the state vector x includes the deterministic parameters y and the stochastic parameters p , of which the former are constant at all epochs such as ambiguities and the latter vary over epochs such as satellite and receiver clocks. After measurement update, Eq. (13) can be expressed as Combining Eq. (19) and the state equations Eq. (7) the information matrix for time update reads: Applying the Householder transformation, we obtain the following equation: Obviously, the outdated parameters p i−1 can be removed from Eq. (21) and the state vector can be predicted by solving the following equation: After the time update, the state parameters of the current epoch are introduced and Eq. (22) is ready as initial information for the measurement update.

Measurement update
Combining the observation equations at the epoch i of Eq. (6) and the initial information of Eq. (22) The standard full solution with the estimates and their covariance matrix can be derived straightforward as: However, in practice instead of the above full solution of Eq. (26), there are also other options for saving computation time according to the requirement of each individual application. For example, if no estimated parameters are desired in real-time, the solution step can be simply skipped. In the clock estimation, we need the clock estimates every epoch, but their STD which is used as the quality indicator can be updated at a much lower frequency, or can be simply approximated by the number of observation involved in the estimation.
Another reason for solving for parameters and their covariance matrix is to obtain the posterior residuals and their STDs for the online quality control which is essential for real-time data processing. It is worth to point out that the SRIF algorithm makes the posteriori residuals available before solving for the parameters, so that the quality control based on the posteriori residuals can be carried out in any case. This provides the opportunity to speed up the computation and makes SRIF more efficient than KF and LSQ.
For the next epoch i + 1, Eq. (26) is considered as a prior information the same as Eq. (19) for time-update, so that sequential processing can be carried out for real-time processing.

Parameter elimination
The deactivated parameters, such as the ambiguity parameters for previous data arcs and the parameters to mitigate the effect of outliers should be eliminated from SRIF to save computer memory and to reduce computation burden. Assuming that a set of deterministic parameters y 1 will be eliminated and the remaining parameters are included in y 2 , that is, y are divided into the two groups, then by ignoring the epoch index Eq. (25) can be rewritten as with the index b for the SRI before the elimination. In order to remove the parameter y 1 , the parameters vector in Eq. (27) should be re-sorted as (26) By applying Householder transformation on Eq. (28), we get Finally, the parameter vector y 1 can be removed as it is independent of the existing parameters and will never be used afterwards, and we can obtain the following equation for the time update for the next epoch or measurement update:

Real-time quality control of SRIF
In this section, we develop a practical real-time quality control procedure for SRIF based on the study by Shi et al. (2008), dedicating to the processing of huge GNSS networks, and make effort to improve its computational efficiency. The new approach follows the well-known procedure with three sequential steps: detection, identification, and adaptation (Teunissen, 1990(Teunissen, , 1998(Teunissen, , 2018Yang et al., 2010) and is presented in detail hereafter.

Detection
The first step of the SRIF QC is to detect the existence of undetected cycle slips or outliers in the observation vector. Starting with the measurement update model of Eq. (23) including the information after the time update and the observations at the current epoch, Eq. (24) is obtained by applying the orthogonal transformation T . In Eq. (24) e i is the corresponding vector of the posterior residuals at the current epoch which theoretically satisfies e i ∼ N (0, I) , and the variance of unit weight can be calculated as: where e T i e i follows the Chi-square distribution with the degrees of freedom m.
In principle, we can test the individual residuals or the variance of unit weight of Eq. (31) to detect whether there is any outlier or not. In this study, we check both by constructing the following simple hypothetical test: where k 1 and k 2 are the thresholds for the residuals and unit weight STD, respectively, and theoretically they can be set according to their distributions. In this study for the multi-GNSS clock estimation, we use the empirical value of 5.0 and 1.5, respectively, which are based on our long-term operational processing experience at the GFZ IGS real-time analysis center. It is worth to mention that the unit weight STD depends on the a priori STD of the observations, which should be fine-tuned for each station, so that the unit weight STD is on average close to 1. This can affect the selection of the aforesaid thresholds. For example, a too optimistic the a priori STD for observations will results in larger residuals and consequently a large STD of unit weight, resultinng in more false outliers detected. In contrast, with a too large the a priori STD for observations, we may fail to detect the real outliers.
If H 0 is accepted, there is no problematic observation at this epoch; otherwise, usually at least one outlier exists in the observation vector. It should be noted that the rejection of H 0 can also be caused by the inaccurate modeling of the state parameters. However, in this contribution, we focus only on the quality control for observation blunders.

Identification
After H 0 is rejected, the second step of the SRIF quality control is to find out which observations are contaminated by undetected cycle slips or blunders. The basic idea is to extend the observation model to include the possible outliers and then check whether under the extended model the H 0 can be accepted. We first present the approach and then discuss how to select the outlier candidates from thousands of observations.
Assuming that there are n b possible outliers with the observation index ip(k), k = 1,2,…n b , the function model Eq. (23) can be extended by introducing the corresponding outlier parameters with the a priori values of zero and variance matrix D = R −1 R −T as: where θ � is a matrix of m × n b with unit vector for all the columns, and for its k-th column only the ip(k)-th element equals to 1. By applying the same orthogonal transformation as that used to Eq. (23) in the measurement update at epoch i, we can get where S ,i can be considered as the sensitivity matrix of the residual vector e i with respect to outliers . It means that the magnitude of the un-detected cycle slips or blunders is mapped into the a posteriori residual vector through S ,i as: It is clear that the residual vector e i is a combination of outliers and observations noises, so the outlier parameters can be solved from Eq. (36) using the least square adjustment and the solution reads as Then the same hypothetical test H 0 can be conducted with the residuals and STD of the above solution. If the test is passed, the set of outliers should be accepted as the finally identified ones; otherwise, an additional problematic observation should be selected and added, then the same procedure from Eq. (33) to (37) is repeated until H 0 is accepted. For the positioning of a single station, only tens of observations are involved per epoch, and hence it is possible to test all possible outlier combinations to identify the right ones. However, it is not possible for the real-time clock estimation using a global network, as the observation number per epoch reaches up to several thousands. Therefore, it is very important to select the outlier candidates which are critical to balance the product reliability and the computation efficiency.
Considering the situation of a global network with about 100 stations, there are enough redundant observations to detect the outliers which should show the largest residuals in the least squares adjustment. Therefore, we developed an empirical approach to select the potential outlier candidates according to the magnitude of the residuals. The residuals are sorted according to their absolute values and the largest one is selected as the most likely outlier in the associated adjustment. Then the above-mentioned identification procedure, that is, Eq. (33) to Eq. (37), is carried out. The whole identification will stop if the hypothesis H 0 is accepted. Otherwise, the updated residuals of Eq. (37) will be used to find out the next most likely outlier according to their absolute values and added to the already selected ones for the next iteration. The procedure is repeated until the H 0 is accepted or a certain number of observations has been marked as outliers, for example 100, in our operational data processing. It must be pointed out that a sophisticated pre-processing should be implemented to find out as many problematic observations as possible, particularly the large outliers, and only few small outliers can remain undetected to reduce the computation cost of quality control and to avoid numerical instability caused by large outliers. We have designed for each satellite-station pair a channel-filter to identify outliers by checking the length of data gaps, the variation of raw phase and range observations, and especially the Melbourne-Wübbena combination and ionosphere-delay combination for jumps. From our operational results, after the channel-filtering the remaining outliers are quite few which will be presented in the experimental validation.
Another critical issue is the computation of the sensitivity vector S because there are thousands of observations, and each could be an outlier. Theoretically, we can compute and save the sensitivity for all the observations, however it is not only too time-consuming, but also needs a very huge memory. According to Eq. (18) where the matrix multiplication is replaced by vector addition, the sensitivity vector S for a single observation can be calculated by applying the orthogonal transformation to the corresponding vector. Therefore, only the sensitivity vectors for the observations selected as outlier candidates have to be calculated, which is usually very small in quantity, for example in the operational processing the epochs with one outlier account to only about 15.58%, and those with two account to about 5.28%..

Adaptation
Once the outliers are identified, the negative impact of the identified outliers must be removed from the filter.
As it is very difficult to down-weight or weight negatively the corresponding observations, we adapt the filter by extending the function model as is done in the identification step. This is also because the most important effect comes from the cycle slips which can only be mitigated by adding an outlier parameter. The extended function model for the final identified outliers is already achieved with the help of the sensitivity vectors.
Assuming that there are n b outliers found at the last step of the identification, the adapted model with the n b outliers is already available, i.e., Eq. (34) where obviously R i is obtained in the measurement update, R ,i and S ,i are already calculated as the sensitivity vectors in the identification step.
We can simply apply the same type of triangularization to the last n b column to achieve the SRI for further data processing: Furthermore, the outlier parameters must be eliminated sooner. For a range outlier the introduced parameter can be removed directly, whereas the parameter for a phase observation could be either a cycle slip or an outlier. For a cycle slip, it can be replaced by the difference between the ambiguities before and after the cycle slip, then the previous ambiguity will be eliminated and only the new ambiguity is kept.
Assuming the ambiguities are n 1 and n 2 , the n 12 = n 2 − n 1 , we have n 1 = n 2 − n 12 . Replace n 1 in the Eq. (34) by n 2 − n 12 , we only have to triangularize the last n b columns, otherwise the triangularization must start from the column where n 1 is. After the triangularization of the last n b columns, the outlier parameters can be eliminated, and the information equation is ready for deriving an adapted solution and for the processing of the next epoch.

Realization of SRIF processing procedure
We upgraded the SRIF module in the PANDA software package, which was developed originally at Wuhan University, China, with the above-mentioned quality control approach, and the performance of the module is evaluated in an operational run at the GFZ IGS Real-Time analysis center. In this section we present the whole processing procedure of the upgraded SRIF module and its data flowchart. The key aspects related to the software realization and computational efficiency are emphasized.
The flowchart of the SRIF module is illustrated in Fig. 1. After the input of the information about stations and satellites involved in the processing and the corresponding processing setup, the initialization of the filter is carried out accordingly with the initial constraints and the normalized state equations. Then the processing of the epoch-wise observations starts, including reading real-time observations, data pre-processing, time update, measurement update, quality control and parameter estimation. The equations employed are listed on the right side of the processing steps, accordingly. It must be pointed out that the processing of the matrix transformation is realized based on the BLAS/ LAPACK library which will speed up the computation significantly according to the performance of the computer facility in use.
In addition to the processing procedure, the information matrix of the SRIF is also very important, as almost all the steps in the flowchart are working on this array and its change in the procedure is actually the data flowchart of the processing. Figure 2 shows the SRI matrix.
After the initialization, the initial constraints of all the parameters are transferred to R 0 and Z 0 which are put on the upper-right marked with R and Z 1 , respectively. The size for R is dynamically controlled by the total number of the active parameters n tot , then the dimension of R is n tot × (n tot + 1) . The parameters are divided into two groups: stochastic and deterministic parameters. The Fig. 1 Flowchart of the SRIF processing with the new quality control approach stochastic parameters will be eliminated epoch-by-epoch and they are put in the front of the others for saving computation cost of elimination. Worth to mention that the ambiguity parameters are kept at the end of the parameter list as they are maintained according to the actual observed satellites. The rows below R are for observation update and usually should be able to include all observations at one epoch, otherwise the update must be carried out in batches. In order to make the parameter elimination easy, the SRI matrix is extended for the parameters to be removed. While removing a set of parameters, their information is moved into the extend place and the original is prepared ready for the same parameter but for the next epoch, and the parameters can be removed after applying a corresponding triangularization. This can also be utilized for time update since its major processing is to remove the state parameters at the previous epoch after adding the state equations. In the observation update, the u vector with respect to the column-wise transformation matrix is saved in the same column but the low-left part for the computation of the sensitivity vectors used in the quality control.

Operational validation
Network A network of 85 globally distributed multi-GNSS stations is employed in the operational clock estimation at the GFZ real-time analysis center for the validation. All the stations can track both GPS and GLONASS signals, while 72 and 61 stations can track Galileo and BDS-2 signals, respectively. Most of the real-time observation streams are provided by IGS and the Multi-GNSS Experiment (MGEX), and several stations are provided by GFZ and the German Federal Agency for Cartography and Geodesy (BKG). Figure 3 shows the distribution of stations and the corresponding tracking status. It should be mentioned that the real-time data streams are not as stable as the recorded observation files because of the data communication. In general, at least 75 stations are available for clock estimation.

Software and processing strategy
The upgraded PANDA software package with the new quality control procedure is utilized for the validation processing on a Linux computer with eight Intel Xeon 3.5 GHz processors. The data processing strategy for multi-GNSS satellite clock estimation is summarized in Table 1. It should be mentioned that in the operational real-time clock estimation the ambiguities are estimated as float constant and not fixed.

Performance of quality control
The validation processing has been running since 1. June 2021 and there is no unexpected system crash or restart. The performance of the online quality control is investigated using the statistics of the detected outliers and computational cost.
Taking the result on 8 June 2021 as an illustration, the computation time at each epoch either with or without quality control is shown in Fig. 4. The epochs without outliers are marked in blue, and those with outliers in red. We can see that there is no significant increase in the computation time when the outliers are detected, and the average computation time is 2.06 s and 2.11 s for the epochs without and with outliers, respectively. Compared to the result by Fu et al. (2019), in which 0.4 s is additionally needed for quality control to an averaged epoch computation time of 2.4 s, the computation time for quality control of the new approach for SRIF in this study is significantly reduced.
In order to explore the distribution of the number of detected outliers and the relationship with computation time, the frequency of epochs with a certain number of outliers is retrieved from the results in the whole month of June and listed in Table 2. Among all the epochs, 67.63% are free of outliers due to the quality control  applied in the pre-processing, and the average computation time is only 2.07 s. Among the remaining epochs, about half of them are with 1 outlier, and the percentage of epochs with more than 10 outliers is about 0.3%. The computation time increases along with the number of detected outliers, in general the time is less than 3 s even with 9 outliers. For those epochs that suffer from a larger number of outliers, i.e., more than 10, the clock estimation usually can still be completed within 5 s with very few exceptions. However, if more stations and satellites are involved in the processing, it might be a problem to keep the processing time within 5 s to meet the current requirements.

Clock precision
For further evaluation, the operational real-time clock estimates are compared with the GBM products, the MGEX final products provided by GFZ (Deng et al., 2016). To remove the impact of clock datum used in different solutions, we formed the single-differenced clocks in each system between these two products, using G01, R01, E01, and C11 as the reference clock for GPS, GLO-NASS, GALILEO, and BDS, respectively. Moreover, we also removed the radial orbit differences between the different orbits used herein and by GBM. Figure 5 presents the monthly average STD of the clock differences for GPS, GLONASS, GALILEO, and BDS, from top to bottom. Note that we did not present the mean or RMS values of the clock differences as the systematic biases between different solutions are much larger than the STD value, but they have insignificant impact for the positioning users and can be compensated by the ambiguity parameters.
For the GPS satellites, the STD values are generally smaller than 0.1 ns and the average value is 0.06 ns. For the GLONASS satellites, the average STD value is 0.24 ns. The Galileo satellites have a comparable performance to GPS with an average STD 0.06 ns. As for BDS, the average STD value is 0.11 ns and 0.36 ns for the MEO/IGSO and GEO satellites, respectively. The relatively larger STD values of BDS GEO satellites are caused by the poor tracking geometry with regional network, and the poor orbit quality which has a meter-level agreement between different ACs. As shown in the figure, the real-time clock estimates have a good agreement compared to the postprocessing product, indicating the quality control procedure works well in the operational processing.

Conclusions
To meet the requirements in the high-rate update of precise satellite clock estimation for real-time precise positioning service, the SRIF estimator is suggested because of its high numerical stability and computational efficiency.
To improve the processing efficiency, we developed a new quality control procedure in SRIF which is capable of processing networks with many stations and observations. In the procedure, the most-like outliers are selected according to the posterior residuals one by one and are tested together based on the sensitivity vectors of the observations as candidate outliers. As there are usually only a few outliers at each epoch after the data pre-processing and only the sensitivity vectors of the candidates are calculated, the computation cost for the quality control is neglectable compared to the total epoch processing time.  The new quality control procedure is implemented into the SRIF module of the PANDA software package. An experimental validation has been running operationally at GFZ since 1 June 2021 for the estimation of satellite clocks using the data streams from a global network of 85 GNSS stations. For the integrated processing of GPS, GLONASS, Galileo and BDS, on average about 75 active stations and 72 satellites are processed. The computation time for each epoch is about 2.07 s when there are no outliers. There are about 67.63% of the total epochs without any outlier, and 15.58% with one outlier, and 0.3% with more than 10 outliers. The computation time increases along with the number of detected outliers but still stays within 3 s for all the epochs with less than 10 outliers. Among the other 0.3%, there are only very few epochs with rather large number of detected outliers and the processing time could be longer than 5 s, which should be further investigated.
The operational real-time satellite clocks show good agreement with the GBM final products, with an average of 0.06 ns, 0.24 ns, 0.06 ns, 0.11 and 0.36 ns for GPS, GLONASS, Galileo, BDS MEO&IGSO, and BDS GEO satellites, respectively. This result further confirms the efficiency of the SRIF processing and the new quality control method.