Improved Strategies for BeiDou Ultrarapid Satellites’ Clock Bias Prediction Using BDS-2 and BDS-3 Integrated Processing

GNSS ultrarapid clock biases are key inputs of rapid high-accuracy applications, especially for its prediction parts. With the fast development of the BeiDou system (BDS), the system performances are mainly represented by orbit and clock products. However, it is suggested that the BDS-predicted clock biases cannot meet the requirement of real-time or near real-time services. In this research, the BDS satellite-predicted ultrarapid clock bias products are optimized with three methods, namely, one-step strategy, intersatellite correlation, and variogrammodel, using the combined estimation of BDS-2 and BDS-3 satellites. Firstly, considering the traditional two-step strategy for modelling clock bias prediction, we take all terms (including trend and periodic terms) into one-step solution of model estimation based on the sparse modelling in machine learning. Secondly, because of the much more stable on-board atomic clock of BDS-3 satellites, the intersatellite correlations between BDS-2 and BDS-3 are utilized to enhance the solution of model coefficients. -irdly, to further improve the model, the temporal correlations in model residuals are used to reconstruct the stochastic function obtained by variogram. In addition, to verify the proposed improved strategies, 12 schemes of BDS clock bias prediction experiments are designed and analyzed with different conditions. According to the results of predicted clock biases, it is indicted that (1) the stability of BDS-3 on-board clocks is more optimal compared with BDS-2, which can be used to strengthen the solution of the clock bias prediction model; (2) the one-step estimation of the clock bias model by sparse modelling can slightly increase the accuracy of prediction results; (3) both BDS-2and BDS-3-predicted clock biases benefited each other by inserting the intersatellite correlations into the weight matrix, in which the accuracy of 18-hour period with one-step strategy can be improved by 28.6% and 27.2% for BDS-2 and BDS-3, respectively; and (4) after the introduction of the variogram model in updating the weight matrix, the clock bias predictionmodel is further corrected by 8.0% and 11.1% for BDS-2 and BDS-3. In summary, improved strategies for BDS ultrarapid satellites’ clock bias prediction using BDS-2 and BDS-3 integrated processing are meaningful for the current BDS ultrarapid satellites’ clock bias prediction products.


Introduction
China BeiDou satellite system (BDS) is constructed with three steps, called BDS-1 (the demonstration system), BDS-2 (the regional service system), and BDS-3 (the global cover system) [1,2]. Currently, more than 40 BDS satellites are in orbit, including 16 BDS-2 and 25 BDS-3. With the flourishing development of the BDS innovative applications, the precision and accuracy of satellite orbit and clock products impose a more important role for different navigation satellite systems [3]. In addition, compared with GPS, the precise products (especially for orbit and clock) of BDS should be paid more attention [4], which represents the main parts of system service performances. Among orbit and clock bias products, the satellite on-board atomic clocks, as one of the main loads, are the references for precise orbit determination and high-accuracy services [5,6]. Currently, considering different accuracy, precision, and time delays, six types of clock bias products are issued by GNSS analysis centers (ACs) of the International GNSS Service (IGS) [7].
However, regarding the real-time or near real-time centimeter-level services, the ultrarapid clock products should be selected and inserted into applications [8][9][10].
According to the products issued by IGS, it is shown that the accuracy of observed and predicted parts of GPS clock biases are about 0.15 ns and 3.00 ns, respectively. However, the values are increased to 0.60 ns and 6.00 ns for the international GNSS Monitoring and Assessment Service (iGMAS). Moreover, for BDS ultrarapid clock products, the median errors obtained from iGMAS ACs are between 0.60 ns and 1.50 ns and 2.50 ns and 4.50 ns for the observed and predicted products, which is far behind the needs of real-time centimeter-lever applications [8]. Hence, more attention should be paid on the improvement of accuracy and precision for BDS satellites' ultrarapid clock biases, especially the prediction parts, in aspect of the availability, stability, and reliability for GNSS products.
According to research studies about satellite clock bias prediction, main points are focused on three contents: the preprocessing strategies for modelling satellites' clock biases [11,12]; analyzing and improving the prediction models [7,8,11,[13][14][15]; and testing and discussing the impact of environments and clock characteristics on modelling satellites' clock biases [16,17]. Furthermore, for algorithms and strategies of clock bias prediction, researchers mainly highlight the investigation of long-term clock products and its optimal solutions, such as the expanded state modelling [18], the introduction of artificial neural network algorithm [19], the combination of multiple periodic terms with trend terms based on an iterative method [13], and the sidereal filtering algorithm including the terms of subdaily clock bias variations [7]. Overall, the widely used clock bias prediction models can be summarized with two types, namely, the polynomial model and the nonlinear model [13,20]. Nevertheless, compared with the predicted strategy proposed by Mao et al. [5], all methods for satellite clock bias predictions rule out the impact of intersatellite correlation on the estimation of model parameters, while it was verified to reduce the accuracy of BDS-2 satellites' clock bias prediction. Moreover, some researchers present that the correlation can be excavated as a constraint condition to improve the estimated parameters, such as satellite orbit determination [21,22]. erefore, the intersatellite correlation should be considered as an alternative method to conduct the prediction of BDS-2 and BDS-3 satellites' integrated ultrarapid clock biases.
Currently, for BDS ultrarapid clock bias predictions, the widely used strategies, combining the trend (polynomial model) with periodic terms [11,23], take the methods of GPS as references. Additionally, an optimized strategy was selected to conduct the prediction of BDS ultrarapid satellites' clock biases, in which the nonlinear terms were analyzed and modelled by a backpropagation neural network (BPNN) algorithm [8]. At present, more than 40 BDS satellites, including BDS-2 and BDS-3, jointly provide the high-accuracy positioning, navigation, and timing services (PNT) [24]. For BDS-3 satellite on-board clock, the advantages and characteristics of BDS-3 experimental satellite clock bias prediction with the aid of intersatellite linking observations were discussed [25]. Furthermore, it is indicated that the stability of BDS-3 satellite on-board atomic clock frequency is better than that of BDS-2 with more than an order of magnitude [26]. In our previous research studies, the preliminary results of BDS-2 and BDS-3 satellites' integrated prediction for ultrarapid products were outputted and analyzed by the introduction of the intersatellite correlation [6], and its impact on ultrarapid estimation was also investigated [24]. According to one-month period of experiments, the results indicated that BDS ultrarapid satellitepredicted clock bias products can be optimized with the intersatellite correlation, which increased the accuracy of 18hour prediction by 30.7%-47.3% and 49.9%-59.3% for BDS-2 and BDS-3, respectively. erefore, to further improve the quality of BDS ultrarapid clock bias products, the integrated processing of BDS-2 and BDS-3 clock bias series should be considered in the following research studies.
Compared with the widely used strategy of modelling clock bias prediction [5,8], two main issues should be discussed and optimized in the estimation of model parameters: (1) the steps of modelling are usually divided into two parts, namely, trend and periodic terms, in which the correlations among unknown parameters are ignored [6]; (2) only the significant periodic terms (two in BDS-2 and three in BDS-3) in modelling were considered [24], which inevitably results in significant fitting residuals. erefore, to obtain a more accurate and reliable clock bias prediction model, the strategy should be further improved and optimized, such as the optimal clock bias model parameters. Shi et al. [27] discussed the optimal dimensions for different on-board atomic clocks, in which different types of clocks demonstrated different characteristics. In this research, the parameters of the clock bias model will be estimated in onestep solution with a unified model, which must take the accuracy of the prediction model and the numbers and types of parameters (numbers of periodic and trend terms) into consideration. Moreover, to select an optimal model for clock bias prediction, the machine learning will be referred as the wealth achievements in sparse modelling [28].
In this contribution, the prediction strategy of BDS ultrarapid satellite clock biases is investigated with the integrated processing of BDS-2 and BDS-3 clock bias series in aspect of the disadvantages in traditional methods. Firstly, an optimal model selection method is proposed by the sparse modelling to improve the two-step strategy for clock product predictions. en, we design an integrated processing strategy to fully use the intersatellite correlations between BDS-2 and BDS-3 satellites. Finally, considering the spatial and temporal variations of prediction model residuals, an empirical model will be established to adjust the clock bias model from updating the weight matrix.

Methods
In modelling of clock bias series prediction, the BDS-3 satellites are equipped with high-quality atomic clocks compared with BDS-2. In order to typically obtain the benefit information of BDS-3 satellites and improve the solution of BDS-2 and BDS-3 integrated estimation, this study proposes a method in which intersatellite correlation is extracted in calculating the model parameters. Moreover, according to the disadvantages of strategy for modelling clock bias, such as two-step estimation, all parameters of the clock bias model 2 Mathematical Problems in Engineering (trend and periodic terms) are estimated in one step. Because of different characteristics for different types of satellite clocks, the optimal model should be selected for each satellite, which will be obtained by the method of sparse modelling. Moreover, the residuals of clock biases will be modelled with the empirical function by variogram.

Modelling BDS-2 and BDS-3 Satellites' Integrated
Processing. After rigorously preprocessing original clock bias series, such as gross detection and repair and denoising [6], the clock bias series of the estimated parts can be presented as L. Referring to the relevant research studies [5,6,8], we can assume the clock bias series model as where t i and k are the t i -th epoch time and the k-th satellite, respectively; a 0 , a 1 , and a 2 are the polynomial coefficients of trend terms; n and j are the total numbers for periodic terms and the j-th period; A j, B j , and T j represent the amplitudes of the periodic terms and their corresponding period; and (t i ) are the model residuals of equation (1). In the estimation of periodic terms T j , the fast Fourier transform (FFT) is typically chosen to calculate the parameters of T j based on the residuals by removing trend terms. erefore, the satellite clock bias error equations are denoted as where V k is the clock bias series error. In the widely used strategy of modelling clock bias series [6,24], the two-step method for the solution of model coefficients, trend and periodic terms, is divided into two independent parts, which should be theoretically obtained in one-step estimation. Moreover, due to the differences among on-board atomic clocks, the number and types of model coefficients cannot be expressed as same forms for all satellites. erefore, to acquire a more accurate and precise clock bias model, it is necessary to further improve the parameter estimation and selection. Firstly, in equation (2), all parameters are considered into one-step estimation. en, equation (2) can be expressed as b satellites and s epochs (all satellites and epochs in one solution), whose matrix form is as follows: Mathematical Problems in Engineering where (t i ) can be presented as In equation (4), A, X, and P are the coefficient matrices of the clock bias model, the model parameters, and the weight matrix of each epoch, respectively. In equation (5), δ bb can be calculated from δ bb � STD b , which denotes the standard deviation of the residuals for the b-th satellite. en, we can take the covariance of the k-th and the b-th satellite asδ kb � |r kb | · ������ � δ kk · δ bb . It should be noted that r kb is the intersatellite correlation factor, which should be deduced from the correlation function between BDS-2 and BDS-3 clock bias series. AssumeL k as the average values of satellites' clock bias series; then, In the estimation of unknown parameters of clock bias models, the weight matrix must be adjusted by epoch time: where Δt denotes the interval of adjacent epoch. From equations (2)-(7), all coefficients of clock bias model are estimated in one-step solution, which is augmented by extracting the intersatellite correlations of BDS-2 and BDS-3 integrated processing.

Sparse Modelling for Clock Bias Prediction.
Based on Section 2.1, the estimation of model parameters takes the intersatellite correlation into consideration to indirectly improve the solutions. However, two potential flaws can cause the failure of parameter estimation: (1) obviously, all terms (trend and periodic terms) considered can describe the clock bias model with more accuracy; because of the limited observations (observed clock bias series) and model overfitting, the stability and accuracy of coefficients are inevitably influenced; (2) different satellites with different types of on-board atomic clocks will present different characteristics [27], which cannot be described by an uniform model (same parameters for all satellites). erefore, the strategy for solving the modelling of clock biases should be further optimized. In this study, the sparse modelling of machine learning is referred to accurately construct the clock bias model. In theory, the introduction of all independent variables for modelling can efficiently improve the results of prediction, which also increases the complexity and overfitting for parameter estimation. In practice, to improve the interpretation and accuracy of the prediction model, the independent variables should be used, namely, model selection (or variable selection) [28,29]. In general, the model selection can only compare and select a limited range of alternative models. However, the sparse modelling can effectively choose the basic data and reduce the number of independent variables, so as to optimize the accuracy and computational efficiency of prediction.
To conduct the sparse modelling of clock bias series, the Lasso (least absolute shrinkage and selection operator) algorithm is used, which was one of the variable selection methods proposed by Tibshirani [30]. In Lasso algorithm, the parameter estimation and model selection are simultaneously solved by minimizing the target function of multivariable linear regression model and imposing the constraint on the absolute value of regression coefficients. According to the solution of Lasso algorithm, some variables are compressed and excluded. e algorithm implementation of Lasso is obtained by L 1norm based on the least squares (LS): where ω is a tuning parameter. Assume the solution of LS estimation as d ; then, set ω 0 � 3+2n d�1 |X d | ≤ ω; the model selection is obtained as the value of ω causing shrinkage of the solutions towards zero; thus, some coefficients of X can be exactly equal to zero. According to the Lagrange multiplier algorithm, the target function for the Lasso can be set as where λ denotes the Lagrange multiplier, which will be selected by adjusted R-square method from an alternative dataset. Equation (9) can be explained as the combination of sum of errors squares and the regularization function. Compared with Tikhonov regularization algorithm (a L 2norm algorithm), equation (9) (a L 1 -norm regularization) can automatically select the most optimal model by sparse modelling without loosing the convexity of target function. To obtain the numerical solution of equation (9), the iteratively reweighted least squares (IRLS), one of the most effective and simplest methods, is used to solve the Lasso problem [31]. Equation (9) can be read as Set X (m) as the result of the m-th iteration; then, the approximate solution of equation (10) is outputted as where U (m) is a diagonal weight matrix, in which the element is determined based on the results of m − 1-th step solution.
For example, the f-th diagonal element is the f-th value inX (m− 1) . So, From the matrix inversion lemma, the solution of X (m) is During the iteration processing, the initial values of model coefficients can be set asX (0) � (A T PA) −1 A T PL; the initial weight matrix can be denoted as an identity matrix.

Modelling Residuals by the Variogram.
According to the solution of clock bias model, the prediction model can be well described with the estimated coefficients. However, due to model errors and noise of satellite on-board clock, the residuals of prediction model should be further improved to increase the accuracy and precision of modelling clock biases. In previous research studies, several strategies were mainly focused on the neural network to analyze the part of residuals in the prediction model [5,8]. Moreover, the combination of the partial least squares (PLS) and the BPNN algorithms was proposed and discussed to conduct the ultrarapid clock bias prediction [6,24]. However, due to the lack of accuracy information for satellite on-board clock, it is difficult to model and predict the model residuals with a desired method, such as the abnormal values of residuals occurring in certain epochs [24]. Because of the characteristics of the clock bias series, the model residuals are inevitably presented with the spatiotemporal correlation for each satellite. In this contribution, to fully use the potential spatial and temporal correlation between BDS-2 and BDS-3 satellites' clock bias series, the variogram of residuals is taken as an alternative method. More details concerning variogram are discussed and summarized in several relevant research studies [27,[32][33][34]. In fact, we can use a modified equation to express the variogram of residuals as [27,32] 2c where2c(h)is the estimated variogram value; (t) represents the residuals of clock bias model in equation (1); t l -t i � h; and |N(h)|is the number of elements of N(h). According to equation (14), the variogram can be further expressed as where C is the covariance operation. C(t i − t l ) can be obtained by where cov is the covariance operation. e temporal characteristic of model residuals will be outputted by the quantity variogram. After presenting the quantitative estimation of the residuals for each satellite, an empirical variogram model will be summarized, such as spherical, exponential, and Gaussian. Based on the model variogram (equations (17) and (18)), as h ⟶ ∞, 2c(h) ⟶ 2C(0); thus, from equation (15), the covariance C(h) can be determined by the model variogram. Hence, the weight matrix of parameter estimation will be constructed based on equation (16), which will be taken into equation (3) to calculate the model coefficients again.
From the strategies of intersatellite correlations, sparse modelling, and variogram, the clock bias prediction model of BDS satellites will be obtained. Moreover, to illustrate and show more details of the improved strategies by the integrated processing of BDS-2 and BDS-3 ultrarapid clock bias predictions, we listed all steps of data processing in Figure 1.

Experiments and Results
In this section, the improved strategies for BDS-2 and BDS-3 integrated clock bias prediction are tested and analyzed with 12 schemes from the intersatellite correlations, sparse modelling, and variogram model. In experiments, one month (DOY 41-70, 2019) of BDS-2 and BDS-3 integrated ultrarapid orbit determinations was firstly implemented to prepare the original ultrarapid clock bias products. Details of strategy for orbit determination mainly refer to our previous research [35]. Additionally, observations of MGEX and iGMAS are collected for orbit determination, in which 27 stations can track BDS-3 satellites with B1I and B3I signals marked by different colors in Figure 2. As described above, experiments are classified as three groups from intersatellite correlations, sparse modelling, and the variogram model. Because of the low accuracy of BDS GEO satellite clock biases, we exclude it from all experiments.
Before the clock bias prediction, to illustrate the feasibility of improving the BDS prediction model by advantages of Mathematical Problems in Engineering BDS-3 satellites, the frequency stability of BDS-2 and BDS-3 on-board clock is taken as an example. We select the Allan deviation to describe the differences of on-board clock between BDS-2 and BDS-3 satellites. In Figure 3, results of DOY 41, 2019, are listed with different satellite types. Moreover, the average daily RMS values of one month are calculated in Tables 1 and 2 for BDS-2 and BDS-3, respectively.
From Figure 3 and Tables 1 and 2, the results of experiments indicate that the performance of frequency stability of BDS-3 satellites' clocks is better than that of BDS-2. Due to the undesired quality of clock bias series based on the limited tracking stations (B1I and B3I observations), the experimental results about the order of magnitude for the differences between BDS-2 and BDS-3 satellites are smaller than the widely discussed ones. In summary, because of the more accurate information of BDS-3 satellites on-board clock, it is suggested that the intersatellite correlations between BDS-2 and BDS-3 satellites can be used to improve the estimation of BDS satellite clock bias models based on the combined processing of BDS-2 and BDS-3.
en, three groups of experiments are conducted to verify the proposed methods of BDS clock bias prediction.   e first group of experiments includes four schemes to analyze the methods of one-step strategy with the aid of sparse modelling.
Scheme 1: based on trend and periodic terms, the clock bias prediction model is established. In this scheme, the coefficients of trend terms are obtained by fitting the clock bias series, while periodic terms are firstly determined by FFT based on the residuals of models. It should be noted that the trend terms are modelled by three coefficients (a 0 , a 1 , a 2 ) of polynomial function for all satellites, and the significant periods are chosen as discussed in Wang et al. [24]. Scheme 2: similar to Scheme 1, the periodic terms of BDS-2 satellites are taken from Huang et al. [8] and Mao et al. [5] as references. ree significant periodic terms of BDS-3 are used to conduct the clock bias prediction model. Scheme 3: compared with Schemes 1 and 2, FFT is used to analyze the periodic terms by removing trend terms. en, all periodic terms are considered into modelling clock bias series. It should be noted that the trend terms (a 0 , a 1 , a 2 ) and the amplitudes of all periodic terms are included in the unknown coefficients. Scheme 4: to test the sparse modelling, the Lasso algorithm is used to estimate the coefficients of the clock bias model by one-step strategy based on Scheme 3.
According to above four schemes, the BDS-2 and BDS-3 satellites' integrated clock biases are predicted with one-day period. To assess the accuracy of prediction model, the rapid BDS satellites' clock biases issued by WHU (Wuhan University, iGMAS AC) are taken as references. e residuals of the prediction model and accuracy of predicted clock biases for BDS-2 and BDS-3 satellites are shown in Figures 4 and 5, respectively. To analyze the prediction model, the average daily RMS values of fitting residuals within 30-day experiments are calculated in Table 3. Furthermore, in Table 4, the       According to four schemes of experiments, the results for prediction model can be summarized as follows: (1) in Figure 4 and Table 3, the fitting residuals of BDS-3 satellites are significantly smaller than that of BDS-2, which can reflect that the performance of BDS-3 on-board clocks is better; (2) Schemes 1 and 2 present the similar fitting residuals, while Scheme 3 shows the minimum values with all periodic terms being considered; (3) to analyze the characteristic of sparse modelling, two time periods, 24 h and 18-24 h, are tested by the average daily RMS values. It is found that the model residuals decreased with the increase of time period beyond 6 h, which is led by time factors in weight matrix of equation (7). e RMS values are smaller for 18-24 h period based on sparse modelling, compared with Schemes 1 and 2. (4) From the accuracy of clock bias prediction in Figure 5, Schemes 1 and 2 almost output the same results for BDS-2 as the previous studies about BDS-2 satellite clock bias prediction. Moreover, it is revealed that Scheme 3 is more accurate than the traditional methods (Schemes 1 and 2), while Scheme 4 shows slightly improved results than Scheme 3 (except C13). (5) For the results of BDS-3 in Figure 5, the most obvious difference is Scheme 3 becoming worst as all terms are estimated by the one-step solution. However, after the sparse modelling (Scheme 4), the accuracy of predicted clock bias is further improved with Schemes 1 and 2.
Furthermore, regardless of the prediction accuracy of clock biases, the repeatability of the model coefficients for the trend terms can indirectly reflect the performance of the BDS clock prediction model. In Table 4, the STD of BDS-2 model coefficients is presented as an example to discuss the differences of four schemes. In general, the prediction model of Scheme 4 with the sparse modelling gives more optimal results. It should be noted that the prediction accuracy is about 6.0 ns and 5.0 ns for BDS-2 and BDS-3 satellite clock biases within 18-hour period, respectively. Further improvement should be introduced in the process of model coefficients.
erefore, the second group of experiments includes two schemes for the verification of the introduction of the intersatellite correlation in the estimation of model coefficients.
Scheme 5: referring Scheme 1, the intersatellite correlations are inserted into the solution of model coefficients. More details of experiments and results were typically described in our previous research studies [6,24]. Scheme 6: similarly, based on Scheme 4, the intersatellite correlation is obtained and used to estimate the model coefficients for each satellite.
In experiments, before analyzing the accuracy of predicted clock biases, the values of intersatellite correlation factor are firstly discussed. e intersatellite correlation on DOY 41 is shown in Table 5, from which the mean values are about 0.31-0.45. As mentioned in our previous research studies, the intersatellite correlation values represent the low-correlation degree (0.3-0.5). Moreover, the results of intersatellite correlations between BDS-2 and BDS-3 are slightly smaller than the intrasatellite correlation of BDS-2 and BDS-3 satellites. In this study, equipped with more stable on-board clock of BDS-3, we take the correlation of BDS satellites into consideration to improve the solution of clock bias prediction.
In addition, the ten-day accuracy of BDS-predicted clock biases for C14 and C24 is listed in Figure 6 (average daily values) as an example to explain the variation of four schemes. Moreover, to highlight the differences of results for different schemes, Table 6 illustrates the average daily RMS values of one-month experiments and their improvements after the introduction of intersatellite correlation. According to results shown in Figure 6 and Table 6, the intersatellite correlation can impose a significant improvement on the BDS satellite clock bias prediction, which can reach 32.9% and 16.9% for BDS-2 and BDS-3 satellites, compared to traditional method (Scheme 1). Moreover, for Scheme 4, the improvement of 27.2% and 28.6% can be realized for BDS-2 and BDS-3 after using the intersatellite correlation.
Overall, getting rid of the trend and periodic terms, the model residuals are still contained in the clock biases, which can reduce the accuracy of clock bias prediction. erefore, a new strategy which takes the spatial and temporal variations of model residuals into consideration is proposed to further improve the prediction model. e third group of experiments is designed to verify the prediction model of ultrarapid clock bias by using an empirical variogram function to update the weight matrix in the estimation of model coefficients. Before conducting experiments of prediction, the quantitative measure of the temporal correlation in residuals for each satellite is presented in terms of Scheme 4. In Figure 7, the dashed lines (different colors for different satellites) represent the experimental variogram with onemonth residuals. According to the experimental variogram, an empirical model will be assumed as the parametric variogram model. We select the spherical model to describe variogram as equations (17) and (18).
From the above analyses, the model variogram can be summarized as equations (17) and (18) for BDS-2 and BDS-3 satellites, respectively.
Scheme 7: according to BDS-2 and BDS-3 satellites' integrated ultrarapid clock biases, the clock products are predicted by Scheme 1, in which the intersatellite correlations are considered. Scheme 8: according to BDS-2 and BDS-3 integrated ultrarapid clock biases, the clock products are predicted by Scheme 2, in which the intersatellite correlations are considered. Scheme 9: according to BDS-2 and BDS-3 integrated ultrarapid clock products, the clock biases are predicted from Scheme 6, in which the time period is set as 18 h. Moreover, the residuals of prediction model are taken into BPNN to predict the nonlinear parts in clock biases. Scheme 10: similar to Scheme 9, the residuals of prediction model are processed by the gray model. e prediction of clock bias will be constructed by the combination of two parts. Scheme 11: based on Scheme 6, the strategy for processing model residuals discussed in Wang et al. [24] is adopted, namely, the PLS plus BPNN, to optimize the predicted residuals. Scheme 12: based on Scheme 6, the empirical variogram function is used to adjust the weight matrix of equation (4), in which the nondiagonal elements of P are replaced by the covariance of residual series based on equations (17) and (18).
e residuals of prediction model are also processed by PLS plus BPNN.
To fully analyze the prediction model, the ultrarapid clock biases acquired by the integrated estimation are typically compared with 12 schemes. Similarly, the rapid clock bias products of WHU are also assumed as references. To assess the model accuracy based on different experiments, three schemes are presented in Figure 8 with the model residuals of C14 and C24 on DOY 41, 2019. It can be easily concluded that the prediction model improved with sparse modelling, intersatellite correlation, and PLS plus BPNN for residuals (Scheme 11) can obtain a desired clock bias model as the minimum model residuals.
Moreover, details about the accuracy of clock bias prediction for Schemes 7-12 are shown in Figures 9 and 10. Because the ultrarapid clock products are updated with a 3 h latency and 6 h interval period, the accuracy of 18 h in one day is analyzed as discussed in Hu et al. [6].
Besides, the details of one-day accuracy, the average daily RMS values within 30 days for Schemes 7-9 and Schemes 10-12, and its corresponding improvements are presented in Tables 7 and 8, respectively. From Figures 8-10 and Tables 7  and 8, the results can be summarized as follows: (1) compared with the traditional methods (Schemes 1 and 2), the use of intersatellite correlation can slightly reduce the model residuals by BDS-2/BDS-3 combined estimation; (2) it is suggested that the period terms selected by the widely used methods caused the decrease of prediction accuracy by −4.3% and −21.1% for BDS-2 and BDS-3, respectively; (3) according to the processing of the gray model, the accuracy of BDS-2 and BDS-3 clock bias prediction can be improved by 23.3% and 16.9%, respectively. However, it should be mentioned that Scheme 9 is optimized from the traditional method, which needs to be further compensated with the model errors. erefore, the BPNN and the PLS plus BPNN are tested in Schemes 10 and 11.
In Table 8, we compare different methods in aspect of model residuals based on our proposed strategy (Scheme 6).
e PLS plus BPNN can achieve slight improvements, 0.5% and 6.2% for BDS-2 and BDS-3, than traditional BPNN only on modelling residuals. Moreover, due to the existence of spatial and temporal correlation in residuals series, the variogram function is modelled to update the weight matrix in estimation of coefficients. Compared with Scheme 11, the predicted clock biases can be improved by 8.0% and 11.1% for BDS-2 and BDS-3, respectively. It should be noted that the combination of PLS and BPNN is also utilized to predict the residual series.

Discussion
In this contribution, we optimize the strategy for BDS-2 and BDS-3 integrated ultrarapid clock bias prediction based on our previous research studies. One-month experiments with    12 schemes are designed to fully verify the improved methods.
According to the results of three groups of experiments, conclusions can be summarized and analyzed as follows.
Firstly, compared with the traditional method (a twostep solution), the one-step strategy is used to estimate the coefficients of BDS clock bias models. Moreover, the model is selected by sparse modelling and solved by Lasso algorithm. As shown in Schemes 1-4, results of one-step solution can output the minimum model residuals, while the sparse modelling based on one-step solution acquires the most accurate prediction of BDS clock biases. Because of the low quality of BDS clock bias (estimated by B1I and B3I observation), random noise and wrong periodic terms are included in the estimation of model coefficients. erefore, on the one hand, the gross detection and repair and denoising algorithm, such as Tikhonov regularization, should be inserted into the processing of clock bias prediction; on the other hand, the adaptive model optimization like machine learning is needed to select the optimal numbers, types, and arc length for clock bias prediction. In proposed strategy, the sparse modelling gives a better result than that of widely used models. In addition, compared with Schemes 1 and 2, the results of BDS-2 are almost unchanged, while it is more significant for BDS-3. e periodic terms and relationship in BDS-3 satellites' clock and orbit should be paid more attention in the following research studies.
Secondly, before conducting the clock bias prediction, we analyzed the frequency stability of BDS-2 and BDS-3 one-month clock products. It is indicted that BDS-3 satellite on-board clocks exhibit better performances than that of BDS-3. erefore, the advantage of BDS-3 information can be utilized to enhance the solution of combined estimation for BDS-2/BDS-3. In Schemes 5 and 6, the intersatellite correlations are tested based on Schemes 1 and 4, respectively. It is found that the prediction of clock biases is significantly improved, especially for the traditional method (Scheme 1). In fact, the combined estimation of BDS-2/BDS-3 clock biases is intended to overcome the limited observations (27 tracking stations with B1I and B3I signals) and strengthen the solution of model coefficients. Moreover, the intersatellite correlation about BDS-2 and BDS-3 integrated clock bias prediction was already fully analyzed in our previous research studies [5,6,24]. Furthermore, it must be noted that all experiments depended on the B1I and B3I observations. In aspect of different observations for the estimation of BDS-2 and BDS-3 integrated clock bias series, the differences of signals and satellite types must be fully analyzed, such as intersatellite bias.
irdly, after improving the clock bias prediction model by one-step strategy with sparse modelling and intersatellite correlation, the accuracy of prediction can be increased by more than 30.0% for BDS-2 and BDS-3 integrated clock biases. But the model residuals cannot be avoided in our strategy. Some prediction methods for model residuals are widely used to improve the prediction models, such as the gray model, BPNN, and PLS plus BPNN model. For the temporal correlation, the residuals are modelled by variogram to update the weight matrix in the estimation of model coefficients. According to the third group of experiments, it is shown that the improvements are almost less than 10.0% after the introduction of variogram. In the following work, a more accurate variogram model is a main point to further analyze, such as exponential, Gaussian, power, and nugget models.
Consequently, as discussed above, the proposed strategy can effectively optimize the BDS-2 and BDS-3 integrated ultrarapid clock bias prediction. However, because of the low accuracy of clock bias series or estimation models, the prediction model cannot fully use the advantages in BDS-3 clock bias series. erefore, all results in this research are initial tests, which should be further verified with more data and tested in calculating BDS ordinary clock products, especially for BDS-3 satellites. Additionally, due to different types of BDS-3 satellite manufacturers, CASC (China Aerospace Science and Technology Corporation) and SECM (Shanghai Engineering Center for Microsatellites), the differences of orbit dynamics are inevitable. Moreover, the clock biases are strongly correlated with the orbit radial component. In the following studies, we will take the satellite types, orbit errors, and the intersatellite-type correlations (GEO, IGSO, and MEO) into consideration to further optimize the BDS clock bias prediction products.

Conclusions
BDS ultrarapid clock biases are the key inputs of real-time or near real-time services for BDS rapid and high-accuracy services. However, it is found that the predicted parts of BDS ultrarapid satellites' clock biases have some defects, which inevitably reduces the accuracy of clock bias prediction. In this study, the traditional strategy is improved with three methods to optimize the prediction models.
Due to the inaccurate two-step solution of clock bias model coefficients, the one-step strategy is proposed to estimate the models (trend and periodic terms). e sparse modelling of machine learning is used to improve the estimation of models. According to the results of clock bias prediction based on rapid clock bias products, the one-step strategy can slightly increase the accuracy of the predicted clock biases. Moreover, considering the differences in BDS satellites (BDS-2 and BDS-3) and its characteristics of onboard clock biases, the intersatellite correlations are assumed as an indirect method to enhance the solution of model coefficients. It is concluded that the accuracy of BDS-2-and BDS-3-predicted clock biases can be improved by 27.2% and 28.6% by one-step strategy, respectively. In addition, the model residuals should be further analyzed in modelling clock bias prediction. We output the residuals of prediction model and construct the empirical variogram function to extract the temporal correlation in model residuals, which is taken into updated weight matrix in estimating model coefficients. It is suggested that the accuracy of clock bias model can be improved by 8.0% and 11.1%, compared with the strategy proposed in recent research studies [6,24]. In summary, according to the results of predicted clock biases, it is indicted that the stability of BDS-3 on-board clocks is more optimal compared with BDS-2, which can be used to strengthen the solution of clock bias prediction model; the one-step estimation of the clock bias model by sparse modelling can slightly increase the accuracy of prediction results; both BDS-2-and BDS-3-predicted clock biases benefited each other by inserting the intersatellite correlations into the weight matrix, in which the accuracy of 18-hour period with one-step strategy can be improved by 28.6% and 27.2% for BDS-2 and BDS-3, respectively; after the introduction of the variogram model in updating the weight matrix, the clock bias prediction model is further corrected by 8.0% and 11.1% for BDS-2 and BDS-3. erefore, improved strategies for BDS ultrarapid satellites' clock bias prediction using BDS-2 and BDS-3 integrated processing are meaningful for the current BDS ultrarapid satellites' clock bias prediction products.
Overall, the proposed strategy for BDS-2 and BDS-3 integrated ultrarapid clock bias prediction can be used to improve the accuracy for iGMAS ACs. However, the improved strategies are just tested based on one-month experiments.
e availability and accuracy of the proposed methods for BDS-2 and BDS-3 integrated clock products should be further investigated.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
Zhongyuan Wang and Chao Hu conceived and defined the research scheme. Chao Hu verified the feasibility of the method and implemented the software algorithm. Zhongyuan Wang checked the data processing results and manuscript. Chao Hu wrote the manuscript.