**1. Introduction**

The precise satellite clock product is indispensable to obtain high-precision and reliable precise point positioning (PPP) services [1], which can be divided into the real-time, ultra-rapid, rapid and final satellite clock products [2,3]. Compared with the other satellite clock products, the ultra-rapid satellite clock product offers better performance in terms of stability and latency [4,5]. Therefore, the ultra-rapid satellite clock product can be beneficial for the PPP with the relatively low real-time requirement [6]. However, the accuracy of the ultra-rapid satellite clock product needs to be improved to meet the centimeter-level PPP application.

The ultra-rapid satellite clock product typically includes the observation and prediction sessions [7]. The observation session is batch estimated based on the observed receiver–satellite range [8]. In the prediction session, the satellite clock is obtained by the fitting and extrapolation based on the observation session and an accurate satellite clock prediction model [9]. It is noted that the satellite clock batch estimation accuracy determines the quality of the ultra-rapid satellite clock product in the observation and prediction sessions. For the PPP application based on the ultra-rapid satellite clock product, the PPP performance is affected by the satellite clock batch estimation accuracy [3]. Therefore, clarifying the error propagation of the satellite clock batch estimation and its effect on PPP can improve the accuracy of the ultra-rapid satellite clock product and the performance for the PPP application.

The batch observation model needs to be constructed for the satellite clock batch estimation. The traditional observation models for the satellite clock estimation generally include undifferenced, epoch-differenced and mixed-differenced models [10,11]. The

undifferenced and mixed-differenced observation models for the satellite clock estimation are analytically equivalent [10], whilst the epoch-differenced model suffers from the quality of the initial satellite clock, resulting in significant satellite clock biases [11]. For the generation of the ultra-rapid satellite clock product in the observation session, the undifferenced observation model is extensively applied due to the lower observation noise and smaller satellite clock biases than the others. Given the batch observation model, the batch-estimated satellite clock assimilates the ranging errors due to the model strength. Hence, the satellite clock bias is introduced so that the accuracy of the batch-estimated satellite clock is affected. The ranging errors absorbed by the batch-estimated satellite clock mainly include the code hardware delay for satellites and the fixed satellite orbit error. Considering the characteristic of the hardware delay over time, the time-invariant and time-variant hardware delays affect the satellite clock estimation accuracy [10]. For the satellite clock batch estimation within observation arcs, the time-variant hardware delay can be roughly ignored due to its stability [12]. Additionally, the impacts of the satellite orbit error on the real-time satellite clock estimation were already analyzed and verified preliminarily in previous studies [13,14]. The radial and tangential satellite orbit errors can be absorbed to various extents by the estimated satellite clock. This absorption increases as the size of the station network used for the satellite clock estimation decreases, especially for the regional distributed network [13]. In view of the ultra-rapid satellite clock product for the global PPP service, the satellite clock is batch estimated based on the global distributed network [8]. The impacts of the satellite orbit error on the satellite clock batch estimation based on the global distributed network need to be further analyzed. Therefore, the satellite clock batch estimation observation model can be constructed in the observation session by considering the time-invariant hardware delays and the satellite orbit error.

Based on the constructed batch estimation observation model, the satellite clock bias is induced by the assimilated ranging errors and the model accuracy. For the real-time satellite clock estimation based on the undifferenced observation model, the satellite clock biases induced by the model are grouped into the initial clock bias and time-dependent bias [10,15]. The initial clock bias is the satellite clock error at the initial epoch, which depends on the accuracy of the initial satellite clock and the bias introduced by the pseudorange observation [16,17]. The time-dependent bias is determined by phase observation [10]. This classification can be used as the reference for the satellite clock batch estimation. However, for the satellite clock batch estimation, the initial clock bias was predetermined with relatively high precision. The impacts of the initial satellite clock accuracy on the initial clock bias can be ignored, which will not be further analyzed in this contribution. Moreover, the accuracy of the fixed satellite orbit product is not precise enough. The satellite orbit error needs to be considered in the classification of the batch-estimated satellite clock bias. Therefore, referring to the typical classification of the satellite clock bias in [10], we clarify the source of the batch-estimated satellite clock biases to analyze the impacts of these biases on the satellite clock batch estimation accuracy.

In addition to satellite clock biases, observation redundancy is another critical factor affecting the satellite clock batch estimation accuracy. Generally, the parameter estimation accuracy depends on the observation redundancy [18,19]. Since the satellite clock estimation accuracy becomes approximately saturable when the observation redundancy has adequately increased. The benefit of increasing more observation redundancy is negligible for the accuracy improvement of the satellite clock estimation [20]. Furthermore, the excessive observation redundancy will burden the efficiency and may increase the update delay, which is not conducive for the generation of the ultra-rapid satellite clock product with low latency [4,5]. The lower latency for the batch-estimated satellite clock brings the better performance for the satellite clock prediction and the PPP based on the predicted satellite clock [21,22]. The update interval and the update rate simultaneously determine the latency of the batch-estimated satellite clock, which can be improved by the efficient sliding satellite clock batch estimation with a short window and the satellite clock prediction with a short period [23,24]. In order to achieve the counterbalance between efficiency and estimation accuracy, many studies have shown the number of stations suitable for the real-time satellite clock estimation [11,20]. Such a number of stations may not be suitable for satellite clock batch estimation because the observation redundancy has already been improved for the batch estimation model by combining observations within observation arcs. In addition, the variations in the number of stations may affect the station distribution. Since the ultra-rapid satellite clock product for the global PPP users is based on the global distributed network, the impacts of the baseline distance among stations on the satellite clock batch estimation can be ignored and will not be considered in this contribution [25]. Therefore, the number of stations suitable for the satellite clock batch estimation needs to be determined to achieve the counterbalance between the efficiency and the batch estimation accuracy.

For the PPP based on the estimated satellite clock, the impacts of the estimated satellite clock on PPP have already been studied in many studies [10,26]. The absolute accuracy, latency and sampling interval of the satellite clock product are the main factors affecting PPP performance, including the convergence time and the positioning accuracy after convergence [21,26]. Furthermore, for the PPP based on the real-time satellite clock, the impacts of the classified satellite clock biases on PPP are analyzed [10]. Two such types of satellite clock biases, i.e., the initial clock bias and the time-dependent bias, affect the pseudorange and phase observations in the typical PPP observation model, which determine the convergence time and the positioning accuracy after convergence, respectively [15,27]. However, few studies have performed the comprehensive analysis for the impacts of the batch-estimated satellite clock on PPP. Therefore, we clarify the error propagation of the batch-estimated satellite clock biases in the PPP observation model and analyze their impacts on PPP.

We begin by constructing the batch observation model for the satellite clock batch estimation and clarify the error propagation of the satellite clock batch estimation. Then, the batch-estimated satellite clock biases are classified according to their impacts on the satellite clock batch estimation accuracy. Furthermore, we derive the analytic relationship between the observation redundancy and the satellite clock batch estimation accuracy to determine the number of stations suitable for satellite clock batch estimation. Moreover, we clarify the error propagation of the batch-estimated satellite clock biases in the PPP observation model. In the experimental verification, we present the data collection and validation strategies for the satellite clock batch estimation and PPP. The effect analysis of the satellite clock batch estimation and the PPP based on the batch-estimated satellite clock are verified. Finally, we give the discussion and summarize the main conclusions of this study.

#### **2. Methodology**

We first introduce the typical observation model of the satellite clock batch estimation. The source of the satellite clock bias induced by the batch observation model is clarified and classified. Then, based on the constructed batch observation model, the observation redundancy suitable for the satellite clock batch estimation is determined by deriving the covariance–variance matrix. Finally, the impacts of the classified satellite clock biases from the batch estimation on PPP are illustrated.

#### *2.1. Sources of Batch-Estimated Satellite Clock Biases*

The undifferenced observation model based on the batch estimation is adopted due to the lower observation noise and the smaller satellite clock biases than other traditional satellite clock estimation models [10]. Limited by the model strength, the batch observation model causes satellite clock biases in the batch-estimated satellite clock [10,28]. According to the impacts of the resulting satellite clock biases on the satellite clock batch estimation accuracy, they are categorized to clarify the error propagation in the satellite clock batch estimation.

The dual-frequency ionosphere-free (IF) combination is widely used for the undifferenced observation model to eliminate the higher-order ionospheric delay. The tropospheric delay is usually corrected for its dry component with the a priori model, and the wet

component is estimated with the zenith tropospheric delay and its corresponding mapping function [29]. The station coordinate and the satellite orbit are usually fixed in the IF observation model [30]. It is assumed that the station coordinate is precisely known, which is normally fixed to the PPP weekly solution, while the error exists in the fixed ultra-rapid satellite orbit product due to the different processing strategies for the satellite orbit determination such as the satellite attitude model and the solar radiation pressure model [31,32]. For the batch observation model, all the observations and parameters in the observation session are constructed, which can be expressed as

$$
\begin{bmatrix} \mathbf{v}\_{P,IF,r}^{i} \\ \mathbf{v}\_{L,IF,r}^{i} \end{bmatrix} = \begin{bmatrix} \mathbf{G}\_{r} & \mathbf{G}^{i} & \mathbf{G}\_{\mathcal{Z}\mathcal{T}\mathcal{D}} & \mathbf{G}\_{r} & \mathbf{G}^{i} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{G}\_{r} & \mathbf{G}^{i} & \mathbf{G}\_{\mathcal{Z}\mathcal{T}\mathcal{D}} & \mathbf{G}\_{r} & \mathbf{G}^{i} & \mathbf{G}\_{r} & \mathbf{G}^{i} & \mathbf{G}\_{N} \end{bmatrix} \begin{bmatrix} \delta t\_{r} \\ \delta t^{i} \\ \mathbf{Z}\mathbf{T}\mathbf{D}\_{r} \\ \mathbf{b}\_{r} \\ \mathbf{b}^{i} \\ \mathbf{B}\_{r} \\ \mathbf{B}\_{r} \\ \mathbf{N}\_{r}^{i} \end{bmatrix} - \begin{bmatrix} \mathbf{f}\_{P,IF,r}^{i} + \mathbf{G}^{i} \cdot \mathbf{\mu}\_{r}^{i} \delta \mathbf{x}^{i} \\ -\left[\mathbf{f}\_{P,IF,r}^{i} + \mathbf{G}^{i} \cdot \mathbf{\mu}\_{r}^{i} \delta \mathbf{x}^{i}\right] \\ \mathbf{f}\_{r}^{i} \end{bmatrix} \tag{1}
$$

where *<sup>G</sup><sup>r</sup>* <sup>=</sup> *<sup>I</sup><sup>m</sup>* <sup>⊗</sup> *<sup>e</sup>n*, *<sup>G</sup><sup>s</sup>* <sup>=</sup> <sup>−</sup>*e<sup>m</sup>* <sup>⊗</sup> *<sup>I</sup>n*, *<sup>G</sup>ZTD* <sup>=</sup> *<sup>I</sup><sup>m</sup>* <sup>⊗</sup> *<sup>M</sup>n*, *<sup>G</sup><sup>N</sup>* <sup>=</sup> *<sup>I</sup><sup>m</sup>* <sup>⊗</sup> *<sup>e</sup>n*, *<sup>G</sup>* is the designed matrix for the satellite clock batch estimation model, *<sup>s</sup>* is the satellite, *<sup>r</sup>* is the receiver, *<sup>I</sup><sup>i</sup>* is an *<sup>i</sup>* <sup>×</sup> *<sup>i</sup>* identity matrix, *<sup>e</sup><sup>i</sup>* is an *<sup>i</sup>* <sup>×</sup> 1 vector with all elements equal to 1, *<sup>M</sup><sup>i</sup>* is an *<sup>i</sup>* <sup>×</sup> <sup>1</sup> vector containing the tropospheric mapping function for each receiver, *m* is the number of the observations, *n* is the number of the observed satellites at different stations, ⊗ is the Kronecker product operation [33], *P* is the pseudorange observation, *L* is the phase observation, *IF* is the IF combination, *v* is the vector of the a posterior observation residual for each observation, *μ<sup>s</sup> <sup>r</sup>* is the unit vector from the satellite to the receiver, *<sup>δ</sup>x<sup>s</sup>* is the vector of the satellite orbit error for each satellite, *<sup>δ</sup>t<sup>r</sup>* is the vector of the receiver clock for each receiver, *<sup>δ</sup>t<sup>s</sup>* is the vector of the satellite clock for each satellite, *ZTD* is the vector of the zenith tropospheric delay for each receiver, *N* is the vector of the ambiguity for each satellite and each receiver, *<sup>b</sup><sup>r</sup>* is the vector of the code hardware delay for each receiver, *<sup>b</sup><sup>s</sup>* is the vector of the code hardware delay for each satellite, *<sup>B</sup><sup>r</sup>* is the vector of the phase hardware delay for each receiver, *B<sup>s</sup>* is the vector of the phase hardware delay for each satellite and *l* is the vector of the a priori observation residual for each observation. In addition to the ranging errors in (1), the satellite and receiver antenna phase center offsets and variations, relativity, tidal loadings and phase windup should be modeled and precisely corrected in the a priori observation residual [30].

It can be seen from (1) that all the unknown parameters except the zenith tropospheric delay are linearly dependent. The satellite clock, receiver clock, zenith tropospheric delay and ambiguity parameters are the parameters to be estimated. Hence, other unknown parameters are absorbed by the estimated parameters or remain in the observation residual. Such absorption should satisfy the requirement of minimizing the weighted sum of the squares of residuals. According to the correlation between the parameters, and the consistency of the batch observation equation, *<sup>δ</sup>x<sup>s</sup>* , *br*, *b<sup>s</sup>* , *<sup>B</sup><sup>s</sup>* and *<sup>B</sup><sup>r</sup>* are absorbed by *<sup>δ</sup>t<sup>s</sup>* , *<sup>δ</sup>tr*, *<sup>δ</sup>t<sup>s</sup>* , *Ns <sup>r</sup>* and *<sup>N</sup><sup>s</sup> <sup>r</sup>*, respectively.

We focus on the ranging errors absorbed by the batch-estimated satellite clock, i.e., the code hardware delay for satellites and the satellite orbit error. Since the hardware delay can be treated as the constant within the observation session, the absorbed code hardware delay for satellites causes the time-invariant bias dependent with satellites. The timeinvariant bias of the batch-estimated satellite clock will not affect the standard deviation (STD). For the satellite orbit error, not all of the satellite orbit errors are absorbed by the batch-estimated satellite clock. The assimilation ability depends on the station network distribution used for the satellite clock batch estimation [13]. The pseudorange and phase observations of the global distributed network are used for generating the worldwide satellite clock product. Since the satellite can be approximately regarded always directly above the global distributed network, the projection directions of the satellite orbit radial and tangential errors for each satellite in the signal propagation direction are identical and opposite, respectively. It is well known that the angles of the satellite parallax for the global

distributed stations are different, resulting in different projection amounts of the satellite orbit errors. Therefore, the common satellite orbit radial errors observed by different stations in the signal propagation direction can be absorbed by the batch-estimated satellite clock. On the contrary, the satellite orbit tangential error cannot be absorbed due to the opposite projection direction. Therefore, the common satellite orbit error can be absorbed by the batch-estimated satellite clock epoch-wise, which belongs to the time-variant bias and affects the STD of the satellite clock batch estimation. Furthermore, the satellite orbit error not absorbed by the batch-estimated satellite clock is absorbed by the observation residual. The STD of the satellite clock batch estimation is still affected.

It can be seen from (1) that the rank of the coefficient matrix presents deficiency due to the linear dependency between the satellite clock and receiver clock. The clock constraint should be introduced to avoid the rank deficient [34,35]. Therefore, the batch-estimated receiver and satellite clocks are aligned as the clock relative to the reference clock. The reference clock strategies include selecting one satellite clock as zero, one receiver clock as zero and the mean value of satellite clocks as zero, which are equivalent to each other [36]. In order to ensure the stability of the satellite clock batch estimation, we select the zeromean condition as the reference clock [37]. Inevitably, the batch-estimated satellite clock contains the timescale difference caused by the clock bias of the selected reference clock [15]. The timescale difference belongs to the time-invariant bias independent with satellites, which does not affect the satellite clock batch estimation accuracy when the accuracy of the reference clock is better than 10−<sup>6</sup> s [15].

In summary, the reparametrized satellite clock batch estimation model due to the parameter assimilation can be redefined 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}\_{ZTD} & \mathbf{0} \\ \mathbf{G}\_{r} & \mathbf{G}^{s} & \mathbf{G}\_{ZTD} & \mathbf{G}\_{N} \end{bmatrix} \begin{bmatrix} \delta \mathbf{f}\_{r} \\ \delta \overline{\mathbf{f}}^{s} \\ \mathbf{Z} \mathbf{T} \mathbf{D}\_{r} \\ \overline{\mathbf{N}}\_{r}^{s} \end{bmatrix} - \begin{bmatrix} \mathbf{f}\_{P,IF,r}^{s} \\ \mathbf{f}\_{L,IF,r}^{s} \end{bmatrix} \tag{2}
$$

where *<sup>δ</sup>tr*, *<sup>δ</sup>t s* , *Ns 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, satellite clock, ambiguity, a priori pseudorange observation residual and a priori phase observation residual, respectively. They can be rewritten as

$$\begin{cases} \delta \overline{t}\_{\mathcal{I}} = \delta t\_{\mathcal{I}} + \mathsf{b}\_{\mathcal{I}} \\ \delta \overline{\mathcal{F}}^{s} = \delta t^{s} + \alpha \mu\_{r}^{s} \delta \mathbf{x}^{s,R} + \beta \mu\_{r}^{s} \delta \mathbf{x}^{s,T} + \mathsf{b}^{s} + \mathsf{D}^{s} \\ \end{cases} \tag{3}$$
 
$$\begin{cases} \overline{\mathcal{V}}\_{r}^{s} = \mathcal{N}\_{r}^{s} + \mathsf{B}\_{r} - \mathcal{B}^{s} + \mathsf{b}^{s} - \mathsf{b}\_{r} \\ \overline{\mathcal{I}}\_{P,IF,r}^{s} = \mathcal{I}\_{P,IF,r}^{s} - \mathcal{G}^{s} \cdot (1 - \alpha) \mu\_{r}^{s} \delta \mathbf{x}^{s,R} - \mathcal{G}^{s} \cdot (1 - \beta) \mu\_{r}^{s} \delta \mathbf{x}^{s,t} \\ \overline{\mathcal{I}}\_{L,IF,r}^{s} = \mathcal{I}\_{L,IF,r}^{s} - \mathcal{G}^{s} \cdot (1 - \alpha) \mu\_{r}^{s} \delta \mathbf{x}^{s,R} - \mathcal{G}^{s} \cdot (1 - \beta) \mu\_{r}^{s} \delta \mathbf{x}^{s,T} \end{cases} \tag{3}$$

where *D<sup>s</sup>* is the timescale difference, *<sup>δ</sup>xs*,*<sup>R</sup>* and *<sup>δ</sup>xs*,*<sup>T</sup>* are the satellite orbit radial and tangential errors, *α* and *β* are the assimilation proportion of the batch-estimated satellite clock to the satellite orbit radial and tangential errors. The absorption ability can be calculated by using the range of the satellite parallax angle during satellite regular motion [13]. Taking GPS as an example, the estimated satellite clock can absorb 97.1% of the satellite orbit error in the radial direction at least, and 24.0% in the tangential direction at most. The specific assimilation proportion depends on the station network distribution used for the satellite clock batch estimation.

The impact of the observation noise needs to be accounted for the satellite clock batch estimation. Influenced by the ambiguity in the phase observation, the variation and the absolute bias of the batch-estimated satellite clock are determined by the phase and pseudorange observation, respectively [38]. Furthermore, any biases in the ambiguity or the undetected cycle slips in the data preprocessing will cause the time-invariant bias in the batch-estimated satellite clock. Moreover, for the real-time satellite clock estimation, the time-invariant bias will also be introduced by the initial satellite clock with poor accuracy, which arises from the broadcast clock or the predicted satellite clock [10]. However, the impacts of the initial satellite clock on the satellite clock batch estimation can be ignored. This is because the initial value of the satellite clock can be predetermined with high precision in the batch estimation.

In order to categorically analyze the impacts of the satellite clock biases on the accuracy of the batch-estimated satellite clock, the time-invariant and time-variant bias induced by the batch observation model are denoted as the initial clock bias (ICB) and time-dependent bias (TDB), respectively. The ICB and TDB are both estimated satellite clock biases, which cause the satellite clock error compared with the satellite clock reference product. The only difference between them is that the ICB and TDB contribute to the constant and variable components of the batch-estimated satellite clock errors, respectively. The ICBs caused by the code hardware delay assimilation, the reference clock selection and the pseudorange observation noise will not affect the STD of the batch-estimated satellite clock, but will affect the root mean square (RMS). The TDB induced by the satellite orbit error and the phase observation noise will affect the STD and RMS.
