## 3 Robust Helmert Variance Component Estimation (RVCE)  
In Equation (5), the weight matrix is the combined weight of the three GNSSs. As the observation qualities of individual GNSSs are not the same due to the different signal structures, frequencies, accuracies of satellite orbit and clock products, signal-to-noise ratios and signal transmission paths, neither individual GNSSs nor the individual satellites should be treated equally [15; 19]. Although the existing LS-VCE method can be employed to allocate the weight of different kinds of observation [25], it has lack of robust function. Moreover, in the area of precise time and frequency transfer based on multi-GNSSs, outliers inevitably exist in observations [26; 27]. These outliers will further distort the results of receiver clock estimation, and do not reasonably represent the characteristics of the time and frequency references. Therefore, a robust approach based on the Helmert VCE is introduced in this study to allocate the appropriate weights and to control for the influence of outliers, which are generally referred to as the posterior measurement residuals. The error equation of the three combined GNSSs for time and frequency can be written as:  
\[V=A\cdot\hat{X}-L=\left[ \begin{array}{c}V_{G}\\ V_{C}\\ V_{E}\end{array} \right]=\left[ \begin{array}{c}A_{G}\\ A_{C}\\ A_{E}\end{array} \right]\cdot\hat{X}-\left[ \begin{array}{c}L_{G}\\ L_{C}\\ L_{E}\end{array} \right] \tag{7}\]  
where, \(V\) denotes the combination posterior measurement residuals, which can be further written as the residuals of the individual systems (i.e., GPS, BDS and Galileo) and \(\hat{X}\) is the estimated vector of unknown parameters. Helmert VCE (i.e., posterior variance stochastic model estimation) is based on the residuals derived from the initial inappropriate weights (i.e., prior equal weights) for the observations. In the practical application of combining the three GNSSs for time and frequency transfer, the range of the GNSS residual values may differ due to their different systematic errors. The individual residuals therefore must be unified to a common range using their corresponding variances and cofactors, the formula for this unification is given as:  
\[\widetilde{V}(i)=\frac{V(i)}{\sigma_{0}^{2}\sqrt{Q(i)}} \tag{8}\]  
where \(\widetilde{V}\) denotes the standard residual matrix. The unit weight variance is \(\sigma_{0}^{2}\) and \(Q\) is the corresponding cofactor, which can be determined using Equations (9) and (10), respectively:  
\[\sigma_{0}^{2}=\frac{V^{T}PV}{n-m} \tag{9}\]  
\[Q=P^{-1}-AN^{-1}A^{T} \tag{10}\]  
where, \(n\) is the number of available GNSS satellites, \(m\) denotes the estimated parameter number and \(N\) is the matrix of the normal equation. Therefore, according to the Helmert VCE theory, the estimated variance components can be determined using the following formulae:  
\[\left[ \begin{array}{c}\hat{\sigma}_{0,G}^{2}\\ \hat{\sigma}_{0,C}^{2}\\ \hat{\sigma}_{0,E}^{2}\end{array} \right]=\left[ \begin{array}{ccc}s_{G,G}&s _{G,C}&s_{G,E}\\ s_{G,C}&s_{C,C}&s_{C,E}\\ s_{G,E}&s_{C,E}&s_{E,E}\end{array} \right]^{-1}\left[ \begin{array}{c}\widetilde {V}_{G}^{T}P_{G}\widetilde{V}_{G}\\ \widetilde{V}_{C}^{T}P_{C}\widetilde{V}_{C}\\ \widetilde{V}_{E}^{T}P_{E}\widetilde{V}_{E}\end{array} \right] \tag{11}\]  
\[s_{k,k}=n_{k}-2 tr\left(N^{-1}N_{k}\right)+tr\left(N^{-1}N_{k}N^{-1}N_{k}\right)\] \[s_{k,j}=n_{k}-tr\left(N^{-1}N_{k}N^{-1}N_{j}\right)\] \[N=A^{T}PA \tag{12}\]  
where, \(tr\) is the trace operation for matrix. The initial weights for three GNSSs are typically inappropriate, and the corresponding variance components \(\hat{\sigma}_{0,G}^{2}\), \(\hat{\sigma}_{0,C}^{2}\), and \(\hat{\sigma}_{0,E}^{2}\) are also unequal.