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

Abstract

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-specific 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 offsets 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 offsets 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 five stations in the Asia-Pacific region. According to the practical engineering experience, 2 dm and 5 dm are defined for static and kinematic PPP, respectively, as a convergence threshold. Then, in the static PPP, the improved model is demonstrated to be effective, 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-Pacific 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.

1. 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, the standard positioning accuracy achieves five meters in both horizontal and vertical directions for Asia-Pacific users and 10 meters for global users [2].

With the development of business requirements and the processes of technological innovation, higher and higher accuracies of positioning and navigation are needed. Generally, the accuracies of satellite orbit and clock products directly decide BDS positioning performance [3,4]. With the great efforts of international GNSS service (IGS) and international GNSS monitoring and assessment service (iGMAS), several types of precise orbit and clock products for BDS satellites are provided to satisfy various positioning demands with different accuracies [5,6,7,8]. Among these products, the ultra-rapid product is generated for possible real-time applications by calculating orbits and clocks over the next 24 h. Although the orbit product with an about five centimeter uncertainty can be neglected, the precision of predicted clocks cannot satisfy the requirement of real-time precise point positioning (RT-PPP) [6,9]. Therefore, improving the accuracies of predicted clocks is of considerable interest in order to make RT-PPP workable. In this paper, we will focus on an improved clock prediction method for all BDS-2 satellites and most BDS-3 MEOs, since the MEO satellite is the backbone type of BDS-3, which has better orbit and clock products than GEOs and IGSOs, due to its favorable observation geometry. Another reason is that the qualities of the BDS-3 GEOs and IGSOs products are very poor.

Due to the solution and broadcast delay, clock prediction is indispensable for positioning services with the prediction time varying from several seconds to a few days based on different application scenarios. For the RT-PPP calculation, critical real-time and accuracy of the predicted clock must be met at the same time. Since the intrinsic physical characteristics of satellite clocks are very complex and instable, as well as varying space environment [10,11], it is more challenging to achieve the actual prediction process, and a simple linear function describing the clock offset series cannot meet the abovementioned demands. Despite such difficulties, numerous researches are carried out and many models are given for clock prediction, e.g., quadratic polynomial model (QPM), grey model (GM), spectrum analysis model (SAM), Kalman filter (KF), and neural network (NN) [12,13,14].

Usually, the first- and second-order polynomial models, as well as the GM, are widely used for the global positioning system(GPS) satellite clocks prediction. However, the limitations are obvious. It is certificated that the clock errors of QPM are accumulated rapidly with the prediction time increasing [11]. As for the traditional GM, its exponential coefficients are closely related to the number of elements, and this model requires plenty of sampling data. Although many efforts have been made, the accuracy of clock prediction is still not sufficient for real-time and high-precision applications [15]. Since the unmodeled dynamics in the residuals will be absorbed by clock parameters, including clock bias, drift, and drift rates, which leads to misestimation, based on the fast Fourier transform (FFT) method, Senior (2008) achieved four main periodic items for the GPS satellite clocks and analyzed the relationship between its amplitude and orbit period [16]. They show the incorporation of fixed-period harmonics, especially for the GPS satellite clocks, which provides a very accurate predicting model. Even though the SAM is an enhancement of the QPM with periodic terms included, it is pointed out that the periodic functions have to be determined reasonably based on a very long clock time series [17,18]. Meanwhile, researches indicate that clock parameters can be estimated by employing the KF method with a small amount of recent data. Epstein (2003) achieves an accuracy of 8–9 ns by predicting GPS clocks over six hours with the KF approach [19]. Davis (2012) developed an improved model based on the KF method along with a stochastic model, and a sub-nanosecond clock is obtained with a centimeter-level positioning accuracy [13]. However, the estimation is not reliable if a frequency jump happens. Moreover, artificial NNs are employed by scientists to predict clock offsets [11,20,21,22], e.g., back-propagation neural network (BPNN) and wavelet neural network (WNN), which have the ability to estimate non-linear time series by sample training. However, the topology structure of the WNN is hard to set up, while the BPNN involves an iterative procedure to determine appropriate weights, which requires considerable time consumption.

While most of the abovementioned achievements focus on GPS clocks, BDS satellites attract more and more attention since their dramatically development. It is reported that ultra-rapid products with 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.

2. 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. The types of BeiDou navigation satellite system (BDS) satellite clocks.

System	Satellite Type	Clock Type	Satellite Pseudo Random Number (PRN)
BDS-2	Geostationary earth orbit (GEO)	Rubidium(Rb)	C01, C02, C03, C04, C05
	Inclined geosynchronous orbit (IGSO)	Rb	C06, C07, C08, C09, C10, C13, C16
	Medium Earth orbit (MEO)	Rb	C11, C12, C14
BDS-3	GEO	Rb	C59
	IGSO	Rb	C38, C39
	MEO	Rb-II	C19, C20, C21, C22, C23, C24, C32, C33, C36, C37
	MEO	Passive hydrogen maser (PHM)	C25, C26, C27, C28, C29, C30, C34, C35

). 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.

2.1. 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 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:

$$clk\left(t_i\right)=a_0+a_1\left(t_i-t_0\right)+a_2{\left(t_i-t_0\right)}^2+{\epsilon }_i,$$

(1)

where a₀ is the clock bias correction coefficient, a₁ is the clock drift correction coefficient, and a₂ is clock drift rate correction coefficient; t₀ 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.

2.2. 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]:

$$clk\left(t_i\right)=a_0+a_1\left(t_i-t_0\right)+a_2{\left(t_i-t_0\right)}^2+{\displaystyle {\sum }_r^l{\mathsf{\Lambda }}_r}\sin \left[2\pi f_r\left(t_i-t_0\right)+{\varphi }_r\right]+{\epsilon }_i,$$

(2)

where a₀ is the clock bias correction coefficient, a₁ is the clock drift correction coefficient, and a₂ is clock drift rate correction coefficient; t₀ 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; ${\mathsf{\Lambda }}_r$ and f_r are the amplitude and frequency of periodic terms, respectively; ${\varphi }_r$ is its corresponding initial phase; and ε_i is the residuals of the prediction model.

2.3. 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 ${\left\{x_i,y_i\right\}}_{i=1}^N$, $x_i\in R^n$. Then, the linear regression function in a high-dimensional eigenspace is shown as follow:

$$y_i=w^T\phi \left(x_i\right)+b,$$

(3)

where $\phi \left(x\right)$ 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:

$$\{\begin{array}{l}\min J\left(w,b,e\right)=\frac{1}{2}{\Vert w\Vert }^2+\frac{\gamma }{2}{\displaystyle \sum _{i=1}^N{e}_i^2}\\ \mathrm{s}.\mathrm{t}.y_i=w^T\phi \left(x_i\right)+b+e_i\end{array},$$

(4)

where $e_i=y_i-\left(w^T\phi (x_i)+b\right)$, and $\gamma$ is the regularization parameter. Then, the solution of the LS-SVM regressor can be obtained with the aid of constructing a Lagrangian function:

$$L\left(w,b,e,\alpha \right)=\frac{1}{2}{\Vert w\Vert }^2+\frac{\gamma }{2}{\displaystyle \sum _{i=1}^N{e}_i^2-{\displaystyle \sum _{i=1}^N{\alpha }_i\left(w^T\phi \left(x_i\right)+b+e_i-y_i\right)}},$$

(5)

where ${\alpha }_i\in R^n$ are the Lagrange multipliers. The conditions for optimality are shown as follows:

$$\{\begin{array}{l}\frac{\partial L}{\partial w}=0\to w={\displaystyle \sum _{i=1}^N{\alpha }_i\phi \left(x_i\right)}\\ \frac{\partial L}{\partial b}=0\to {\displaystyle \sum _{i=1}^N{\alpha }_i}=0\\ \frac{\partial L}{\partial e_i}=0\to {\alpha }_i=\gamma e_i,i=1,\dots ,N\\ \frac{\partial L}{\partial {\alpha }_i}=0\to y_i=w^T\phi \left(x_i\right)+b+e_i,i=1,\dots ,N\end{array}.$$

(6)

Elimination of w and e will yield a linear system instead of a quadratic programming problem [24], which is described as:

$$\left[\begin{array}{cc}0& 1_N^T\\ 1_N& \mathsf{\Omega }+{\gamma }^{-1}{I}_N\end{array}\right]\left[\begin{array}{c}b\\ \alpha \end{array}\right]=\left[\begin{array}{l}0\\ Y\end{array}\right],$$

(7)

with $Y={\left[y_1,\cdots ,y_N\right]}^T$, $1_N={\left[1,\cdots ,1\right]}^T$, and $\alpha ={\left[{\alpha }_1,\dots ,{\alpha }_N\right]}^T$; $I_N$ is an $N\times N$ identity matrix, and $\mathsf{\Omega }\in R^{N\times N}$ is the kernel matrix defined by ${\mathsf{\Omega }}_{ij}=\phi {\left(x_i\right)}^T\phi \left(x_j\right)=K\left(x_i,x_j\right)$, 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:

$$K\left(x,x_i\right)=\exp \left(-{\Vert x-x_i\Vert }^2/{\sigma }^2\right),$$

(8)

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].

2.4. 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 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:

$$clk\left(t_i\right)=a_0+a_1\left(t_i-t_0\right)+a_2{\left(t_i-t_0\right)}^2+{\displaystyle {\sum }_r^l{\mathsf{\Lambda }}_r}\sin \left[2\pi f_r\left(t_i-t_0\right)+{\varphi }_r\right]+f_{LS-SVM}+{\epsilon }_i,$$

(9)

where clk(t_i), a₀, a₁, a₂, t_i, t₀, l, ${\mathsf{\Lambda }}_r$, f_r, ${\varphi }_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.

3. Realization of the Improved Model

3.1. 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 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:

$$f_i=\frac{c_{i+1}-c_i}{\tau },$$

(10)

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:

$$MAD=\mathrm{Median}\left\{\left|f_i-m\right|/0.6745\right\},$$

(11)

where m is the median of the epoch-differenced clock series, namely $m=\mathrm{Median}\left\{f_i\right\}$. Then, the abnormal clock offset will be deleted when $\left|f_i\right|>\left(m+n\cdot MAD\right)$ 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.

3.2. 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. BDS satellites clock periodic terms analysis.

Type	PRN	Main Period (h)	Secondary Period (h)	Type	PRN	Main Period (h)	Secondary Period (h)
BDS-2 GEO	C01	23.7449	11.8725	BDS-3 MEO Rb	C19	12.8502	16.3025
	C02	23.7449	12.0029		C20	12.7008	27.3067
	C03	23.7449	11.8725		C21	12.7008	16.3025
	C04	14.1853	19.8594		C22	12.7008	27.3067
	C05	7.9727	12.0029		C23	12.5548	15.6038
BDS-2 IGSO	C06	14.3719	23.2397		C24	12.7008	27.3067
	C07	23.7449	11.8725		C32	12.7008	6.5016
	C08	11.8725	23.7449		C33	12.7008	6.4631
	C09	11.8725	14.9626	BDS-3 MEO PHM	C25	8.6688	14.0034
	C10	11.8725	23.2397		C26	13.4848	15.8300
	C13	12.0029	24.2726		C27	13.0032	8.6688
	C16	12.0029	5.9201		C28	6.3875	4.2834
BDS-2 MEO Rb	C11	12.8502	1.2497		C29	12.7008	10.8145
	C12	12.8502	8.8086		C30	12.7008	8.5333
	C14	14.9626	9.6661		C34	6.4251	13.003

.

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.

3.3. 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. 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
1	x1, x2, …, xm	xm+1
2	x2, x3,…, xm+1	xm+2
…	…	…
n–m	xn-m, xn-m+1, …, xn-1	xn

. 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.

3.4. 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¹⁰], 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. The optimal tuple <γ,σ> and the corresponding prediction RMS error values with 95% error bands.

BDS-2	γ	σ2	RMS Error (ns)	BDS-3	γ	σ2	RMS Error (ns)
C01	55.120	10.341	0.126	C19	124.904	5.348	0.087
C02	16.663	18.447	0.209	C20	128.731	12.291	0.028
C03	228.961	15.132	0.053	C21	319.386	14.623	0.181
C04	163.164	2.645	0.395	C22	20.842	5.6403	0.085
C05	337.170	10.487	0.077	C23	13.284	2.4277	0.085
C06	35.025	11.708	0.266	C24	14.664	13.409	0.031
C07	637.809	17.668	0.188	C25	32.892	10.882	0.084
C08	239.127	13.622	0.107	C26	49.435	7.298	0.036
C09	91.678	5.928	0.148	C27	7.9522	14.315	0.037
C10	617.356	8.934	0.121	C28	71.917	6.116	0.042
C11	429.478	10.040	0.261	C29	56.923	16.666	0.126
C12	957.869	8.763	0.113	C30	18.375	4.1953	0.043
C13	186.411	12.851	0.206	C32	163.368	13.679	0.064
C14	39.434	5.137	0.237	C33	125.343	11.789	0.084
C16	11.499	8.921	0.129	C34	20.009	6.583	0.092

.

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.

4.1. 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.

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 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.

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. The predicted RMS errors for different clock types with different prediction times (ns). IM denotes the improved model.

System	Prediction time	1 h			2 h			3 h
	Type	QPM	SAM	IM	QPM	SAM	IM	QPM	SAM	IM
BDS-2	GEO	1.652	1.476	0.769	2.401	2.231	0.996	3.162	2.851	1.581
	IGSO	1.381	1.035	0.740	2.292	1.939	0.937	3.022	2.437	1.223
	MEO	1.124	0.994	0.497	1.970	1.650	0.573	2.950	2.173	0.965
BDS-3	MEO, Rb-II	0.536	0.463	0.413	0.738	0.591	0.362	1.291	0.962	0.663
	MEO, PHM	0.313	0.282	0.203	0.632	0.463	0.260	1.114	0.767	0.654
System	Prediction time	6 h			12 h			24 h
	Type	QPM	SAM	IM	QPM	SAM	IM	QPM	SAM	IM
BDS-2	GEO	3.629	3.104	2.126	4.729	4.204	3.696	8.229	7.304	6.436
	IGSO	3.440	3.053	2.072	4.240	3.953	3.337	8.026	6.738	5.877
	MEO	3.195	2.656	1.373	4.007	2.895	2.173	6.197	4.795	4.073
BDS-3	MEO, Rb-II	1.913	1.439	1.137	3.953	2.769	1.967	5.713	4.569	3.962
	MEO, PHM	1.873	1.426	1.022	3.873	2.566	1.932	5.473	4.266	3.410

, 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.

4.2. 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.

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. Strategies and parameters for PPP validation.

Solutions	Setting
Stations	10 stations from the Multi-GNSS Experiment (MGEX)
Observations	Ionospheric-free combination BDS-2: B1/B2; BDS-3: B1/B3
Sample rate	30 second
Elevation cutoff angle	7 degree
Weight	Elevation-dependent
Trop zenith delay	Initial model and random walk
Phase windup	Phase windup correction
Earth tides correction	Solid
Antenna phase left	From igs14_2062.atx
Orbit	GFZ ultra-rapid orbit
Earth rotation parameters	Products by GFZ
Differential code biases	Products by CODE
Receiver clock	White noise parameter

. 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.

4.2.1. Static PPP

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. The positioning RMS error statistics of global stations and Asia-Pacific stations in east (E), north (N), and up (U) directions of the daily static PPP solutions using a 3 h clock prediction time.

Area	Direction	QPM	SAM	IM	Enhancement
Global stations	E	0.168	0.117	0.103	38.7%
	N	0.155	0.124	0.110	29.0%
	U	0.288	0.260	0.228	20.8%
Asia-Pacific stations	E	0.138	0.110	0.091	34.1%
	N	0.130	0.105	0.094	27.6%
	U	0.231	0.204	0.188	18.6%

. 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.

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. The convergence time comparisons of the static PPP solutions (unit: min). IM denotes the improved model.

Method	BSHM	HARB	KZN2	TONG	YEL2	CEDU	HKSL	KARR	LHAZ	ULAB	Average
QPM	42	43	58	38	35	46	34	39	41	62	43.8
SAM	36	29	36	27	33	35	32	28	39	45	34.0
IM	25	27	31	26	28	27	28	24	27	36	27.9

, 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.

4.2.2. 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 positioning RMS errors of the IGS tracking stations in north (N), east (E), and up (U) directions using the kinematic PPP solutions after the initialization phase with a 5 min sample interval (unit: m). IM denotes the improved model.

Component	Method	BSHM	HARB	KZN2	TONG	YEL2	CEDU	HKSL	KARR	LHAZ	ULAB	Average
E	QPM	0.352	0.315	0.425	0.385	0.418	0.331	0.395	0.322	0.408	0.464	0.381
	SAM	0.311	0.301	0.397	0.352	0.379	0.299	0.327	0.284	0.387	0.440	0.347
	IM	0.294	0.288	0.355	0.326	0.347	0.287	0.304	0.246	0.365	0.391	0.320
N	QPM	0.216	0.197	0.225	0.208	0.182	0.143	0.135	0.184	0.163	0.249	0.190
	SAM	0.195	0.177	0.202	0.183	0.168	0.137	0.128	0.179	0.158	0.240	0.176
	IM	0.189	0.173	0.199	0.172	0.162	0.135	0.126	0.170	0.143	0.226	0.169
U	QPM	0.462	0.487	0.526	0.518	0.548	0.447	0.427	0.522	0.450	0.542	0.492
	SAM	0.458	0.469	0.494	0.472	0.501	0.424	0.402	0.496	0.410	0.482	0.460
	IM	0.377	0.433	0.401	0.436	0.428	0.303	0.384	0.457	0.382	0.397	0.399

.

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. The convergence time comparisons of kinematic PPP solutions (unit: min). IM denotes the improved model.

Method	BSHM	HARB	KZN2	TONG	YEL2	CEDU	HKSL	KARR	LHAZ	ULAB	Average
QPM	55	52	68	50	48	54	46	49	45	72	53.9
SAM	35	34	47	39	41	38	37	36	39	55	40.1
IM	27	29	37	28	31	29	24	32	29	32	29.8

, 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.

5. 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.