Improving Clock Prediction Algorithm for BDS-2 / 3 Satellites Based on LS-SVM Method

: The satellite clock prediction is crucial to support real-time global satellite precise positioning services. Currently, the clock prediction for the Chinese BeiDou navigation satellite system (BDS) is still challenging to satisfy the precise positioning applications. Based on the exploration of existing prediction models, an improved model combing the spectrum analysis model (SAM) and the least-squares support-vector machine (LS-SVM) is proposed especially for BDS-2 / 3 satellites. Considering satellite-speciﬁc characteristics, the parameters of the LS-SVM method are optimized satellite by satellite, including input length, regularization and kernel parameters. The improved model is evaluated by comparing the predicted clocks of existing methods and the improved model. The bias of the predicted clock o ﬀ sets are within ± 1.0 ns for most medium Earth orbits (MEOs) over three hours employing the improved model, which is better than that of the existing methods and can be applied for several real-time precise positioning applications. The predicted clock o ﬀ sets are further evaluated by applying clock corrections to precise point positioning (PPP) in both static and kinematic modes for 10 international GNSS service (IGS) Multi-GNSS Experiment (MGEX) stations, including ﬁve stations in the Asia-Paciﬁc region. According to the practical engineering experience, 2 dm and 5 dm are deﬁned for static and kinematic PPP, respectively, as a convergence threshold. Then, in the static PPP, the improved model is demonstrated to be e ﬀ ective, and positioning accuracies of some stations obtain more than 15% improvements on average for each direction, which enables them to get sub-decimeter positioning, especially in the Asia-Paciﬁc region. In the kinematic PPP, the improved model performs much better than the others in terms of both the convergence time and the positioning accuracy. The convergence time can be shortened from 1.0 h to below 0.5 h, while the positioning accuracies are enhanced by 16.3%, 10.8%, and 18.9% on average in east, north, and up direction, respectively.


Introduction
The Chinese BeiDou navigation satellite system, abbreviated as BDS [1], has finished its regional system BDS-2 consisting of 14 satellites in 2012, which provides positioning, navigation, timing, and short message communication services for Asia-Pacific users. During the past several years, it has developed dramatically for the purpose of constructing a global navigation system as BDS-3. By the end of 2018, a total of 19 BDS-3 satellites had been launched to achieve a primary system for global services. With the help of altogether 33 operational satellites, i.e., five geostationary earth orbits (GEOs), seven inclined geosynchronous orbits (IGSOs), and 21 medium Earth orbits (MEOs) in orbit, an accuracy of about 2-6 ns are one type of the most important products in providing decimeter-level real-time positioning services for BDS. However, evident nonlinear system errors, for example periodic signals, are detected in ultra-rapid clock products [11,14]. To satisfy the requirements of high-precision application scenarios, such as landslide early warning and hydrographic surveying, an improved model is required to further enhance its prediction accuracy. Thanks to the dramatic development of artificial intelligence, some new ideas for satellite clock prediction were proposed [14]. Among these algorithms, least-squares support-vector machines (LS-SVM) is a machine learning method based on strict mathematical principles, which was developed by Suykens et al. [23]. Since it has high generalization ability, which is also benefit for low-and high-dimensional datasets, as well as fast learning speeds and low computational cost, LS-SVM is widely applied in the areas of pattern recognition, classification, and nonlinear function fitting. In this paper, an improved model based on the SAM and LS-SVM is developed to deal with BDS satellite clock prediction.
In the following sections, the related typical models for satellite clock prediction are summarized in Section 2, including QPM, SAM, and LS-SVM. Then, an improved model for BDS satellite clock prediction is proposed, and satellite-specific optimal parameters including the sample input length, regularization parameter, and kernel parameter are determined by experiments in Section 3. After preprocessing of clock offsets, clock accuracy, and static and kinematic precise point positioning (PPP) are validated and analyzed in Section 4. Finally, Section 5 concludes with a discussion.

The Improved Model
The current BDS constellation consists of GEO, IGSO, and MEO satellites, which are equipped with a rubidium atomic frequency standard (RAFS) for BDS-2, and a new-generation high-quality rubidium clock and passive hydrogen masers (PHMs) for BDS-3 (Table 1). Since the GEOs and IGSOs of BDS-3 are the test satellites without public ultra-rapid clock products, our discussion in this paper is limited to MEOs of BDS-2 (C01 to C16) and MEOs of BDS-3 (C19 to C34). Additionally, satellites C35-C37 are quite new with few observations, which are also excluded because of the poor data quality. For a better understanding of the improved model for clock prediction, traditional methods including the QPM and SAM approaches are introduced firstly. Then, the improved model is presented as a combination of the SAM and LS-SVM methods in detail.

QPM
According to the performance of satellite clocks, the systematic parts of clock offsets can be described by polynomial models with different orders. For BDS satellite clock offsets, normally a Remote Sens. 2019, 11, 2554 4 of 20 second-order polynomial model is adopted as a basic model, which is proved effectively in many contributions [11,14]. The basic QPM method is presented as: where a 0 is the clock bias correction coefficient, a 1 is the clock drift correction coefficient, and a 2 is clock drift rate correction coefficient; t 0 is the reference epoch of clock offset series; t i and clk( i ) denote the i-th epoch and its corresponding clock offset; and ε i is the residuals of the prediction model.

SAM
For GPS clock offsets, the QPM method can precisely fit most features of the satellite clock in most cases [11]. However, since the stability of BDS satellite clock is not as high as that of GPS, periodic terms have to be considered during clock offsets fitting and predicting. Taking the QPM method as a basic function, two significant periodic terms are employed to describe BDS clock offsets, called the SAM method. Hence, the function with periodic terms is expressed as [11]: where a 0 is the clock bias correction coefficient, a 1 is the clock drift correction coefficient, and a 2 is clock drift rate correction coefficient; t 0 is the reference epoch of clock offset series; t i and clk(t i ) denote the i-th epoch and its corresponding clock offset; l is the number of preiodic terms; Λ r and f r are the amplitude and frequency of periodic terms, respectively; φ r is its corresponding initial phase; and ε i is the residuals of the prediction model.

LS-SVM Model
Generally, during the LS-SVM processing, nonlinear mapping is applied to map datasets from the original space into a high-dimensional feature space, so that the nonlinear fitting problem in the original space becomes a linear fitting problem in the high-dimensional feature space. Thereby, the basic theory of LS-SVM can be described as follows.
Suppose the sample datasets are x i , y i N i=1 , x i ∈ R n . Then, the linear regression function in a high-dimensional eigenspace is shown as follow: where ϕ(x) is the nonlinear function mapping datasets from the original space into a high-dimensional feature space; w and b are the regression parameters to be estimated from the sample datasets. According to the minimal structural risk principle, the LS-SVM estimation can be transformed into an optimization problem, expressed as follows: where e i = y i − w T ϕ(x i ) + b , and γ is the regularization parameter. Then, the solution of the LS-SVM regressor can be obtained with the aid of constructing a Lagrangian function: where α i ∈ R n are the Lagrange multipliers. The conditions for optimality are shown as follows: Elimination of w and e will yield a linear system instead of a quadratic programming problem [24], which is described as: with Y = [y 1 , · · · , y N ] T , 1 N = [1, · · · , 1] T , and α = [α 1 , . . . , α N ] T ; I N is an N × N identity matrix, and Ω ∈ R N×N is the kernel matrix defined by Ω ij = ϕ(x i ) T ϕ x j = K x i , x j , which satisfies the Mercer's conditions [25].
In order to realize the LS-SVM method, the kernel function has to be determined firstly. The commonly used kernel functions include linear functions, polynomial functions, and radial basis functions (RBFs). Generally, the kernel function is selected by experience and the RBF is reported as the first choice in kinds of issues. The reason is that the RBF is suitable for both small and large quantities of datasets, especially for the data with no prior information; meanwhile, it has wider applications and possesses higher precision. Comparatively, the linear function has a linearly separable condition, and the polynomial function may lead to more parameters to be estimated when the order of polynomial is high.
The RBF kernel function can be described as follows: where σ is the kernel parameter, which determines the generalization behavior of the RBF kernel approach. Once the RBF kernel function is fixed, the potential factors that influence prediction accuracy include input sample length m, regularization parameter γ, and kernel parameter σ. Firstly, data preprocessing is employed, and a training set and a testing set are constructed based on the input sample length m. Secondly, regularization parameter γ and kernel parameter σ are determined based on the grid search algorithm [26], that is, after setting the search range by experience, an iterative computation in which k-fold cross verification is used as an evaluation criterion is done to find the tuple <γ,σ>, corresponding to the minimum squared error in the regression sample dataset. Usually, more than one tuple of <γ,σ> are obtained with a similar high accuracy level. Since a large γ value may result in overfitting and poor predictions, an optimal parameter pattern with a small γ value is the best choice [27].

The Improved Model
Since the stochastic characteristic of BDS satellite clock is sophisticated and non-stationary, which cannot be described well using the abovementioned existing models, an improved model combing the SAM and LS-SVM methods is proposed for clock prediction. It is known that LS-SVM has high non-linear mapping ability and strong robust property, which is better for approaching the remaining clock residuals after the removal of trend and periodic terms. Because the LS-SVM Remote Sens. 2019, 11, 2554 6 of 20 expression is very complex, the part of LS-SVM is depicted as a function of f LS-SVM . Hence, the improved model is proposed as: where clk(t i ), a 0 , a 1 , a 2 , t i , t 0 , l, Λ r , f r , φ r and ε i parameters are the same as those in Equation (2). Further, the process diagram of the improved model is given in Figure 1. Firstly, preprocessing for raw clock offset series is executed to get rid of outliers, and a preliminary processed clock offset series is obtained. Then, the QPM method is employed to systematically eliminate trends of clock offsets. After that, the SAM method is used to weaken periodic characteristics, and the residuals are considered as an input dataset for LS-SVM. After choosing the RBF kernel function by experience, an input sample length is determined by experiments, and the parameter pattern <γ,σ> are optimized satellite by satellite to consider satellite-specific characteristics for BDS satellite clocks. Finally, the well-trained improved model is used for clock prediction.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 21 remaining clock residuals after the removal of trend and periodic terms. Because the LS-SVM expression is very complex, the part of LS-SVM is depicted as a function of fLS-SVM. Hence, the improved model is proposed as: where clk(ti), a0, a1, a2, ti, t0, l , r Λ , fr, r φ and εi parameters are the same as those in Equation (2).
Further, the process diagram of the improved model is given in Figure 1. Firstly, preprocessing for raw clock offset series is executed to get rid of outliers, and a preliminary processed clock offset series is obtained. Then, the QPM method is employed to systematically eliminate trends of clock offsets. After that, the SAM method is used to weaken periodic characteristics, and the residuals are considered as an input dataset for LS-SVM. After choosing the RBF kernel function by experience, an input sample length is determined by experiments, and the parameter pattern <γ,σ> are optimized satellite by satellite to consider satellite-specific characteristics for BDS satellite clocks. Finally, the well-trained improved model is used for clock prediction.  Additionally, in an actual system, the machine learning algorithm uses historical data to learn the complex, nonlinear relationships between the past and the future. In the training process, a subset of the historical data (70%) are used for training and the rest of the historical data are used for testing. Then, the well-trained model is running online, and provides real-time clock offset predictions with the help of sliding window strategy. In the online operation process, the filter does not need to train itself with every dataset. When the filter fails in some cases, the system is degraded (for example, a polynomial model would succeed), but still works. Meanwhile, a new model is trained offline. Additionally, in an actual system, the machine learning algorithm uses historical data to learn the complex, nonlinear relationships between the past and the future. In the training process, a subset of the historical data (70%) are used for training and the rest of the historical data are used for testing. Then, the well-trained model is running online, and provides real-time clock offset predictions with the help of sliding window strategy. In the online operation process, the filter does not need to train itself with every dataset. When the filter fails in some cases, the system is degraded (for example, a polynomial model would succeed), but still works. Meanwhile, a new model is trained offline.

Preprocessing for Satellite Clock Offsets
It was reported that the ultra-rapid clock products of BDS-2/3 are made available by several GNSS analysis centers, for example the German Research Centre for Geosciences (GFZ) and Wuhan Remote Sens. 2019, 11, 2554 7 of 20 University(WHU), and can be used as an input in data processing. However, the raw clock time series may conceal abnormal data, such as gross errors and blunders, which should be detected before executing clock prediction. The epoch-differenced clocks are calculated in which outliers can be shown more apparently and then the outliers are detected through the median absolute deviation (MAD), which is a robust but simple tool [11,14]. Here, the epoch-differenced clock series are obtained as follows: where f i is the epoch-differenced clock series with index i, c i is the raw clock series with index i, and τ is the sample interval of the raw clock series. Meanwhile, the MAD can be expressed as follows: where m is the median of the epoch-differenced clock series, namely m = Median f i . Then, the abnormal clock offset will be deleted when f i > (m + n · MAD) is valid, where n is an integer empirically set as 3 in this processing. It follows the important principle not to delete outliers more than 5% of total data, which avoids losing a great deal of useful information [22]. Currently, this preprocessing approach is used for each BDS satellite employing the finally products of GFZ with a 30 s sample interval.
The statistic results show that more than 95% of outliers (approximately 0.27%-0.63% of total data) are detected and eliminated for each satellite.

Determination of Periodic Terms
Significant periodic characteristics are found in BDS satellite clock residuals after quadratic polynomial fitting and periodic terms changes from satellite to satellite due to different system configurations [11]. In order to identify the potential periods for each BDS satellite clock, the FFT approach is employed to analyze the clock residuals after quadratic polynomial fitting using 60-day data from day of year (DOY) 020 to 079 in 2019. The relationships between amplitude and frequency are identified for each type of satellites. Considering satellite clock stability is inversely proportional to its operation time, the newest launched satellites C02, C13, C14, C33, and C34, belonging to GEOs, IGSOs, and MEOs of BDS-2, as well as BDS-3 MEOs equipped with Rubidium (Rb) and PHM, respectively, are chosen as examples, and results are given in Figure 2. Additionally, a special periodic term of about 1.3 h has been identified for satellite C11, and the result is also shown in Figure 2. Further, the two outstanding periodic terms of each BDS satellite clock are given in Table 2. on the prediction accuracy. Therefore, the periodic terms of the clock offsets must be considered carefully.    For most of GEOs and IGSOs, the two main periodic terms, 24 h and 12 h, are identified, and the periodic terms are caused by the orbital operation and general relativistic effect. Results also show that 12 h is the most significant periodic term for most of MEOs equipped with Rb with different amplitudes. As for BDS-2, a 12 h harmonic with a smaller amplitude than for BDS-3 is noticeable. Comparatively, the clock periods of BDS-3 MEOs equipped with PHM are not apparent. The thermal insulation materials for satellites are not good enough at an early time, which leads to a changing environment for satellite clocks. As the angle of the sun above the orbital plane changes over time depending on the orbital position, the temperature of satellite clock is also varying, related to the orbital position. This is the main exterior reason why the clock offsets exhibit periodic oscillation features. Hence, the most pronounced periods of BDS satellite clocks are related to their corresponding orbit periods. Interestingly, a special periodic term of about 1.3 h has been identified with a noticeable amplitude for satellite C11, and it is different from the normal period of the MEO orbital geometric configuration. This obvious period may be caused by the hardware of the atomic clock itself, which is the main interior reason for this special periodic term and has a large impact on the prediction accuracy. Therefore, the periodic terms of the clock offsets must be considered carefully.

Input Length Determination
During the clock prediction processing, the sample dataset is usually divided into two parts, with one being a training sample dataset and the other being a testing sample dataset. Meanwhile, the performance of LS-SVM can be measured in terms of its plasticity and generalization ability, that is, the plasticity is a good fitting of the trained model with the training sample dataset, and the generalization ability is checked by the testing sample dataset, which means a good prediction. An LS-SVM with better performance should have better plasticity and generalization ability at the same time.
In order to build a suitable LS-SVM model, a sample dataset, the remaining clock residual after the removal of trend and periodic terms, is divided into a training sample dataset and a testing sample dataset in an appropriate way. Normally, 70% of sample data are employed for training and the rest are used for testing [28]. In the training processing, the training sample dataset is further divided into a number of successive multiple inputs, single output (MISO) groups, as shown in Table 3. Herein, the input data length affects LS-SVM performance significantly and should be determined by experiments firstly. Table 3. The structure of the training sample dataset, in which n is the number of training sample data and m is the dimension number of the input vector at a training iteration.

No.
Multiple Input Single Output x n-m , x n-m+1 , . . . , x n-1 x n Then, we can judge which input length is the best for each BDS satellite by checking the agreement of true (GFZ precise clock products) and prediction with different lengths of input data. In the LS-SVM construction phase, the length of input data successively increases from 5 to 200 with 5 as an interval to investigate how the length of input data affects LS-SVM performance. Figure 3 gives predicted root mean square (RMS) errors of each satellite type, and they are on the sub-nanosecond level when the prediction time is less than 3 h, regardless of the length of input data rising from 5 to 200. Regarding different satellite types, the RMS error results of GEOs and IGSOs are larger than those of MEOs. The prediction RMS values for all MEOs are within 0.4 ns for most cases, with the best one at a length of 120 as about 0.1 ns. Although MEOs of both BDS-2 and BDS-3 have a similar performance, the MEOs belonging to BDS-3 achieve a slightly better performance than those of BDS-2, especially for BDS-3 MEOs equipped with PHMs. Considering the prediction time of this experiment is three hours, the following inference is usually reasonable and the prediction accuracy will be more precise if the prediction time is not less than 3 h. Based on many experiments, it can be identified that, when the prediction time is less than 3 h, the recommended length of input data is 120 (1 h) for most satellites as a best compromise between prediction accuracy and computational complexity. the rest are used for testing [28]. In the training processing, the training sample dataset is further divided into a number of successive multiple inputs, single output (MISO) groups, as shown in Table 3. Herein, the input data length affects LS-SVM performance significantly and should be determined by experiments firstly. Then, we can judge which input length is the best for each BDS satellite by checking the agreement of true (GFZ precise clock products) and prediction with different lengths of input data. In the LS-SVM construction phase, the length of input data successively increases from 5 to 200 with 5 as an interval to investigate how the length of input data affects LS-SVM performance. Figure 3 gives predicted root mean square (RMS) errors of each satellite type, and they are on the sub-nanosecond level when the prediction time is less than 3 h, regardless of the length of input data rising from 5 to 200. Regarding different satellite types, the RMS error results of GEOs and IGSOs are larger than those of MEOs. The prediction RMS values for all MEOs are within 0.4 ns for most cases, with the best one at a length of 120 as about 0.1 ns. Although MEOs of both BDS-2 and BDS-3 have a similar performance, the MEOs belonging to BDS-3 achieve a slightly better performance than those of BDS-2, especially for BDS-3 MEOs equipped with PHMs. Considering the prediction time of this experiment is three hours, the following inference is usually reasonable and the prediction accuracy will be more precise if the prediction time is not less than 3 h. Based on many experiments, it can be identified that, when the prediction time is less than 3 h, the recommended length of input data is 120 (1 h) for most satellites as a best compromise between prediction accuracy and computational complexity.

Optimization of LS-SVM Parameters
The regularization parameter and kernel parameter are the two key parameters in the LS-SVM method. The kernel function used for BDS clock offsets is the RBF function with σ as the kernel parameter and is crucial for a better prediction. When the value of σ is small enough, i.e., σ→0, each sample data is classified into a single group, leading to over-fitting. When the value of σ is large enough, i.e., σ→∞, the distance between any points is zero, which means that all the points are divided into one group, resulting in under-fitting. Therefore, it is worth investigating the kernel parameter selection for a better balance. Meanwhile, when the regularization parameter γ is large, the training sample dataset will be fully fitted, introducing more complexity and inefficiency. When the value of γ is small, the fitting function is as simple as possible, which may not fully describe characteristic properties of the training sample dataset.
It is important to note that the regularization parameter γ and the kernel parameter σ are coupled. Therefore, the grid search algorithm with a cross-validation method is applied to find an optimal tuple <γ,σ> [29]. Firstly, the coupled simulated annealing (CSA) algorithm determines suitable starting points within the search range [e -10 ,e 10 ], based on experience. Secondly, these starting points are then given to the LS-SVM model to get the highest accuracy of prediction without over-fitting or under-fitting. The search range is gradually narrowed from a coarse mesh to a fine mesh and the threshold value is gradually decreased. Finally, the optimal tuple <γ,σ> are chosen under the principle that the regularization parameter γ is the smallest in order to decrease the computing complexity.
With the aid of optimal input length as 120 (1 h) determined in Section 3.3, the optimal tuple <γ,σ> and its corresponding prediction RMS errors for each BDS satellite are presented in Table 4.

Assessment of the Predicted Clocks
Generally, all of the satellite clock prediction accuracy can be evaluated by the deviation of prediction from its true value. However, the actual value of in-orbit clock is unknown. Therefore, a reference benchmark must be given. On the other hand, the predicted clock can be checked together with its related orbit by means of PPP. In this paper, both evaluations are carried out and the latter one is more important, as it will put forward a preliminary assessment of the positioning accuracy, which is the main goal of this research.

Clock Comparison
First of all, GFZ precise clock products are chosen as the reference datum to evaluate the accuracy of clock prediction. To demonstrate the effectiveness of the improved model (combing the SAM and LS-SVM methods), two other predicted clocks are used for comparison, i.e., clocks using the QPM method and clocks using the SAM method. Meanwhile, the latest launched satellites of each type are chosen as examples to compare their predicted RMS errors, as shown in Figure 4.
First of all, GFZ precise clock products are chosen as the reference datum to evaluate the accuracy of clock prediction. To demonstrate the effectiveness of the improved model (combing the SAM and LS-SVM methods), two other predicted clocks are used for comparison, i.e., clocks using the QPM method and clocks using the SAM method. Meanwhile, the latest launched satellites of each type are chosen as examples to compare their predicted RMS errors, as shown in Figure 4.  As shown in Figure 4, the QPM and SAM predictions have very similar trends and the SAM method performs better. Overall, their prediction errors are all within ±2.5 ns, when the prediction time is smaller than 3 h. Comparatively, the performance enhancement of the improved model is noticeable and the prediction accuracy has been further increased. No matter whether the  Figure 4, the QPM and SAM predictions have very similar trends and the SAM method performs better. Overall, their prediction errors are all within ±2.5 ns, when the prediction time is smaller than 3 h. Comparatively, the performance enhancement of the improved model is noticeable and the prediction accuracy has been further increased. No matter whether the prediction time is 1 h or 3 h, the prediction errors of all satellite clocks are kept within ±1.0 ns. Particularly, when the prediction time is smaller than 1.5 h, most of the prediction errors are within ±0.5 ns.

As shown in
In order to further evaluate the prediction accuracy and stability of the improved model, 60-day clock predictions are analyzed statistically. As abovementioned, the predicted clock errors increase when the prediction time increases. In practice, the daily clock prediction, i.e., a 24 h prediction, may not meet the reality application requirements from time to time, so that we evaluated the prediction precision with different prediction times. In this study, six typical prediction times, i.e., 1, 2, 3, 6, 12, and 24 h, are considered. The average clock differences between the predictions and GFZ precise clock products at different prediction times are calculated and listed in Table 5, which also demonstrates that the improved model is evidently superior to the others.  Table 5 shows the predicted overall RMS errors for different clock types with different prediction times. No matter whether the QPM, SAM methods, or the improved model are applied, the MEOs achieve a relatively good accuracy across all prediction times, and GEOs and IGSOs both have poor prediction accuracies, since they have poor observed structures and observations. GEO satellites have the poorest clock predictability, which may be mainly affected by their inaccurate orbit products. Among all MEO satellites, the Rb clocks belonging to BDS-2 and BDS-3 MEOs have a similar accuracy, and the PHM clocks belonging to BDS-3 MEOs have the best prediction precision, regardless of the prediction time being 1 h or 24 h. Table 5 also shows that the clock prediction strategy of the QPM method is suitable for short-term clock prediction, but its 24 h clock prediction is not as good as the short-term product. Compared with the QPM method, the SAM method achieves better accuracy for most cases. With prediction times of 1 h, 2h, and 3 h, prediction errors of the SAM method decrease, on average, by 14.1%, 17.1%, and 22.4%, respectively. Similarly, the improved model achieves better prediction precision than the SAM method for all different clock types and different prediction times. Generally, the improved model achieves an accuracy of better than 1.5 ns for the 3 h prediction time, that of 2 ns for the 6 h prediction time, that of 3.5 ns for the 12 h prediction, and that of 6 ns for the 24 h prediction time. The prediction precision for GEO satellites may be lower sometimes, but we believe the clock prediction error will be improved with a better precise orbit product. Additionally, the improved model seems more suitable for the MEO satellites, and the maximum improvements are increased by 70.9% and 50.9% for BDS-2 and BDS-3 MEO Rb clocks, respectively, for the 2 h prediction time. For the PHM clocks belonging to BDS-3 MEOs, the accuracy is still enhanced, more than 37.7% after the 24 h prediction time, when the improved model is employed.

PPP Validation
The previous section compared the predicted clocks with the reference clock products using different methods. The correctness and efficiency of the predicted clock is further validated by applying PPP processing. Considering the BDS is a newly built global system, 10 IGS Multi-GNSS Experiment (MGEX) stations, i.e., BSHM, CEDU, HARB, HKSL, KARR, KZN2, LHAZ, TONG, ULAB, and YEL2, are selected, including five stations in the Asia-Pacific region. The distribution of the stations is illustrated in Figure 5. All of these stations are capable of tracking BDS-2 and BDS-3 signals, and the observation data are publicly accessible on the Crustal Dynamics Data Information System (CDDIS) serves [12]. Additionally, the GBU orbit is used to investigate the potential positioning accuracy of PPP, considering the compatibility between orbit and clock products.  Here, it should be noted that there are various error sources that could contribute to the PPP solutions. However, as listed in Table 6, these errors can be eliminated or reduced by different methods. For example, hardware delays are the systematic errors between two GNSS code observations at the same or different frequencies, and can be corrected by the biases provided by the analysis centers. Meanwhile, the ionosphere delay is eliminated by forming an ionosphere-free combination observation, and the tropospheric zenith delay is estimated using a random walk process. Further, the IGS tracking stations are fixed based on the strict environmental criteria, and the same observations are employed to do the comparison, so the multipath effects are almost the same. In addition, phase windup correction, earth tides correction, and antenna phase center correction are also used. Therefore, we believe that the improvements observed in the experiments are from the factor of the LS-SVM method.  The observation data used in this study are collected from DOY 031 to 060 in 2019 with a 30 s interval. As stated, the precision of the 24 h clock predictions is considered very poor for the PPP processing, so the 3 h clock prediction results are applied in the PPP processing. The open-source software package RTKLIB is employed to run this experiment, and the main processing strategies and parameters for both static and kinematic PPP are listed in Table 6. In both static and kinematic PPP, only BDS observations are employed for the validation. The accuracy of PPP is assessed by comparing the positioning results with the station coordinates provided by IGS in east, north, and up directions.
Here, it should be noted that there are various error sources that could contribute to the PPP solutions. However, as listed in Table 6, these errors can be eliminated or reduced by different methods. For example, hardware delays are the systematic errors between two GNSS code observations at the same or different frequencies, and can be corrected by the biases provided by the analysis centers. Meanwhile, the ionosphere delay is eliminated by forming an ionosphere-free combination observation, and the tropospheric zenith delay is estimated using a random walk process. Further, the IGS tracking stations are fixed based on the strict environmental criteria, and the same observations are employed to do the comparison, so the multipath effects are almost the same. In addition, phase windup correction, earth tides correction, and antenna phase center correction are also used. Therefore, we believe that the improvements observed in the experiments are from the factor of the LS-SVM method. The positioning RMS errors of static PPP solutions in east, north, and up directions based on the IGS coordinates for the 10 stations from DOY 031 to 060 in 2019 are presented in Figure 6, using the QPM, SAM methods, and the improved model, respectively.
Compared with QPM, all stations are enhanced using the SAM method in the east direction and further optimized by the improving model with the largest improvement of 48.2% at BSHM station. Similar conclusions can be drawn in the north direction with the best improvement of 37.8% at LHAZ station using the improved model. Meanwhile, accuracies of all stations coordinates in the up direction have been increased, and the maximum improvement of 26.3% at LHAZ and the minimum improvement of 8.8% at HKSL are obtained when the improved model is employed. In general, the QPM method can achieve an accuracy on a decimeter level for most stations, using static PPP, while an accuracy on a sub-decimeter level is achieved for some special directions with the SAM method. When the improved model is employed, accuracies of all stations coordinates are enhanced further, and more than 70% stations obtain sub-decimeter-level positioning accuracy in horizontal direction.
The overall precisions of the static PPP are compared, and the results are presented in Table 7. It is indicated that static PPP precision with the aid of the clock predicted by the improved model have significantly enhanced, especially in the east and north directions. The overall precision of the static PPP with a 3 h clock prediction time reaches about a 10 cm horizontal precision and a vertical precision of above 23 cm. The RMS errors of PPP with the clock predicted by the improved model are enhanced by more than 20% compared to the results computed by the QPM method. With a 3 h clock prediction time, the clock error reaches 1.0 ns for most MEOs, and it can meet numerous real-time applications. With a shorter prediction length, a better static PPP precision can be expected. The positioning RMS errors of static PPP solutions in east, north, and up directions based on the IGS coordinates for the 10 stations from DOY 031 to 060 in 2019 are presented in Figure 6, using the QPM, SAM methods, and the improved model, respectively. Compared with QPM, all stations are enhanced using the SAM method in the east direction and further optimized by the improving model with the largest improvement of 48.2% at BSHM station. Similar conclusions can be drawn in the north direction with the best improvement of 37.8% at LHAZ station using the improved model. Meanwhile, accuracies of all stations coordinates in the up direction have been increased, and the maximum improvement of 26.3% at LHAZ and the  Meanwhile, in order to investigate the effect of station distribution on PPP, the stations are divided into two groups, namely global stations and Asia-Pacific stations. Compared with QPM, the improved model achieves 34.1%, 27.6%, and 18.6% improvements on average in east, north, and up components, respectively, for Asia-Pacific stations, while the respective improvements for global stations are 38.7%, 29.0%, and 20.8%. Although the enhancement rates of global stations are larger than those of Asia-Pacific stations, the absolute accuracy of positioning in the Asia-Pacific region is better. The reason is that the Asia-Pacific stations track more BDS satellites with better observations. Additionally, the convergence times of static PPP solutions are presented in Table 8, using the QPM and SAM methods, as well as the improved model. Generally, when the QPM method is employed, the convergence times of static PPP solutions for most stations are maintained at about 45 min. Compared with QPM, the convergence times of all stations are enhanced using the SAM method and further optimized by the improved model. On average, the convergence time of the improved model in the static PPP mode is less than 30 min.

Kinematic PPP
The performance of the clock predicted by the improved model is further verified using the kinematic PPP processing. In this experiment, a 3 h clock prediction product is adopted to mitigate the accumulated prediction error. Two stations from MGEX network, named CEDU and ULAB, are used, and the data are collected on DOY 045 in 2019. The positioning results of the stations are compared to the IGS reference coordinates, and the positioning error of each station is presented in Figure 7 and the corresponding statistics of the two stations are listed in Table 9. The performance of the clock predicted by the improved model is further verified using the kinematic PPP processing. In this experiment, a 3 h clock prediction product is adopted to mitigate the accumulated prediction error. Two stations from MGEX network, named CEDU and ULAB, are used, and the data are collected on DOY 045 in 2019. The positioning results of the stations are compared to the IGS reference coordinates, and the positioning error of each station is presented in Figure 7 and the corresponding statistics of the two stations are listed in Table 9.      Figure 7 intuitively demonstrates that the improved model provides a better accuracy and a shorter convergence time. As for the QPM method, position residuals converge after 1 h for CEDU and 1.2 h for ULAB. Comparatively speaking, the convergence time of the SAM method is shorter than that of the QPM method. Figure 7 also shows the improved model achieves more stable position residuals for all three components after convergence; the position error in the up direction may reach four decimeters, and the precisions of the east and north components stay at two decimeters. Meanwhile, position residuals obtained by the improved model converge rapidly within half an hour for CEDU and ULAB and are kept within an accuracy of ±0.2 m in the rest of time. Table 9 shows the RMS error values in the east, north, and up directions at each IGS tracking stations, which are calculated after convergence. It shows that the clock predicted by the improved model achieves better PPP results than that of the QPM and SAM methods. When the SAM method is used, the positioning precision improvement reaches 8.8%, 7.1%, and 6.5% on average in the east, north, and up directions, respectively. When the improved model is employed, the further improvement reaches 16.3%, 10.8%, and 18.9%, on average, in the east, north, and up directions, respectively, compared with that using the QPM method. Statistical results also show that the 3-dimension (3D) precision of the kinematic PPP results using the improved model is less than 0.4 m on average.
Further, the convergence times of kinematic PPP solutions are presented in Table 10, using the QPM and SAM methods, as well as the improved model. Generally, when the QPM method is employed, the convergence times of kinematic PPP solutions for most stations are more than 50 min. Compared with QPM, the convergence times of all stations are enhanced using the SAM method with the largest improvement of 36.4% at BSHM station. When the improved model is used, the convergence times of kinematic PPP solutions decrease significantly, and most stations only require 70% of SAM's convergence times.

Conclusions
The predicted satellite clock offsets are crucial to support real-time global satellite precise positioning applications. Although there are already several real-time precise positioning service providers, unfortunately, not all users can use the correction information due to either cost of the service and limitation of their equipment or out of the service coverage. An alternative way is to enhance the accuracy of the predicted satellite clocks for real-time precise positioning. Based on the study of existing prediction models, an improved model combing the SAM method and the LS-SVM model is proposed; by further considering satellite-specific characteristics, the parameters of the LS-SVM method are optimized satellite by satellite, including input length, regularization, and kernel parameters.
The results indicate that GEOs and IGSOs have typical periodic terms corresponding to the orbital periods of 12 h and 24 h. It also show that there exist some special periods, e.g., 1.3 h for C11, beyond our understanding that periodic terms are related to the orbit period, which may be caused by the hardware of the atomic clock itself. Experiments show that 120 (1 h data) is the best input length achieving the highest prediction accuracy for most satellites. Additionally, the optimal tuple <γ,σ> is searched by the cross-validation method based on a prediction precision comparison. Further, the improved model with satellite-specific parameters is validated, and results show that the clock accuracy of the improved model is generally better than those of the QPM and SAM methods. The overall accuracy of predicted clocks is better than ±1.0 ns for most MEOs within 3 h, which is significantly enhanced with respect to the ultra-rapid products, and the accuracy is good enough for several real-time precise applications.
The predicted clock offsets are further evaluated by applying them to PPP in both static and kinematic solutions for 10 IGS MGEX stations, including five stations in the Asia-Pacific region. The positioning accuracies obtained with the QPM and SAM methods, as well as the improved model, are compared and their convergence times are also analyzed. In the static PPP, the improved model is demonstrated to be effective, and helps to achieve more than 15.0% improvements on average for each direction, which enables sub-decimeter positioning, especially in the Asia-Pacific region. In the kinematic PPP, the improved model performs much better than the others in terms of both convergence time and positioning accuracy. The convergence time can be shortened from 1.0 h to below 0.5 h, while the positioning accuracy is enhanced by 16.3%, 10.8%, and 18.9% on average in east, north, and up directions, respectively.