A Practical Adaptive Clock O ﬀ set Prediction Model for the Beidou-2 System

: The predicted navigation satellite clock o ﬀ sets are crucial to support real-time global navigation satellite system (GNSS) precise positioning applications, especially for those applications di ﬃ cult to access the real-time data stream, such as the low earth orbit (LEO) autonomous precise orbit determination. Currently, the clock prediction for the Chinese BeiDou system is still challenging to meet the precise positioning requirement. The onboard clocks of BeiDou satellites are provided by di ﬀ erent manufacturers, and the clocks’ switch events are more frequent. Considering the satellite-speciﬁed and temporal variation of the BeiDou clocks characteristics, we intend to use an adaptive model for BeiDou clock prediction. During clock prediction, we identify di ﬀ erent models for BeiDou clocks’ characteristics, and then address the optimal model with a cross-validation procedure. The model achieving the minimum variance in the cross-validation procedure is used for the ﬁnal clock prediction. We compared the prediction results of our method with two well-recognized BeiDou ultra-rapid clock products, named GBU-P and ISU-P, respectively. The comparison results indicate that the adaptive model achieves about 1-ns precision for 3-h prediction, which corresponds to 47.3% and 32.1% precision improvement compared to the GBU-P and ISU-P products, respectively. The e ﬃ ciency of the predicted clocks is further validated with the precise point positioning (PPP) data processing. The results indicate that the static PPP solution precision is improved by 21.6%–30.0% compared to the current predicted clock product. The precision improvement in kinematic PPP is even more signiﬁcant, which reaches 46.7%–53.9% with respect to these GBU-P and ISU-P products. Therefore, the proposed adaptive model is a practical and an e ﬃ cient way to improve the BeiDou clock prediction.


Introduction
Precise point positioning (PPP) is a popular satellite positioning technique providing a global scale, centimeter-level accuracy user position, thus becoming increasingly popular in GNSS communities [1][2][3].PPP data processing requires a precise handling of various error sources, including the satellite orbit and the satellite clock offset.The post-processing PPP relies upon the precise orbit and clock products computed from the globally-distributed GNSS tracking network on a daily basis.Meanwhile, the real-time precise point positioning (RTPPP) service has been widely used in natural hazard early warning systems, landslide monitoring and atmospheric water vapor sensing [4], which is supported by the real-time or predicted ephemeris and clock products [5][6][7].The real-time products can be estimated and distributed to users with real-time data stream access [8,9].For many cases, short-term prediction of the orbit and clock with the historical data is more practical and computationally efficient for RTPPP users.
For example, the low earth orbiters' (LEO) autonomous precise orbit determination (POD) can only use the predicted products due to a lack of stable data connection between the satellite and the ground facilities.The real-time POD capacity will improve the LEO satellites' performance, such as cutting the radio occultation data processing latency [10,11], reducing the dependency of the ground control point for the remote sensing satellite imagery data processing, and supporting the LEO navigation augmentation applications [12][13][14].The Luojia-1A satellite has demonstrated the potential benefit of navigation signal augmentation from the LEO satellite [15,16].
Supporting the real-time PPP by the predicted clocks attracts much attention in the GNSS communities.The international GNSS service (IGS) released its clock prediction product (known as the ultra-rapid products) to support the RTPPP applications.Currently, 24-h orbit prediction achieves a few centimeters to one decimeter precision [17,18], which is good enough for most precise positioning applications, whereas the predicted clock products only obtained 3 ns precision, which is still challenging for real-time centimeter-level PPP.The 3 ns clock bias corresponds to about 90 cm ranging error in line-of-sight direction, which will deteriorate the positioning precision.
The satellite clock prediction model for the GPS system has been well studied with long-term observations [19,20], since the precision of the clock offset prediction is governed by the characteristics of the onboard clocks.The IGS ultra-rapid clock products (IGU-P) use a polynomial model with a sinusoid term for GPS clock offset prediction, which achieves about 1-4 ns precision for a 24-h prediction [21].A number of clock prediction methods have been proposed to further improve the IGU-P products, such as the polynomial model with more periodical terms, the autoregressive integrated moving average (ARIMA), the gray model, the wavelet neural network (WNN) model and the fusion-model [19,[21][22][23][24][25][26].Most of these models focus upon handling the stochastic terms in the clock.Currently, centimeter-to decimeter-level PPP results have been achieved with different GPS clock prediction products [6,27].Hence, the predicted clock offset is a practical way of supporting such real-time applications.The BeiDou navigation system (BDS), an emerging global navigation satellite system deployed by China, plans to provide worldwide positioning and navigation services by 2020 [28,29].At the current stage, the prediction algorithm for the BeiDou clocks is inherited from GPS clock prediction, but the prediction precision is currently still not satisfied [30][31][32].Presently, the 24-h prediction of the BeiDou clock offset only achieves about 6-7 ns precision, which still cannot meet the need for prediction accuracy.There are a few researchers attempting to improve the precision of BeiDou ultra-rapid products by considering the non-linear effects in BeiDou clock time series.Huang et al. [33] attempted to mitigate the non-linear effects in clock offset by the backward propagation neural network (BPNN) model, and achieved promising results.He et al. [34] consider the non-linear effects by the generalized regression neural network (GRNN) and the prediction results are enhanced.Zhou, et al. uses the ARIMA method to improve short-term BeiDou clock prediction [35].However, it still far from fully understands the characteristics of BeiDou clocks.
According to extensive data analysis, we find that the key to improving BeiDou clock prediction is not the non-linear effects, since there is no substantial difference between the GPS and BeiDou clock in terms of stability.According to [31], the Allan variance for GPS clocks and BeiDou clocks are at the same level.The key issue of BeiDou clock prediction is the satellite-specified clock characteristics.Unlike the GPS system, the BeiDou employs a hybrid constellation, which consists of different satellite types, and the onboard clocks also come from different manufacturers.Therefore, the behaviors of the BeiDou clocks are highly satellite-specified.Each BeiDou satellite is equipped with several atomic clocks, and these clocks also have different in-orbit performance.In addition, the BeiDou clock presents frequency clock steering or phase adjustment in operation, which causes many jumps and turning points [36] in the clock time series.Another BeiDou clock characteristic is that the optimal clock model for a particular satellite may change over time.These satellite-specified and temporal variations of the clock characteristics have a dramatic impact upon clock performance, so it is unreasonable to describe all of the BeiDou clocks with one unified and static model.Instead, building a 'smart' model for BeiDou clocks is the major challenge of the clock prediction issue.In order to solve this problem, we propose an adaptive model for BeiDou clock offset prediction.
In this study, we built different models for different BeiDou clock types and then employed a cross-validation algorithm to select the optimal model for each clock.In this way, the satellite-specified and temporal variation of the satellite clock behaviors is properly considered in the clock model.The optimal prediction model is always driven by the observed clock data, so that the model can automatically adapt different clock types and even clock behavior changes.
The remaining sections of the paper are organized as follows: The selection of the clock models for the BeiDou satellites is discussed in Section 2, and the procedure of the adaptive model is described in Section 3. The performance of the adaptive clock prediction model is compared with the GBU-P product and the ISU-P products in Section 4, and the benefit of improving the clock prediction error in the precise point positioning (PPP) is investigated in Section 5.The conclusion comes in Section 6.

Modeling the BeiDou Clocks
The current BeiDou constellation consists of Geosynchronous Orbit (GEO) satellites, Inclined Geosynchronous Satellite Orbit (IGSO) satellites, and Medium Earth Orbit (MEO) satellites.All of these BeiDou-2 satellites are equipped with the rubidium atomic frequency standard, while the onboard BDS-2 satellites and corresponding clock types are listed in Table 1.Since the Beidou-3 satellites are still in the evaluation phase, our discussion in this study is limited to BDS-2 satellites clocks.

Satellites Type
Clock Type PRN The most popular clock model is the quadratic polynomial combining with periodic terms, which has been used in ultra-rapid clock products for BDS-2 clocks.The model can be expressed as follows [22,33]: where Y(t) is the clock offset phase observation, t 0 is the reference epoch, t i represents the observation epoch; a 0 , a 1 and a 2 are the clock offset, clock velocity and clock drift, respectively; and l is the number of periodic terms Λ i , f r are the amplitude and frequency of periodic terms, respectively, ϕ r is the phase, and ∆ i represents the residuals of the prediction model.Generally, the non-periodical terms are caused by the intrinsic property of the clocks.The clock velocity is caused by the constant bias in the clock frequency.The clock drift is also known as the frequency drift, which is caused by the linear variation of the clock frequency.The frequency drift depends upon the clock age, the space environment, etc.If the clock frequency drift is negligible over 24-h, the clock can be modeled as the linear model (2), otherwise, the quadratic model (3) should be used.Since the onboard clocks are all of a high-performance frequency standard, higher order, non-linear effects in clock frequency are usually negligible.
The periodical terms are caused by the general relativistic effect, or coupled with the orbit error.
From the long-term analysis, there is a diurnal or semi-diurnal period presence in the clock offset time series.The procedure of the satellite clock prediction includes a model fitting procedure and a prediction procedure.In the model fitting procedure, the clock model is often described with an empirical clock model, and this empirical model is normally a reduced version of the model described in Equation (1).In order to achieve good prediction accuracy, the model coefficient is fitted with a period of historical clock data.Then the fitted model is used to predict the clock of the next period.Some methods also involve a time series model or a neural network model to handle the small, un-modeled errors.Currently, the BeiDou clocks are often predicted by the quadratic model with two periodical terms.The model is fitted with 24-h clock data.
In this section, we use the final BeiDou clock products released by Wuhan University (known as WHM) for the clock prediction model.About 5-week precise clock data (DOY 203-237, 2018) are used to reveal the characteristics of the BeiDou clocks.The clock prediction model is discussed from three aspects: The optimal fitting period length, the optimal polynomial order, and the optimal periodical terms.

The Optimal Fitting Period Length
How to choose the fitting period is the key issue in the model fitting step.A too short period causes an inaccurate model coefficient, and a too long period may introduce nuisance non-linear biases.Currently, the 24-h period is the most popular choice in BeiDou clock prediction, but whether it is the optimal choice has not been systematically examined.The root mean squares (RMS) method is used to analyze the prediction result, which can be expressed as: where y i and ŷi are the predicted and the observed clock offset at the i th epoch, respectively.In this study, we tried to fit the BeiDou clocks with the linear model and the quadratic model with different fitting periods, and the results are presented in Figure 1.The figure presents the root mean squares (RMS) of the 24-h clock prediction error with different clock models and different fitting periods, using the five weeks of clock data.In this study, the fitting period varies from 12-h to 72-h.The figure shows that a longer fitting period means a larger fitting error for the linear model.The fitting error increases dramatically as the fitting period increases, which means that the linear model should not use a long fitting period.In contrast, the fitting error of the quadratic model presents a 'U' shape, and the optimal fitting period is 48-h.Both a longer and shorter fitting period will decrease the fitting precision.A longer fitting period decreases the prediction error, which may be caused by the accumulated non-linear clock frequency drift effect.A shorter fitting period may cause a poor prediction precision due to a poor observability of the frequency drift term.The figure also indicates that the linear model outperforms the quadratic model with a fitting period shorter than 36 h, whereas the quadratic model outperforms the linear model when the fitting period is longer than 36-h.Therefore, the selection of the fitting period depends upon the underlying clock model.Different clock models correspond to the different optimal fitting period length.In this study, a 24-h fitting period for the linear model, and a 48-h fitting period for the quadratic model, is selected, as their optimal fitting period length, respectively, since they minimize the prediction error.The optimal fitting period length depends upon the specified clock prediction model.

The Optimal Polynomial Order
The second issue of clock modeling is the optimal polynomial order for BeiDou clock prediction.The polynomial order of the satellite clocks depends upon the frequency drifting characteristics, which is considered as an intrinsic property of the clock.The ultra-rapid clock prediction products (IGU-P) employ a linear model or a quadratic model, according to the satellite and atomic clock type for GPS.Ray suggests a quadratic plus sinusoid model for the GPS Block II/IIA Rubidium clocks and a linear plus sinusoid model for the remaining clocks [21].Huang et al. [33] employ a quadratic model for BeiDou clock prediction empirically.Hence the optimal polynomial order for BeiDou clock prediction is worth investigating.We compared the clock prediction error with both the linear model and the quadratic model, and the results are presented in Figure 2. Figure 2 indicates that none of these two models can fit all of the BeiDou satellite clocks.For some satellites, the linear model outperforms the quadratic model, whereas the quadratic model outperforms for the rest of the satellites.The outperformed satellites for the two models are nearly fifty-fifty, so there is no method that has significant advantages compared to the other.The difference of the mean prediction RMS reaches 2 ns with different clock models, which means that the clock prediction error can be further mitigated by optimally selecting the underlying clock model instead of using a particular model to all satellites.

The Optimal Polynomial Order
The second issue of clock modeling is the optimal polynomial order for BeiDou clock prediction.The polynomial order of the satellite clocks depends upon the frequency drifting characteristics, which is considered as an intrinsic property of the clock.The ultra-rapid clock prediction products (IGU-P) employ a linear model or a quadratic model, according to the satellite and atomic clock type for GPS.Ray suggests a quadratic plus sinusoid model for the GPS Block II/IIA Rubidium clocks and a linear plus sinusoid model for the remaining clocks [21].Huang et al. [33] employ a quadratic model for BeiDou clock prediction empirically.Hence the optimal polynomial order for BeiDou clock prediction is worth investigating.We compared the clock prediction error with both the linear model and the quadratic model, and the results are presented in Figure 2.

The Optimal Polynomial Order
The second issue of clock modeling is the optimal polynomial order for BeiDou clock prediction.The polynomial order of the satellite clocks depends upon the frequency drifting characteristics, which is considered as an intrinsic property of the clock.The ultra-rapid clock prediction products (IGU-P) employ a linear model or a quadratic model, according to the satellite and atomic clock type for GPS.Ray suggests a quadratic plus sinusoid model for the GPS Block II/IIA Rubidium clocks and a linear plus sinusoid model for the remaining clocks [21].Huang et al. [33] employ a quadratic model for BeiDou clock prediction empirically.Hence the optimal polynomial order for BeiDou clock prediction is worth investigating.We compared the clock prediction error with both the linear model and the quadratic model, and the results are presented in Figure 2. Figure 2 indicates that none of these two models can fit all of the BeiDou satellite clocks.For some satellites, the linear model outperforms the quadratic model, whereas the quadratic model outperforms for the rest of the satellites.The outperformed satellites for the two models are nearly fifty-fifty, so there is no method that has significant advantages compared to the other.The difference of the mean prediction RMS reaches 2 ns with different clock models, which means that the clock prediction error can be further mitigated by optimally selecting the underlying clock model instead of using a particular model to all satellites.Figure 2 indicates that none of these two models can fit all of the BeiDou satellite clocks.For some satellites, the linear model outperforms the quadratic model, whereas the quadratic model outperforms for the rest of the satellites.The outperformed satellites for the two models are nearly fifty-fifty, so there is no method that has significant advantages compared to the other.The difference of the mean prediction RMS reaches 2 ns with different clock models, which means that the clock prediction error can be further mitigated by optimally selecting the underlying clock model instead of using a particular model to all satellites.
A closer analysis of the clock behavior indicates that the optimal clock model for a particular satellite may also change over time.Two examples are presented in Figure 3.We selected C06 and C08 satellites and analyzed the clock prediction error with different models at different times.The figure indicates that the optimal clock model may change over time.The panel (a) and (b) shows that the optimal prediction model for C06 is the linear model at DOY 121, but it switched to the quadratic model at DOY A similar phenomenon also presents in the C08 satellite.Although the clock switch and steering information of C06 and C08 is not publicly accessible, it can conclude that the optimal prediction model switching over time is not an occasional phenomenon.Hence, it is impossible to optimally select a static prediction model for all of the BeiDou satellites, since the optimal prediction model may change over time.Therefore, the optimal clock prediction model should be determined according to the clock data, rather than using the empirically fixed model.A closer analysis of the clock behavior indicates that the optimal clock model for a particular satellite may also change over time.Two examples are presented in Figure 3.We selected C06 and C08 satellites and analyzed the clock prediction error with different models at different times.The figure indicates that the optimal clock model may change over time.The panel (a) and (b) shows that the optimal prediction model for C06 is the linear model at DOY 121, but it switched to the quadratic model at DOY A similar phenomenon also presents in the C08 satellite.Although the clock switch and steering information of C06 and C08 is not publicly accessible, it can conclude that the optimal prediction model switching over time is not an occasional phenomenon.Hence, it is impossible to optimally select a static prediction model for all of the BeiDou satellites, since the optimal prediction model may change over time.Therefore, the optimal clock prediction model should be determined according to the clock data, rather than using the empirically fixed model.

The Optimal Periodic Terms
Introducing the periodical term can significantly reduce the GPS clock prediction error [21,38], and thus it has been applied in BeiDou clock modeling [33].The characteristics of the GPS clocks have been analyzed, and the periodic terms were detected from the long-term clock offset time series [37,38].Wang et al. [30] investigated the periodic terms in BDS satellites clocks variations by spectrum analysis, indicating the evident 1 cycle per revolution (CPR) and 2 CPR periodicities.In order to analyze how the periodic terms impact prediction accuracy, we performed a long-term spectrum analysis with one-year BeiDou clock final products and addressed their first two major periods.The analysis indicates that different BeiDou satellites have different periods, and their periods are listed in Table 2.

The Optimal Periodic Terms
Introducing the periodical term can significantly reduce the GPS clock prediction error [21,37], and thus it has been applied in BeiDou clock modeling [33].The characteristics of the GPS clocks have been analyzed, and the periodic terms were detected from the long-term clock offset time series [37,38].Wang et al. [30] investigated the periodic terms in BDS satellites clocks variations by spectrum analysis, indicating the evident 1 cycle per revolution (CPR) and 2 CPR periodicities.In order to analyze how the periodic terms impact prediction accuracy, we performed a long-term spectrum analysis with one-year BeiDou clock final products and addressed their first two major periods.The analysis indicates that different BeiDou satellites have different periods, and their periods are listed in Table 2.In order to investigate the impact of periodic terms, we compared the clock fitting accuracy and the 24-h clock prediction accuracy with one periodic term, two periodic terms and without a periodic term model.The RMS of the BeiDou clock fitting residual with different periodic terms are presented in Figure 4.The two panels show the fitting error of the linear model and the quadratic model, combined with different periodic terms, respectively.The figure indicates that adding periodic terms indeed improves the clock fitting accuracy.Adding two periodic terms dramatically improves the fitting accuracy of the GEO satellites.Adding one periodic term significantly improves these IGSO satellites.The MEO also benefit from the additional periodic terms, but less significantly.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 25 In order to investigate the impact of periodic terms, we compared the clock fitting accuracy and the 24-h clock prediction accuracy with one periodic term, two periodic terms and without a periodic term model.The RMS of the BeiDou clock fitting residual with different periodic terms are presented in Figure 4.The two panels show the fitting error of the linear model and the quadratic model, combined with different periodic terms, respectively.The figure indicates that adding periodic terms indeed improves the clock fitting accuracy.Adding two periodic terms dramatically improves the fitting accuracy of the GEO satellites.Adding one periodic term significantly improves these IGSO satellites.The MEO also benefit from the additional periodic terms, but less significantly.
(a) (b) Adding periodic terms indeed improves the clock fitting accuracy, but by not necessarily also improving the clock prediction.We further investigated the impact of periodic terms on the clock prediction accuracy.The accuracy improvement of adding periodic terms, compared to this no periodic terms case, is presented in Figure 5.The polynomial model is used as the reference, and the positive value means the prediction precision improved by adding periodic terms.The figure shows that the impact of the periodic terms on clock prediction is not the same as its impact on clock fitting.The periodic terms have the most significant impact on GEO clock prediction, and have little impact on the MEO clock prediction.For the linear model, adding two periodic terms achieves a 52 percent improvement at maximum, but it has a negative impact on the quadratic model.Introducing one or two periodic terms has a less significant impact on the IGSO satellites.The impact of periodic terms on MEO satellites is almost negligible.The periodic terms are only obvious in long-term time series analysis.For one or two days, the periodic terms are not obvious.Adding periodic terms might lead to over fitting, so that it produces a negative impact on the clock prediction.Overall, whether the periodic terms improve clock prediction depends on the polynomial terms and the satellite type.The linear model combined with the two periodic terms model can further improve the clock prediction precision than the linear model, but the quadratic model does not present obvious precision improvement by adding additional period terms.Adding periodic terms indeed improves the clock fitting accuracy, but by not necessarily also improving the clock prediction.We further investigated the impact of periodic terms on the clock prediction accuracy.The accuracy improvement of adding periodic terms, compared to this no periodic terms case, is presented in Figure 5.The polynomial model is used as the reference, and the positive value means the prediction precision improved by adding periodic terms.The figure shows that the impact of the periodic terms on clock prediction is not the same as its impact on clock fitting.The periodic terms have the most significant impact on GEO clock prediction, and have little impact on the MEO clock prediction.For the linear model, adding two periodic terms achieves a 52 percent improvement at maximum, but it has a negative impact on the quadratic model.Introducing one or two periodic terms has a less significant impact on the IGSO satellites.The impact of periodic terms on MEO satellites is almost negligible.The periodic terms are only obvious in long-term time series analysis.For one or two days, the periodic terms are not obvious.Adding periodic terms might lead to over fitting, so that it produces a negative impact on the clock prediction.Overall, whether the periodic terms improve clock prediction depends on the polynomial terms and the satellite type.The linear model combined with the two periodic terms model can further improve the clock prediction precision than the linear model, but the quadratic model does not present obvious precision improvement by adding additional period terms.Overall, we addressed two optimal models for the BeiDou clocks.The analysis results indicate that the optimal clock model should be either one of the models listed in Table 3, but which is the optimal one, is satellite dependent.The clock analysis indicates that the BeiDou clocks behave in different manners, so that it is difficult to optimally predict with one uniform model.Determining the clock prediction model according to the satellite type may be sound for the GPS satellites, but BeiDou clocks seem not to follow this rule.The optimal clock model for BeiDou is satellite specified, and may change over time.The exhaustive clock models for different BeiDou status have been built, but how to reasonably select the model has not been addressed.
In this study, we employ an adaptive model to reasonably select the optimal clock model for each BeiDou satellite.The adaptive model employs a cross-validation procedure to discriminate the optimal clock model.During the cross-validation procedure, a small part of the observation period is reserved for the validation purpose.The model achieves minimum variance in the cross-validation section, and is thus considered as the optimal model for the particular satellite.This cross-validation procedure is based upon the assumption that the clock frequency drift property is the same for the validation section and the prediction period.This assumption generally holds, except when the clock is steered or switched.If such clock events happen, the satellite clock becomes unpredictable.
The whole adaptive model consists of three steps: Data preprocessing, cross-validation and clock prediction.The flowchart of the proposed adaptive model is illustrated in Figure 6.The outliers and daily boundary problems in the raw precise clock product are handled in the preprocessing stage.Then the optimal clock model is selected by the cross-validation procedure, after which the clock is predicted with the model identified by the cross-validation procedure.The cross-validation stage can be further divided into three sub-steps: Data segmentation, multi-model clock fitting and model selection, which are illustrated in Figure 6.Overall, we addressed two optimal models for the BeiDou clocks.The analysis results indicate that the optimal clock model should be either one of the models listed in Table 3, but which is the optimal one, is satellite dependent.

Adaptive Model for BDS Clocks Offsets
The clock analysis indicates that the BeiDou clocks behave in different manners, so that it is difficult to optimally predict with one uniform model.Determining the clock prediction model according to the satellite type may be sound for the GPS satellites, but BeiDou clocks seem not to follow this rule.The optimal clock model for BeiDou is satellite specified, and may change over time.The exhaustive clock models for different BeiDou status have been built, but how to reasonably select the model has not been addressed.
In this study, we employ an adaptive model to reasonably select the optimal clock model for each BeiDou satellite.The adaptive model employs a cross-validation procedure to discriminate the optimal clock model.During the cross-validation procedure, a small part of the observation period is reserved for the validation purpose.The model achieves minimum variance in the cross-validation section, and is thus considered as the optimal model for the particular satellite.This cross-validation procedure is based upon the assumption that the clock frequency drift property is the same for the validation section and the prediction period.This assumption generally holds, except when the clock is steered or switched.If such clock events happen, the satellite clock becomes unpredictable.
The whole adaptive model consists of three steps: Data preprocessing, cross-validation and clock prediction.The flowchart of the proposed adaptive model is illustrated in Figure 6.The outliers and daily boundary problems in the raw precise clock product are handled in the preprocessing stage.Then the optimal clock model is selected by the cross-validation procedure, after which the clock is predicted with the model identified by the cross-validation procedure.The cross-validation stage can be further divided into three sub-steps: Data segmentation, multi-model clock fitting and model selection, which are illustrated in Figure 6.

Data Preprocessing
Data preprocessing is a vital procedure for the clock prediction, which removes the outliers and handles the daily boundaries to produce clean data for the following data processing.The data with gross error or outliers cannot reflect the nature of clocks, thus should not be involved in clock modeling and prediction.The popular median absolute deviation (MAD) method is used to detect the gross error.This method has a high breakdown point and is robust to multiple outliers, which can be expressed as follows [39]: where the  presents the clock frequency sequence,  is the median of the clock frequency sequence, m = median{ }.If the clock frequency satisfies  >  +  *  or  <  −  * , the clock offset will be regarded as a gross error.In this case,  is set as 3, and the detected clock offset gross error is excluded in the data processing.
Another problem in preprocessing is the daily boundary discontinuity (DBD) problem.Since the final orbit and clock products are estimated on a daily basis, then the DBD is introduced into the clock offset due to the datum difference between different orbit arcs [40].Although the presence of DBD does not affect the positioning performance, it has an impact on the prediction precision assessment and clock modeling over multiple days [41].In this study, the DBD is estimated first, and then its effect is removed.The DBD is estimated by comparing the difference between the first epoch of the prediction section and the first epoch of real-estimated clock offset.We only use the last hour data of the observed section to fit the clock model and predict clock for DBD estimation in order to ensure the prediction accuracy.After that, the clock offset of the first day is aligned to the second day by adding the estimated DBD value.

Cross-Validation for Model Selection
The key step for the adaptively optimal model selection is the cross-validation procedure.This procedure includes three sub-steps: Data segmentation, multi-model clock fitting and model selection.
The major drawback of the current clock prediction strategy is that the prediction procedure is an open-loop procedure.There is no feedback on the quality of clock prediction, thus it is difficult to ensure the quality when the clock model becomes complex.In this study, we form a closed-loop procedure for the model selection, so that the optimal decision can be made based upon the feedback information.

Data Preprocessing
Data preprocessing is a vital procedure for the clock prediction, which removes the outliers and handles the daily boundaries to produce clean data for the following data processing.The data with gross error or outliers cannot reflect the nature of clocks, thus should not be involved in clock modeling and prediction.The popular median absolute deviation (MAD) method is used to detect the gross error.This method has a high breakdown point and is robust to multiple outliers, which can be expressed as follows [39]: where the y i presents the clock frequency sequence, m is the median of the clock frequency sequence, m = median y i .If the clock frequency satisfies y i > m + n * MAD or y i < m − n * MAD the clock offset will be regarded as a gross error.In this case, n is set as 3, and the detected clock offset gross error is excluded in the data processing.
Another problem in preprocessing is the daily boundary discontinuity (DBD) problem.Since the final orbit and clock products are estimated on a daily basis, then the DBD is introduced into the clock offset due to the datum difference between different orbit arcs [40].Although the presence of DBD does not affect the positioning performance, it has an impact on the prediction precision assessment and clock modeling over multiple days [41].In this study, the DBD is estimated first, and then its effect is removed.The DBD is estimated by comparing the difference between the first epoch of the prediction section and the first epoch of real-estimated clock offset.We only use the last hour data of the observed section to fit the clock model and predict clock for DBD estimation in order to ensure the prediction accuracy.After that, the clock offset of the first day is aligned to the second day by adding the estimated DBD value.

Cross-Validation for Model Selection
The key step for the adaptively optimal model selection is the cross-validation procedure.This procedure includes three sub-steps: Data segmentation, multi-model clock fitting and model selection.
The major drawback of the current clock prediction strategy is that the prediction procedure is an open-loop procedure.There is no feedback on the quality of clock prediction, thus it is difficult to ensure the quality when the clock model becomes complex.In this study, we form a closed-loop procedure for the model selection, so that the optimal decision can be made based upon the feedback information.
In the clock prediction procedure, the precise satellite clock estimated with the GNSS signals is denoted as the observation section, which is used for clock modeling.The period for clock prediction is known as the prediction section.It is clear that no prior clock information is available before the clock prediction procedure.In order to capture the local clock characteristics for modeling, the observation section is usually adjacent to the prediction section.A demonstration of the data segmentation strategy is illustrated in Figure 7.In the cross-validation procedure, the observation section is further divided into the fitting part and the validation part.The validation part comes after the fitting part, which is used to check the model performance.How to reasonability divide the data depends upon the user choice.A longer validation part may affect the model fitting and a shorter validation part may affect the validation accuracy.In this study, we use a four-hour window for the validation purpose.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 25 In the clock prediction procedure, the precise satellite clock estimated with the GNSS signals is denoted as the observation section, which is used for clock modeling.The period for clock prediction is known as the prediction section.It is clear that no prior clock information is available before the clock prediction procedure.In order to capture the local clock characteristics for modeling, the observation section is usually adjacent to the prediction section.A demonstration of the data segmentation strategy is illustrated in Figure 7.In the cross-validation procedure, the observation section is further divided into the fitting part and the validation part.The validation part comes after the fitting part, which is used to check the model performance.How to reasonability divide the data depends upon the user choice.A longer validation part may affect the model fitting and a shorter validation part may affect the validation accuracy.In this study, we use a four-hour window for the validation purpose.In the cross-validation procedure, different clock models are fitted with the fitting part, one by one.The model coefficient estimated from the fitting part is then used to predict the clock of the validation part.The prediction error can be obtained by comparing the predicted clock and the observations in the validation part.Then the prediction error is used as the feedback information to assist the model selection.According to Table 3, we consider the two models which achieve the best BeiDou clock prediction accuracy in this study, so both of them are used in the multi-model clock fitting stage.Since the fitting part is not enough to estimate the periodic terms precisely due to the data segment, we use the model selection method to choose the linear model or the quadratic model.According to the previous analysis, the length of the observation section is model dependent.The 24-h observation section is adapted for the linear model and the 48-h observation section is used for the quadratic model.
The cross-validation procedure is used to address the optimal clock model from the prediction error.The root mean square (RMS) is used as the overall prediction error indicator, which is defined in Equation (5).In this study, we only discriminate the optimal clock model from two candidates, hence the unified clock prediction model can be expressed as: where  is the adaptive factor controlling the model selection. and  are the RMS of the prediction error, subject to the linear model and quadratic model in the validation part.When  = 0, Equation ( 6) becomes the linear model combined with two periodic terms and 24-h observation data fitting is used for this model.When  = 1, Equation ( 6) becomes the quadratic model and 48-h observation data is preferred in data fitting.If data outage occurs in the 48-h data, 24-h data is also acceptable for the quadratic model fitting to ensure the prediction data availability.Hence, the outcome of the cross validation is the adaptive factor .A different adaptive factor means the final clock prediction procedure uses a different clock model.In the cross-validation procedure, different clock models are fitted with the fitting part, one by one.The model coefficient estimated from the fitting part is then used to predict the clock of the validation part.The prediction error can be obtained by comparing the predicted clock and the observations in the validation part.Then the prediction error is used as the feedback information to assist the model selection.According to Table 3, we consider the two models which achieve the best BeiDou clock prediction accuracy in this study, so both of them are used in the multi-model clock fitting stage.Since the fitting part is not enough to estimate the periodic terms precisely due to the data segment, we use the model selection method to choose the linear model or the quadratic model.According to the previous analysis, the length of the observation section is model dependent.The 24-h observation section is adapted for the linear model and the 48-h observation section is used for the quadratic model.
The cross-validation procedure is used to address the optimal clock model from the prediction error.The root mean square (RMS) is used as the overall prediction error indicator, which is defined in Equation (5).In this study, we only discriminate the optimal clock model from two candidates, hence the unified clock prediction model can be expressed as: where β is the adaptive factor controlling the model selection.RMS linear and RMS quadratic are the RMS of the prediction error, subject to the linear model and quadratic model in the validation part.When β = 0, Equation ( 6) becomes the linear model combined with two periodic terms and 24-h observation data fitting is used for this model.When β = 1, Equation ( 6) becomes the quadratic model and 48-h observation data is preferred in data fitting.If data outage occurs in the 48-h data, 24-h data is also acceptable for the quadratic model fitting to ensure the prediction data availability.Hence, the outcome of the cross validation is the adaptive factor β. A different adaptive factor means the final clock prediction procedure uses a different clock model.

Clock Offset Prediction
The final step is the final clock prediction, which follows the normal procedure.The difference of the adaptive model is that we chose the optimal clock model for each BeiDou satellite according to the cross-validation results.In this step, the clock model will be fitted based on the whole observation section of the clock offset.Then the fitted model is used to predict the next 24-h clock offset.The adaptive model introduces the cross-validation procedure to form a close-loop during the model selection, so that it allows the optimal selection of the clock model based on the feedbacks of the validation.As a result, the algorithm is 'smarter' to choose the right model.In case of the clock characteristics changed over time, the adaptive can also respond according to the feedbacks and still give the optimal clock model.

Accuracy Assessment of the Adaptive Model
In order to evaluate the performance of the adaptive model, we compared the predicted clock with the officially released ultra-rapid clock products.In this study, the ultra-rapid clock products provided by the German Research Centre for Geosciences (GFZ) and the International GNSS Monitoring & Assessment System (iGMAS) are used for verification purposes and are denoted as GBU and ISU, respectively.The ultra-rapid clock products provide a 24-h clock prediction with about a 3 h delay and 6 h updating interval.Both GBU and ISU products include a 24-h observation section and a 24-h prediction section.These two organizations provide stable, regularly released ultra-rapid clock products to support various real-time applications.The detailed strategies for clock prediction of these two institutes are not published, but their products are trustful and world recognized.Due to different clock estimation strategies and clock datum, the clock products from different institutes cannot be mixed up directly, so we use the precise clock product from each institute and perform clock prediction respectively.The prediction part of GBU and ISU are denoted as GBU-P and ISU-P, respectively.The predicted clock offsets by the adaptive model based upon GBU and ISU products are denoted as GBU-A and ISU-A for short.

Data Description
In this study, the GBU and ISU ultra-rapid products from DOY 1 to DOY 259 in 2018 are used, and these products are publicly accessible from their websites (ftp://ftp.gfz-potsdam.de,http://www.igmas.org/)[42].Before the evaluation, we exclude the satellites with clock jumps or with excessive outliers.If the outlier rate of the clock data is higher than 10% within 24-h, the data are then excluded during the analysis.The data feasibility of the GBU and ISU products during the examination period is presented in Figure 8.The figure indicates that the ISU products have better data integrity and quality than the GBU products.The data quality of GEO satellites is not as good as the IGSO and MEO satellites for both GBU and ISU products.Particularly, it is noticed that the C04 satellite has very poor data quality in both products.

Clock Offset Prediction
The final step is the final clock prediction, which follows the normal procedure.The difference of the adaptive model is that we chose the optimal clock model for each BeiDou satellite according to the cross-validation results.In this step, the clock model will be fitted based on the whole observation section of the clock offset.Then the fitted model is used to predict the next 24-h clock offset.The adaptive model introduces the cross-validation procedure to form a close-loop during the model selection, so that it allows the optimal selection of the clock model based on the feedbacks of the validation.As a result, the algorithm is 'smarter' to choose the right model.In case of the clock characteristics changed over time, the adaptive can also respond according to the feedbacks and still give the optimal clock model.

Accuracy Assessment of the Adaptive Model
In order to evaluate the performance of the adaptive model, we compared the predicted clock with the officially released ultra-rapid clock products.In this study, the ultra-rapid clock products provided by the German Research Centre for Geosciences (GFZ) and the International GNSS Monitoring & Assessment System (iGMAS) are used for verification purposes and are denoted as GBU and ISU, respectively.The ultra-rapid clock products provide a 24-h clock prediction with about a 3 h delay and 6 h updating interval.Both GBU and ISU products include a 24-h observation section and a 24-h prediction section.These two organizations provide stable, regularly released ultra-rapid clock products to support various real-time applications.The detailed strategies for clock prediction of these two institutes are not published, but their products are trustful and world recognized.Due to different clock estimation strategies and clock datum, the clock products from different institutes cannot be mixed up directly, so we use the precise clock product from each institute and perform clock prediction respectively.The prediction part of GBU and ISU are denoted as GBU-P and ISU-P, respectively.The predicted clock offsets by the adaptive model based upon GBU and ISU products are denoted as GBU-A and ISU-A for short.
In this study, the GBU and ISU ultra-rapid products from DOY 1 to DOY 259 in 2018 are used, and these products are publicly accessible from their websites (ftp://ftp.gfz-potsdam.de,http://www.igmas.org/)[42].Before the evaluation, we exclude the satellites with clock jumps or with excessive outliers.If the outlier rate of the clock data is higher than 10% within 24-h, the data are then excluded during the analysis.The data feasibility of the GBU and ISU products during the examination period is presented in Figure 8.The figure indicates that the ISU products have better data integrity and quality than the GBU products.The data quality of GEO satellites is not as good as the IGSO and MEO satellites for both GBU and ISU products.Particularly, it is noticed that the C04 satellite has very poor data quality in both products.

Accuracy Assessment of Clock Predicted by the Adaptive Model
In order to evaluate the prediction precision of the adaptive model, we compared the BeiDou clock predicted by the adaptive model and the GBU-P/ISU-P products.The clock prediction process is performed on a daily basis.The predicted clock products using the adaptive model is denoted as the GBU-A and ISU-A, respectively, according to the different precise clock sources.The prediction clock error is evaluated by comparing the predicted clock to the observed part of the ultra-rapid clock products.The RMS is used as the overall prediction error indicator.
To demonstrate the effect of the adaptive strategy, we compared the clock prediction results with the adaptive strategy and without the adaptive strategy.The clock prediction without, using the selected optimal model in the validation period, are denoted as GBU-W and ISU-W.The average RMS of the 24 h prediction results is shown in Table 4, which clearly shows that the adaptive strategy is effective.In most cases, the GBU-W and ISU-W are better than GBU-P and ISU-P, which is because the two model candidates are carefully selected and analyzed in Section 2. The selected model candidates are effective.Furthermore, the comparison between GBU-A/ISU-A and GBU-W/ISU-W also shows that the model selection strategy is effective.With proper model selection, the prediction error can be further improved.Therefore, the following section will further compare the details of GBU-P/ISU-P and GBU-A/ISU-A.The overall RMS of 24-h clock prediction subject to the GBU-P product and ISU-P product are presented in Figure 9. Generally, the adaptive model achieves higher clock prediction accuracy than the GBU-P and ISU-P products.The left panel shows the overall RMS of a 24-h clock predicted by GBU-P and GBU-A, respectively.The figure indicates that the adaptive model achieves a comparable or higher prediction accuracy than the GBU-P products.Particularly, the adaptive model achieves a significant prediction accuracy improvement for C03, C05, C06, C12 and C14, and the maximum accuracy improvement reaches 5 ns.The right panel shows the overall RMS comparison between the ISU-P products and ISU-A products.The figure shows that the adaptive model improves the clock prediction accuracy for most BeiDou satellites.The results indicate that GFZ and iGMAS may adopt different clock prediction strategy.The figure also indicates that the adaptive model can dramatically decrease the prediction error for some satellites, such as the C14 for the GBU-P product and C05 for the ISU-P product.The precision improvement for these satellites exceeds 50%.It means that the original algorithms may not fit the characteristics of these satellites.Generally, the precision improvement is not homogeneous.The adaptive model achieves a more significant improvement for some satellites, and less significant for the others.It also confirms that the BeiDou clocks should not be predicted with a single uniform model.The adaptive model can properly handle most satellite clocks, but there are still a few BeiDou satellites having poor predictability, and all three methods cannot achieve good predict precision, such as the C04 satellite.As stated, the C04 satellite has poor data quality, and how to improve it still needs further investigation.The RMS comparison can only give an overview comparison.We further studied the relationship between the prediction error and the prediction length.In order to reflect the prediction error variation with multiple days' data, we defined an epoch-wise RMS which reflects the RMS of the prediction error for the same prediction epoch over multiple days.The RMS for the i epoch over multiple days is defined as:

GBU-P GBU-A GBU-W ISU-P ISU-
where  means the predicted clock offset at the  epoch on the  day, and  presents the estimated clock offset in the observation part of ultra-rapid products at this  epoch on this  day, and  is the number of days.The epoch-wise  is computed from the prediction error of the  epoch over multiple days, so that it gives a more reliable quality indicator on the prediction error on each prediction epoch.Figures 10 and 11 illustrate the prediction error change along the prediction length.Figure 10 describes the characteristic of prediction error for the BeiDou satellites C05, C12 and C14 in GBU-P and GBU-A products.The top row shows the prediction error of each day and the epoch-wise RMS for the GBU-P products.The middle row shows the same for the GBU-A products and the bottom row compared the epoch-wise RMS for the GBU-P and GBU-A products.The figure indicates that the prediction error generally increases as the prediction length increases.The figure indicates that the adaptive model achieving a better prediction precision is not an occasional case, since the prediction precision improvement is indeed caused by the employing more reasonable clock model for each satellite in the adaptive model.
A similar comparison between the ISU-P and ISU-A products is also performed and the results are presented in Figure 11.The figure shows the prediction error for the satellites C02, C08, and C13.The figure shows that the adaptive model achieves similar results to the GBU products as Figure 10.The figure demonstrates that the adaptive model is not only suitable for GBU products, but is also valid for the ISU products.The figure also indicates that the epoch-wise RMS varies between 10-25 ns at the last epoch of 24-h prediction for the GBU-P and ISU-P products, whereas the adaptive model can reduce to 5-10 ns.The RMS comparison can only give an overview comparison.We further studied the relationship between the prediction error and the prediction length.In order to reflect the prediction error variation with multiple days' data, we defined an epoch-wise RMS which reflects the RMS of the prediction error for the same prediction epoch over multiple days.The RMS for the i th epoch over multiple days is defined as: where y ji means the predicted clock offset at the i epoch on the j day, and ŷji presents the estimated clock offset in the observation part of ultra-rapid products at this i epoch on this j day, and d is the number of days.The epoch-wise RMS i is computed from the prediction error of the i epoch over multiple days, so that it gives a more reliable quality indicator on the prediction error on each prediction epoch.Figures 10 and 11 illustrate the prediction error change along the prediction length.Figure 10 describes the characteristic of prediction error for the BeiDou satellites C05, C12 and C14 in GBU-P and GBU-A products.The top row shows the prediction error of each day and the epoch-wise RMS for the GBU-P products.The middle row shows the same for the GBU-A products and the bottom row compared the epoch-wise RMS for the GBU-P and GBU-A products.The figure indicates that the prediction error generally increases as the prediction length increases.The figure indicates that the adaptive model achieving a better prediction precision is not an occasional case, since the prediction precision improvement is indeed caused by the employing more reasonable clock model for each satellite in the adaptive model.
A similar comparison between the ISU-P and ISU-A products is also performed and the results are presented in Figure 11.The figure shows the prediction error for the satellites C02, C08, and C13.The figure shows that the adaptive model achieves similar results to the GBU products as Figure 10.The figure demonstrates that the adaptive model is not only suitable for GBU products, but is also valid for the ISU products.The figure also indicates that the epoch-wise RMS varies between 10-25 ns at the last epoch of 24-h prediction for the GBU-P and ISU-P products, whereas the adaptive model can reduce to 5-10 ns.The improvement of the adaptive model can also be reflected by checking the consistency between the satellites.Figure 12 shows the prediction error of all BeiDou satellites on DOY 121, The figure indicates that the prediction error of the adaptive model is more concentrated around zero than the GBU-P and ISU-P products.The GBU-P model can predict clock offset well for most BeiDou satellites, but there are a few satellites which possess extremely poor prediction precision.The ISU-P model has a poorer prediction precision than the GBU-P products for most satellites, but its worst case is better than the GBU-P model.Particularly, the adaptive model reduces the maximum prediction error for both the GBU-P and ISU-P products.It is because the adaptive model can automatically select the best model for each satellite, so that it has better adaption on these satellites with different characteristics.The figure also indicates that the observed section has a strong impact on the clock prediction error.Both GBU-A and ISU-A employs the same prediction method, but the prediction error is completely different, which is because they employ a different observation section in the clock modeling phase.The improvement of the adaptive model can also be reflected by checking the consistency between the satellites.Figure 12 shows the prediction error of all BeiDou satellites on DOY 121, The figure indicates that the prediction error of the adaptive model is more concentrated around zero than the GBU-P and ISU-P products.The GBU-P model can predict clock offset well for most BeiDou satellites, but there are a few satellites which possess extremely poor prediction precision.The ISU-P model has a poorer prediction precision than the GBU-P products for most satellites, but its worst case is better than the GBU-P model.Particularly, the adaptive model reduces the maximum prediction error for both the GBU-P and ISU-P products.It is because the adaptive model can automatically select the best model for each satellite, so that it has better adaption on these satellites with different characteristics.The figure also indicates that the observed section has a strong impact on the clock prediction error.Both GBU-A and ISU-A employs the same prediction method, but the prediction error is completely different, which is because they employ a different observation section in the clock modeling phase.

Accuracy Comparison with Different Prediction Length
As aforementioned, the clock prediction error increases as the prediction length is also increasing.In practice, a 24-h clock prediction may not meet the accuracy requirements sometimes, so that we also evaluated the prediction precision with different prediction length.In this study, four typical prediction lengths, say 3-h, 6-h, 12-h and 24-h, are considered.We calculated the overall RMS of prediction for different satellite types with all 259-day data.The results are presented in Table 5 and Table 6.
Table 5 shows the overall RMS comparison of the GBU-P and GBU-A products with different prediction lengths.It shows that the GBU-A model achieves more than a 1 ns precision improvement for most cases.The improvement is extremely significant for short-term prediction.The GBU-A model contributes to 47.3% precision improvement for 3-h prediction.For 24-h prediction, the value decreases to 14.2%, but the absolute RMS improvement is quite similar.Regarding different satellite types, the GEO satellites have poor prediction accuracy, especially the satellite C04.The MEO satellites achieve the best clock prediction precision.The overall prediction error for the GBU-P products is smaller than 3 ns for 6-h prediction, whereas the adaptive model extends it to 12-h.
Similarly, the comparison results between the ISU-P and ISU-A products are listed in Table 6.Similar to Table 5, the adaptive model achieves better prediction precision than the ISU-P for all different satellite types and different prediction lengths.For 24-h prediction, the precision improvement of ISU-A for the GEO, IGSO and MEO satellites reaches 35.7%, 34.7% and 17.3%, respectively.The table indicates that the ISU-P achieves a very good prediction precision for these MEO satellites, but not for the GEO and IGSO satellites.With regard to the GEO and IGSO satellite, the ISU-P products outperform the GBU-P for short-term prediction, and achieve poorer precision for long-term prediction.It concludes that the prediction strategy of the ISU-P products is extremely suitable for MEO satellite prediction, but may not be the best choice for the GEO and IGSO satellites.The adaptive model takes advantage of its better adaption on different clock types, and achieves high prediction accuracy for all satellite types.

Accuracy Comparison with Different Prediction Length
As aforementioned, the clock prediction error increases as the prediction length is also increasing.In practice, a 24-h clock prediction may not meet the accuracy requirements sometimes, so that we also evaluated the prediction precision with different prediction length.In this study, four typical prediction lengths, say 3-h, 6-h, 12-h and 24-h, are considered.We calculated the overall RMS of prediction for different satellite types with all 259-day data.The results are presented in Tables 5 and 6.
Table 5 shows the overall RMS comparison of the GBU-P and GBU-A products with different prediction lengths.It shows that the GBU-A model achieves more than a 1 ns precision improvement for most cases.The improvement is extremely significant for short-term prediction.The GBU-A model contributes to 47.3% precision improvement for 3-h prediction.For 24-h prediction, the value decreases to 14.2%, but the absolute RMS improvement is quite similar.Regarding different satellite types, the GEO satellites have poor prediction accuracy, especially the satellite C04.The MEO satellites achieve the best clock prediction precision.The overall prediction error for the GBU-P products is smaller than 3 ns for 6-h prediction, whereas the adaptive model extends it to 12-h.
Similarly, the comparison results between the ISU-P and ISU-A products are listed in Table 6.Similar to Table 5, the adaptive model achieves better prediction precision than the ISU-P for all different satellite types and different prediction lengths.For 24-h prediction, the precision improvement of ISU-A for the GEO, IGSO and MEO satellites reaches 35.7%, 34.7% and 17.3%, respectively.The table indicates that the ISU-P achieves a very good prediction precision for these MEO satellites, but not for the GEO and IGSO satellites.With regard to the GEO and IGSO satellite, the ISU-P products outperform the GBU-P for short-term prediction, and achieve poorer precision for long-term prediction.It concludes that the prediction strategy of the ISU-P products is extremely suitable for MEO satellite prediction, but may not be the best choice for the GEO and IGSO satellites.The adaptive model takes advantage of its better adaption on different clock types, and achieves high prediction accuracy for all satellite types.
It also concludes that the clock offset prediction precision is also affected by the prediction length, satellite type, precise clock source and the prediction algorithms.Generally, longer prediction length presents a larger prediction error.
MEO satellites present smaller clock prediction error than the IGSO and GEO satellites.GEO satellites have the poorest clock predictability, which may be affected by their inaccurate orbit product.Currently, the precise orbit product for BeiDou GEO is significantly worse than the other satellites, due to poor observability, and the remaining orbit error may be assimilated into the precise clock product automatically.The clock prediction strategy of ISU-P is more suitable for short-term clock prediction, but its 24-h clock prediction is not as good as the GBU-P product.The adaptive model can significantly improve the clock prediction accuracy for different clock products and different prediction lengths.The GBU-A products improve clock prediction error by over 40% for short-term prediction and 14% for 24-h prediction, compared to the GBU-P product.The ISU-A product improves over 30% clock prediction error compared to the ISU-P product for all different prediction lengths.Generally, the precision of BeiDou clock prediction with the adaptive model achieves better than 1.5 ns for 3-h prediction, 2 ns for 6-h prediction, 3 ns for 12-h prediction and 6 ns for 24-h prediction.The prediction precision for BeiDou GEO satellites may be lower sometimes, but we believe the clock prediction error will be improved with a better precise orbit product.

Accuracy Assessment of PPP with the Predicted Clock
The previous section compared the predicted clock with the reference clock products.The correctness and efficiency of the predicted clock is further tested by applying them to PPP data processing.Considering the Beidou-2 system only provides the navigation service in the Asia-Pacific region, we select six stations from the multi-GNSS experiment network (MGEX) to test the PPP performance by applying the predicted clock product [43].The distribution of the stations is illustrated in Figure 13.All of these stations are capable of tracking Beidou-2 signals, and the observation data are publicly accessible at ftp://igs.ign.fr/.
In the PPP data processing, both the predicted orbit and clock are used to investigate the potential positioning accuracy for real-time PPP data processing.The data are processed with the GBU-P, GBU-A, ISU-P and ISU-A clock products, respectively.Considering the compatibility between the orbit and the clock products, the GBU-P and ISU-P orbit products are used when the GBU-A and ISU-A clocks are applied respectively in PPP data processing.
In the PPP data processing, both the predicted orbit and clock are used to investigate the potential positioning accuracy for real-time PPP data processing.The data are processed with the GBU-P, GBU-A, ISU-P and ISU-A clock products, respectively.Considering the compatibility between the orbit and the clock products, the GBU-P and ISU-P orbit products are used when the GBU-A and ISU-A clocks are applied respectively in PPP data processing.

Static Precise Point Positioning
The observation data used in this study is collected on DOY 205, 2018 with a 30-s interval.As stated, the precision of the 24-h clock prediction is considered too poor for PPP data processing, so the 12-h clock prediction result is used in PPP data processing.The open-source software package RTKLIB is used for PPP data processing [24] and the process strategy is listed in Table 7.In the static PPP data processing, only BeiDou observation is used and the bidirectional filtering technique is applied to eliminate the convergence effect.The positioning accuracy of PPP is assessed by comparing the positioning results to the reference station coordinates provided by the international GNSS service (IGS).The daily RMS of the static PPP positioning error for the six stations is illustrated in Figure 14, which indicates that the solution with the clock offset predicted by the adaptive model achieves a better precision than the solution computed with the GBU-P and ISU-P products, especially in the east and north direction.The RMS values of the adaptive model are less than 0.2 m in the east and north direction, but many stations have their precision poorer than 0.2 m if the GBU-P and ISU-P products are used.In the up direction, the improvement of the adaptive model is not as significant as the horizontal direction.

Static Precise Point Positioning
The observation data used in this study is collected on DOY 205, 2018 with a 30-s interval.As stated, the precision of the 24-h clock prediction is considered too poor for PPP data processing, so the 12-h clock prediction result is used in PPP data processing.The open-source software package RTKLIB is used for PPP data processing [24] and the process strategy is listed in Table 7.In the static PPP data processing, only BeiDou observation is used and the bidirectional filtering technique is applied to eliminate the convergence effect.The positioning accuracy of PPP is assessed by comparing the positioning results to the reference station coordinates provided by the international GNSS service (IGS).The daily RMS of the static PPP positioning error for the six stations is illustrated in Figure 14, which indicates that the solution with the clock offset predicted by the adaptive model achieves a better precision than the solution computed with the GBU-P and ISU-P products, especially in the east and north direction.The RMS values of the adaptive model are less than 0.2 m in the east and north direction, but many stations have their precision poorer than 0.2 m if the GBU-P and ISU-P products are used.In the up direction, the improvement of the adaptive model is not as significant as the horizontal direction.The overall precision of the static PPP results is compared, and the results are presented in Table 8.The solution with better precision is highlighted in the table.The table indicates that static PPP results with the clock predicted by the adaptive model have significantly improved the positioning precision, especially in the east and north direction.Since the precision of the predicted clock is not as good as the final products, so corresponding PPP precision is also degraded.The overall precision of the PPP results with 12-h clock prediction reaches about 10 cm horizontal precision and 20-30 cm vertical precision.The adaptive model can improve the PPP precision by adopting both GBU and ISU products.The average PPP precision improvement for the GBU products reaches 54.7%, 32.6% and 7.6% in the east, north and up directions, and the improvement for the ISU products reaches 47.9%, 54.2% and 13.7%, respectively.The 3D RMS of PPP with the clock predicted by the adaptive model is improved by 21.6% and 30.0%compared to the results computed with GBU-P and ISU-P products.With 12-h clock prediction, the clock prediction error reaches 2-3 ns, which corresponds to a few decimeters ranging error on line-of-sight direction, but it still concludes that the 12-h clock prediction is capable of supporting about a 1-2 decimeter level horizontal precision PPP positioning, which can meet the requirement of many real-time precise positioning applications.With shorter prediction length, a better static PPP precision can be expected.The overall precision of the static PPP results is compared, and the results are presented in Table 8.The solution with better precision is highlighted in the table.The table indicates that static PPP results with the clock predicted by the adaptive model have significantly improved the positioning precision, especially in the east and north direction.Since the precision of the predicted clock is not as good as the final products, so corresponding PPP precision is also degraded.The overall precision of the PPP results with 12-h clock prediction reaches about 10 cm horizontal precision and 20-30 cm vertical precision.The adaptive model can improve the PPP precision by adopting both GBU and ISU products.The average PPP precision improvement for the GBU products reaches 54.7%, 32.6% and 7.6% in the east, north and up directions, and the improvement for the ISU products reaches 47.9%, 54.2% and 13.7%, respectively.The 3D RMS of PPP with the clock predicted by the adaptive model is improved by 21.6% and 30.0%compared to the results computed with GBU-P and ISU-P products.With 12-h clock prediction, the clock prediction error reaches 2-3 ns, which corresponds to a few decimeters ranging error on line-of-sight direction, but it still concludes that the 12-h clock prediction is capable of supporting about a 1-2 decimeter level horizontal precision PPP positioning, which can meet the requirement of many real-time precise positioning applications.With shorter prediction length, a better static PPP precision can be expected.

Kinematic Precise Point Positioning
The performance of the clock product predicted with the adaptive model is verified with the kinematic PPP data processing.In this test, we use a 3-h clock prediction product to mitigate the accumulated prediction error.In this test, four stations from the MGEX network named KARR, MRO1, PERT and PNGM are used, and the data was collected on DOY 121.Both GPS and BeiDou observations are employed in this data processing, while the other process strategies are the same as the static PPP experiment.The positioning results of the stations is compared to the IGS reference coordinates, and the positioning error of each station is presented in Figures 15 and 16. Figure 15 which shows that there is no significant difference in convergence time between GBU-P and GBU-A results.However, when the GBU-P product is utilized, the precision gradually reduces after 2 h, and the error in the up direction may reach to 5 decimeters, which mainly caused by the increasing clock prediction error.Also, the precision of the east and north components is declined significantly.On the other hand, the kinematic PPP precision more stable after convergence with the GBU-A products, and the precision stays in the one decimeter level.The kinematic PPP results by adopting ISU-P and ISU-A products are illustrated in Figure 16, which indicates that the error of up direction is increasing significantly with the ISU-P product.The ISU-P product can lead a 5 decimeters positioning error in the up direction, whereas the ISU-A product can achieve even better than 2 decimeters kinematic PPP precision.

Kinematic Precise Point Positioning
The performance of the clock product predicted with the adaptive model is verified with the kinematic PPP data processing.In this test, we use a 3-h clock prediction product to mitigate the accumulated prediction error.In this test, four stations from the MGEX network named KARR, MRO1, PERT and PNGM are used, and the data was collected on DOY 121.Both GPS and BeiDou observations are employed in this data processing, while the other process strategies are the same as the static PPP experiment.The positioning results of the stations is compared to the IGS reference coordinates, and the positioning error of each station is presented in Figure 15 and Figure 16. Figure 15 which shows that there is no significant difference in convergence time between GBU-P and GBU-A results.However, when the GBU-P product is utilized, the precision gradually reduces after 2 h, and the error in the up direction may reach to 5 decimeters, which mainly caused by the increasing clock prediction error.Also, the precision of the east and north components is declined significantly.On the other hand, the kinematic PPP precision is more stable after convergence with the GBU-A products, and the precision stays in the one decimeter level.The kinematic PPP results by adopting ISU-P and ISU-A products are illustrated in Figure 16, which indicates that the error of up direction is increasing significantly with the ISU-P product.The ISU-P product can lead a 5 decimeters positioning error in the up direction, whereas the ISU-A product can achieve even better than 2 decimeters kinematic PPP precision.

Conclusions
The clock prediction is a promising method to support real-time PPP applications.The GPS clock prediction methods have been well studied, but the same method does not achieve satisfactory precision in BeiDou clock prediction.This study points out the uniqueness of the BeiDou clocks is that the clock characteristics are highly satellite-specified, and the characteristics may change by clock operations or clock events.Hence it is difficult to achieve good clock prediction results with a unified and static model for all BeiDou satellites.By analyzing the characteristics of BeiDou clocks, we addressed two optimal clock models for different BeiDou clock status and proposed a new adaptive method for BeiDou clock prediction.In the adaptive model, a cross-validation procedure is introduced to optimally select the clock models.The cross-validation procedure tests the RMS of the clock prediction error with different clock models with a small portion of data, and then selects the optimal model with the test results.The closed-loop design increases the flexibility of the adaptive model, so that it achieves the best prediction performance for different types of clocks in the BeiDou constellation.The performance of the proposed method is validated by the ultra-rapid clock products released by two different organizations (GBU and ISU), and the comparison based on 259-day data indicates that the adaptive model improves the GBU-P prediction precision by 47.3%, 41.0%, 28.8% and 14.2% for 3-h, 6-h, 12-h and 24-h prediction, respectively.The corresponding precision improvements comparing to the ISU-P products are 32.1%,34.07%, 36.4% and 33.2%, respectively.The adaptive model can achieve 1 ns precision for a 3-h prediction, and better than 2 ns for a 6-h prediction.Moreover, the clock prediction precision for the MEO satellites is better than that the GEO and IGSO satellites.
The performance of the adaptive model is further validated by the static PPP and the kinematic PPP data processing.As for the static PPP, observation data from six stations in the MGEX network are processed with different predicted clock products, and the results indicate that the clock predicted with the adaptive model improves the precision remarkably.The 3D RMS of static PPP with the clock offset predicted by the adaptive model is improved by 21.6% and 30.0%compared to use the GBU-P and ISU-P products.The north direction achieves the most significant precision improvement, which even reaches around 50%.In the kinematic PPP mode, the improvement of the clock predicted by the adaptive model is even more significant.Comparing with the GBU-P product, the adaptive model achieves a 29.0%, 75.7% and 47.0% positioning precision improvement in east, north, and up directions, respectively.The positioning precision by the adaptive model improves 52.3%, 19.9% and 64.7% in the three components compared to the ISU-P product.Therefore, the adaptive model can significantly improve the BeiDou clock prediction precision and the PPP positioning precision by considering the satellite-specified clock characteristics, it concludes that the adaptive model is a practical and efficient clock prediction method for BeiDou clock prediction.

Figure 1 .
Figure 1.Average root mean square (RMS) of 24-h clock prediction with different fitting length.

Figure 2 .
Figure 2. BeiDou clock offset prediction error, comparing the linear model with the quadratic model.

Figure 1 .
Figure 1.Average root mean square (RMS) of 24-h clock prediction with different fitting length.

25 Figure 1 .
Figure 1.Average root mean square (RMS) of 24-h clock prediction with different fitting length.

Figure 2 .
Figure 2. BeiDou clock offset prediction error, comparing the linear model with the quadratic model.

Figure 2 .
Figure 2. BeiDou clock offset prediction error, comparing the linear model with the quadratic model.

Figure 4 .
Figure 4. RMS of a clock offset fitting residual with polynomial model combined periodic terms.(a) Fitting residuals of the linear model plus different periodic terms; (b) fitting residuals of the quadratic model plus different periodic terms.

Figure 4 .
Figure 4. RMS of a clock offset fitting residual with polynomial model combined periodic terms.(a) Fitting residuals of the linear model plus different periodic terms; (b) fitting residuals of the quadratic model plus different periodic terms.

Figure 5 .
Figure 5. Improvement of the clock prediction result by adding different periodic terms.(a) Improvement of the linear model with periodic terms, comparing with the basic linear model; (b) Improvement of the quadratic model with periodic terms, comparing with the basic quadratic model.

Figure 5 .
Figure 5. Improvement of the clock prediction result by adding different periodic terms.(a) Improvement of the linear model with periodic terms, comparing with the basic linear model; (b) Improvement of the quadratic model with periodic terms, comparing with the basic quadratic model.

Figure 6 .
Figure 6.Flowchart of the adaptive BeiDou clock prediction model.

Figure 6 .
Figure 6.Flowchart of the adaptive BeiDou clock prediction model.

Figure 7 .
Figure 7. Demonstration of data segment strategy.

Figure 7 .
Figure 7. Demonstration of data segment strategy.

Figure 8 .
Figure 8. Data feasibility.The red line, blue line and green line are presenting the Geosynchronous Orbit (GEO), Inclined Geosynchronous Satellite Orbit (IGSO) and Medium Earth Orbit (MEO) satellites, respectively.(a) The GBU products; (b) the ISU products.

Figure 10 .
Figure 10.Comparison of clock prediction precision between the adaptive model and the GBU-P product for long-term data.(a) C05 prediction error with GBU-P; (b) C12 prediction error with GBU-P; (c) C14 prediction error with GBU-P; (d) C05 prediction error with GBU-A; (e) C12 prediction error with GBU-A; (f) C14 prediction error with GBU-A; (g) RMS of C05 prediction error; (h) RMS of C12 prediction error; (i) RMS of C14 prediction error.

Figure 10 .
Figure 10.Comparison of clock prediction precision between the adaptive model and the GBU-P product for long-term data.(a) C05 prediction error with GBU-P; (b) C12 prediction error with GBU-P; (c) C14 prediction error with GBU-P; (d) C05 prediction error with GBU-A; (e) C12 prediction error with GBU-A; (f) C14 prediction error with GBU-A; (g) RMS of C05 prediction error; (h) RMS of C12 prediction error; (i) RMS of C14 prediction error.

Figure 11 .
Figure 11.Comparison of clock prediction precision between the adaptive model and ISU-P product for long-term data.(a) C02 prediction error with ISU-P; (b) C08 prediction error with ISU-P; (c) C13 prediction error with ISU-P; (d) C02 prediction error with ISU-A; (e) C08 prediction error with ISU-A; (f) C13 prediction error with ISU-A; (g) RMS of C02 prediction error; (h) RMS of C08 prediction error; (i) RMS of C13 prediction error.

Figure 11 .
Figure 11.Comparison of clock prediction precision between the adaptive model and ISU-P product for long-term data.(a) C02 prediction error with ISU-P; (b) C08 prediction error with ISU-P; (c) C13 prediction error with ISU-P; (d) C02 prediction error with ISU-A; (e) C08 prediction error with ISU-A; (f) C13 prediction error with ISU-A; (g) RMS of C02 prediction error; (h) RMS of C08 prediction error; (i) RMS of C13 prediction error.

Figure 12 .
Figure 12.Comparison of the prediction error between the adaptive model and GBU-P, ISU-P products.(a) prediction error with GBU-P; (b) prediction error with GBU-A; (c) prediction error with ISU-P; (d) prediction error with ISU-A.

Figure 12 .
Figure 12.Comparison of the prediction error between the adaptive model and GBU-P, ISU-P products.(a) prediction error with GBU-P; (b) prediction error with GBU-A; (c) prediction error with ISU-P; (d) prediction error with ISU-A.

Figure 14 .
Figure 14.Comparison of daily RMS of precise point positioning (PPP) results from six MGEX stations using GBU-P, ISU-P and the adaptive clock product.(a) RMS in the east direction with GBU; (b) RMS in the east direction with ISU; (c) RMS in the north direction with GBU; (d) RMS in the north direction with ISU; (e) RMS in the up direction with GBU; (f) RMS in the up direction with ISU.

Figure 14 .
Figure 14.Comparison of daily RMS of precise point positioning (PPP) results from six MGEX stations using GBU-P, ISU-P and the adaptive clock product.(a) RMS in the east direction with GBU; (b) RMS in the east direction with ISU; (c) RMS in the north direction with GBU; (d) RMS in the north direction with ISU; (e) RMS in the up direction with GBU; (f) RMS in the up direction with ISU.

Figure 15 .
Figure 15.The kinematic PPP positioning results of four MGEX stations shown at (a) GBU-P products; (b) GBU-A products.Figure 15.The kinematic PPP positioning results of four MGEX stations shown at (a) GBU-P products; (b) GBU-A products.

Figure 15 .
Figure 15.The kinematic PPP positioning results of four MGEX stations shown at (a) GBU-P products; (b) GBU-A products.Figure 15.The kinematic PPP positioning results of four MGEX stations shown at (a) GBU-P products; (b) GBU-A products.

Table 3 .
The optimal clock prediction model for BeiDou satellites.

Table 3 .
The optimal clock prediction model for BeiDou satellites.

Table 4 .
The average RMS of 24 h prediction errors of GBU and ISU products (Unit: Nanosecond).

Table 5 .
The prediction error between GBU-P and GBU-A products with different prediction lengths (Unit: Nanosecond).
1IMP presents the improvement of the adaptive model (Unit: %).

Table 6 .
The prediction error between ISU-P and ISU-A products with different prediction lengths (Unit: Nanosecond).

Table 7 .
Data strategy for static PPP.

Table 7 .
Data strategy for static PPP.

Table 8 .
Comparison of overall RMS of PPP results based on GBU-P, ISU-P and the adaptive clock product (Unit: Meter).

Table 8 .
Comparison of overall RMS of PPP results based on GBU-P, ISU-P and the adaptive clock product (Unit: Meter).

Table 9 .
The RMS values of four MGEX stations in the North(N), East(E), and Up(U) components of kinematic PPP solutions (Unit: Meter).