1. INTRODUCTION
As a major component of satellite navigation technology, it is important to analyse the characteristics and bias of satellites’ on board atomic clocks. The clocks may experience relativity and gravitational redshift effects due to height and relative motion (Cacciapuoti et al., Reference Cacciapuoti, Dimarcq, Santarelli, Laurent, Lemonde, Clairon, Berthoud, Jornod, Reina, Feltham and Salomon2007). In addition, compared with ground atomic clocks, the environment of an on board atomic clock is more complex. The atomic clock is sensitive to various outside and intrinsic factors, such as temperature, humidity, collision, relative motion, etc. Therefore, atomic clocks could suffer large deviations without effective corrections (Tang et al., Reference Tang, Xu, Cao, Liu and Wang2015). It is well known that the fundamental principle of satellite navigation and positioning is the measurement of time. Therefore, high stability and accuracy of on board atomic clocks is a key issue for precise applications such as accurate navigation and positioning (Shi et al., Reference Shi, Xu, Li and Gao2015), long baseline navigation (Batista, Reference Batista2015), precise point positioning (Li and Zhao, Reference Li and Zhao2012), precise timing (Griggs et al., Reference Griggs, Kursinski and Akos2014) and so on.
Satellite clock monitoring and bias prediction for satellite navigation systems such as the Global Positioning System (GPS), Beidou, GLONASS and Galileo has been achieved by capturing real-time clock information from observing stations from wide areas. Analysis and prediction results are released after several days or weeks of processing. For example, the International Global Navigation Satellite System Service (IGS) provides clock bias and precise ephemeris products for GPS and GLONASS periodically. The International GNSS Monitoring and Assessment System (iGMAS) can provide precise orbit and clock bias products for GPS, Beidou, GLONASS and Galileo. However, these results and products are generally time-late and may not satisfy some real-time precise applications when satellite signals are lost. In addition, the prediction length is not long enough for some long-term precise applications. To date, a number of Satellite Clock Bias (SCB) analysis and prediction models or algorithms, such as the Quadratic Polynomial (QP) model (Huang et al., Reference Huang, Zhang and Xu2014; Wang et al., Reference Wang, Lv, Wang, Li and Gong2016), Spectrum Analysis (SA) (Heo et al., Reference Heo, Cho and Heo2010; Senior et al., Reference Senior, Ray and Beard2008), Grey system Model (GM) (Yuan et al., Reference Yuan, Wang, Dong, Wu and Qu2008; Zheng et al., Reference Zheng, Lu and Chen2008), Auto Regressive Integrated Moving Average (ARIMA) model (Xi et al., Reference Xi, Cai, Li, Li, Li and Deng2014; Stein and Evans, Reference Stein and Evans1990), Kalman Model (Galleani and Tavella, Reference Galleani and Tavella2010; Davis et al., Reference Davis, Bhattarai and Ziebart2012; Pratt et al., Reference Pratt, Axelrad, Larson, Lesage, Gerren and DiOrio2013), Least Square Support Vector Machine (LS-SVM) model (Liu et al., Reference Liu, Chen, Liu and Sun2013), Wavelet Neural Network (WNN) model (Ai et al., Reference Ai, Xu, Li and Xiong2016; Wang et al., Reference Wang, Lv, Qu, Li and Wang2017), etc. In addition, combination prediction methods were also proposed by combining typical multiple models (Lei et al., Reference Lei, Hu and Zhao2014; Wang et al., Reference Wang, Hu, He, Wu, Ma and Wang2011; Xu et al., Reference Xu, Zhou, Shi, Hu and Zhou2016). However, existing models or algorithms generally do not take fully into account the influences of factors such as periodic items, sampling length, stochastic items and so on, which result in low precision and usability. SCB is caused by multiple factors, which result in the various characteristics. Therefore, in order to research and design more precise and reasonable prediction methods or models for particular applications, the influences of periodic and stochastic effects as well as sampling length should be considered.
For more precise and stable prediction of SCB, we propose a fusion-based prediction method by considering various characteristics and fitted residue. First, an Empirical Mode Decomposition (EMD) algorithm is used to pre-process and decompose the observed SCB series into multiple components with various periodic characteristics. Then, existing typical prediction models are used to fit each component with different characteristics and sampling length. The fitting residue is analysed and the model with best fitting results is selected to reconstruct the SCB. Next, the synthetic residue between the original and reconstructed data (referred to as primary fitted error) is computed and fitted by multiple prediction models, and the model with best fitting results is selected. Finally, the selected models are fused to predict future SCB in the short-, mid- and long-term. The proposed fusion-based prediction method is tested and analysed using the GPS SCB from the IGS.
2. RELATED WORKS
In this section, some typical SCB prediction models (QP, SA, GM and ARIMA) will be briefly explained, using their mathematical expressions. Their performance will also be briefly summarised.
2.1. Typical prediction model
2.1.1. QP model
The multinomial model is one of the most popular models in SCB analysis. Under this scheme, QP is a more typical model, which can be expressed as
$$x_{i} =\alpha_{0} +\alpha_{1} \Delta t_{i} +\alpha_{2} \Delta t_{i}^{2} +\int_{t_{0} }^{t_{i} } {f(t)dt}$$
where x i is the observed SCB at time t i, for the SCB dealt with differential processing, it could be the first or higher order difference of SCB. α0, α1 and α2 are the fitting parameters to be computed. Δti = ti – t 0 is the difference between observed time t i and reference time t 0. 
$\int_{t_{0}}^{t_{i} } {f(t)dt}$ is the stochastic error with a normal distribution, which is denoted as e i. Generally, the stochastic term is the clock noise, which is typically composed of five noise processes, that is White Phase Modulation (WPM), Flicker Phase Modulation (FPM), White Frequency Modulation (WFM), Flicker Frequency Modulation (FFM) and Random Walk Frequency Modulation (RWFM). Thus, the noise is non- stationary.
Denote
$${\bi \alpha}=\left[\matrix{ \alpha_{0} \cr \alpha_{1} \cr \alpha_{2} }\right],\quad {\bf A}=\left[\matrix{ 1 &\Delta t_{1} &\Delta t_{1}^{2} \cr 1 &\Delta t_{2} &\Delta t_{2}^{2} \cr \cdots &\cdots &\cdots \cr 1 &\Delta t_{n} &\Delta t_{n}^{2} }\right],\quad {\bf X}=\left[\matrix{ x_{1} \cr x_{2} \cr \cdots \cr x_{n} }\right],\quad {\bf e}=\left[\matrix{ e_{1} \cr e_{2} \cr \cdots \cr e_{n} }\right],$$then
$${\bf X}={\bf A}\alpha +{\bf e}$$
where 
$\hat{\bi \alpha}$ can be estimated by the least squares method as
$$\hat{\bi {\alpha}}=({\bf A}^{{\rm T}}{\bf A})^{-1}{\bf A}^{{\rm T}}{\bf X}.$$With the estimated parameters, we can use
$$\hat{x}_{j} =\hat{\alpha}_{0} +\hat{\alpha}_{1} (t_{j} -t_{0} )+\hat{\alpha}_{2} (t_{j} -t_{0} )^{2}$$to obtain the fitted value of the observed data (j ≤ n) and predict the future data (j > n).
2.1.2. SA model
For SCB with periodic characteristics, Spectrum Analysis (SA) can describe it better. The model can be described as
$$x_{i} =\alpha_{0} +\beta_{0} \Delta t_{i} +\sum\limits_{k=1}^m {A_{k} \sin (2\pi f_{k} \Delta t_{i} +\phi_{k} )} +e_{i}$$
where α0 and β0 are constant and coefficient terms during long-term variations, respectively. m is the number of the main periodic characteristics. f k is the frequency corresponding to the periodic characteristics and A k and ϕk are amplitude and phase, respectively. e i is the residue of x i. m and f k can be determined by fast Fourier transform. Denote 
$\alpha_{k} =A_{k} \cos \phi_{k} $, 
$\beta_{k} =A_{k} \sin \phi_{k}$, we then have
$$x_{i} =\alpha_{0} +\beta_{0} \Delta t_{i} +\sum\limits_{k=1}^m \lpar \alpha_{k} \sin (2\pi f_{k} \Delta t_{i} )+\beta_{k} \cos (2\pi f_{k} \Delta t_{i})\rpar +e_{i}$$Denote
$${\bf H}=\left[\matrix{ 1 &\Delta t_{1} &\sin (2\pi f_{1} \Delta t_{1}) &\cos (2\pi f_{1} \Delta t_{1}) && &\sin (2\pi f_{m} \Delta t_{1}) &\cos (2\pi f_{m} \Delta t_{1} ) \cr 1 &\Delta t_{2} &\sin (2\pi f_{1} \Delta t_{2}) &\cos (2\pi f_{1} \Delta t_{2}) && &\sin (2\pi f_{m} \Delta t_{2}) &\cos (2\pi f_{m} \Delta t_{2} ) \cr \cdots &\cdots &\cdots &\cdots &\cdots &\cdots &\cdots &\cdots \cr 1 &\Delta t_{n} &\sin (2\pi f_{1} \Delta t_{n}) &\cos (2\pi f_{1} \Delta t_{n}) &&&\sin (2\pi f_{m} \Delta t_{n}) &\cos (2\pi f_{m} \Delta t_{n}) }\right]^{{\rm T}},$$
${\bf a}=\left[ \matrix{ \alpha_{0} &\beta_{0} &\alpha_{1} &\beta_{1} &\cdots &\cdots &\alpha_{k} &\beta_{k} }\right]^{{\rm T}}$, 
${\bf e}=\left[ \matrix{ e_{1} &e_{2} &\cdots &e_{n} }\right]^{{\rm T}}$, then, there is
$${\bf X}={\bf Ha}+{\bf e}$$where â can be estimated by the least squares method as
$$\hat{{\bf a}}=({\bf H}^{{\rm T}}{\bf H})^{-1}{\bf H}^{{\rm T}}{\bf X}.$$With the estimated parameters, we can use
$$\hat{x}_{j} =\hat{\alpha}_{0} +\hat{\beta}_{0} \Delta t_{j} +\sum\limits_{k=1}^m \lpar \hat{\alpha}_{k} \sin (2\pi f_{k} \Delta t_{j} )+\hat{\beta}_{k} \cos (2\pi f_{k} \Delta t_{j})\rpar $$to obtain the fitted value of the observed data (j ≤ n) and predict the future data (j > n).
2.1.3. Grey system Model (GM)
The Grey system Model (GM) is a prediction system with incomplete definite information, that is, part of the information is known and part unknown, and some relationships between various factors are generally uncertain. The principle is: first, process the original data series with accumulated additive or subtractive factors to generate a new data series. This process decreases the randomisation of the original data and results in some regular characteristics. Then, establish a model based on differential equations to predict the future data. On the basis of this principle, multiple types of GM models have been presented by setting various parameters. Among these models, GM(1, 1), which is constructed by a first order differential equation with only one variable, is the most popular.
We denote the observed data series as 
$\{x_{i}^{(0)} ,i=1,2,\ldots n\}$ and the corresponding times are 
$\{t_{i}^{(0)} ,i=1,2,\ldots n\}$. Next we calculate the cumulative sum of the observed data series and the resulting series is 
$\{x_{i}^{(1)} ,i=1,2,\ldots n\}$. Then, for first order GM 
${dx^{(1)} \over dt}+ax^{(1)}=u$, the differential equation after derivative discretisation can be expressed as follows
$${\bf X}={\bf AU}$$
where 
${\bf X}=\left[\matrix{ x_{2}^{(0)} \cr x_{3}^{(0)} \cr \cr x_{n}^{(0)} }\right], {\bf A}=\left[\matrix{ -[x_{1}^{(1)} +x_{2}^{(1)} ]/2 \cr -[x_{2}^{(1)} +x_{3}^{(1)} ]/2 \cr \cdots \cr -[x_{n-1}^{(1)} +x_{n}^{()} ]/2 }\quad \matrix{ 1 \cr 1 \cr \cdots \cr 1 }\right], {\bf U}=\left[\matrix{ a \cr u \cr }\right]$.
Then, the value Û can be estimated by the least squares method as
$$\hat{{\bf U}}=({\bf A}^{{\rm T}}{\bf A})^{-1}({\bf A}^{{\rm T}}{\bf X})=\left[\matrix{ \hat{a} \cr \hat{u} }\right].$$With the estimated parameters, we can use
$$\hat{x}_{j}^{(0)} =\lpar 1-e^{\hat{a}}\rpar \left[x_{1}^{(0)} -{\hat{u} \over \hat{a}}\right]e^{-\hat{a}(j-1)}$$to obtain the fitted value of the observed data (j ≤ n) and predict the future data (j > n).
2.1.4. Auto Regressive Integrated Moving Average (ARIMA) model
ARIMA is derived from the Auto Regression Moving Average (ARMA) model by pre-processing the original data series with d-th order difference, and can be denoted as ARIMA(p, d, q). When d = 0, ARIMA(p, d, q) is ARMA(p, q). ARMA is a popular model used to fit stable sequences and includes Auto Regression (AR), Moving Average (MA) and ARMA. For most non-stable sequences such as SCB series, it always becomes stable after dealing with first or higher order differences. Therefore, ARIMA can be used to fit and predict SCB.
Denote the SCB series after d-th order difference as 
$\{x_{i}^{(-d)} ,i=1,2,\ldots n\}$, then, ARIMA(p, d, q) can be expressed as
$$x_{i}^{(-d)} =\phi_{1} x_{i-1}^{(-d)} + \cdots \phi_{p} x_{i-p}^{(-d)} +\varepsilon_{i} -\theta_{1} \varepsilon_{i-1} -\cdots -\theta_{q} \varepsilon_{i-q}$$
where ϕ1,…,ϕp are parameters of auto regression, θ1,…,θq are parameters of a moving average and p and q are orders of auto regression and a moving average, respectively. εi ,εi-1 ,… ,εi-q are white noise with normal distribution 
$\lpar 0,\sigma_{\varepsilon }^{2}\rpar $. When p = 0 or q = 0, the model is MA(q) or AR(p), respectively.
The construction and application of ARIMA includes two parts: estimation of parameters and determination of model orders. The former part is a complex procedure. The popular methods used to estimate parameters are moment estimation, maximum likelihood estimation, or least squares estimation. Moment estimation is a rough method with lower computation complexity, but the accuracy is lower than the other two methods. The least squares method is a special case of maximum likelihood method with some specific conditions. Performance of the least squares method is the best, but the procedure is more complex. For a data series, there may be multiple groups of p and q satisfying the requirement. The model should be optimised by determining the best p and q. At present, Akaike Information Criterion (AIC) norm and Bayesian Information Criterion (BIC) norm are two widely used methods for parameter optimisation.
The performance of the above four typical models have been researched and analysed (Wang et al., Reference Wang, Lv, Wang, Li and Gong2016; Senior et al., Reference Senior, Ray and Beard2008; Yuan et al., Reference Yuan, Wang, Dong, Wu and Qu2008; Xi et al., Reference Xi, Cai, Li, Li, Li and Deng2014; Wang et al., Reference Wang, Lv, Gong, Zhou and Wang2015), and can be described as follows.
(1) For short-term prediction, all of the four models perform well. However, SA considers the periodic variation, while the exponential coefficients of GM(1, 1) are related to the number of elements in the observed SCB series. Sufficient sampling data should be provided, otherwise, prediction precision cannot be guaranteed.
(2) For mid- and long-term prediction, SA and GM(1, 1) will perform well with sufficient sampling data. QP and ARIMA(p, 1, q) are not suitable for mid- and long-term prediction because of error accumulation.
(3) For periodic characteristics analysis, SA is more suitable. However, SA performs better only when there is sufficient sampling data.
(4) As to satellite clock types, ARMA(p, 1, q) is more easily affected by clock type and prediction condition. The results from Rubidium atomic clocks will be more accurate than those of Caesium atomic clocks (Heo et al., 2010; Wang et al., Reference Wang, Lv, Gong, Zhou and Wang2015).
2.2. Existing combination prediction method
For a particular prediction model, the performance can be influenced by various factors, such as sampling length, periodic characteristics, prediction length and type of atomic clock. For example, some models are appropriate for short-term prediction, while some are appropriate for mid- and long-term prediction. Some models are appropriate for periodic items, while some are appropriate for stochastic ones. Therefore, for particular applications, a combination prediction model that is superior to a single model can be obtained by combining multiple models, making the best use of the advantages of each. Existing combination prediction models can be divided into two categories: one is to combine the results of each model by weighted linear or nonlinear combination, such as the algorithms proposed by Wang et al. (Reference Wang, Hu, He, Wu, Ma and Wang2011) and Wang (Reference Wang2010). The other category is to select models for different components, and then superimpose them as the final prediction results. For example, Xu et al. (Reference Xu, Zhou, Shi, Hu and Zhou2016) combined QP and AR. Xu et al. (Reference Xu, Wang and Yang2013a; Reference Xu, Wang and Yang2013b) combined QP and a functional network. Lei et al. (Reference Lei, Hu and Zhao2014) combined QP and a least squares support vector machine. The principle of these two categories are illustrated in Figure 1, the left panel denotes the former category, and the right panel denotes the latter.
Principle of existing combination prediction model for SCB.
It can be seen from Figure 1 and existing research (Xu et al., Reference Xu, Zhou, Shi, Hu and Zhou2016; Reference Xu, Wang and Yang2013a; Reference Xu, Wang and Yang2013b; Lei et al., Reference Lei, Hu and Zhao2014) that the combination methods generally did not take fully into account the characteristics of the single models under different conditions, as well as various characteristics of observed SCB. Therefore, it has been difficult to obtain a combination prediction model with a superior and more stable performance under a variety of conditions.
3. FUSION-BASED SCB PREDICTION METHOD
First, this section analyses the characteristics of SCB. Then, an appropriate pre-process method will be presented to decompose an SCB series into multiple components with various characteristics. Finally, a fusion-based prediction method will be presented.
3.1. Characteristic analysis of the SCB
Due to influences of intrinsic physical characteristics and outside environments, the SCB series will show non-stable characteristics (Heo et al., Reference Heo, Cho and Heo2010). In addition, the on board atomic clocks are generally caesium or rubidium atomic clocks, and the frequency drift over an extended period of time is usually nonlinear (Galleani et al., Reference Galleani, Sacerdote, Tavella and Zucca2003; Zucca and Tavella, Reference Zucca and Tavella2015). The factors above will influence stability and accuracy of satellite clocks, and lead to clock offset. Generally, the offset consists of systematic and stochastic offsets.
The systematic offset is caused by initial phase offset, frequency offset and frequency drift after a long period of working. Initial phase offset is a constant. The intrinsic frequencies are different for different atomic clocks. The difference will increase linearly as time elapses, which results in a frequency offset. Frequency drift is the intrinsic characteristic of all standard atomic frequencies and will increase non-linearly as time elapses. Frequency offset and drift display as trend characteristics in SCB, which leads to nonlinear characteristics after a long period. In addition, on board atomic clocks are influenced by multiple factors such as temperature, humidity, vibration, radiation and gravitation. Accompanied with the rotation and revolution of the Earth, and periodic variation of the orbit, these factors lead to the periodic variation characteristics of atomic clocks. For bias with this characteristic, the periodic variation function is generally used for better fitting.
The stochastic offset is caused by noise inside the electronic equipment. Although the physical mechanism of the noise is not presently clear, it can be described as multiple independent energy spectrum noises. There are five modulations of noise, namely random walk frequency modulation, flicker frequency modulation, white frequency modulation, flicker phase modulation and white phase modulation. In addition, GPS rubidium atomic clocks are also influenced by very low frequency noise flicker walk frequency modulation and random run frequency modulation. The superposition of these noises gives a stochastic model of SCB. Research (Guo, Reference Guo2006) shows that the energy spectrum of most stochastic noise emanates at a low Fourier frequency. In other words, the noise will show non-stationary characteristics after long-term working, which leads to non-stationary SCB.
In summary, SCB is generally the superposition of multiple types of data such as periodic data, trend data, and stochastic data. At the same time, SCB shows non-stationary and nonlinear characteristics.
3.2. EMD and the usability for SCB analysis
EMD was proposed in 1998 by Huang et al. (Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998) of the National Aeronautics and Space Administration (NASA). EMD aims to provide a new adaptive processing method for time-frequency signals. The principle is firstly to decompose a time series signal as finite number of Intrinsic Mode Functions (IMFs) with different time-scales, followed by processing each IMF with a Hilbert transform to obtain the time-frequency spectrogram. The signal can be reconstructed by superimposing the IMFs.
The IMFs must satisfy the following requirements:
(1) During the whole time range, the number of the local extreme points must be equal to that of the zero-crossing point for the function, or differ by one only.
(2) At any time, an average of the local maximum value envelope (upper envelope solid) and the local minimum value envelope (lower envelope solid) must be zero.
EMD is excellent in processing nonlinear and non-stationary signals and for satellite clocks, can reflect the characteristics of SCB more exactly.
First, the SCB series is pre-processed, and decomposed by EMD into multiple IMFs with different characteristics. Then, each IMF is analysed and the related future IMFs are predicted. Finally, the residue is analysed and the predicted SCB are obtained by reconstruction. For the IMFs, when the periodicity is not obvious, the data may vary like a stochastic sequence. Then, IMFs with this property will be added as one component for processing. For example, the minimum satellite period is about 1.5 hours, and there will be six clock biases when the sampling interval is 15 minutes. If an IMF fluctuates frequently within six consecutive datasets, then, it may be more appropriate to consider it as a stochastic sequence. Note that this paper aims at selecting an appropriate prediction model for SCB with various characteristics. Therefore, when there are two or more IMFs performing stochastic characteristics, they will be added as one component for further analysis and prediction. Other IMFs stay unchanged.
The procedure of SCB decomposition and reconstruction based on EMD is illustrated in Figure 2. Firstly, select n prediction models {fj (x), j ∈ [1,n]} for the fitting and prediction of n IMFs {imfi, i ∈ [1,n]}. Then, select an appropriate model to fit and predict the residue e 0 between reconstructed signal 
$\hat{x}^{\prime}_{0} $ and the original signal. The estimated original signal can be obtained by adding the results ê 0 to 
$\hat{{x}}^{\prime}_{0}$.
Procedure of fusion and EMD-based SCB prediction.
3.3. Construction of fusion-based prediction model and SCB prediction
SCB is a superposition of data series with various characteristics. At the same time, performance of different prediction models varies for data with different characteristics, and the accuracy and stability of a single model are usually poor. Therefore, for SCB prediction, an appropriate data analysis algorithm such as EMD should first be used to decompose the observed data into multiple components (such as IMFs). Then, the best prediction model should be selected to fit and estimate each component. Finally, the most accurate and stable prediction results can be obtained by fusing each component predicted by the selected models. According to this principle, a fusion-based SCB prediction method is proposed by considering SCB characteristics and fitted residue. The proposed method includes two stages: the first is construction of the fusion-based prediction model, and the second is prediction of future SCB. The procedure of model construction is illustrated in Figure 3.
Framework of the proposed fusion-based SCB prediction method.
The framework in Figure 3 is described as follows:
(1) Determine the length of observed SCB according to future SCB to be predicted. If the length of SCB to be predicted is p, then the previous observed SCB with length 1, 2, …, p will be used for analysis. The i-th data is denoted as x i.
(2) Characteristic analysis and data decomposition. For the previous observed SCB with length 1, 2, …, p, pre-process it with first order difference and decompose it by methods such as EMD. For observed data of each length, a set {xi,j, j ∈ [1,n]} with n components such as IMFs can be obtained, where
$x_{i} =\sum_{j=1}^{n} x_{i,j}$.(3) Data component fitting and prediction. Analyse the characteristics of each component, such as periodic, trend and stochastic characteristics. Select the prediction model set for each component, and denote them as {fi,j, j ∈ [1,m 1]}, {fi,j, j ∈ [1,m 2 ]}, …, {fi,j, j ∈ [1,mn]}, respectively. Then, fit the previously observed SCB with length 1, 2, …, p and predict the future corresponding components.
(4) Prediction model selection and primary fitted residue computation. According to the above fitting results, take the prediction model with the minimum fitted variance as the selected model for the specific length of the observed data. The fitted components are denoted as
$\hat{x}_{i,n} = \hat{x}_{i,n}^{(f_{n,j})}\vert _{\min \lpar var\lpar x_{i,n} -\hat{x}_{i,n}^{(f_{n,j})}\rpar \rpar }$. The primary fitted results of the observed data are obtained by superimposing the fitted results of each component, that is $\hat{x}_{i}^{\prime} =\sum_{j=1}^{n} \hat{x}_{i,j}$. The primary fitted residue is calculated with $\hat{e}_{i}^{\prime} =x_{i} -x_{i}^{\prime}$.(5) Residue prediction and fusion model determination. For the primary fitted residue of observed data of each length, analyse the characteristics and select the prediction model with minimum fitted variance. The estimation result of the primary fitted residue is
$\hat{e}_{i} = \hat{e}^{\prime}_{i}\vert _{\min \lpar var\lpar e^{\prime}_{i} -\hat{e}_{i}^{\prime (f_{j})}\rpar \rpar }$. Finally, take the observed data of the length with the minimum synthetic fitted variance as the reference for SCB prediction.
After the length of observed data is determined, the fusion-based prediction model can be determined. This consists of the selected prediction models for each component and the primary fitted residue. Then, the expected future SCB is calculated by 
$\hat{x}_{i} =\hat{x}^{\prime}_{i} +\hat{e}_{i}$.
4. APPLICATION EXAMPLE AND EXPERIMENTAL RESULTS ANALYSIS
This section will present an applied example of the proposed prediction method. The method will then be tested on an actual SCB series to analyse its performance.
4.1. Fusion-based prediction model application example
According to the framework illustrated in the last section, four typical models (QP, SA, GM(1,1) and ARIMA(p,1,q)) will be taken as examples to test the proposed fusion-based prediction method (referred to as MergM). Two existing combination prediction methods mentioned in Section 2.2 will be compared with the proposed method. One is the combination of QP and AR proposed by Xu et al. (Reference Xu, Zhou, Shi, Hu and Zhou2016) (referred to as LS+AR), the other one is the weight-based combination model proposed by Wang (Reference Wang2010) (referred to as ComT).
The actual SCB series is from GPS "final" precise SCB data published by IGS, which were sampled at 15 minute intervals from 2 November 2014 to 1 March 2015. There are a total of 120 days from the first day of GPS week 1817 to the first day of GPS week 1834 without SCB hops or interruption. The first half of the data is used for construction and parameters determination of the fusion model, and the other half is used for prediction and performance analysis. For coverage of various types of satellites and on board atomic clocks, the satellites be analysed are BLOCK IIA PRN10, BLOCK IIR PRN02, BLOCK IIR-M PRN05 and BLOCK IIF PRN27. At the same time, the corresponding clock types of PRN02, PRN05, PRN10 and PRN27 are IIR Rb (built in 2004), IIR-M Rb (built in 2009), IIA Cs (built in 1996) and IIF Cs (built in 2013), respectively.
The Root Meat Square (RMS) and maximum absolute difference (Range) of prediction error are taken as criteria to evaluate the performance. RMS is calculated by 
$\sqrt{{1 \over n}\sum_{i=1}^{n} e_{i}^{2}}$, where n is the number of elements in the predicted SCB series, e i is the prediction error. Range is calculated by 
$e_{\max } -e_{\min } $, where e max and e min are the maximum and minimum prediction errors, respectively.
4.2. Prediction results and comparison analysis
We use four typical models (QP, SA, GM(1,1) and ARMA(p,1,q)), two existing combination models (LS+AR and ComT) and the proposed fusion-based prediction method MergM to predict the future SCB: one day (1d, short-term), three days (3d, short-term), 15 days (15d, mid-term), 30 days (30d, mid-term) and 60 days (60d, long-term), respectively. For existing models, the length of used observed data is equal to the prediction length, that is, we use the previous 1 d, 3 d, 15 d, 30 d and 60 d observed data to predict the future 1 d, 3 d, 15 d, 30 d and 60 d SCB, respectively. RMS and Range of the results are shown in Table 1, where the results in italic and bold are the best two. Taking 3 d, 30 d and 60 d predictions of PRN05 and PRN10 as examples, the prediction errors of each method are plotted in Figure 4 for detailed comparison.
Comparison of SCB prediction errors for PRN05 in the left panels and PRN10 in the right panels.
Clock bias prediction results comparison (ns).
The comparison results in Table 1 show that, compared with existing models, the proposed fusion-based prediction model possess better stability and accuracy. The results of the proposed method are the best under almost all conditions. In addition, most existing methods will accumulate errors over time, which is directly reflected in Figure 4. In comparison, the prediction results of the proposed method are more accurate and stable for all short-term, mid-term and long-term predictions. For example, the RMS for 30 d and 60 d predictions can be maintained under 36 ns and 90 ns, respectively. The results show that the proposed method possesses a better synthetic performance.
In addition, the prediction results in Table 1 for IIR Rb of PRN02 perform the best, and that for IIA Cs of PRN10 perform the worst. When the prediction length is no more than 15 days, the prediction results for rubidium are better than caesium. When prediction lengths are 30 d and 60 d, the prediction results for the modern atomic clock IIF caesium are better than for IIR-M rubidium. The results and analysis show that the prediction results of the proposed method perform better for rubidium and modern caesium atomic clocks, but perform worse for caesium atomic clocks which have been working for a comparatively long period. The reason may be that there is a larger level of uncertainty about the performance of aging atomic clocks, which indicate that the atomic clock should be replaced or upgraded. The facts are also consistent with these inferences; the caesium atomic clock IIA of PRN10 terminated its service on 16 July 2015, and was replaced by a modern caesium atomic clock in a IIF satellite on 31 October 2015.
4.3. Determination of parameters during prediction
The parameters involved in the fusion-based prediction method were recorded and analysed. The parameters include the quantity of components obtained by EMD, the selected models for each component, the primary fitted error and the length of observed data needed for optimum prediction. The time consumed from data decomposition to prediction was also recorded, as well as that by the combination prediction methods LS+AR and ComT. The experiment was conducted on an ordinary laptop (Windows 7 64bit, Intel(R) Core(TM) i5-4210U @1·70 GHz 2·40 GHz, RAM: 4 GB). The results are listed in Table 2, where f i is the prediction model determined for the i–th component.
Parameters used by the proposed method and the consumed time.
It can be seen from the results in Table 2 that, for the prediction of various lengths of future SCB, the length of observed data and quantity of components obtained by EMD are different for different satellites. This indicates that there are differences between characteristics of SCB generated by different types of on board atomic clocks. For example, for the prediction of PRN27 1 d, 3 d, 15 d, 30 d and 60 d into the future, all the lengths of observed data needed for optimum prediction are a constant 7 d, and the quantities of EMD components are a constant 7 d. This indicates that the characteristics of SCB of PRN27 can be described adequately with 7 d observed data. In comparison, for PRN05, the length of observed data needed for predicting future 1 d, 3 d, 15 d and 30 d are 8 d, while for 60 d it is 59 d. These results could provide a better reference to determine the length of observed data required when designing a prediction method.
It can also be seen from the selected models that, although the quantity of components for various satellites are different, only SA and ARIMA are selected for these components. What is more, SA is mostly selected for the last component. This is because the last component is the main one, which generally shows periodic characteristics and is less affected by stochastic factors. Therefore, SA is more appropriate than ARIMA, which coincides with the results pointed out in Section 2.1. For the primary fitted error, GM11 is generally selected for short-term and mid-term prediction, while SA is generally selected for long-term prediction. These results show that it is important to decompose and predict SCB according to the periodic or stochastic characteristics.
In addition, for future prediction of SCB no longer than 30 d, apart from the time consumed for the future 30 d prediction of PRN10, which is 66·95 seconds, all the other calculations take less than 60 seconds (1 minute). For the future 60 d prediction for all the satellites, the times consumed are all within 90 seconds (1·5 minutes). Compared with the two existing combination methods, the time consumed by the proposed method is longer. The main reason is that most of the IMFs are determined by using the ARIMA model, which is more complex and time consuming. This indicates that the complexity is influenced by associated single models. Therefore, one way to improve the proposed method and reduce the complexity is to select more time saving models with higher performance. For practical SCB prediction applications using observed data sampled at 15 minutes’ interval, the proposed fusion-based prediction method is a more appropriate selection for higher precision and stability.
In summary, the experiments covered almost all types of GPS satellites and current operational on board atomic clocks, such as BLOCK IIR Rb, BLOCK IIR-M Rb, BLOCK IIA Cs and BLOCK IIF Cs. The prediction results of the proposed method perform well for all of them, and are better for modern types of satellites and atomic clocks. The above experimental results and analysis indicate that the proposed method has good practical usability and outstanding performance.
5. CONCLUSIONS
For SCB prediction, a fusion-based framework is proposed through the analysis of data characteristics and fitted residue. The proposed framework also considers issues such as length and pre-processing of the observed data, and fusion of multiple models. We take four typical models as examples, and a prediction algorithm following the proposed framework is presented. Experimental results show that the prediction accuracy and stability of the proposed approach performs better compared with single models and existing combination models, which indicates the validity and usability of the proposed method.
Performance of the prediction model by fusing multiple models will be influenced by the single models involved to a certain extent. As the quantity and performance of the involved models increase, the performance of the proposed framework could also be improved further. The proposed method is validated by experimental analysis, and the models for components are selected by simple variance analysis. In future work, the principle of each single model will be theoretically analysed, and then, the fusion-based method with better performance will be presented by optimising the procedure of model selection and parameter setting based on rigorous theoretical analysis.
ACKNOWLEDGMENTS
This work was supported in part by grants from the National Natural Science Foundation of China (nos. 61602508 and 61379151).
