**1. Introduction**

Precise satellite orbits and clock offsets are essential for obtaining precise tempo-spatial information in precise point positioning (PPP), which plays an important role in precise orbit determination (POD) for low Earth orbit (LEO) satellites, atmospheric retrieval, precise positioning, timing, and so on. The International GNSS Service (IGS) has officially provided

**Citation:** Jiao, G.; Song, S. High-Rate One-Hourly Updated Ultra-Rapid Multi-GNSS Satellite Clock Offsets Estimation and Its Application in Real-Time Precise Point Positioning. *Remote Sens.* **2022**, *14*, 1257. https:// doi.org/10.3390/rs14051257

Academic Editor: Yunbin Yuan

Received: 30 December 2021 Accepted: 2 March 2022 Published: 4 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

final, rapid, and ultra-rapid products since 5 November 2000 [1]. As summarized in IGS product introductions (https://igs.org/products/#about, accessed on 1 March 2022), the IGS final products with highest precision are delayed for 12~18 days; the rapid products are delayed for 17~41 h [2]. With the development of global navigation satellite systems (GNSSs) and the increasing requirements for positioning accuracy and low latency, the final and rapid products can hardly meet the real-time and near-real-time user needs. The ultra-rapid products with low latency and high accuracy are of great interest in GNSS real-time and near-real-time fields. IGS integrated ultra-rapid GPS-only products from eight analysis centers (AC), Center for Orbit Determination in Europe (CODE), Natural Resources Canada (NRCan), and so on, which are updated every six hours [3,4]. At present, the latency and time resolution of the ultra-rapid products released by most ACs are 6 h and 15 min, respectively. The ultra-rapid orbits and clock offsets consist of observation and prediction sessions. The high latency means that we need to predict the satellite orbit and clock offsets for a longer time. With the increase in prediction time, the accuracy of satellite orbits and clock offsets will be reduced correspondingly. Especially for satellite clock offsets prediction, the accuracy of satellite clock offsets will be reduced with the longer prediction time. Therefore, with the six-hourly updated ultra-rapid GPS-only products, it is difficult to meet the high-precision real-time applications due to the low sampling and accuracy loss of clock offsets prediction.

To support better real-time services, IGS established the Real Time Working Group (RTWG) and provided the real-time service (RTS) to the GNSS community [5]. The Center National d'Etudes Spatiales (CNES) adopted the Kalman Filter (KF) estimator for satellite precise clock estimation (PCE) based on the raw observations [6]. The EPOS-RT and PANDA software used the Square Root Information Filter (SRIF) to obtain the real-time clock offsets based on the mixed-differenced model [7–9]. Whether in the raw observation receiving or real-time product broadcasting, a good network is essential for the real-time PCE. However, the available observations for real-time product solution will be affected by the network environment, and its corresponding product accuracy will be decreased [5]. In addition, it is difficult for users with poor network environment and communication ability to receive real-time precise products. On the contrary, the ultra-rapid products estimation is not expected to require high network conditions. It is not necessary for server and client to receive data continuously. All the same, the high latency and low samplings for ultra-rapid products are the problems that cannot be ignored. Therefore, it is important to provide the ultra-rapid orbit and clock products with low latency, high precision, and high samplings. To reduce the latency, Li et al. [10] used multi-thread technology and decreased the sampling rate (5 to 10 min) to estimate one-hourly updated multi-GNSS products based on undifferenced (UD) ionosphere-free (IF) models. Similarly, Chen et al. [11] updated data processing strategy and optimized software to obtain one-hourly updated multi-GNSS products with 5 min sampling. Wuhan University (WHU) GNSS AC released one-hourly updated multi-GNSS products with 5 min sampling [12,13]. Although, the latency and the samplings of ultra-rapid clock offsets has been improved from 6 to 1 h and 15 to 5~10 min, respectively. However, the samplings are not adequate for some applications. The time resolution of satellite clock offsets will affect the PPP in terms of convergence time and positioning accuracy, especially for the kinematic PPP applications [7,14,15]. To meet real-time needs and improve the positioning performance, it is meaningful to develop the high-rate one-hourly ultra-rapid PCE.

Furthermore, the single GPS is gradually developing to multi-constellation GNSS [16]. Multi-GNSS observations can improve the positioning accuracy, convergence time, and the stability and reliability of GNSS services [17]. However, the development of multi-GNSS has brought severe challenges to GNSS data processing. Figure 1 depicts the number of estimated parameters for the different PCE models with the increase in the number of the tracking stations. Due to the large number of estimation parameters (including ambiguities) in multi-GNSS, the computation is time-consuming. Bock et al. [18] proposed the epoch-difference (ED) PCE model to achieve high-rate GPS-only PCE for final products. In addition, most GNSS ACs used this model to obtain high-rate rapid and final products. Although the ED model greatly reduces the time-consuming, the current ED model and data processing strategy are difficult to realize the high-rate one-hourly updated ultra-rapid PCE. The research on high-rate one-hourly updated ultra-rapid multi-GNSS PCE is slightly insufficient. To meet the real-time and near-real-time requirements for multi-GNSS PPP and ensure the stability and reliability of GNSS precise service, it is necessary to investigate the multi-GNSS high-rate one-hourly updated ultra-rapid multi-GNSS PCE.

**Figure 1.** The number of estimated parameters for the different PCE models with the increase in the number of the tracking stations, in which UD PCE, ED PCE, and MED PCE denote undifferenced ionosphere-free PCE model, epoch-difference PCE model, and modified epoch-difference PCE, respectively. Suppose a station can track six satellites for each constellation in one epoch.

With this background, we introduce an efficient model and design a new framework for high-rate one-hourly updated ultra-rapid PCE. We investigate the influence of tropospheric delay on ED clock offsets based on the traditional ED model and further propose the modified ED PCE model to reduce the latency and improve calculation efficiency. Furthermore, Open Multi-Processing (OpenMP) and Intel Math Kernel Library (Intel MKL) technologies are used to improve the PCE processing efficiency and improve the time resolution. First, we present the mathematical models of multi-GNSS modified ED PCE and the methods of restoring ED clock offsets to absolute clock offsets in Section 2. Then, the data processing strategies of multi-GNSS modified ED PCE and new data processing framework are depicted in Section 3. In the following section, we verify the accuracy of the estimated clock offsets. Furthermore, the real-time PPP (RT-PPP) performance using real-time products and high-rate one-hourly updated products are also compared. Finally, we summarize the contributions and give some recommendations.

#### **2. Methods**

This section begins with the single-frequency observation models. Then, the modified ED models using GPS, BDS-2, BDS-3, GLONASS, and Galileo observations are developed in detail. Additionally, the methods of restoring ED clock offsets to absolute clock offsets are introduced.

The linearized equations of single-frequency observation equation can be expressed as [19]:

$$\begin{cases} P\_{r\_{\vec{\gamma}}}^{\varepsilon} = R\_r^{\varepsilon} + \mathfrak{c} \cdot \left( \delta t\_{\mathbf{r}} + b\_{\mathbf{j}}^{\varepsilon} + \widetilde{b\_{\mathbf{j}}^{\varepsilon}} \right) - \mathfrak{c} \cdot \left( \delta t^{\mathbf{f}} + b\_{\mathbf{j}}^{\varepsilon} + \widetilde{b\_{\mathbf{j}}^{\varepsilon}} \right) + g\_{\mathbf{w}} \cdot \tau\_{\mathbf{w}} + I\_{r\_{\vec{\gamma}}}^{\varepsilon} + \varepsilon p\_{\mathbf{j}} \\ \Phi\_{r\_{\vec{\gamma}}}^{\varepsilon} = R\_r^{\varepsilon} + \mathfrak{c} \cdot \left( \delta t\_{\mathbf{r}} + B\_{\mathbf{j}}^{\varepsilon} + \widetilde{B}\_{\mathbf{j}}^{\varepsilon} \right) - \mathfrak{c} \cdot \left( \delta t^{\mathbf{s}} + B\_{\mathbf{j}}^{\mathbf{s}} + \widetilde{B}\_{\mathbf{j}}^{\mathbf{s}} \right) + g\_{\mathbf{w}} \cdot \tau\_{\mathbf{w}} - I\_{r\_{\vec{\gamma}}}^{\mathbf{s}} + \lambda\_{\mathbf{j}} \cdot N\_{r\_{\vec{\gamma}}}^{\mathbf{s}} + \varepsilon\_{\Phi\_{\rangle}} \end{cases} (1)$$

where *r*, *s*, and *j* depict the receiver, satellite, and the frequency band, respectively; *P* and Φ are the range and phase observations; *R<sup>s</sup> <sup>r</sup>* is the geometrical range between satellite to receiver; *c* is the lightspeed; *δtr* and *δt <sup>s</sup>* denote the receiver and satellite clock offsets; *b<sup>r</sup> j* and *b<sup>s</sup> <sup>j</sup>* represent the receiver and satellite constant time-invariant uncalibrated code delays (UCDs); *B<sup>r</sup> <sup>j</sup>* and *<sup>B</sup><sup>s</sup> <sup>j</sup>* depict the corresponding constant time-invariant uncalibrated phase delays (UPDs); !*b<sup>r</sup> j* , !*bs j* , *B*!*r <sup>j</sup>* , and *<sup>B</sup>*!*<sup>s</sup> <sup>j</sup>* represent time-varying part; *τ<sup>w</sup>* and *gw* depict the zenith wet delay (ZWD) and the corresponding wet mapping function; *I<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the slant ionospheric delay; *N<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* denote the integer ambiguity with its wavelength *λj*; *εPj* and *ε*Φ*<sup>j</sup>* are the range and phase observation noises containing multipath and unmodeled error.

Most IGS and International GNSS Monitoring and Assessment System (iGMAS) ACs adopted undifferenced IF model to obtain the satellite clock offsets [1], which can be expressed as:

$$\begin{cases} P\_{r,IF\_{\hat{i},\hat{j}}}^{s} = R\_r^s + \boldsymbol{\c} \cdot \delta \mathbf{t}\_{IF\_{\hat{i},\hat{j}}}^r - \boldsymbol{\c} \cdot \delta \mathbf{t}\_{IF\_{\hat{i},\hat{j}}}^s + \boldsymbol{g}\_w \cdot \boldsymbol{\pi}\_w + \boldsymbol{\varepsilon}\_{P\_{IF\_{\hat{i},\hat{j}}}}\\ \Phi\_{r,IF\_{\hat{i},\hat{j}}}^s = R\_r^s + \boldsymbol{\c} \cdot \delta \mathbf{t}\_{IF\_{\hat{i},\hat{j}}}^r - \boldsymbol{\c} \cdot \delta \mathbf{t}\_{IF\_{\hat{i},\hat{j}}}^s + \boldsymbol{g}\_w \cdot \boldsymbol{\pi}\_w + \lambda\_{IF\_{\hat{i},\hat{j}}} \cdot \mathbf{N}\_{r,IF\_{\hat{i},\hat{j}}}^s + \boldsymbol{\varepsilon}\_{\Phi\_{I\hat{i},\hat{j}}} \end{cases} \tag{2}$$

The IF combination for observations and hardware delays can be described as [19]:

$$\begin{cases} \mathbf{a}\_{i,j} = \frac{f\_i^2}{f\_i^2 - f\_j^2}; \boldsymbol{\beta}\_{i,j} = \frac{-f\_i^2}{f\_i^2 - f\_j^2} \\ (\cdot)\_{I\mathbf{F}\_{i,j}} = \mathbf{a}\_{i,j} \cdot (\cdot)\_i + \boldsymbol{\beta}\_{i,j} \cdot (\cdot)\_j \\ (\cdot) = \mathbf{P}\_r^s, \Phi\_r^s, \mathbf{b}^r, \mathbf{b}^s, \mathbf{B}^r, \mathbf{B}^s, \overline{\mathbf{b}}^r, \overline{\mathbf{b}}^s, \overline{\mathbf{b}}^r, \overline{\mathbf{b}}^s \\ \delta \mathbf{I}\_{I\mathbf{F}\_{i,j}} = \delta t\_r + \mathbf{b}^s\_{I\mathbf{F}\_{i,j}} + \overline{\mathbf{B}}\_{I\mathbf{F}\_{i,j}}^s \\ \delta t\_{I\mathbf{F}\_{i,j}}^s = \delta t^s + \mathbf{b}^s\_{I\mathbf{F}\_{i,j}} + \overline{\mathbf{B}}\_{I\mathbf{F}\_{i,j}}^s \\ \lambda\_{I\mathbf{F}\_{i,j}} \cdot \mathbf{N}\_{r,I\mathbf{F}\_{i,j}}^s = \left( \mathbf{B}\_{I\mathbf{F}\_{i,j}}^s - \mathbf{b}\_{I\mathbf{F}\_{i,j}}^s \right) - \left( \mathbf{B}\_{I\mathbf{F}\_{i,j}}^s - \mathbf{b}\_{I\mathbf{F}\_{i,j}}^s \right) + \mathbf{a}\_{i,j} \cdot \lambda\_j \cdot \mathbf{N}\_{r,j}^s + \beta\_{i,j} \cdot \lambda\_i \cdot \mathbf{N}\_{r,i}^s \end{cases} (3)$$

where *αi*,*<sup>j</sup>* and *βi*,*<sup>j</sup>* are frequency factors; *fi* and *fj* depict the *i*th and *j*th frequency; *P<sup>s</sup> r*,*IFi*,*<sup>j</sup>* , Φ*s r*,*IFi*,*<sup>j</sup>* , *b<sup>r</sup> IFi*,*<sup>j</sup>* , *b<sup>s</sup> IFi*,*<sup>j</sup>* , *B<sup>r</sup> IFi*,*<sup>j</sup>* , *B<sup>s</sup> IFi*,*<sup>j</sup>* , !*br IFi*,*<sup>j</sup>* , !*bs IFi*,*<sup>j</sup>* , *<sup>B</sup>*!*<sup>r</sup> IFi*,*<sup>j</sup>* , and *<sup>B</sup>*!*<sup>s</sup> IFi*,*<sup>j</sup>* are IF combinations for the corresponding observations and hardware delays, respectively; *NIF* denote the float ambiguity; *<sup>ε</sup>PIFi*,*<sup>j</sup>* and *<sup>ε</sup>*Φ*IFi*,*<sup>j</sup>* are IF observations noises for pseudorange and carrier phase, and *<sup>ε</sup>PIFi*,*<sup>j</sup>* will absorb !*b<sup>r</sup> IFi*,*<sup>j</sup>* − *<sup>B</sup>*!*<sup>r</sup> IFi*,*<sup>j</sup>* + *<sup>B</sup>*!*<sup>s</sup> IFi*,*<sup>j</sup>* − !*b<sup>s</sup> IFi*,*<sup>j</sup>* .

As mentioned above, it is difficult to achieve the high-rate one-hourly updated PCE because of many estimated parameters for the UD model. To ensure the accuracy, stability, and computational efficiency of one-hourly ultra-rapid PCE, the ED model is applied to perform PCE, which can be expressed as:

$$\begin{cases} \Delta \mathbf{f}\_{r, \mathrm{IF}\_{i\bar{j}}}^{\mathrm{g}}(t\_{k+1}, t\_{k}) = \boldsymbol{\varepsilon} \cdot \Delta \delta \mathbf{f}\_{\mathrm{IF}\_{i\bar{j}}}^{\mathrm{r}}(t\_{k+1}, t\_{k}) - \boldsymbol{\varepsilon} \cdot \Delta \delta \mathbf{f}\_{\mathrm{IF}\_{i\bar{j}}}^{\mathrm{s}}(t\_{k+1}, t\_{k}) + \Delta g\_{\mathrm{w}}(t\_{k+1}, t\_{k}) \cdot \boldsymbol{\pi}\_{\mathrm{w}}(t\_{k+1}) + \boldsymbol{\varepsilon}\_{\Delta \mathrm{IF}\_{\bar{i}, \bar{j}}} \\ \Delta \Phi^{s}\_{r, \mathrm{IF}\_{i\bar{j}}}(t\_{k+1}, t\_{k}) = \boldsymbol{\varepsilon} \cdot \Delta \delta \mathbf{f}\_{\mathrm{IF}\_{i\bar{j}}}^{\mathrm{r}}(t\_{k+1}, t\_{k}) - \boldsymbol{\varepsilon} \cdot \Delta \delta \mathbf{f}\_{\mathrm{IF}\_{i\bar{j}}}^{\mathrm{s}}(t\_{k+1}, t\_{k}) + \Delta g\_{\mathrm{w}}(t\_{k+1}, t\_{k}) \cdot \boldsymbol{\pi}\_{\mathrm{w}}(t\_{k+1}) + \boldsymbol{\varepsilon}\_{\Delta \mathrm{M}\_{\bar{I},i}} \end{cases} \tag{4}$$

where Δ is the ED operator; *s* = *GPS*, *BDS*, *GLONASS*, *Galileo*; for example, the ED clock offsets are Δ*δt r IFi*,*<sup>j</sup>* (*tk*+1, *tk*) = *δt r IFi*,*<sup>j</sup>* (*tk*<sup>+</sup>1) − *δt r IFi*,*<sup>j</sup>* (*tk*). The ambiguity can be removed by ED operation when there are no cycle slips in carrier phase observations. The inter-system bias is quite stable over a day and can be eliminated by ED [20]. Therefore, only ZWDs and the ED clock offsets of receiver and satellite remain. One receiver equipped high-precision atomic clock is selected as reference clock to avoid the singularity between receiver and satellite clock offsets.

Although the traditional ED model and the corresponding data processing strategy are difficult to achieve the high-rate one-hourly updated ultra-rapid PCE, the ED model greatly reduces the number of estimated parameters and reduces the consuming time. Based on the characteristic of traditional ED model and the corresponding data processing strategy, we investigated the influence of tropospheric delays on ED PCE. Generally, the tropospheric delays cannot be ignored in the ED PCE. The systematic bias caused by mapping function and systematic variation of tropospheric delays will be absorbed by absolute clock offsets [7,21]. The ZWD are generally solved once an hour. If ZWD can be eliminated by ED, as shown in Figure 1, the estimated parameters will be further reduced. In addition, there is no correlation between epochs on the premise of estimating ED clock offsets as white noise. The parallel processing between epochs can be realized to further improve the solution efficiency. The modified ED model without ZWD estimations can be expressed as:

$$
\begin{bmatrix}
\Delta P\_{r,IF\_{i\bar{j}}}^{s} \\
\Delta \Phi\_{r,IF\_{\bar{i}\bar{j}}}^{s} \\
 ck\_{ref}
\end{bmatrix} = \begin{bmatrix}
e\_n^T \otimes e\_2 \otimes e\_m & e\_k^T \otimes e\_1 \otimes e\_m \\
e\_1^T \otimes e\_1 \otimes e\_l & 
\end{bmatrix} \begin{bmatrix}
\Delta \delta t\_{IF\_{\bar{i},\bar{j}}}^r \\
\Delta \delta t\_{IF\_{\bar{i},\bar{j}}}^s
\end{bmatrix} + \begin{bmatrix}
\varepsilon\_{\Delta P\_{IF\_{\bar{i},\bar{j}}}} \\
\varepsilon\_{\Delta \Phi\_{IF\_{\bar{i},\bar{j}}}}
\end{bmatrix} \tag{5}
$$

where *n* is the number of stations, *m* represents the number of observations, *k* depicts the number of the satellites, *l* is the number of observations at reference station. *en*, *em*, *ek*, and *el* are *n*-row, *m*-row, *k*-row, and *l*-row vector in which all values are 1; *clkref* is the reference clock. For simplicity, SHA1 and SHA2 denote the modified ED PCE model (ignoring ZWD) and the traditional ED PCE model (estimating ZWD), respectively.

The estimated ED clock offsets must be restored to the absolute clock offsets before it can be used normally. The low-rate clock offsets are used as control points to obtain absolute clock offsets, which can be written as [18]:

$$
\begin{bmatrix}
\tilde{\delta}t\_{fix}(t\_1) \\
\Delta\delta(t\_2, t\_1) \\
\Delta\delta(t\_3, t\_2) \\
\vdots \\
\Delta\delta(t\_i, t\_{i-1}) \\
\tilde{\delta}t\_{fix}(t\_i) \\
\Delta\delta(t\_{i+1}, t\_i) \\
\vdots \\
\Delta\delta(t\_n, t\_{n-1}) \\
\tilde{\delta}t\_{fix}(t\_n)
\end{bmatrix} = \begin{bmatrix}
\delta t(t\_1) \\
& -\delta t(t\_2) & \delta t(t\_3) \\
& & \ddots & \vdots \\
& & & -\delta t(t\_{i-1}) & \delta t(t\_i) \\
& & & & \delta t(t\_i) \\
& & & & -\delta t(t\_i) & \delta t(t\_{i+1}) \\
& & & & & \ddots & \vdots \\
& & & & & & \ddots & \vdots \\
& & & & & & -\delta t(t\_n - 1) & -\delta t(t\_n)
\end{bmatrix} \tag{6}
$$

where *δ* !*tfix*(*t*1), *δ* !*tfix*(*ti*), and *δ* !*tfix*(*tn*) are the low-rate clock offsets; Δ*δt*(*ti*+1, *ti*) is the estimated ED clock offsets; *n* denotes the number of epochs of the low-rate clock offsets.

#### **3. Datasets and Processing Strategies**

#### *3.1. Datasets*

Approximately 120 stations are used to estimate the high-rate one-hourly ultra-rapid clock offsets. The hourly updated observations from IGS MGEX and iGMAS are used to perform the PCE. The selected GNSS stations are shown in Figure 2. To verify the performance of high-rate hourly updated ultra-rapid clock offsets, the GNSS stations without PCE are selected for RT-PPP. The stations marked red are used for high-rate PCE, and the other GNSS stations marked blue are used for RT-PPP.

**Figure 2.** Distribution of the selected GNSS tracking stations for PCE and PPP, in which the stations marked red are used for PCE and the stations marked blue are used for RT-PPP.

#### *3.2. Processing Strategies*

The one-hourly updated ultra-rapid multi-GNSS PCE processing strategies are illustrated in Table 1. The dual-frequency observations of L1 (1575.42 MHz) and L2 (1227.60 MHz) for GPS, G1 (1602+k\*9/16 MHz) and G2 (1246+k\*7/16 MHz) for GLONASS, E1 (1575.42 MHz) and E5a (1176.45 MHz) for Galileo, B1I (1561.098 MHz) and B3I (1268.52 MHz) for BDS-2 and BDS-3 are used to generate high-rate hourly updated ultra-rapid multi-GNSS clock offsets. The satellite orbit and station coordinates are fixed to the SHAO one-hourly updated ultra-rapid solutions [11]. The reference frame and time system of SHAO (Shanghai Astronomical Observatory) precise products are ITRF 2014 and GPS time (GPST), and the satellite orbit and clock offsets is 5 min sampling [11]. As for the noteworthy clock, one receiver equipped high-accuracy atomic clock is selected as the reference clock to avoid the singularity between receiver and satellite ED clock offsets. Because the low-rate clock offsets are used as control points to obtain absolute clock offsets, the reference clock is consistent with SHAO to avoid additional bias. The dry tropospheric delay is corrected using the modified Hopfield model and Global Pressure and Temperature 3 (GPT3), and the Vienna mapping functions 3 (VMF3) is used to obtain the mapping functions of both dry and wet parts [22]. Regarding the wet part, the two schemes mentioned above are proposed to investigate the effect of troposphere estimation on PCE. As for phase center offset (PCO) and phase center variations (PCV), GPS, GLONASS, Galileo satellite antennas are corrected using the antenna file data provided by IGS MGEX [23,24]. BDS-2 and BDS-3 are corrected using the BeiDou official data (http://www.beidou.gov.cn/, accessed on 1 March 2022).


**Table 1.** The data processing strategy for precise clock estimation (PCE).


**Table 1.** *Cont.*

Figure 3 illustrates the slide window for one-hourly updated ultra-rapid orbit and clock offsets products generation. To achieve one-hourly updated ultra-rapid orbit and clock offsets, the computation time of POD and high-rate PCE is within one hour. SHAO can solve low-rate ultra-rapid multi-GNSS orbit and clock offsets in 35 to 40 min based on the Lenovo SR650 server (4\*Inter (R) Xeon (R) Gold 6244 CPU @3.60 GHz), which means that the high-rate PCE needs to be solved in about 20 min. Moreover, the latency of one-hourly ultra-rapid multi-GNSS satellite products is one hour, we need to predict it to meet the real-time applications. The orbit prediction method is integrating the prediction orbit by the fitted orbit, from 24 h, and adjacent afore 24 h orbits (observation session) [11]. The clock offsets prediction model is modified Auto Regressive Integrated Moving Average (ARIMA) [27].

**Figure 3.** The slide window for one-hourly updated ultra-rapid orbit and clock offsets products generation.

According to the characteristics of the modified ED PCE model, we propose a corresponding new data processing flow, as shown in Figure 4, to achieve high-rate one-hourly updated ultra-rapid PCE. First, the IGS real-time stream and one-hourly updated observations obtained from IGS MGEX and iGMAS are merged and preprocessed. After that, the POD is performed to obtain one-hourly updated ultra-rapid precise products including orbit, low-rate clock offsets, Earth rotation parameters (ERP), and station coordinates. Then, the one-hourly updated ultra-rapid precise products and observations are imputed into high-rate PCE system marked red in Figure 4. In the processing of calculating correction information, there is no correlation between stations, the OpenMP can be used for parallel

processing of multiple stations. In the processing of calculating ED clock offsets, the Intel MKL is used to accelerate the matrix solution. Because the SHA1 model does not need to estimate ZWD, the OpenMP can be used for parallel processing of multiple epochs. However, the OpenMP cannot be used for SHA2 parallel processing of multiple epochs. In the process of obtaining absolute clock offsets, the OpenMP can be used for parallel processing of multiple satellites and Intel MKL can be used to accelerate the matrix solution. This high-rate PCE system has strong flexibility, and it can be connected to the ultra-rapid, rapid, final, and real-time POD system to obtain high-rate clock offsets.

**Figure 4.** The data processing flow for one-hourly updated ultra-rapid orbit and clock offsets products generation.

### **4. Validation and Results**

This part begins with the validation of the estimated clock offsets. Because the highrate one-hourly updated ultra-rapid products are mainly oriented to real-time and near-realtime users, we compared the estimated high-rate clock offsets and real-time clock offsets

provided by CNES (CLK93) in terms of the clock offsets accuracy and PPP performances. Finally, the performance of the estimated clock offsets is discussed.

### *4.1. The ED Clock Offsets*

The time series of estimated ED clock offsets and its comparison with GFZ rapid (GBM) and CODE final (COD) multi-GNSS products are depicted in Figure 5. To show ED clock offsets more clearly, the ED clock offsets are shifted by the same amount (1 ns) to avoid overlapping. Figure 6 shows their differences. Because BDS constellations are composed of Geostationary Earth Orbit (GEO), Inclined Geosynchronous Satellite Orbit (IGSO), and Medium Earth Orbit (MEO) satellites, different orbit types of BDS satellites have different accuracy in POD [28,29]. Therefore, the corresponding clock offsets have different accuracy for different orbit types. In order to analyze BDS-2 objectively, we discuss it separately. From Figures 5 and 6, we can find that the difference between GBM and the estimated ED clock offsets for SHA1 and SHA2 is very small, and they are basically distributed on the diagonal with slope of 1. From Figures 5 and 7, we can draw a similar conclusion that the estimated ED clock offsets are basically consistent with COD. This means that the performance of the estimated ED clock is basically consistent with GBM and COD. Certainly, we can also find that the performance of BDS-2 is slightly worse than other systems, especially for GEO satellites. There are three reasons for the difference. First, BDS-2 is a regional navigation system, and the stations tracked BDS-2 signals are basically located in Asia and Europe. There are less observations for BDS-2 POD and PCE because of the unevenly distributed stations. Secondly, BDS-2 contains five GEO and seven IGSO satellites. The orbit accuracy of BDS-2 GEO and IGSO satellites is poor due to the force models and the observation geometry. Therefore, the estimated ED clock offsets calculated by fixed GEO orbit is slightly poor. Last, but most important, there are some differences between GFZ and SHAO in BDS data processing strategies in terms of the precise attitude model, the solar radiation pressure model, and DCB corrections [11,30], which leads to the poor consistency of the results. In order to more accurately and objectively explain the accuracy of the estimated ED clock offsets compared with GBM and COD, the mean, STD, and RMS of the estimated ED clock offsets between GBM, COD and SHA1, SHA2 for multi-GNSS satellites are listed in Table 2. Our calculation is basically consistent with the results of GBM and COD. The difference in ED clock offsets between GBM, COD and SHA1, SHA2 for GPS and Galileo satellites range from 6 to 7 ps, and the difference for BDS and Galileo range from 10 to 16 ps. The results in Table 2 confirm the findings of Figures 5–7.


**Table 2.** The statistics of means, STD, and RMS for GPS, BDS-2 GEO, BDS-2 IGSO, BDS-2 MEO, BDS-3, GLONASS, and Galileo satellites ED clock (unit: ps).

**Figure 5.** The time series of estimated ED clock offsets and its comparison with GFZ rapid and COD final clock products, in which CC denotes the correlation coefficient.

Now, turn to the influence of tropospheric delays on ED clock offsets estimation. Figures 5–7 not only show the difference between GBM, COD, and the estimated ED clock offsets but also depict the difference between SHA1 and SHA2. As shown in Figure 5, SHA1 and SHA2 are highly overlapped, and it is difficult to distinguish the difference. From Figures 6 and 7 and Table 2, we can see that the difference between SHA1 and SHA2 is very small, ranging from 0 to 0.35 ps in terms of the statistics of mean, STD, and RMS for GPS, BDS-2 GEO, BDS-2 IGSO, BDS-2 MEO, BDS-3, GLONASS, and Galileo satellites, which can be ignored. To prove this point more favorably and reflect the situation of each satellite clearly, the relationship between SHA1 and SHA2 in terms of the mean, STD, and RMS of each satellite is depicted in Figure 8. The mean of each satellite for SHA2 is closer

to 0 than SHA1, but the difference is very tiny. In addition, the STD and RMS of each satellite are basically distributed on the diagonal with slope of 1. In general, the ED clock of SHA1 and SHA2 are basically consistent both in mean, STD, and RMS. However, from a calculation time-consuming point of view, adding the ZWD parameter will increase the computation complexity. The corresponding time-consuming will increase. As shown in Figure 9, the mean time-consuming of SHA1 is about 16.5 min, and that of SHA2 is about 23.5 min, based on Lenovo SR650 server (4\*Inter (R) Xeon (R) Gold 6244 CPU @3.60GHz) using OpenMP and Intel MKL. In case of poor data quality, the time-consuming of SHA1 is about 20 to 23 min. In normal conditions, SHA1 takes less than 20 min.

**Figure 6.** The difference between GFZ rapid clock offsets and estimated ED clock offsets, in which (a, b, c) denotes (mean, STD, and RMS).

**Figure 7.** The difference between CODE final clock offsets and estimated ED clock offsets, in which (a, b, c) denotes (mean, STD, and RMS).

As mentioned above, the time-consuming of POD and high-rate PCE must be less than 1 h to achieve one-hourly updated orbits and clocks. Because the difference between SHA1 and SHA2 is very tiny, and its difference is less than 0.35 ps in terms of mean, STD, and RMS. As stated earlier, SHA1 can use the OpenMP for parallel processing of multiple stations, multiple epochs, and multiple satellites, while SHA2 can only use the OpenMP for parallel processing of multiple stations and multiple satellites. The computational efficiency of SHA1 is significantly higher than that of SHA2. Considering the calculation time and accuracy, SHA1 as the optimal strategy is applied to obtain one-hourly updated ultra-rapid clock.

**Figure 8.** The mean, STD, and RMS of each satellite for GPS, BDS-2, BDS-3, GLONASS, and Galileo using two tropospheric processing strategies.

#### *4.2. The Absolute Clock Offsets*

The ED clock offsets calculated by the SHA1 strategy are restored to the absolute clock offsets by using Equation (6). The double-difference (DD) method of selecting one satellite to eliminate the clock datum [31] is used to assess the absolute clock. GPS PRN01 (G01), BDS-2 PRN07 (C07), BDS-3 PRN19 (C19), GLONASS PRN24 (R24), and Galileo PRN01 (E01) are selected as the reference satellite to assess the clock of other GNSS satellites. Figures 10 and 11 show the STD and RMS of high-rate and low-rate one-hourly updated GNSS ultra-rapid clock for observation. In addition, the ultra-rapid products are more targeted to real-time and near real-time users; we added the real-time stream products CLK93 released by CNES to discuss the precision of estimated clock and CLK93. When users use the low-rate clock, they need to use mathematical interpolation algorithm to interpolate the low-rate clock, which causes the loss of clock accuracy [14]. In order to reflect the actual applications of low-rate clock, the clock interpolated by the ninth order Lagrange interpolation algorithm was also evaluated.

**Figure 9.** The computation time for 1-day high-rate multi-GNSS PCE using 120 GNSS tracking stations, in which Non, OMP, and MKL denote no acceleration method, OpenMP, and Intel MKL, respectively.

It can be seen from Figure 10 that the STD of the estimated high-rate clock (SHA30s (PCE)), low-rate clock (SHA300s), and clock offsets interpolated by mathematical interpolation algorithm (SHA30s (Math)) for GPS, BDS-2 GEO, BDS-2 non-GEO, BDS-3, GLONASS, and Galileo satellites are basically better than 0.2, 1.2, 0.3, 0.3, 0.4, and 0.2 ns, respectively. From SHA30s (PCE), SHA300s, and SHA30s (Math), we can find that SHA30s (PCE) has a certain improvement compared with SHA30s (Math), but the improvement is not obvious compared with SHA300s. There are three reasons for the little increase in SHA30s (PCE) compared with SHA300s. Firstly, the weight of the control points is higher than that of the ED clock offsets when the ED clock offsets are restored to the absolute clock offsets by using low-rate clock offsets as control points. Secondly, the accuracy of SHA300s is obtained by evaluating the clock offsets on time node (300s intervals), and there is no precision loss caused by the interpolation algorithm. Finally, the accuracy of SHA300s provide by SHAO is already good.

It is worth mentioning that the STD of the estimated clock offsets is slightly better than CLK93. There are two possible explanations: the clock offsets precision is inevitably disturbed by the stability of real-time observation streams, more commonly, the available observation data for real-time product solution are less than ultra-rapid solutions. In other words, one-hourly updated ultra-rapid orbit and clock offsets products can ensure better accuracy in the case of poor network environment. Whether it is used as backup data or practical data, one-hourly updated ultra-rapid orbit and clock offsets products are very reliable.

**Figure 10.** The STD of multi-GNSS ultra-rapid estimated clock for observation session, in which SHA30s (PCE), SHA300s, SHA30s (Math), and CLK93 depict the estimated high-rate clock, lowrate clock, clock interpolated by mathematical interpolation algorithm, and CNES real-time clock, respectively. The Ref.SATE means the reference satellite.

Note that RMS, the RMS of SHA30s (PCE), SHA300s, and SHA30s (Math) for GPS, BDS-2, BDS-3, GLONASS, and Galileo satellites are basically better than 0.5, 3.0, 0.9, 1.2, and 0.4 ns. The RMS of satellites such as G04, G14, G18, G23, C11, C12, and C14 are especially large compared with other satellites. Because the different ACs adopt different processing strategies such as the precise attitude model, the solar radiation pressure model, and DCB corrections [11,30,32], these inconsistencies lead to the differences between SHAO and GBM on some satellites. Moreover, the RMS of GLONASS satellites is worse than other GNSSs, which is due to the possible bias caused by the float ambiguity solutions in POD and frequency division multiple access (FDMA) signal transmit mechanism. From the comparison between SHA30s (PCE) and CLK93, we can find that the RMS of SHA30s (PCE) is better than that of CLK93. As commented earlier, the main reasons are as follows: on the one hand, different ACs have differences in data processing strategies, on the other hand, the available observation data for real-time product solution is not as stable as that of ultra-rapid solutions.

**Figure 11.** The RMS of high-rate and low-rate one-hourly updated GNSS ultra-rapid clock for observation session, in which SHA30s (PCE), SHA300s, SHA30s (Math), and CLK93 depict the estimated high-rate clock, low-rate clock, clock interpolated by mathematical interpolation algorithm, and CNES real-time clock, respectively.

Now, pay attention to the mean STD and RMS for the observation session and prediction session. Combined with Figures 10–12, we can obtain four key findings. In terms of the precision of the estimated high-rate clock offsets and SHAO low-rate clock offsets, GPS and Galileo are the best, BDS-3 is the second, and GLONASS and BDS-2 are relatively worse. Secondly, SHA30s (PCE) is optimal in both the observation and prediction session. The absolute clock of the prediction and observation sessions has been modified by the improvement of the time resolution. The STD and RMS are improved by (0.97 to 9.09% and 0.12 to 5.56%) in the observation session, (2.82 to 23.08% and 0.95 to 9.09%) in the first hour of the prediction session, and (0 to 3.85% and 0.12 to 4.19%) in the second hour of the prediction session, respectively. Thirdly, the precision of SHA30 (PCE) is slightly better than CLK93 in the observation section and basically the same as CLK93 in the prediction

section. Finally, there are some differences between SHAO, GFZ, and CNES in the data processing strategy, and the clock offsets of different ACs show differences in RMS.

**Figure 12.** Average STD and RMS of high-rate and low-rate one-hourly updated GNSS ultra-rapid clock for observation session and prediction session.

#### *4.3. Real-Time PPP Validation*

The RT-PPP experiments will be carried out to verify the practicability and reliability of the estimated clocks. The PPP solutions using GBM products (30 s interval) is taken as the reference. Figure 13 depicts the PPP convergence at HUEG for GPS-only, BDSonly, GLONASS-only, Galileo-only, and multi-GNSS (GPS + BDS + GLONASS + Galileo). To show the PPP convergence clearly, only the first 2 h are shown in Figure 13. From Figure 13, we can see that the convergence time for SHA30s (PCE) is obviously improved, and its performance is basically same as that of GBM products. Figure 14 shows the average convergence time. During the observation session, SHA30s (PCE) can improve the convergence time by (48.16, 31.55, 37.06%) for GPS-only PPP, (32.88, 21.96, 22.60%) for BDS-only PPP, (40.78, 31.15, 35.05%) for GLONASS-only PPP, (40.49, 23.49, 25.57%) for Galileo-only PPP, and (64.75, 60.28, 55.07%) for multi-GNSS PPP, compared with SHA300s in north, east, and up directions, respectively. In the first hour of the prediction session, SHA30s (PCE) can improve the convergence time by (49.84, 32.47, 37.63%) for GPS-only PPP, (33.62, 24.53, 28.12%) for BDS-only PPP, (44.16, 27.69, 28.69%) for GLONASS-only PPP, (44.82, 25.03, 23.95%) for Galileo-only PPP, and (65.75, 60.98, 56.21%) for multi-GNSS PPP, respectively. In the second hour of the prediction session, SHA30s (PCE) can improve the convergence time by (57.27, 34.85, 36.03%) for GPS-only PPP, (39.76, 27.80, 31.31%) for BDS-only PPP, (43.99, 28.29, 26.30%) for GLONASS-only PPP, (45.74, 26.30, 28.29%) for Galileo-only PPP, and (62.12, 55.51, 57.14%) for multi-GNSS PPP, respectively. The estimated high-rate one-hourly ultra-rapid clock greatly shortens the PPP convergence

time, which is basically the same as GBM and slightly better than CLK93 in the observation session. In the prediction session, SHA30s (PCE) performs good, and the convergence time of SHA30s (PCE) is at the same level as CLK93.

**Figure 13.** PPP convergence at HUEG for GPS-only, BDS-only, GLONASS-only, Galileo-only, and multi-GNSS (GPS + BDS + GLONASS + Galileo) from 0:00 to 2:00 on DOY 121, 2021.

As for positioning accuracy, the final convergence accuracy of SHA30s (PCE) and SHA300s is basically the same in the observation section. However, the positioning accuracy of SHA300s decreases slightly in the prediction session, while the performance of SHA30s (PCE) is still good. From the Figures 15 and 16, it is obvious that the positioning accuracy will decrease with the increase in prediction time. However, using high-rate clock offsets can appropriately avoid the decrease in positioning accuracy. The positioning accuracy of SHA30s (PCE) for GPS-only, BDS-only, GLONASS-only, Galileo-only, and multi-GNSS PPP is improved by 14.90, 21.68, 13.38, 9.06, and 25.74% in the observation session, 12.34, 1.52, 13.37, 8.66, and 12.99% in the first hour of the prediction session, and 7.97, 9.86, 7.76, 22.45, and 11.65% in the second hour of the prediction session, respectively. Similarly, SHA30s (PCE) improves the PPP positioning accuracy, which is slightly better than CLK93 in the observation session. In the prediction session, the positioning accuracy of SHA30s (PCE) is at the same level as CLK93. Even in the case of poor network environment, the estimated high-rate one-hourly ultra-rapid clock can ensure the better PPP performance.

**Figure 14.** The average convergence time for GPS-only, BDS-only, GLONASS-only, Galileo-only, and multi-GNSS PPP from DOY 121 to 130, 2021.

**Figure 15.** PPP positioning error at HUEG for GPS-only, BDS-only, GLONASS-only, Galileo-only, and multi-GNSS from 20:00 on DOY 121 to 2:00 on DOY 122, 2021.

**Figure 16.** The average RMS for GPS-only, BDS-only, GLONASS-only, Galileo-only, and multi-GNSS PPP from DOY 121 to 130, 2021.

In summary, the estimated high-rate one-hourly ultra-rapid clock can improve the PPP performance in terms of positioning accuracy and convergence time and can be well applied to real-time PPP applications.

#### **5. Conclusions**

This contribution focused on the high-rate one-hourly updated ultra-rapid PCE and its application in RT-PPP. We modified the ED PCE model by investigating the influence of tropospheric delays on ED PCE and proposed a new framework to obtain high-rate one-hourly updated ultra-rapid clock. The OpenMP and Intel MKL are used to construct a parallel processing system to improve computation efficiency. Through the analysis and discussion, we can draw the following conclusions.

As for the influence of tropospheric delays on ED PCE, the difference between SHA1 and SHA2 is very tiny, ranging from 0 to 0.35 ps in terms of the statistics of mean, STD, and RMS for multi-GNSS satellites, which can be ignored. From a calculation time-consuming point of view, adding the ZWD parameter will increase the computation complexity. The corresponding time-consuming will increase. The time-consuming of SHA1 is about 16.5 min, and that of SHA2 is about 23.5 min based on the Lenovo SR650 server (4\*Inter (R) Xeon (R) Gold 6244 CPU @3.60GHz). Therefore, SHA1 as the optimal strategy with low time-consuming and high accuracy is recommended to obtain one-hourly updated ultra-rapid clock offsets. From the precision of the clock offsets point of view, improving the time resolution can improve the precision of the absolute clock offsets in both prediction and observation sessions. The STD and RMS are improved by (0.97 to 9.09% and 0.12 to 5.56%) in the observation session, (2.82 to 23.08% and 0.95 to 9.09%) in the first hour of the prediction session, and (0.11 to 3.85% and 0.12 to 4.19%) in the second hour of the prediction session, respectively.

Similarly, the estimated clock offsets perform well in RT-PPP. The estimated clock offsets improve the PPP positioning accuracy and greatly shorten the convergence time. The positioning accuracy can be improved by 9.06~25.74% in the observation session, 1.52~13.37% in the first hour of the prediction session, and 7.76~22.45% in the second hour

of the prediction session, respectively. The convergence time can be significantly improved by 21.96~64.75, 23.95~65.75, and 26.30~62.12% for the observation session, the first hour of the prediction session, and the second hour of the prediction session compared with low-rate products, respectively. The RT-PPP performance of SHA30s (PCE) is basically the same as GFZ rapid products and slightly better than CLK93. Whether in the prediction session or in the observation session, the estimated clock shows good performance. The estimated high-rate one-hourly ultra-rapid clock can ensure the better PPP performance in the case of poor network environment.

Synthesizing the analysis and discussion above, the estimated high-rate one-hourly updated ultra-rapid precise clock offsets (URL: ftp://igsdepot.ign.fr/pub/igs/products/ mgex/, accessed on 1 March 2022) have excellent performance, which can be well applied to real-time or near real-time applications and research.

**Author Contributions:** Conceptualization, G.J. and S.S.; methodology, G.J.; software, G.J.; validation, G.J. and S.S.; formal analysis, G.J.; investigation, G.J. and S.S.; resources, S.S.; data curation, G.J. and S.S.; writing—original draft preparation, G.J.; writing—review and editing, G.J. and S.S.; visualization, G.J.; supervision, S.S.; project administration, S.S.; funding acquisition, S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China grant number No. 41730109 and No. 12073063.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The datasets analyzed during this study are available in the IGS MGEX repository (ftp://cddis.gsfc.nasa.gov/pub, accessed on 1 March 2022). The multi-GNSS precise orbit and clock products are available in the GFZ (ftp://ftp.gfz-potsdam.de/GNSS/products/mgex/, accessed on 1 March 2022), CODE (ftp://igs.ign.fr/pub/igs/, accessed on 1 March 2022), and CNES (http://www.ppp-wizard.net/products/REAL\_TIME/, accessed on 1 March 2022) ACs.

**Acknowledgments:** The authors would like to thank Ke Su for his valuable suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

