Rapid Estimation of Undifferenced Multi-GNSS Real-Time Satellite Clock Offset Using Partial Observations

Abstract

Real-time satellite clock offset is a crucial element for real-time precise point positioning (RT-PPP). However, the elapsed time for undifferenced (UD) multi-global navigation satellite system (GNSS) real-time satellite clock offset estimation at each epoch is increased with the growth of stations, which may fall short of real-time application requirements. Therefore, a rapid estimation method for UD multi-GNSS real-time satellite clock offset is proposed to improve the computation efficiency, in which both the dimension of the normal equation (NEQ) and the number of redundant observations are calculated before adjustment; if these two values are larger than the predefined thresholds, the elevation mask is gradually increased until they are less than the predefined thresholds. Then, the clock offset estimation is conducted; this method is called clock offset estimation using partial observations. Totals of 50, 60, 70 and 80 stations are applied to perform experiments. Compared to clock offset estimation using all observations, the elapsed times of clock offset estimation using partial observations can be reduced from 6.80 to 3.10 s, 7.93 to 2.97 s, 12.04 to 3.14 s for 60, 70 and 80 stations, respectively. By using the proposed method, the elapsed time of the clock offset estimation at each epoch is less than 5 s. The estimated clock offset accuracy for GPS, BDS-3, Galileo and GLONASS satellites are better than 0.04, 0.05, 0.03 and 0.16 ns when using the partial observations to estimate clock offset with 50, 60, 70 and 80 stations, respectively. For the multi-GNSS kinematic PPP using the estimated clock offset from 50, 60, 70 and 80 stations with partial observations, the positioning accuracy at 95% confidence level in the east, north and up direction are better than 2.70, 2.20 and 5.60 cm, respectively.

1. Introduction

The prerequisite of RT-PPP [1] is real-time satellite orbit and clock offset. RT-PPP is an efficient technology that can be used in many areas, such as precise agriculture [2], time transfer [3] and weather monitoring [4]. The accuracy of RT-PPP significantly relies on the performance of satellite orbit and clock offset, in which real-time satellite orbit can be predicted with high-precision up to several hours [5], whereas the satellite clock offset cannot be predicted precisely for several hours due to its complicated variation [6]. Therefore, it is necessary to estimate the high-performance and real-time satellite clock offset in real-time.

In 2001, the International global navigation satellite system (GNSS) Service (IGS) real-time working group was established, aiming to provide real-time products for real-time applications [7]; these real-time products (such as satellite orbit, clock offset and observation-specific biases) have been broadcast to global users since 1 April 2013. These products are also broadcast by many analysis centers [8]. The update rate for satellite clock offset from many analysis centers [9] is 5 s. Therefore, rapid estimation of satellite clock offset is necessary and it should be updated within 5 s.

There are two types of research for the rapid estimation of satellite clock offset: For the first, the clock offset estimation uses high-performance computer equipment. A linear algebra package (LAPACK) is used to conduct a matrix operation; compared to the matrix operation without LAPACK, the computation efficiency can be improved by 20.0% and 56.5% when 28 and 75 stations are applied to estimate BDS and GPS satellite clock offset, respectively [10]. The open multi-processing and Intel math kernel library technology were employed to estimate clock offset, and the improvement of computation efficiency was 76.6% [11]. Tao employed the graphics processing unit to accelerate the GPS/GLONASS/Galileo/BDS satellite clock offset estimation, and its time-consuming is about 0.86 s with 90 stations [12]. On the other hand, the rapid algorithm is investigated to accelerate clock offset estimation. Cao et al. used a single observation and recursive Kalman filter to avoid the complex inverse of the matrix: when the dimensions for a normal equation (NEQ) matrix is about 1300, the clock offset estimation time is about 0.5 s at each epoch [13]. Zuo et al. used Square Root Information Filter (SRIF) to estimate multi-GNSS real-time satellite clock offset; its computation time is 2.11 s at each epoch [14]. Based on this, Gong et al. conducted QR factorization, then applied SRIF to estimate clock offset: the elapsed time is about 0.5 s for each epoch [15]. Furthermore, the clock offset estimation also applies the dual-threads method, and the elapsed time is about 0.2 and 2 s for single and multi-GNSS at each epoch, respectively [16,17]. Some scholars employed the carrier range to estimate clock offset [18,19], and the estimation time is shorter than 1 and 2 s for single and multi-GNSS, respectively. There are limitations for these two categories: For the first one, high-performance processing equipment is required, which cannot be obtained by the common user. In terms of the latter, all selected observations are used for clock offset estimation; thus, with the increase of stations, data processing and quality control [20,21] become more and more challenging, and the elapsed time may also be longer. Therefore, it is necessary to investigate the method of rapid clock offset estimation.

We focus on the rapid estimation of UD multi-GNSS real-time satellite clock offset in this contribution. A rapid estimation method is proposed, in which the thresholds for the dimension of NEQ and the number of redundant observations are set, respectively. When the dimension of NEQ and the number of redundant observations is larger than the predefined threshold, the elevation mask is increased to remove some observations until the dimension of NEQ and the number of redundant observations is smaller than the thresholds. Thereafter, the adjustment is conducted. This method is called clock offset estimation using partial observations. The estimated clock offset is validated in terms of computation efficiency, clock offset accuracy and kinematic PPP. The satellite clock offset estimation model and clock offset estimation using partial observations are introduced in Section 2. Then, the experiment is designed and conducted, the computation efficiency, clock offset accuracy, kinematic PPP performance are analyzed in Section 3. Finally, the discussions and conclusion are given in Section 4 and Section 5, respectively.

2. Methods

The multi-GNSS satellite clock offset estimation model will be introduced first, then the proposed and in-depth description of the satellite clock offset estimate method using partial observations.

2.1. Multi-GNSS Satellite Clock Offset Estimation Model

The GNSS raw code $P$ and carrier phase $L$ observation can be expressed as:

$$P_{r,f}^{s,g}={\rho }_r^{s,g}+c\left(t_r^g-t^{s,g}\right)+c\left(d_{r,f}^g-d^{s,g}\right)+{\gamma }_f{I}_{r,1}^s+m_r^{s,g}{T}_r+{\epsilon }_{r,f}^{s,g}$$

(1)

$$L_{r,f}^{s,g}={\rho }_r^{s,g}+c\left(t_r^g-t^{s,g}\right)+{\lambda }_f^g\left(N_{r,f}^{s,g}+b_{r,f}^g-b^{s,g}\right)-{\gamma }_f{I}_{r,1}^s+m_r^{s,g}{T}_r+e_{r,f}^{s,g}$$

(2)

where $s$ and $g$ are the satellite, navigation satellite system, respectively. $r$ and $f$ indicate receiver and frequency, respectively. ${\rho }_r^{s,g}$ is the geometric distance with the unit of meters, $c$ indicates light speed. $t_r^g$ denotes the receiver clock, $t^{s,g}$ is the satellite clock offset. $d_{r,f}^g$ denotes the receiver code hardware delay, $d^{s,g}$ indicates the satellite code hardware delay, their units are seconds. $b_{r,f}^g$ is the receiver phase hardware delay, $b^{s,g}$ denotes the satellite phase hardware delay. ${\gamma }_f$ is the coefficient of ionospheric delay at frequency $f$. $I_{r,1}^s$ indicate the ionosphere delay at the first frequency. $T_r$ and $m_r^{s,g}$ are the zenith wet delay and its mapping function for the troposphere, respectively. $N_{r,f}^{s,g}$ is the ambiguity of carrier phase. ${\epsilon }_{r,f}^{s,g}$ and $e_{r,f}^{s,g}$ indicate the measurement noise and other errors of code and phase observations, respectively. To eliminate the impact of the first-order ionosphere delay, the dual-frequency ionosphere-free (IF) combination is used. The satellite phase center offset (PCO) and variation (PCV), relativistic effect and tide are corrected using the existing model [22]. Then, the error model is expressed as:

$$v_{r,pc}^{s,g}=c\left({\overline{t}}_r^g-{\overline{t}}^{s,g}\right)+m_r^{s,g}{T}_r+{\epsilon }_{r,pc}^{s,g}$$

(3)

$$v_{r,lc}^{s,g}=c\left({\overline{t}}_r^g-{\overline{t}}^{s,g}\right)+{\lambda }_{IF}^g{\overline{N}}_{r,IF}^{s,g}+m_r^{s,g}{T}_r+e_{r,lc}^{s,g}$$

(4)

where $v_{pc,r}^{s,g}$ and $v_{lc,r}^{s,g}$ are “observed values minus computed values“ (OMC) for pseudorange and carrier phase observation of each station–satellite pair, respectively. ${\overline{t}}_r^g$ and ${\overline{t}}^{s,g}$ are the reparametrized receiver and satellite clock offset, respectively. ${\overline{N}}_{r,IF}^{s,g}$ is the reparametrized carrier phase ambiguity [23,24].

The hardware delay and clock offset datum are different for different GNSS; this bias is called inter-system bias (ISB) [25], and should be considered. Therefore, the reparametrized UD error equation can be expressed as:

$$\left\{\begin{array}{c}\begin{array}{c}{v}_{r,pc}^{s,G}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,G}\right)+{m_r^{s,G}T}_r^s+{\epsilon }_{r,pc}^{s,G}\\ v_{r,pc}^{s,R}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,R}\right)+{ISB}_r^{R-G}+m_r^{s,R}{T}_r^s+{\epsilon }_{r,pc}^{s,R}\end{array}\\ \begin{array}{c}{v}_{r,pc}^{s,C3}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,C}\right)+{ISB}_r^{C-G}+m_r^{s,C}{T}_r^s+{\epsilon }_{r,pc}^{s,C}\\ v_{r,pc}^{s,E}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,E}\right)+{ISB}_r^{E-G}+m_r^{s,E}{T}_r^s+{\epsilon }_{r,pc}^{s,E}\end{array}\end{array}\right.$$

(5)

$$\left\{\begin{array}{c}\begin{array}{c}{v}_{r,lc}^{s,G}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,G}\right)+{\lambda }_f^G{\overline{N}}_{r,IF}^{s,G}+{m_r^{s,G}T}_r^s+e_{r,lc}^{s,G}\\ v_{r,lc}^{s,R}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,R}\right)+{{\lambda }_f^R\overline{N}}_{r,IF}^{s,R}+{ISB}_r^{R-G}+m_r^{s,R}{T}_r^s+e_{r,lc}^{s,R}\end{array}\\ \begin{array}{c}{v}_{r,lc}^{s,C3}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,C}\right)+{\lambda }_f^C{\overline{N}}_{r,IF}^{s,C}+{ISB}_r^{C-G}+m_r^{s,C}{T}_r^s+e_{r,lc}^{s,C}\\ v_{r,lc}^{s,E}=c\left({\overline{t}}_r^G-{\overline{t}}^{s,E}\right)+{\lambda }_f^E{\overline{N}}_{r,IF}^{s,E}+{ISB}_r^{E-G}+m_r^{s,E}{T}_r^s+e_{r,lc}^{s,E}\end{array}\end{array}\right.$$

(6)

where the symbols ‘G, R, C, E’ represent the GPS, GLONASS, BDS and Galileo, respectively. $IS{B}_r^{R-G}$, $IS{B}_r^{C-G}$ and $IS{B}_r^{E-G}$ are the ISB between GLONASS, BDS, Galileo and GPS, respectively. To reduce computation burden, the GLONASS inter-frequency bias is not considered [20]. Two rank deficiencies exist and should be eliminated. The linear relationship between receiver and satellite clock offset is the first one, a receiver clock is selected to eliminate it [26]; the second is the linear-dependency between satellite clock offset and ISB, the zero-mean condition [27] is applied to each system for stations, it can be expressed:

$$\left\{\begin{array}{c}{\displaystyle {\displaystyle \sum }_{r=1}^n}IS{B}_i^{R-G}=0\\ {\displaystyle {\displaystyle \sum }_{r=1}^n}IS{B}_i^{C-G}=0\\ {\displaystyle {\displaystyle \sum }_{r=1}^n}IS{B}_i^{E-G}=0\end{array}\right.$$

(7)

After imposing the constraints, the rank deficiency can be removed and the clock offset is estimable.

2.2. Satellite Clock Offset Estimation Using Partial Observations

With the development of the IGS, more and more stations can be used to perform GNSS precise data processing. The number of stations versus the year is presented in Figure 1: more than 500 IGS stations can be applied as of the end of 2021 (https://www.igs.org/network/#station-map-list, accessed on 11 October 2022), which provides abundant data resources for the GNSS data processing. However, it also brings a computation burden when many stations are applied to GNSS data processing. For the UD satellite clock offset estimation model, its computation efficiency depends on the elapsed time of the inverse of the NEQ. The elapsed time of the inverse of the NEQ for different dimensions is presented in Figure 2; the elapsed time increases exponentially as the increase of dimensions of the NEQ because many ambiguities should be estimated in UD satellite clock offset estimation. If the dimension of the NEQ is about 2000, the elapsed time is less than 2.5 s, while its elapsed time will be longer than 5 s when the dimension of an NEQ is larger than 2500. However, the update rate for real-time satellite clock offset is less than 5 s. Therefore, when using many stations to estimate clock offset, it cannot be promptly updated.

When estimating clock offset using 60, 80 and 100 stations, the clock offset accuracy difference is marginal [28]. If more stations are utilized to estimate clock offset, the estimated parameters are increased; with one observation of a station–satellite pair, i.e., adding one code and phase observation, an ambiguity parameter will be added. In UD satellite clock offset estimation, most estimated parameters are ambiguities, whereas ambiguities are byproducts. Therefore, GNSS satellite clock offset estimation using the partial observations method is proposed to improve the computation efficiency. In this method, many stations are selected, then, the thresholds for the dimension of the NEQ and the number of redundant observations is set. When these two indicators are larger than predefined thresholds, some observations are removed by increasing the elevation mask, and the remaining observations with higher elevations are applied to conduct clock offset estimation. It is anticipated the clock offset can be estimated as quickly as possible.

Figure 3 presents the data processing of the clock offset estimation using partial observations. First, the data is read (such as GNSS observations, satellite orbit, earth rotation paraments (ERP)). The elevation mask is set to 7° to remove the observations with low elevation. Next, the dimension of NEQ and the number of redundant observations is calculated, respectively (‘nobs’ and ‘ndim’ represent the number of observations and the dimension of the NEQ in the figure, respectively). From Figure 2, when the dimension of NEQ is about 2000, its computation time for clock offset estimation is about 2.5 s, which is less than 5 s. Therefore, the first condition is that the dimension of an NEQ is less than 2000. The second condition is that the difference between ‘nobs’ and ‘ndim’ is larger than 600, which is called the number of redundant observations; if this value is smaller than 0, the clock offset estimation cannot be estimated, whereas when it is larger, it will lead to time consumption. When these two conditions can be met, the clock offset is estimated; otherwise, the elevation mask is increased by 3° every time, some observations are removed, and this loop is iteratively conducted until two judgment conditions are met. Finally, the clock offset estimation is performed, outputting clock offset products.

3. Results

The experiment setup is introduced first. Next, the computation efficiency is evaluated and analyzed. Thereafter, the clock offset performance is evaluated. Finally, PPP is used to validate the estimated clock offset.

3.1. Experiment Setup

To validate the proposed method, multi-GNSS Experiments (MGEX) stations [29] from a global scale are used to experiment. The station distribution is shown in Figure 4, in which the purple, blue, green and red circles represent the 50, 60, 70 and 80 stations, respectively. The experiment is conducted from the day of the year (DOY) 009 to DOY 015 in 2022, the multi-GNSS (GPS/BDS-3/Galileo/GLONASS) satellite clock offset is estimated using dual-frequency IF combination with L1/L2, B1I/B3I, E1/E5a and G1/G2 frequency, respectively. The detailed processing strategies are shown in

Table 1. Data process for satellite clock offset estimation.

Items	Strategies
Observations	GPS: L1/L2; Galileo: E1/E5a; BDS-3: B1I/B3I; GLONASS: G1/G2
Weight	E ≤ 30°, P = 2sinE; E > 30°, P = 1;
Satellite orbit/ERP	Center for Orbit Determination in European (CODE) [30]
Relativistic effect	Corrected [31]
Wind-up	Corrected [32]
Satellite DCB	CAS [33]
Satellite PCO/PCV	igs14.atx
Station coordinates	Fixed to IGS weekly solutions
Station displacement	Solid tide, ocean tide, pole tide: IERS Convention 2010 [34]
Receiver PCO/PCV	igs14.atx
Estimator	Sequential least square adjustment [20]
Satellite and receiver clock offset	White noise
ISB	White noise
Troposphere	Saastamoinen model [35] + GMF [36], estimated as a piece-wise constant every hour
Ambiguity	Float solution, estimated as a constant for a station-satellite pair if no cycle slip

. The experiment schemes are as follows: multi-GNSS satellite clock offset estimation using all observations with 50, 60, 70 and 80 stations, respectively; multi-GNSS satellite clock offset estimation with 50, 60, 70 and 80 stations using partial observations, respectively. The experiment is conducted on a laptop with an Intel i7-11370h processor.

3.2. Computation Efficiency

The computation efficiency depends on the dimension of the NEQ. Figure 5 shows the dimension of the NEQ at each epoch for clock offset estimation using all and partial observations with 50, 60, 70 and 80 stations, respectively. When satellite clock offset estimation using all observations, the dimension of the NEQ increases as the number of stations increases, because one carrier phase is added, an ambiguity parameter will be added. The number of dimensions is between 1815 and 1929, 2147 and 2297, 2490 and 2656, 2816 and 3020 for 50, 60, 70 and 80 stations, respectively. For clock offset estimation using partial observations, when employing 50 stations, the dimension of the NEQ is less than 2000 for each epoch, which is less than the predefined threshold, so there is no need to delete observations. When 60, 70 and 80 stations are applied to clock offset estimation, the dimension of the NEQ is still less than 2000, because the predefined thresholds for the dimension of the NEQ and the number of redundant observations is set; by increasing the elevation mask, some observations are removed, which is beneficial for improving computation efficiency. The mean dimension of the NEQ is presented in

Table 2. The mean number of dimensions of satellite clock offset estimation using all and partial observations with 50, 60, 70 and 80 stations.

Number of Stations	All Observations	Partial Observations	Reduction Rate
50	1871	1871	0%
60	2221	1932	13.0%
70	2575	1901	26.2%
80	2919	1890	35.3%

. For clock offset estimation using all and partial observations, the dimensions of NEQ for 60, 70 and 80 stations are 2221 and 1932, 2575 and 1901, 2919 and 1890, respectively. Compared to clock offset estimation using all observations, the decreased rate of the mean dimension of the NEQ is 13.0%, 26.2% and 35.3%, respectively.

The number of stations that can be tracked by each GNSS satellite using all and partial observations with 50, 60, 70 and 80 stations are shown in Figure 6. For clock offset estimation using all observations, the tracked stations by each satellite are increased as the growth of number of the stations. When 50 stations are applied, about 15 stations can be tracked for most satellites, while nearly 25 stations can be tracked by each satellite when using 80 stations. The tracked stations for some BDS-3 satellites are fewer than that for other systems, which may be caused by the launch date of these satellites being later; the receiver of stations is not being updated in time and they cannot receive the BDS-3 signals. For some satellites, the number of tracked stations is fewer compared to other satellites (such as G11, G22, G28, E14, E18, R11). For clock offset estimation using partial observations, about 15 stations can be tracked by each GNSS satellite; the number of tracked stations is almost the same whether 50, 60, 70 or 80 stations are applied. Therefore, it is anticipated that the elapsed time is almost the same for clock offset estimation using 50, 60, 70 and 80 stations.

The elevation masks for clock offset estimation using partial observations with 50, 60, 70 and 80 stations are shown in Figure 7. When using 50 stations, both the dimensions of the NEQ and the number of redundant observations were less than the predefined thresholds. Therefore, all observations were retained, and the elevation mask is 7 degrees for each epoch. When 60, 70 and 80 stations are applied, the elevation mask is between 13 and 16, 19 and 22, 25 and 28 degrees, and the mean is 14.62, 21.92 and 26.92 degrees, respectively.

The computation time for clock offset estimation using all and partial observations with 50, 60, 70 and 80 stations at each epoch is shown in Figure 8. For clock offset estimation using all observations, because the dimension of the NEQ is less than 2000 for 50 stations, it can be estimated within 5 s. When 60 to 80 stations are applied, the elapsed time is gradually increased; it is between 5 and 10, 5 and 12.5, 10 and 15 s at each epoch for 60, 70 and 80 stations, respectively—the clock offset cannot be updated within 5 s. In the case of clock offset estimation using partial observations, the dimensions of the NEQ are less than 2000 for 60, 70 and 80 stations; therefore, the elapsed time of each epoch is still less than 5 s, and the real-time applications can be met. The mean elapsed time of clock offset estimation for each epoch using all and partial observations with different stations are presented in

Table 3. Mean elapsed time of clock offset estimation for each epoch using all and partial observations with different stations (unit: second).

Number of Stations	All Observations	Partial Observations	Reduction Rate
50	2.86	2.86	0%
60	6.80	3.10	54.4%
70	7.93	2.97	62.6%
80	12.04	3.14	73.9%

. Compared to clock offset estimation using all observations, the mean elapsed time of clock offset estimation using partial observations can be reduced from 6.80 to 3.10 s, 7.93 to 2.97 s, 12.04 to 3.14 s for 60, 70 and 80 stations, respectively. The reduction rate is 54.4%, 62.6% and 73.9%, respectively. The computation time and improvements show that the proposed method is effective and feasible.

3.3. Clock Offset Accuracy

The clock offset accuracy of each satellite using all and partial observations with 50, 60, 70 and 80 stations is presented in Figure 9. When evaluating the clock offset, the CODE final clock offset products are selected, the clock offset accuracy evaluation method is referred to [6]. For different schemes, the clock offset accuracy for each GPS, BDS and Galileo satellite outperform 0.08, 0.09 and 0.05 ns, respectively. Most GLONASS satellite clock offset accuracies are better than 0.20 ns. Compared to other systems, the GLONASS satellite clock offset performance is poorer because the code inter-frequency bias is neglected [20] and the frequency stability of the cesium atomic clock equipped on GLONASS satellites is poorer than that of other GNSS [37]. The accuracy difference of each Galileo satellite is the smallest among the four GNSS due to their superior atomic clock performance, while the accuracy difference of each GLONASS satellite is the largest. For clock offset estimation using all observations, the accuracy difference is small when using 50, 60, 70 and 80 stations. The improvement of accuracy is marginal with the increase in stations. In terms of partial observations, except for the R09 satellite clock, when 50, 60, 70 or 80 stations are applied, the clock offset accuracy for other GLONASS satellites is almost the same. Furthermore, compared to clock offset estimation using all and partial observations, the clock offset accuracy difference is small whether 50, 60, 70, or 80 stations are employed.

Figure 10 and

Table 4. Mean clock offset accuracy for each system using different schemes (unit: ns).

Number of Stations	GPS		BDS-3		Galileo		GLONASS
	All	Partial	All	Partial	All	Partial	All	Partial
50	0.0360	0.0360	0.0381	0.0381	0.0285	0.0285	0.1215	0.1215
60	0.0350	0.0351	0.0375	0.0378	0.0262	0.0265	0.1173	0.1202
70	0.0338	0.0358	0.0376	0.0416	0.0261	0.0275	0.1157	0.1255
80	0.0322	0.0338	0.0376	0.0455	0.0256	0.0278	0.1155	0.1512

show the estimated clock offset accuracy of each system and the mean, respectively. For all observations of GPS, BDS-3 and Galileo satellites, their clock offset accuracy is better than 0.05 ns, the value is between 0.03 and 0.04, 0.03 and 0.04, 0.02 and 0.03 ns, respectively, which shows fairly good performance. Because GLONASS code inter-frequency bias is neglected and the frequency stability of cesium atomic clocks is inferior, compared to other systems, their satellite clock offset accuracy is inferior. When different stations are applied, the clock offset accuracy is between 0.11 and 0.13 ns. With the increase of stations, the clock offset accuracy improvement is very marginal. For partial observations, the GPS, BDS-3 and Galileo satellite clock offset accuracy is between 0.03 and 0.04, 0.03 and 0.05, 0.02 and 0.03, 0.12 and 0.15 ns, respectively. The clock offset accuracy is slightly poorer when applying partial observations compared to the clock offset estimated from all observations, which can be attributed to the fact that when partial observations are applied, the elevation mask is larger, the observation arc of each station–satellite pair is shorter, the ambiguity initialization is more frequent than clock offset estimation using all observations. Nevertheless, the clock offset accuracy difference for each system is marginal, and this difference is acceptable.

3.4. PPP Validation

The multi-GNSS kinematic PPP experiment is performed using satellite clock offset from all and partial observations with 50, 60, 70 and 80 stations. A total of 85 stations are used to conduct PPP; these stations are not applied to clock offset estimation, the geographical distribution of which is shown in Figure 11. The experiments are conducted from DOY 010 to DOY 015, 2022. The down-weight is conducted for GLONASS observations because of the poorer clock offset accuracy and the neglect of code IFB. The estimated parameters are station coordinates, ZWD, receiver clock offset, ISB, ambiguities. The station coordinates are compared to IGS weekly solutions.

The mean positioning accuracy for each station using clock offset products from different schemes is presented in Figure 12. For PPP using clock offset products with all observations, the positioning accuracy difference is similar whether the estimated clock offset from 50, 60, 70, or 80 stations. The positioning accuracy in east, north and up directions is between 1 and 2 cm, 0 and 2 cm, 1 and 6 cm for most stations, respectively. Although the clock offset is estimated using different stations, the positioning accuracy is almost the same because the clock offset accuracy difference is marginal when 50, 60, 70 and 80 stations are applied to clock offset estimation. For PPP using clock offset products with partial observations, compared to the positioning accuracy using the clock offset from 50 and 60 stations, the positioning accuracy using the clock offset from 70 and 80 stations is slightly worse. This is caused by the slightly worse clock offset accuracy estimated from 70 and 80 stations, especially for GLONASS satellite clock offset, their clock offset accuracy is gradually increased with the increase of stations. Compared to the clock offset from all observations, the positioning accuracy is slightly degraded when clock offset estimation using partial observations with 60, 70 and 80 stations, especially for the up component, because the clock offset is related to the up direction in PPP.

The positioning accuracy using the clock offset from different schemes at the 95% confidence level is shown in

Table 5. Positioning accuracy using the clock offset from different schemes at the 95% confidence level (unit: cm).

Number of Stations	All Observations			Partial Observations
	E	N	U	E	N	U
50	2.30	1.90	5.35	2.30	1.90	5.35
60	2.32	1.87	5.35	2.28	1.88	5.44
70	2.40	1.87	5.24	2.60	1.99	5.52
80	2.29	1.91	5.24	2.68	2.15	5.58

. For all observations, the positioning accuracy using the clock offset estimated from 50, 60, 70 and 80 stations in east, north and up directions are between 2.29 and 2.40, 1.87 and 1.91, 5.24 and 5.35 cm, respectively. The positioning accuracy difference is at the millimeter-level for different stations. Therefore, the positioning accuracy improvement is limited with the growth of stations for clock offset estimation. For partial observations, with the increase of stations, the positioning accuracy is gradually decreased for the up component, the difference is also at the millimeter level, which is acceptable for the kinematic PPP of centimeter-level.

4. Discussions

The computation time is still one of the problems for UD real-time clock offset estimation because the inverse of NEQ is time-consuming, which hinders real-time applications such as RT-PPP since the real-time satellite clock offset should be updated within 5 s. It has been proved that when more and more stations are used for clock offset estimation, the improvement of clock offset accuracy is limited. In this contribution, we investigated the rapid estimation of multi-GNSS real-time satellite clock offset. The thresholds for the number of NEQ and redundant observations are set. Increasing the elevation mask to remove some observations to improve the computation efficiency. The results demonstrated that the elapsed time of clock offset estimation is less than 5 s when 60, 70 and 80 stations are applied. The estimated clock offset accuracy for GPS, BDS-3, Galileo and GLONASS satellites are better than 0.04, 0.05, 0.03 and 0.16 ns.

There are several characteristics of this method. First, although many stations are applied to estimate clock offset, the elapsed time is still within 5 s, which can meet the requirement of the update for clock offset estimation within 5 s. Second, the real-time observations are impacted by various factors such as station environment distribution, network delay, resulting in the observation not being received. When a few stations are applied to estimate the clock offset, if some observations cannot be received, the clock offset cannot be estimated because the number of stations is too few, so the performance of the clock offset will be degraded. If clock offset estimation uses many stations, when all observations can be received, the proposed method can be applied, if some observations cannot be received, the clock offset can still be estimated. However, with the increase of the stations, the elevation mask of observation is larger, resulting in the observation arc being shorter, the initialization of ambiguity is more frequent. Therefore, the clock offset is slightly degraded with the increase of the stations. It is suggested that 70 to 80 stations are applied. When the observations of all stations can be received, the proposed method can be employed. Once the stations are impacted by their environment distribution, network delay, or other factors, some observations cannot be received, estimating the clock offset is still possible.

5. Conclusions

Real-time satellite clock offsets are the prerequisites for RT-PPP. A rapid estimation method of UD multi-GNSS real-time satellite clock offset using partial observations is proposed in this contribution. In this method, the dimension of the NEQ and the number of redundant observations is calculated before clock offset estimation, respectively, when these two values are smaller than the predefined thresholds, then clock offset estimation is performed. Otherwise, the elevation mask is gradually increased until they are smaller than the predefined thresholds. Because the dimension of the NEQ can be reduced, the computation efficiency can be significantly improved compared to clock offset estimation using all observations. Observations from 50, 60, 70 and 80 MGEX stations are applied for validating the proposed method.

The dimensions of NEQ are analyzed. For all and partial observations, the dimensions of NEQ are 2221, 2575, 2919 and 1932, 1901 and 1890 for 60, 70 and 80 stations, respectively. Compared to the clock offset estimation using all observations, the reduced rate of the proposed method is 13.0%, 26.2%, 35.3%, respectively. Compared to clock offset estimation using all observations, the elapsed times of clock offset estimation using partial observations can be reduced from 6.80 to 3.10, 7.93 to 2.97, 12.04 to 3.14 s for 60, 70 and 80 stations, respectively. Therefore, the computation time of the clock offset estimation is less than 5 s. When using the proposed method, the clock offset accuracy of GPS, BDS-3, Galileo and GLONASS satellites are outperform 0.04, 0.05, 0.03 and 0.16 ns. The positioning accuracy of multi-GNSS kinematic PPP at a 95% confidence level is better than 2.70, 2.20 and 5.60 cm for different schemes. The results demonstrate that rapid estimation and reliable clock offset can be achieved using the proposed method.