*2.2. Suitable Observation Redundancy for Satellite Clock Batch Estimation*

In addition to the ICB and TDB, the observation redundancy correlated with the number of stations also affects the satellite clock batch estimation. The accuracy of the batchestimated satellite clock becomes saturable when the observation redundancy adequately increases. The benefit of increased observation redundancy is negligible for the accuracy improvement of the satellite clock batch estimation [20]. Furthermore, excessive observation redundancy will reduce the efficiency, which needs to be considered in the satellite clock batch estimation with low latency. Therefore, based on the clarified batch observation model, we derive the analytic relationship between the observation and the parameter to be estimated, and simultaneously consider the solution efficiency to determine the saturable accuracy of the batch-estimated satellite clock.

In order to reduce the computation burden of estimating all parameters in the normal equation with large dimensions, the parameter pre-elimination and back-substitution are generally adopted in the batch estimation [39]. The parameters to be estimated are divided into the time-variant and time-invariant parameters to perform the parameter pre-elimination and back-substitution [2]. The satellite clock and receiver clock are denoted as the time-variant parameter and estimated as epoch-wise. The zenith tropospheric delay is denoted as the time-invariant parameter and estimated as piece-wise. The ambiguity is generally regarded as the constant in the absence of cycle slips, which is also denoted as the time-invariant parameter. For the processing strategy of the satellite clock batch estimation, the time-variant parameter is generally set up and pre-eliminated epoch-wise due to its potentially large number. The back-substitution is used to obtain the solution of the time-variant parameter. The observation model can be rewritten as

$$
\begin{bmatrix} \mathbf{v}\_{P,IF,r}^{s} \\ \mathbf{v}\_{L,IF,r}^{s} \end{bmatrix} = \begin{bmatrix} \mathbf{G}\_{r} & \mathbf{G}^{s} \\ \mathbf{G}\_{r} & \mathbf{G}^{s} \end{bmatrix} \begin{bmatrix} \delta \mathbf{f}\_{r} \\ \delta \mathbf{f}^{s} \end{bmatrix} + \begin{bmatrix} \mathbf{G}\_{ZTD} & \mathbf{0} \\ \mathbf{G}\_{ZTD} & \mathbf{G}\_{N} \end{bmatrix} \begin{bmatrix} \mathbf{Z} \mathbf{T} \mathbf{D}\_{r} \\ \overline{\mathbf{N}}\_{r}^{s} \end{bmatrix} - \begin{bmatrix} \mathbf{I}\_{P,IF,r}^{s} \\ \mathbf{I}\_{L,IF,r}^{s} \end{bmatrix} \tag{4}
$$

Combining with the pseudorange and phase observations, the observation error is accurately characterized by the corresponding weighting matrix. In order to simply represent the processing of the pre-eliminated parameter and back-substitution, the normal equation of the satellite clock batch estimation can be abbreviated in the matrix form as

$$PV = PAX + PBY - PL\tag{5}$$

where *P* is the weighting matrix related to the observation error, *X* and *Y* are the vectors of the time-variant and time-invariant parameters, *A* and *B* are their corresponding coefficient matrices, *V* and *L* are the vectors of the a posteriori and the a priori observation residuals. According to the principle of the least square parameter estimation based on the parameter pre-elimination and back-substitution, the batch-estimated parameters can be expressed as

$$\begin{cases} \dot{\mathbf{X}} = \mathbf{N}\_{XX}^{-1} \left( \mathbf{U}\_{XL}L - \mathbf{N}\_{XY}\hat{\mathbf{Y}} \right) \\ \dot{\mathbf{Y}} = \left( \mathbf{N}\_{YY} - \mathbf{N}\_{YX}\mathbf{N}\_{XX}^{-1}\mathbf{N}\_{XY} \right)^{-1} \left( \mathbf{U}\_{YL}L - \mathbf{N}\_{YX}\mathbf{N}\_{XX}^{-1}\mathbf{U}\_{XL}L \right) \end{cases} \tag{6}$$

where *<sup>N</sup>XX* <sup>=</sup> *<sup>A</sup>*T*CA*, *<sup>N</sup>XY* <sup>=</sup> *<sup>A</sup>*T*CB*, *<sup>N</sup>YX* <sup>=</sup> *<sup>B</sup>*T*CA*, *<sup>N</sup>YY* <sup>=</sup> *<sup>B</sup>*T*CB*, *<sup>U</sup>XL* <sup>=</sup> *<sup>A</sup>*T*C*, *<sup>U</sup>YL* <sup>=</sup> *<sup>B</sup>*T*C*, *C* <sup>=</sup> *P*T*P*. The diagonal element in the covariance–variance matrix can reflect the STD of the parameter to be estimated. The satellite clock batch estimation accuracy can be obtained by extracting the diagonal elements of the covariance–variance matrix corresponding to the satellite clocks. This evaluated accuracy does not require the additional satellite clock reference product. The clock covariance–variance matrix of the satellite clock batch estimation derived by (6) can be expressed as

$$\text{COV}(\varepsilon\_{\hat{X}'}\varepsilon\_{\hat{X}}) = \mathbf{N}\_{XX}^{-1} + \mathbf{N}\_{XX}^{-1}\mathbf{N}\_{XY}(\mathbf{N}\_{YY} - \mathbf{N}\_{YX}\mathbf{N}\_{XX}^{-1}\mathbf{N}\_{XY})^{-1}\mathbf{N}\_{YX}\mathbf{N}\_{XX}^{-1} \tag{7}$$

It can be found that the clock covariance–variance matrix consists of the coefficient matrix and the weighting matrix for the satellite clock batch estimation model. The satellite clock batch estimation accuracy depends on the number of stations and the observation error, which are closely related to the coefficient matrix dimension and the weighting matrix, respectively. Thus, the analytic relationship between the observation and the parameter to be estimated has clearly been derived. However, the number of stations suitable for the satellite clock batch estimation and the corresponding saturable estimation accuracy still cannot be determined directly from the analytic expression. This is because the coefficient matrix dimension is also affected by the number of visible satellites at different stations and different epochs, which needs to be determined experimentally.

The GPS data from the day of the year (DOY) 045 to 051 in 2021 are used for determining the number of stations suitable for the satellite clock batch estimation and the corresponding saturable estimation accuracy. We ignore the impacts of different stations with identical quantities on the satellite clock batch estimation. This is due to the identical full coverage for visible satellites in the observation session and the little difference in the observation quality between different International GNSS Service (IGS) stations used. The efficiency of the satellite clock batch estimation has limitations due to the low latency requirements, which should be controlled within one hour including the consumption of hourly data collection. We count the processing time to reflect the efficiency of satellite clock batch estimation [40]. The satellite clock batch estimation based on the final satellite orbit product provided by IGS is performed on the PowerEdge R7525 Server (AMD EPYC 7F72 24-Core Processor @3.69 GHz). The STD is extracted from (7) based on the gradually increasing number of the global distributed network stations and the elevation-dependent stochastic model, whose consistency with the STD for the external coincidence is verified. The STD for the external coincidence is conducted by calculating the batch-estimated satellite clock error based on the double difference strategy [10]. Moreover, the criteria for gradually increasing stations serve to maximize the observation, which is mainly relative to the observed satellites. The mean values of the STD for all satellites and the processing time of the satellite clock batch estimation with a gradually increasing number of stations are shown in Figure 1.

**Figure 1.** STD and processing time with a gradually increasing number of stations for one week.

It can be seen from Figure 1 that the satellite clock batch estimation accuracy is gradually improved as the increasing number of stations. The satellite clock batch estimation accuracy is greatly affected by the observation redundancy for less than 40 stations and ranges from 0.03 to 0.10 ns. However, the GPS satellite clock batch estimation accuracy tends to be stable when the number of stations is more than 40. This indicates that the satellite clock batch estimation accuracy becomes saturable. Considering the counterbalance between the efficiency and estimation accuracy, we suggest that the number of stations suitable for GPS hourly updated satellite clock batch estimation is to be 40.

#### *2.3. Impacts of Batch-Estimated Satellite Clock Biases on PPP*

Based on the suitable observation redundancy for the satellite clock batch estimation, the ICB and TDB of the batch-estimated satellite clock affect the PPP performance. We will clarify the error propagation of these satellite clock biases in the PPP observation model. For the conventional dual-frequency IF combined PPP processing, the consistent satellite clock and orbit products are generally fixed [41]. Since the satellite clock batch estimation is performed based on the ultra-rapid satellite orbit product, the error exists in the fixed satellite clock product for the PPP observation model. Due to the consistency between the satellite clock and orbit product used for PPP, the fixed ultra-rapid satellite orbit product is also fixed in PPP. Thus, the satellite clock and orbit errors cannot be ignored in the PPP observation model. For a typical PPP observation model [42,43], the observations and the parameters are epoch-wise constructed, which can be expressed as

$$
\begin{bmatrix} \mathbf{v}\_{P,IF,r}^{s} \\ \mathbf{v}\_{L,IF,r}^{s} \end{bmatrix} = \begin{bmatrix} \mathbf{K}\_{r} & \mathbf{K}\_{r} & \mathbf{K}\_{\mathcal{Z}\text{TD}} & \mathbf{K}\_{r} & \mathbf{K}^{s} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{K}\_{r} & \mathbf{K}\_{r} & \mathbf{K}\_{\mathcal{Z}\text{TD}} & \mathbf{K}\_{r} & \mathbf{K}^{s} & \mathbf{K}^{s} & \mathbf{K}\_{\mathcal{N}} \end{bmatrix} \begin{bmatrix} \mu\_{r}^{s}\delta\mathbf{x}\_{I} \\ \delta t\_{r} \\ \mathbf{Z}\text{TD}\_{r} \\ \mathbf{b}\_{r} \\ \mathbf{b}^{s} \\ \mathbf{B}\_{r} \\ \mathbf{B}^{s} \\ \mathbf{N}\_{r}^{s} \end{bmatrix} - \begin{bmatrix} \mathbf{f}\_{P,I,\mathcal{Z}}^{s} + \mathbf{g}\_{P,\mathcal{Z}}^{s}\delta\mathbf{x}^{s} - \mathbf{K}^{s}\cdot\delta\mathbf{\tilde{F}} \\ \mathbf{f}\_{P,I,\mathcal{Z}}^{s} + \mathbf{K}^{s}\cdot\mu\_{r}^{s}\delta\mathbf{x}^{s} - \mathbf{K}^{s}\cdot\delta\mathbf{\tilde{F}} \\ \mathbf{f}\_{P,\mathcal{Z}}^{s} + \mathbf{K}^{s}\cdot\mu\_{r}^{s}\delta\mathbf{x}^{s} - \mathbf{K}^{s}\cdot\delta\mathbf{\tilde{F}} \\ \mathbf{N}\_{r}^{s} \end{bmatrix} \tag{8}
$$

where *<sup>K</sup><sup>r</sup>* <sup>=</sup> *<sup>I</sup><sup>p</sup>* <sup>⊗</sup> *<sup>e</sup>q*, *<sup>K</sup><sup>s</sup>* <sup>=</sup> <sup>−</sup>*e<sup>p</sup>* <sup>⊗</sup> *<sup>I</sup>q*, *<sup>K</sup>ZTD* <sup>=</sup> *<sup>I</sup><sup>p</sup>* <sup>⊗</sup> *<sup>M</sup>q*, *<sup>K</sup><sup>N</sup>* <sup>=</sup> *<sup>I</sup><sup>p</sup>* <sup>⊗</sup> *<sup>e</sup>q*, *<sup>K</sup>* is the designed matrix for the PPP observation model, *p* is the number of the observations, *q* is the number of observed satellites at different stations and *δxr* is the increment for the a priori receiver position vector. The receiver position increment, the receiver clock, the zenith tropospheric delay and the ambiguity parameters are the parameters to be estimated in the PPP observation model. Given that the hardware delay is time-invariant, the redefined PPP observation model due to the correlation among the unknown parameters can be expressed as

$$
\begin{bmatrix} \boldsymbol{\sigma}\_{P,IF,r}^{s} \\ \boldsymbol{\sigma}\_{L,IF,r}^{s} \end{bmatrix} = \begin{bmatrix} \mathbf{K}\_{r} & \mathbf{K}\_{r} & \mathbf{K}\_{ZTD} & \mathbf{0} \\ \mathbf{K}\_{r} & \mathbf{K}\_{r} & \mathbf{K}\_{ZTD} & \mathbf{K}\_{N} \end{bmatrix} \begin{bmatrix} \boldsymbol{\mu}\_{r}^{s}\boldsymbol{\delta}\mathbf{x}\_{r} \\ \boldsymbol{\delta}\tilde{t}\_{r} \\ \mathbf{Z}\mathbf{T}\mathbf{D}\_{r} \\ \tilde{\mathbf{N}}\_{r}^{s} \end{bmatrix} - \begin{bmatrix} \boldsymbol{\tilde{I}}\_{P,IF,r}^{s} \\ \boldsymbol{\tilde{I}}\_{L,IF,r}^{s} \end{bmatrix} \tag{9}
$$

where *<sup>δ</sup>tr*, *<sup>N</sup><sup>s</sup> r*, *l s <sup>P</sup>*,*IF*,*<sup>r</sup>* and *<sup>l</sup> s <sup>L</sup>*,*IF*,*<sup>r</sup>* are the reparametrized receiver clock, ambiguity, a priori pseudorange observation residual and a priori phase observation residual, respectively. They can be rewritten as

$$\begin{cases} \begin{aligned} \delta \overline{\mathbf{f}}\_{\mathcal{T}} &= \delta \mathbf{f}\_{\mathcal{T}} + \mathbf{b}\_{\mathcal{T}} + \mathbf{K}^{s} \cdot \mathbf{D}^{s} \\ \overline{\mathbf{N}}\_{\mathcal{T}}^{s} &= \mathbf{N}\_{\mathcal{T}}^{s} + \mathbf{B}\_{\mathcal{T}} - \mathbf{b}\_{\mathcal{T}} + \left[ \delta \overline{\mathbf{f}}^{s}(\mathbf{i}\_{0}) - \mathbf{K}^{s} \cdot \mathbf{b}^{s} \right] \\ \end{aligned} \\ \begin{aligned} \overline{\mathbf{I}}\_{\mathcal{P}, \mathrm{IF}, \mathcal{T}}^{s} &= \mathbf{I}\_{\mathcal{P}, \mathrm{IF}, \mathcal{T}}^{s} + \mathbf{K}^{s} \cdot \mu\_{r}^{s} \delta \mathbf{x}^{s, \mathrm{R}} + \mathbf{K}^{s} \cdot \mu\_{r}^{s} \delta \mathbf{x}^{s, \mathrm{T}} - \mathbf{K}^{s} \cdot \sum\_{j=i\_{0}+1}^{i} \Delta \delta \overline{\mathbf{f}}^{s}(\mathbf{i}) \\ \end{aligned} \\ \begin{aligned} \overline{\mathbf{I}}\_{\mathrm{L}, \mathrm{IF}, \mathcal{T}}^{s} &= \mathbf{I}\_{\mathrm{L}, \mathrm{IF}, \mathcal{T}}^{s} + \mathbf{K}^{s} \cdot \mu\_{r}^{s} \delta \mathbf{x}^{s, \mathrm{R}} + \mathbf{K}^{s} \cdot \mu\_{r}^{s} \delta \mathbf{x}^{s, \mathrm{T}} - \mathbf{K}^{s} \cdot \sum\_{j=i\_{0}+1}^{i} \Delta \delta \overline{\mathbf{f}}^{s}(\mathbf{i}) \end{aligned} \end{cases} \end{cases}$$

where *<sup>δ</sup>t s* (*i*0) is the ICB of the fixed batch-estimated satellite clock at the initial epoch *i*0, Δ*δt s* (*i*) is the TDB of the fixed batch-estimated satellite clock at the *i*th epoch.

We focus on the impacts of the satellite clock bias for the fixed batch-estimated satellite clock on the PPP observation model. It can be found from (10) that the timescale difference can be absorbed by the receiver clock in the PPP observation model without affecting the PPP positioning accuracy. The ICB caused by the code hardware delay for satellites can be canceled by the code hardware delay for satellites in the PPP observation model if the identical combination of observation is used for PPP and satellite clock batch estimation, i.e., the dual-frequency IF combination. The TDB induced by the absorbed satellite orbit error can mostly be canceled by the fixed satellite orbit error in the PPP observation model. The reason for the cancellation is that the satellite clock and orbit product used for PPP is consistent. Furthermore, the cancellation depends on the assimilation ability of the batchestimated satellite clock to the satellite orbit error, resulting in incomplete cancellation.

It can be found from (10) that the ICB not caused by the code hardware delay for satellites is absorbed by ambiguity in the PPP observation model. Furthermore, the TDB caused by the unabsorbed satellite orbit error or others is absorbed by the observation residual in the PPP observation model. For the ICB and TDB of the batch-estimated satellite clock caused by these reasons, we investigate their impacts on the PPP positioning performance, including the convergence time and the positioning accuracy after convergence. The convergence time and positioning accuracy after convergence depend on the speed of accurately estimating parameters and the observation noise, respectively. Since the ICB is absorbed by the ambiguity, the estimation of the ambiguity requires a long convergence process and mainly impacts the convergence time of PPP. Furthermore, the inaccurate ambiguity parameter may impact the phase observation, resulting in affecting the positioning accuracy after convergence. The TDB can be treated as impacting the observation noise when estimating the receiver coordinates. The large TDB is equivalent to increasing the observation noise, and mainly impacts the accuracy of the receiver coordinate estimation, i.e., the positioning accuracy after convergence. Meanwhile, the speed of accurately estimating the receiver coordinate, i.e., the convergence time of PPP, may also be affected.

## **3. Results**

The satellite clock batch estimation based on the global distributed network is implemented to verify the impacts of the ICB and TDB. Since the long absence of the satellite clock reference product and the fixed satellite orbit product for G11, the satellite clock batch estimation and the evaluation does not include this satellite. The experimental setup and

processing strategy is shown in Table 1. The station selection of the global distributed network is performed in the preprocessing of satellite clock batch estimation, which depends on the quality of observations at different stations. One example of the station selection is shown in Figure 2. The number of selected stations is 40, which is the suitable observation redundancy for the satellite clock batch estimation verified in the previous experiments. In order to evaluate the batch-estimated satellite clock and the fixed satellite orbit, the GPS final satellite clock and orbit product provided by IGS are introduced as the reference, respectively.


**Table 1.** Experimental setup and processing strategy.

The TDB due to the satellite orbit error assimilation on satellite clock batch estimation IS clarified by the correlation analysis between the satellite clock error and the satellite orbit error. The batch-estimated satellite clock errors are evaluated to verify the impacts of the ICB and TDB on the satellite clock batch estimation. Furthermore, the convergence time and the positioning accuracy after convergence are analyzed to reveal the impacts of the ICB and TDB on PPP.

**Figure 2.** One example of station selection with 40 stations for satellite clock batch estimation. The purple circles represent selected stations.

#### *3.1. Effect Analysis of Satellite Clock Batch Estimation*

The bias and dispersion of the batch-estimated satellite clock error can reveal the ICB and TDB, respectively. We use the double difference strategy to extract the batch-estimated satellite clock error. The reference satellite needs be selected in advance. The reference satellite clock strategies typically include selecting one satellite clock and the mean value of satellite clocks for all satellites. In order to preserve the satellite clock error sequences for all satellites, the mean value of all satellite clocks is selected as the reference in this contribution, which can be expressed as

$$\nabla \Delta \delta \overline{\mathbf{f}}^{\mathbf{s}} = \left[ \delta \overline{\mathbf{f}}^{\mathbf{s}} - \frac{1}{n} \sum\_{i=1}^{n} \delta \overline{\mathbf{f}}^{\mathbf{s}}(i) \right] - \left[ \delta \mathbf{f}\_{ref}^{\mathbf{s}} - \frac{1}{n} \sum\_{i=1}^{n} \delta \mathbf{f}\_{ref}^{\mathbf{s}}(i) \right] \tag{11}$$

where ∇Δ is the double difference operator, *ref* is the satellite clock reference product provided by IGS, *n* is the total number of epochs and *i* is the epoch. In addition, we used the metrics of the RMS and STD to reflect the bias and dispersion for the batch-estimated satellite clock error, respectively.

There are two classifications of biases contributing to the batch-estimated satellite clock, i.e., the ICB and the TDB, in which we will not conduct additional experiments to verify the impacts of the ICB on the satellite clock batch estimation. The ICB is reflected by the RMS of the batch-estimated satellite clock error. The TDB, another classification of the batch-estimated satellite clock, is impacted by the assimilation of the satellite orbit error. The assimilation ability can be indicated by the correlation coefficients between the satellite clock error and the satellite orbit error. The epoch-wise satellite orbit error can be obtained by comparing it with the satellite orbit reference product provided by IGS. The satellite orbit errors in the radial, along-track and cross-track directions can be computed as [44]

$$
\begin{bmatrix}
\Delta e\_R\\ \Delta e\_A\\ \Delta e\_C
\end{bmatrix} = \begin{bmatrix}
\frac{r}{|r|} & \frac{r \times \overline{v}}{|r \times \overline{v}|} \times \frac{r}{|r|} & \frac{r \times \overline{v}}{|r \times \overline{v}|}
\end{bmatrix}^T \begin{bmatrix}
\Delta x\\ \Delta y\\ \Delta z
\end{bmatrix} \tag{12}
$$

where Δ*eR*, Δ*eA* and Δ*eC* are the satellite orbit errors in the radial, along-track and crosstrack directions, respectively. Δ*x*, Δ*y* and Δ*z* are the satellite orbit errors in the Earthcentered Earth-fixed (ECEF) frame. *r* and *v* is the satellite position and the inertial velocity in the ECEF frame.

The satellite clock batch estimation is implemented based on the fixed ultra-rapid satellite orbit products in observation sessions provided by the three different analysis centers including Wuhan University, Geo Forschungs Zentrum (GFZ) and IGS. These three data sources are abbreviated as WHU, GFU and IGU, respectively. The 1DRMS accuracy of the three ultra-rapid satellite orbit products in the observation session for one-year data are 1.86 cm, 1.31 cm and 0.83 cm, respectively. We use the three different ultra-rapid satellite orbit products to investigate the assimilation ability of the batch-estimated satellite clock to different types of the satellite orbit error [45]. The correlation between the batch-estimated satellite clock error and the fixed satellite orbit error is shown in Figure 3.

**Figure 3.** Correlation between batch-estimated satellite clock errors and fixed satellite orbit errors for one year. The row panels from top to bottom show the results based on the fixed satellite orbit product provided by WHU, GFU and IGU, respectively. The column panels from left to right show the results of the satellite orbit error in the radial, along-track and cross-track directions. The dots with different colors represent different satellites.

It can be seen from Figure 3 that the batch-estimated satellite clock error and the fixed satellite orbit error in the radial direction are positively correlated, whilst the correlation in the along-track and cross-track directions can be neglected. This demonstrates that the batch-estimated satellite clock error based on the global distributed network can assimilate the satellite orbit radial error more than the other two directions. Moreover, the correlation coefficient of WHU, GFU and IGU is different, especially that of IGU, which is much lower than WHU and GFU. This is because the station selection for the global distributed network is specific, which depends on the station observation quality correlated with the fixed satellite orbit accuracy. Thus, the common satellite orbit radial error absorbed by the batch-estimated satellite clock is different, resulting in the difference in the assimilation ability and the corresponding correlation coefficient.

In addition to the absorbed satellite orbit error, the unabsorbed satellite orbit error is absorbed by the observation residual and also affects the TDB of the batch-estimated satellite clock. In order to investigate the impacts of the satellite orbit error on TDB, we extracted the batch-estimated satellite clock error based on the fixed satellite orbit products provided by WHU, GFU and IGU. The TDB is indicated by the dispersion of the batchestimated satellite clock error. The bias of the batch-estimated satellite clock error reflects the ICB. Furthermore, the satellite clock batch estimation based on the final satellite orbit product provided by IGS is also performed as the reference for comparison. The batchestimated satellite clock errors based on the fixed four satellite orbit products, named BEST-WHUO, BEST-GFUO, BEST-IGUO and BEST-IGSO, are shown in Figure 4.

**Figure 4.** Batch-estimated satellite clock errors based on fixed WHUO, GFUO, IGUO and IGSO for one year.

It can be seen from Figure 4 that the batch-estimated satellite clock errors based on these four satellite orbit products have significant biases. This is because the batchestimated satellite clock error contains the ICB, which is induced by the assimilation of the code hardware delay for satellites. The dispersion of the batch-estimated satellite clock errors based on the three ultra-rapid satellite orbit products is significant relative to the IGS final satellite orbit product. This is because the satellite orbit error, including the absorbed and unabsorbed components, affects the TDB of the batch-estimated satellite clock. Furthermore, the STD difference among the BEST-WHUO, BEST-GFUO and BEST-IGUO indicates that the TDB of the batch-estimated satellite clock is different, in which the STD of BEST-IGUO is slightly lower than BEST-WHUO and BEST-GFUO. This is because the assimilation of the satellite orbit radial error based on the IGUO is less than that of WHUO and GFUO. Meanwhile, the unabsorbed component of IGUO has better accuracy than WHUO and GFUO.

#### *3.2. Effect Analysis of PPP Based on Batch-Estimated Satellite Clock*

In order to investigate the impacts of the ICB and TDB for the batch-estimated satellite clock on PPP, the PPP in the kinematic mode will be carried out. The observation model of PPP comes from (9). The elevation-dependent weighting model is used for constructing the stochastic model of PPP. The PPP solution is implemented based on the GAMP software [43]. The PPP positioning errors of the east, north and up directions are calculated epoch-wise by comparing with the reference coordinates provided by IGS. We excluded the positioning results from the last 15 min to avoid the impacts of the satellite orbit extrapolation error on PPP. Furthermore, the convergence time and the positioning accuracy after convergence were used as the indicators to evaluate the PPP positioning performance.

We firstly perform the simulation experiment to clarify the impacts of the ICB on PPP by injecting the artificial constant bias into the satellite clock batch estimation result. The injected bias for each satellite is the time-invariant bias in days, which is consistent with the period of PPP solutions. It is noted that the magnitude of the simulated bias is obtained from the normal distribution with the STD of 3 ns for each satellite, which is consistent

with the accuracy of the broadcast clock. One example of the injected bias for each satellite in one day is shown as Figure 5.

**Figure 5.** One example of injected bias for each satellite in one day.

The IGS final satellite orbit product is used for the satellite clock batch estimation and PPP to avoid the impacts of the satellite orbit error on TDB and PPP. The resulting simulated batch-estimated satellite clock is named BEST + ICB. The PPP experiment based on BEST + ICB can illustrate the impacts of the ICB on PPP. Moreover, we set 20 cm as the convergence threshold to determine the convergence time of PPP. The PPP results for BEST and BEST + ICB are shown in Figure 6.

**Figure 6.** PPP results for BEST and BEST + ICB at 20 arbitrary stations for one year. The two subplots on the left show the STD of the convergence time. The two subplots on the right show the 3DRMS of the positioning accuracy after convergence.

The convergence time for each station varies daily. The STD of the convergence time indicates the variation degree for the observation quality from day to day, including the observation accuracy and the satellite distribution. It can be seen from Figure 6 that the STD of the convergence time for BEST + ICB is longer than BEST. This is because the ICB is absorbed by the ambiguity in the PPP observation model and impacts the convergence time of PPP. Furthermore, we can find that the impacts of the ICB on the convergence time across different stations. It can be explained that the convergence time of PPP is affected by the magnitude of the ICB with different satellites. Moreover, the positioning accuracy after convergence for BEST + ICB is slightly worse than BEST, which is attributed to the injected ICB.

In order to explore the impacts of the TDB on PPP, we simulated the random biases epoch-wise and injected the simulated biases into the batch-estimated satellite clock. The injected bias for each satellite is time-variant, which is the main difference from the injected ICBs. It is noted that the STD of the simulated biases is obtained from the normal distribution with the STD of 0.3 ns for all satellites. The value of the simulated biases comes from three times the STD of the satellite clock error. The impacts of the TDBs on PPP can be amplified and verified. Similarly to the injected ICB, the injected random bias for each satellite is also in days. Furthermore, the batch-estimated satellite clock data for simulation is also based on the fixed final satellite orbit product provided by IGS. The resulting simulated batch-estimated satellite clock is named BEST + TDB. The PPP based on BEST + TDB can illustrate the impacts of the TDB on PPP. Moreover, we set the positioning convergence threshold as 50 cm because of the injected simulated random biases. The PPP results for BEST and BEST + TDB are shown in Figure 7.

**Figure 7.** PPP results for the BEST and BEST + TDB at 20 arbitrary stations for one year. The two subplots on the left show the STD of the convergence time. The two subplots on the right show the 3DRMS of the positioning accuracy after convergence.

It can be seen from Figure 7 that the STD of the convergence time is certainly shortened. This is because of the expanding convergence threshold, so that the positioning result for BEST + TDB can be regarded as the convergence. From Figure 7, we can see that the positioning accuracy after convergence for BEST + TDB is worse than BEST. This is because the simulated TDB based on the normal distribution is absorbed by the observation residual, which impacts the positioning accuracy after convergence. Furthermore, the STD of the convergence time for BEST + TDB is slightly longer than BEST, which does not exceed 30 min. This is because the absorbed observation residual impacts the PPP observation model strength.

#### **4. Discussion**

In order to obtain the PPP results at the centimeter level using ultra-rapid satellite products, the satellite clock batch estimation accuracy and its impacts on PPP are analyzed in this contribution.

The traditional observation models for satellite clock estimation generally include undifferenced, epoch-differenced and mixed-differenced models. For satellite clock batch estimation, the undifferenced observation model is typically used. The ICB and TDB are induced by the batch observation model, including the assimilated ranging errors and the model accuracy. The theoretical analysis and experimental results show that the ICB affects the RMS without affecting the STD, whilst the TDB affects the RMS and the STD of the batch-estimated satellite clock. Furthermore, we primarily analyze the impacts of the fixed satellite orbit error on the TDB. The experimental results show that the assimilation of the satellite orbit radial error to the TDB is more than the along-track and the cross-track directions. The assimilation ability depends on the station distribution used for the satellite clock batch estimation.

In addition to the satellite clock biases, the observation redundancy also affects the STD of the batch-estimated satellite clock. Many studies have shown the number of stations suitable for the real-time satellite clock estimation. Such results may not be suitable for the satellite clock batch estimation. We derive the analytic relationship between the observation redundancy and the satellite clock batch estimation accuracy. The experimental results show that the suitable number of stations is suggested to be 40 to achieve the counterbalance between the efficiency and saturable accuracy.

We perform the simulation experiment to clarify the impacts of the ICB and the TDB on PPP. The experiment results show that the convergence time and the positioning accuracy after convergence of PPP are mainly affected by the ICB and TDB of the batch-estimated satellite clock, respectively. Moreover, besides the ICB and the TDB, the impacts of the observation quality on PPP, including the observation accuracy and the satellite distribution, cannot be ignored.
