#### Abstract

High accuracy and reliable predictions of the bias of in-orbit atomic clocks are crucial to the application of satellites, while their clocks cannot transfer time information with the earth stations. It brings forward a new short-term, mid-long-term, and long-term prediction approach with the grey predicting model (GM(1, 1)) improved by the least absolute deviations (GM(1, 1)-LAD) when there are abnormal cases (larger fluctuations, jumps, and/or singular points) in SCBs. Firstly, it introduces the basic GM(1, 1) models. As the parameters of the conventional GM(1, 1) model determined by the least squares method (LSM) is not the best in these cases, leading to magnify the fitting errors at the abnormal points, the least absolute deviations (LAD) is used to optimize the conventional GM(1, 1) model. Since the objective function is a nondifferentiable characteristic, some function transformation is inducted. Then, the linear programming and the simplex method are used to solve it. Moreover, to validate the prediction performances of the improved model, six prediction experiments are performed. Compared with those of the conventional GM(1, 1) model and autoregressive moving average (ARMA) model, results indicate that (1) the improved model is more adaptable to SCBs predictions of the abnormal cases; (2) the root mean square (RMS) improvement for the improved model are 5.7%∼81.7% and 6.6%∼88.3%, respectively; (3) the maximum improvement of the pseudorange errors (PE) and mean absolute errors (MAE) for the improved model could reach up to 88.30%, 89.70%, and 87.20%, 85.30%, respectively. These results suggest that our improved method can enhance the prediction accuracy and PE for these abnormal cases in SCBs significantly and effectively and deliver a valuable insight for satellite navigation.

#### 1. Introduction

Atomic clocks are part of the Global Navigation Satellite System (GNSS) and have been placed on the board of satellites [1, 2]. These atomic clocks are generally the core of sophisticated scientific and technological systems, and their behaviour perform a direct and significant impact on the entire performances of the GNSS [3, 4]. Due to its own performances and the influence of the complex environment in the space, the clock bias needs to be compared and calibrated with the ground station’s atomic clock to maintain more accurate time information. However, it needs to use the existing clock bias information to forecast when the space-borne atomic clock cannot be compared with the ground stations [5, 6]. Also, the accuracy of positioning, navigation, and timing services (PNT) is directly determined by the accuracy of their clock bias prediction. Even a clock bias of 0.1 us can cause pseudorange errors (PE) of 30 metres approximately [7, 8]. The short-term clock bias predicted in advance is necessary to participate in navigation and positioning when communications are interrupted for short periods of time. Meanwhile, in the case of satellite malfunction for a long period of time, the mid-long-term and long-term clock bias predicted beforehand can be used to maintain the basic performances of the whole system. Therefore, prediction accuracy’s improving of satellite clock bias (SCB) is of vital importance for the autonomous navigation capability of satellites and the precision of real-time positioning [9, 10].

Recently, with the in-depth studies on a series of SCBs prediction problem, some short-term, mid-long-term, and long-term SCBs prediction approaches were proposed, e.g., the quadratic polynomial model, grey model, ARMA algorithm, least squares support vector machines (LSSVM) algorithm, and functional network [11–16]. Because the break of the clock information will tremendously affect the real-time precise point positioning (PPP), some approaches were proposed to improve the performances of PPP, such as optimal arc length identifying, diverse polynomial power forecasting, and empirical model made of a sixth-power harmonic function [17, 18]. Furthermore, Wang et al. proposed a predicted approach with a wavelet neural network model and obtained some conclusions through examples [19]. Compared with the IGU-P clock bias products, this method can better improve the prediction accuracy, but this approach did not consider the periodic errors and the frequency shift errors. Panfilo and Tavella used a mathematical method based on stochastic differential equations to forecast the behaviour of space-borne atomic clocks [20]. Concurrently, they considered respective diverse situations with conclusive and stochastic signatures and optimized prediction results of clock bias. The important problem was that the model structure was more complicated. Lu et al. brought forward a fusion-predicted algorithm based on four typical forecast approaches and demonstrated that the method could effectively improve prediction accuracy and stability through the examples; however, the components of the fusion-based method were more difficult to select [21]. Strandjord and Axelrad make use of the repeatable characteristic of clock bias changes and the potential of observed variations to forecast clock bias, and the achievements indicated that it could improve the predicted precision remarkably and effectually, but the more difficult things were that the computational complexities of this method were relatively large [22]. Nevertheless, further studies uncovered some limitations remaining in the classical models when predicting SCBs for abnormal cases (e.g., larger fluctuations, jumps, and/or singular points) that are affected by the phase or frequency of jumps, increasing instability, and/or general unstable tendencies. Such abnormal behaviour has a significant influence on the GNSS and may affect dramatically its entire performances.

A variety of new grey models, such as the grey Verhulst model, grey unbiased GM(1,1) model, and grey N_Verhulst model [23–25], have been proposed. In 2018, based on the grey system theory, a grey unbiased prediction model was developed to forecast the output of China’s shale gas, which showed a relatively good simulation result [24]. Subsequently, by introducing a new nonhomogeneous exponential function in view of the structural defect of the traditional grey Verhulst model, Zeng et al. put forward a new N_Verhulst model and applied to practice and got better results [26]. These new grey models are more suitable for the development law of the nonmonotonic system. As the trend of SCBs usually shows a monotonous change, the conventional GM(1, 1) model and improved GM(1, 1) model have been widely concerned in the field of clock bias prediction. The conventional GM(1, 1) model requires only few data sampling points (generally, more than or equal to four points can be predicted) to make predictions and has a better real-time than artificial intelligence methods, which need more data sampling points to train themselves [27]. Zheng et al. proposed an improved GM(1, 1) model, applied it to the SCB forecasting, and obtained better results than that of the conventional GM(1, 1) model, but they did not consider the SCB forecast in the case of clock bias fluctuations [28]. A discrete GM(1, 1) model based on sequence of stepwise ratio was proposed by Mei et al. [29], and they applied the improved model to predict SCBs of the BDS. They got the prediction precision of geosynchronous Earth orbit (GEO) satellite, inclined geosynchronous satellite orbit (IGSO) satellite, and medium Earth orbit (MEO) satellite increased significantly, compared to the conventional GM(1, 1) model. However, parameters of the conventional GM(1, 1) model are generally determined by the least squares method (LSM) [30], which cannot assure that the parameters are the best in predicting SCBs with larger fluctuations, jumps, and/or singular points. Aiming at these considerable defects, a new method is put forward to predict SCBs using the GM(1, 1) model improved by the least absolute deviations (LAD). During the modelling process of the improved model, the minimum sum of absolute errors is taken as the optimization goal, in order to overcome the problem (the LSM can fall into the local minimum easily). Since the objective function owns nondifferentiable characteristics and the calculation process is quite difficult, the optimization problem is transformed into a linear programming problem through some function transformation, and then, linear programming and the simplex method are used to estimate the parameters of the improved model. All algorithms are demonstrated by examples with larger fluctuations, jumps, and/or singular points in SCBs.

The rest of the work is organized as follows: in Section 2, a description of the prediction mechanism of the conventional GM(1, 1) model is briefly introduced; in Section 3, the conventional GM(1, 1) model improved by LAD is presented, including its mathematical representation, underlying prediction advantage for SCBs with larger fluctuations, jumps, and/or singular points; in Section 4, six separate prediction tests are carried out on GPS clock bias products selecting for SCBs with larger fluctuations, jumps, and/or singular points. Prediction accuracy and pseudorange errors (PE) of the new improved model are then compared with those of conventional GM(1, 1) model and the ARMA model. Simulation pictures are also provided in this section. Finally, some outcomes are drawn in Section 5.

#### 2. Grey Predicting Model

Grey predicting model is widely used in time series forecasting [31, 32]. A new regular sequence was collected from the raw irregular data sequences in order to build the GM(1, 1) model, and then, high-precision predictions of SCBs are made by utilizing the GM(1, 1) model based upon the regularized clock bias data sequences. The GM(1, 1) model could be shown in the following.

The sequences of raw clock bias data are provided aswhere denotes the count of raw clock bias.

The concrete prediction procedures could be described as follows [33–35]:(1)Conduct a cumulative generation operation. A cumulative generation operation to could be delivered as where .(2)Establish a first-power differential equation for the data sequences of generation, as shown in the following equation: where the parameters and are the parameters of equation (3). The parameter is called the development coefficient, which represents the development trend of SCBs sequences. The parameter is regarded as the grey action quantity. Its size reflects the changing relationship of SCBs sequences.(3)Set up the background values and solve for the two parameters and . The sequences can be represented as follows. where . Take the average values of as the background value. First-power differential equation (3) could be described by formula (5) after discretization. Formula (5) can be expressed in the matrix equation form as Then, the two parameters and can be estimated using the LSM as follows. The estimated value of and can be written in the matrix form as follows. where and are the estimated values of and , respectively. Substituting formula (8) into equation (3), we can get The solution of formula (9) is(4)Prediction sequence of the GM(1, 1) model is shown in formula (11).

From formula (11), we can get

From formula (13), when , is the monotone decreasing sequence; when , is the monotone increasing sequence.

From formula (14), when , is the monotone increasing sequence; when , is the monotone decreasing sequence.

Formula (11) is the simulation values and prediction values of the sequences . Therefore, when , the raw clock bias sequence can be approximately replaced with the simulation values .

From the above discussion, it can be seen that as long as the raw data sequence satisfies monotone increasing or monotone decreasing, the grey model based on the LSM can be used for prediction.

#### 3. Grey Predicting Model Improved by the Least Absolute Deviations

Estimation of model parameters is indispensable for calculations and should be determined after the initial recognition of the model structure. It is also a vital procedure closely linked to the predicted results. Traditionally, the estimation of the two essential parameters and is determined by the LSM [36, 37]. Because LAD can better decrease fitting errors of data points than the LSM when there are larger fluctuations, jumps, and/or singular points in SCBs, it is used to determine the two parameters and to enhance the forecast accuracy [38]. It can be seen from formula (7) that the LSM is used to measure the deviations between the real values and when estimating parameters. However, forecasting SCBs using the LSM to estimate the two parameters of the developing coefficient and grey action quantity will vastly amplify fitting errors when the raw SCBs have larger fluctuations, jumps, and/or singular points, likely resulting in poor prediction precision and large forecasting errors.

Therefore, it is not appropriate to employ the LSM to estimate the two parameters and for these abnormal cases in SCBs and then to predict SCBs. The LAD is utilized to measure the deviations between the real values and when estimating the two parameters of the GM(1, 1) model [39]. In other words, LAD is employed to determine optimally and for the GM(1, 1) model, which can increase the predictions accuracy of SCBs. It is superior to other models, such as the GM(1, 1) model based upon the LSM estimation, when there are larger fluctuations, jumps, and/or singular points in SCBs. Moreover, because LAD only considers the first-power of the absolute value of the deviations, the impact of SCBs in the case of larger fluctuations, jumps, and/or singular points are relatively smaller during predictions of SCBs. As a whole, the robustness of LAD is better than the LSM. According to these advantages of LAD, parameters and of the GM(1, 1) model are estimated by using LAD as shown in the following formula:

Taking advantage of LAD to determine the two parameters and and then placing the estimated values and into formula (11), SCB values can be predicted for any future time points. However, the ordinary method cannot be used to solve formula (15) because the optimization problem formula (15) is a nondifferentiable objective function.

To solve the optimization problem formula (15), the function transformation used is shown in the following formula:where .

Therefore, formula (16) can be written as

Consequently, determining the two parameters and using LAD can be attributed to the following linear programming problem:

The linear programming problem formula (18) could be solved using the simplex method, and then, the estimated values and can be obtained. The diagram of the basic process for clock bias forecast with the improved model is revealed in Figure 1. Furthermore, the detailed steps of the improved approach with LAD to forecast clock bias series are outlined as follows. Step 1: make a cumulative generation operation of the raw clock bias time series. Step 2: establish a first-power differential equation for the novel data sequences of generation. Step 3: set up the background values and estimate the development coefficient and the grey action quantity by using LAD. Step 4: introduce the linear programming and the simplex method, and then, search the optimal parameter and corresponding via the comparative analysis and complicated calculation. Step 5: based on the optimal parameters, the model is built to predict the following clock bias time series.

#### 4. Example and Analysis

This section firstly introduces the data selection used for tests and analyses them in detail. Evaluation indicators are then presented. Finally, to fully study the performances of the improved model, six separate prediction tests are carried out with the selected data, and then the prediction results are analysed in detail.

##### 4.1. Experimental Data Selection and Analysis

To study the performances of the improved model in the presence of SCBs with larger fluctuations, jumps, and/or singular points, we take the 1891–1894th GPS week SCBs products, for example. The details of products are summarized by Wang et al. [19]. We analyse the clock bias products of all satellites in orbit during this period. It found that fluctuations, jumps, and/or singular points in SCBs of the PRN24 satellite are larger during this time period, as shown in Figure 2. The right picture (Figure 2(b)) is a partial enlarged view of the left picture (Figure 2(a)). From the figure, we can see that the clock bias is relatively volatile. There are some abnormal clock bias data. Hence, clock bias data from PRN24 satellite are chosen as an example for the prediction tests.

Figure 3(a) is the corresponding frequency change in Figure 2(a). It can be seen from Figure 3(a) that the frequency variation of SCBs is relatively large, and the frequency of abnormal clock bias is also relatively more. Figure 3(b) is the change of the frequency histogram corresponding to Figure 3(a). From this figure, it can be seen that the frequency change of the clock bias does not obey the normal distribution (Figure 3(b)). There are some abnormal clock bias frequency data.

**(a)**

**(b)**

##### 4.2. The Evaluation Criterion

To fully assess the predicted performances of the three models (the ARMA model, the conventional GM(1,1) model, and the improved model), root mean square (RMS) values of the predicted errors and pseudorange errors (PE, PE = RMS × 0.3 m, one nanosecond produces approximately a PE of 0.3 m) values are used to evaluate the prediction accuracy of the three models. Due to the mean absolute errors (MAE), values can better reflect the true situation of the errors of the predicted values, and further evaluations are also conducted with MAE for the predicted performances of the three models [40, 41]. The expression of RMS and MAE can be written aswhere and represent the real and predicted clock bias at time , respectively.

As shown in formulae (19) and (20), values of RMS and MAE are got through using different models. Also, the smaller values of RMS and MAE refer to the better prediction performances of the model.

##### 4.3. Prediction Results and Analysis

In this part, we analyse the predicted performances of the improved model in short-term (12 h and 24 h), mid-long-term (3 d and 5 d), and long-term (14 d and 28 d); six separate predicted tests are carried out in our work. The conventional GM(1, 1) model and the ARMA model are chosen as a contradistinctive algorithm to analyse comparatively the predicted performances of the improved model.

In these short-term, mid-long-term, and long-term prediction tests, the improved model, conventional GM(1,1) model, and ARMA model are established with continuous time to predict the next 12 h, 24 h, 3 d, 5 d, 14 d, and 28 d of SCBs. The predicted errors of the improved model are compared correspondingly with those of the conventional GM(1, 1) model and ARMA model, which are shown in Figure 4. To evaluate the performances of the improved model, RMS, PE, and MAE values of these three models are calculated and listed in Tables 1 and 2. Meanwhile, the percentage improvement of PE and MAE is calculated for the improved model compared with those of the conventional GM(1, 1) model and the ARMA model, which are also shown in Tables 1 and 2. In addition, Figure 5 presents the boxplot of the RMS values of prediction errors and PE values using the conventional GM(1, 1) model, the ARMA model, and this improved method mentioned above. We also display the percentage improvement of MAE and PE using the improved model compared with those of the conventional GM(1, 1) model and ARMA model, which are shown in Figure 6.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(a)**

**(b)**

**(a)**

**(b)**

For convenience, the labels “Method 1,” “Method 2,” and “New Method” stand for the ARMA model, the conventional GM(1, 1) model, and this improved method mentioned above, respectively, in the named system of Figures 4–6.

From Figures 4–6 and Tables 1 and 2, we could find that short-term forecasted performances of the ARMA model and the improved model are roughly in the same precision level, while that of the improved model is better than the ARMA model slightly. But the precision of the conventional GM(1, 1) model in short-term forecast is worse than the improved model and ARMA model. Its prediction accuracy reaches a maximum of 4.16 ns. The MAE reaches 3.37 ns, and the generated PE reaches 1.25 m. Nevertheless, the average forecast errors of the improved model for 24 hours are only 2.28 ns, and the PE can reach as small as 0.68 m (less than 1 m).

As for the mid-long-term predicted performances, the ARMA model and conventional GM(1, 1) model are roughly at the same level, while the predicted performances of the improved model are obviously superior to other two models not only in the improvement of PE accuracy but also in the improvement of MAE accuracy. The PE and MAE of the improved model are 0.97 m (3 d), 0.86 m (5 d) and 2.93 ns (3 d), 2.50 ns (5 d); the accuracy of PE increases by 60.90% (3 d), 83.00% (5 d) and 65.80% (3 d), 77.10% (5 d); and the accuracy of MAE increases by 57.50% (3 d), 81.90% (5 d) and 65.00% (3 d), 78.00% (5 d), compared with those of the ARMA model and conventional GM(1, 1) model, respectively.

With respect to the long-term predicted performances, the forecast accuracy of the improved model is increased significantly both for PE and MAE, compared with those of the ARMA model and conventional GM(1, 1) model. From the statistical results in Tables 1 and 2, the improved model could boost the MAE in 89.50% (14 d), 89.70% (28 d) and 80.70% (14 d), 85.30% (28 d) compared with those of the ARMA model and conventional GM(1, 1) model, respectively. PE of the improved model is about 1.80 m (14 d) and 3.23 m (28 d), and compared with those of the ARMA model and conventional GM(1, 1) model, the percentage improvement of PE is approximately 88.30% (14 d), 87.20% (28 d) and 77.20% (14 d), 81.70% (28 d), respectively.

It could be known from the above quantitative analysis results. The improved model in this paper obtains a better performance than other two models whether in improving PE or MAE. In particular, the mid-long-term and long-term performances have improved significantly for abnormal cases in SCBs.

#### 5. Conclusions

In the conventional GM(1, 1) model prediction process, the parameter estimation is performed using LSM, which cannot guarantee that the parameters are optimal in predicting SCBs with larger fluctuations, jumps, and/or singular points. Aiming at these considerable drawbacks, a novel optimized GM(1, 1) model improved by LAD is put forward to enhance the accuracy of predicted clock bias for these abnormal cases. We analyse the improved model for the short-term, mid-long-term, and long-term SCBs prediction by examples in these cases. Quantitative comparative analysis is also made between the improved model, the conventional GM(1, 1) model, and the ARMA model for different accumulation at the time of SCBs prediction results, respectively.

From the calculation results, we can see that the improved model have a better fitting ability to predict SCBs than the conventional GM(1, 1) model and the ARMA model. For short-term clock bias prediction, the forecast performances of the improved model are approximately similar to that of the ARMA model, while the predicted performances of the conventional GM(1, 1) model is worse than the improved model and the ARMA model. As for the mid-long-term and long-term clock bias forecast, the forecast performances of the improved model show a considerable improvement both for RMS and PE compared with the other two models. PE’s maximum percentage improvement of the improved model can reach up to 81.66% and 88.31% compared with those of the conventional GM(1, 1) model and ARMA model, respectively.

In conclusion, the improved method (GM(1, 1)-LAD model) is more useful and accurate to forecast the short-term, mid-long-term, and long-term SCBs when the SCBs presented with larger fluctuations, jumps, and/or singular points. This study provides a beneficial reference for the theoretical study and practical application of short-term, mid-long-term, and long-term prediction of SCBs. At the same time, this study provides also a novel guidance aspect for satellite clock and clock bias prediction technologies for use in the construction of the GNSS.

#### Abbreviations

IGS: | International GNSS service |

GNSS: | Global Navigation Satellite System |

SCB: | Satellite clock bias |

GM(1, 1): | Grey model |

GM(1, 1)-LAD: | Grey model improved by the least absolute deviations |

LSM: | Least squares method |

RMS: | Root mean square |

ARMA: | Autoregressive moving average |

PE: | Pseudorange errors |

GPS: | Global positioning system |

GLONASS: | GLObal NAvigation Satellite System |

BDS: | BeiDou navigation satellite system |

PNT: | Positioning, navigation, and timing |

LSSVM: | Least squares support vector machines |

PPP: | Precise point positioning |

MAE: | Mean absolute errors |

GEO: | Geosynchronous Earth orbit |

IGSO: | Inclined geosynchronous satellite orbit |

MEO: | Medium Earth orbit. |

#### Data Availability

The data and source code used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors would like to thank Dr. Yulong Ge of National Time Service Center of Chinese Academy of Sciences for his careful revision. This work was supported by the Postdoctoral Science Foundation of China (Grant no. 2019M650801).