#### Abstract

The variation law of satellite clock bias (SCB) can be regarded as a grey system because the spaceborne atomic clock is very sensitive and vulnerable to many factors. GM (1,1) model is the core and foundation of the grey system, which has been highly valued and successfully applied in SCB prediction since its production. However, there are still some problems to be further studied such as the lack of stability of its prediction effect in practical application. In view of this, an improved GM (1,1) model by optimizing the initial condition has been proposed in this paper so as to increase the prediction performance. The new initial condition is obtained by the weighted combination of the latest and oldest components of the original clock bias sequence. And the weight values of these two components are acquired from a method of minimizing the sum of squares of fitting errors. We adopt GPS rapid precision SCB data provided by the International GNSS Service (IGS) for 15 mins, 30 mins, 1 h, 3 h, 6 h, 12 h, and 24 h prediction experiments. The results show that the improved GM (1,1) model is effective and feasible, and its prediction accuracy and stability are significantly better than those of the traditional GM (1,1) model, ARIMA model, and QP model, even for the SCB signal with obvious fluctuation.

#### 1. Introduction

In the global navigation satellite system (GNSS), the clock bias prediction of the spaceborne atomic clock plays an important role in maintaining the time synchronization of the satellite navigation system, optimizing the clock bias parameters of navigation message, meeting the requirements of real-time dynamic precision single point positioning, and providing the prior information required for satellite autonomous navigation [1–3]. Consequently, scholars at home and abroad have carried out a large number of researches on SCB prediction and came up with a variety of prediction models. There are many kinds of prediction models, including linear model, quadratic polynomial (QP) [4, 5], grey system model (GM) [6], auto-regressive integrated moving average model (ARIMA) [7], Kalman filter model [8], support vector machine model (SVM) [9], machine learning [10–13], model designed based on the basic principle of the neural network [14], and combined model [15]. Grey model acts as a significant role in clock bias prediction because of its simple expression and excellent prediction effect with less modelling data [16, 17]. Univariate first-order differential model GM (1,1) is an important part of the grey model, which is widely used in SCB prediction. However, after deeply analysing the modelling mechanism of the GM (1,1) model, we found that it is difficult to strictly approximate the grey differential equation with the fitted differential equation, and the smoothness of the initial sequence involved in the modelling will also affect the prediction accuracy of GM (1,1) model. Aiming at these problems existing in the prediction model of GM (1,1), many scholars have carried out massive studies on the aspects of initial sequence preprocessing, background value reconstruction, time response function optimization, and initial condition optimization.

In terms of initial condition optimization, there are three existing kinds of generation methods. In the first category, a single component of the first-order accumulated generating operation (1-AGO) sequence is used to solve the initial condition. For example, the initial condition of the traditional GM (1,1) model is generated by the oldest component of the 1-AGO sequence. However, some scholars believe that this approach violates the principle of “new information first,” so Li et al. [18] established a kind of GM (1,1) model with the latest component of the 1-AGO sequence as the initial value, which overemphasizes the importance of the latest information and completely ignores the influence of the old information. In addition, Ji and Zhong [19] believe that any component of the 1-AGO sequence can be used to generate initial conditions, and which component to select can be determined by minimizing the function of average relative error. Li et al. [20] have the same idea as Ji, but they use particle swarm optimization algorithm to obtain the optimal initial condition to generate components. The second way uses the linear combination of multiple components of the 1-AGO sequence to produce the initial condition. For instance, Wang et al. [21] proposed a new method of initial condition generation by the latest and oldest components of the 1-AGO sequence and gained the weight coefficients of these two elements by minimizing the sum of square errors. Considering that each component of the 1-AGO sequence will affect the prediction results, Xiong et al. [22], Chen and Li [23], and Ding [24] took the weight combination of each component in the 1-AGO sequence to create initial condition. The third method is to multiply or add coefficients to specific components to get initial conditions. For example, Zhao et al. [25] multiplied the oldest component of 1-AGO by a constant and then solved the constant by minimizing the objective function. And Xie [26] added a constant disturbance component to the latest component of the 1-AGO sequence to generate the initial condition. Furthermore, Madhi and Mohamed [27] gained the estimated value of the initial condition by minimizing the sum of the square errors of the reduced value and the real value.

The above methods are to substitute one or more components of the 1-AGO sequence into the time response function of the whitening equation to solve the initial conditions. Obviously, such ways emphasize the mining, utilization, and weight distribution of information involved in the calculation of initial conditions but ignore the influence of modelling parameters and expression forms on the model itself, which will result in the instability and poor prediction effect of the model. Accordingly, we propose in this paper a new approach to generate the initial condition by using the original sequence. This method focuses on the construction of the prediction model and the influence of parameters on the model. Firstly, the development coefficient and grey action quantity are estimated by 1-AGO and least square method, then the time response function of the whitened equation is restored by the first-order inverse accumulated generating operation (1-IAGO) to obtain the fitting model of the original sequence, and finally, the weighted values of the latest and oldest components of the original sequence are substituted into the fitting model to solve the initial condition to get the GM (1,1) prediction model. We use the GM (1,1) model constructed by this new method to predict the rapid precision SCB published by IGS and verify the effectiveness and superiority of the improved GM (1,1) model in SCB prediction through the comparison with QP model, ARIMA model, and traditional GM (1,1) model.

#### 2. GM (1,1) Prediction Model

GM (1,1) model is used extensively in time series prediction, where the symbol GM (1,1) indicates “first-order grey model in one variable” [27]. GM (1,1) model is a system of the exponential function, which can utilize known sequences before a certain moment as input and output wanton number of sequences behind that certain moment after a series of processing such as 1-AGO, modelling, 1-IAGO, and prediction. Its simple steps and detailed handling process are shown in Figure 1.

There are several explanations for the above flow chart.(1)The signal is a nonnegative sequence of length , represented by . is the 1-AGO sequence of , expressed by , and ， .(2) is called the background value of GM (1,1) model, and its value is solved through grey differential equation and whitened equation. That is, the background value of the model is equal to the area surrounded by and by the interval of -axis, where is the continuous signal corresponding to . However, the expression of is unknown, so the trapezoidal area circled by , and by the interval of -axis is generally used to approximate the representation of , i.e., .(3)The equation is a grey differential equation, and is the whitened equation of the GM (1,1) model. The parameters of and are regarded as development coefficient and grey action quantity, respectively, where represents the development trend of the sequence and the value of reflects the change relationship of [16].(4)If the parameter sequence is expressed as , the grey differential equation can be written in the form of matrix as , and the estimated value of can be obtained through the least square method (LSM) for . The values of and are and .(5)The solution of the whitened equation is also named time response function, and the undetermined coefficient of is called the initial condition of the GM (1,1) model. The value of can be obtained by substituting into the time response function and making . After the time response function is discretised and then 1-IAGO is implemented, the GM (1,1) prediction model is finally given as follows:

In formula (1), when , the value of is the fitted value, while , that of is the predicted value.

#### 3. Improvement of GM (1,1) Prediction Model

In this paper, the GM (1,1) prediction model is improved by adjusting the order of two steps, i.e., modelling grey model and the solution of initial condition. Besides, the data used to solve the initial conditions is the original sequence instead of the 1-AGO sequence. After the signal of is approximated and the parameters and are estimated, the continuous signal of corresponding to the raw sequence is procured with the time response function of the whitened equation, then the value of the initial condition is obtained by one or more components of the raw sequence, and finally, the grey prediction model is achieved by discretising . The improved part is shown in Figure 2.

Use the general solution of the first-order linear differential equation to approximate the time response function, so the function of can be exhibited as

Apply the first derivation of , the continuous form expression of the original sequence is indicated as follows:

Substitute and into formula (3), respectively, and multiply by the corresponding weight to gain the following:

Add formulas (4) and (5) to get the initial condition as follows:

Replace in (3) with the value in formula (6) to obtain the following:

Finally, the output of the GM (1,1) prediction system can be achieved by discretising .

As with the traditional model, we can obtain the predicted value of the improved grey model when through formula (8). However, the weight in the expression is an unknown parameter. Here, we employ the criterion of minimizing the sum of squares of errors [21] to get the value of . Create a function as follows:

Substitute (8) into (9), which gives the following:

Let and get

Obviously, there are parameter and all the components of the original sequence in the expression of and for the improved GM (1,1) model. In this way, the prediction value will not be affected by the grey action quantity for the existent of and will not be introduced by redundant noise due to the 1-AGO operation. Theoretically, this method can reduce the prediction error and improve the prediction accuracy.

Through the above analysis, the general steps of the improved GM (1,1) model for SCB signal prediction are summarized as follows: Step 1: make the value of the original sequence positive, and then draw upon and to calculate the 1-AGO sequence and the background values . Step 2: work out the parameters of and by using vector and matrix . Step 3: solve and then calculate the initial condition by weighting the latest and oldest components of the raw sequence. Step 4: calculate the value of by employing the criterion of minimizing the sum of squares of errors. Step 5: assume and , then substitute the corresponding parameters and sequences into formula (8) to work out the prediction value of the GM (1,1) model.

#### 4. Experimental Results

For the sake of confirming the effectiveness and reliability of the improved GM (1,1) prediction model, we take the SCB data to forecast and compare the results with traditional GM (1,1), QP, and ARIMA models. The root mean square error (RMS) [28] of the predicted value and the true value are used to evaluate the accuracy of the prediction model. The smaller the RMS value is, the higher the accuracy is. The calculation formula of RMS is denoted as follows:where is the predicted value of the clock bias, is the true value of the clock bias, and is the total number of epochs in the forecast duration.

In addition, the confidence interval of the error distribution of the forecast results is utilized to quantitatively analyse the stability of the forecast model. If is used to represent the width of the confidence interval, the smaller the value of is, the more stable the forecast model is. The expressions of and are as follows:

In formula (12), the specific meaning of the parameters such as , , and so on can be seen in literature [29].

##### 4.1. Experimental Data Source

The experimental data used in this article are the rapid precision clock bias from March 6 to April 2, 2022, of the GPS navigation satellite downloaded from the official website of IGS. The sampling interval of this SCB is 5 minutes. In the calculation example, all the satellites of each type with better integrity clock bias data are selected because the prediction results of SCB are closely related to the type of satellite clock [30] and the characteristics of the clock bias signal. During this period, there are 28 satellites with complete clock bias data except PRN02, 11, 20, and 30, whose types and satellite numbers are shown in Table 1. Among them, the clock bias signals of some satellites show an increasing trend and those of some others show a decreasing trend. Furthermore, the clock bias signals of most satellites change monotonically, while those of PRN 05, 08, and 21 fluctuate up and down greatly, especially for PRN 08 satellite. Figure 3 shows the clock bias variation diagram of PRN 08 satellite and its partial magnification. It is worth noting that all the clock bias data have been processed with MAD (median absolute deviation) before prediction.

##### 4.2. Clock Bias Prediction and Result Analysis

In the actual process of the forecast, we use the clock bias data of the first day (March 6, 2022) for modelling to predict 15 mins (minute), 30 mins, 1 h (hour), 3 h, 6 h, 12 h, and 24 h in the future, respectively. In order to further confirm the forecast performance of the improved GM (1,1) model, the SCB data from March 7 to April 2, 2022, are also dealt with like the first day. The prediction results are represented in Tables 2–5 and Figures 4 and 5. For convenience, the labels of “GM” and “IGM” in tables and figures stand for the traditional and the improved GM (1,1) model, respectively.

**(a)**

**(b)**

**(a)**

**(b)**

Tables 2 and 3 list the statistical results of the RMS and RMS improvement percentage of the four models under various forecast durations. The symbols “IGM-QP,” “IGM-ARIMA,” and “IGM-GM” in the tables indicate that IGM is compared with QP, ARIMA, and GM models, respectively. All the data given in the two tables are the average value of 28 satellites on the first day (Table 2) or 27 days (Table 3). Due to the limited space, Figure 3 only shows the RMS values of each satellite predicted the next 1 h and 6 h through the SCB data of the first day.

From Tables 2 and 3 and Figure 4, we can acquire the following information.

###### 4.2.1. Average Accuracy of the First Day

It can be seen from the data given in Table 2 that the prediction accuracy of IGM is the highest, followed by QP, ARIMA, and GM, no matter what the forecast duration is. Compared with the traditional GM (1,1), the accuracy improvement percentage of the improved GM (1,1) model is more than 60% within 1 h prediction time and that of the improved GM (1,1) model is close to 30% even for 24 h forecast duration.

###### 4.2.2. Average Accuracy of the 27 Days

The average prediction accuracy of 27 days of QP, ARIMA, and GM models is lower obviously than those of the first day, while that of IGM model changes very little. Compared with QP, ARIMA, and GM, the prediction accuracy of the IGM model is increased by 30.74%, 44.27%, and 58.31% when the forecast duration is 15 m and that of the IGM model is raised by 15.05%, 17.79%, and 20.90% while the forecast duration is 24 h. The improvement percentages of other forecast durations are between the minimum and maximum.

###### 4.2.3. Influence of Forecast Duration on the Accuracy of Prediction Model

In general, the average prediction accuracy of the four models decreases with the increase in forecast duration, whether for the first day or for 27 days. As shown in Tables 2 and 3, when the forecast duration raises from 15 mins to 24 h, the RMS value increases from 0.215 ns (nanoseconds) to 2.456 ns, which indicates that the accuracy is reduced by 10 times. When the forecast duration is less than 1 h, the average RMS values of IGM and QP are about 0.2 ns and those of ARIMA and GM are slightly higher. When the forecast duration is 24 h, the average RMS of IGM and QP is approximately 1.5 ns and those of ARIMA and GM are all within 2.5 ns.

###### 4.2.4. Influence of Satellite Clock Type on Accuracy of Prediction Model

Among the 28 satellites selected, there are two types of atomic clocks: cesium and rubidium. Due to the poor short-term stability of the cesium atomic clock, the accuracy of the cesium atomic clock predicted by the four models is worse than that of the most rubidium atomic clock. When the forecast duration is longer, this phenomenon is more distinct. Compared with the other three ordinary models, the prediction performance of the improved GM (1,1) model is least affected by the type of satellite clock. For some rubidium atomic clocks such as PRN 3 (IIF Rb), 23 (IIR Rb), and 27 (IIF Rb), the prediction accuracy of GM and IGM is low, while that of QP is very high, which shows that the clock bias of the three satellites displays the tendency of quadratic polynomial rather than that of the exponent. Combined with the satellite types in Table 1 and the RMS values corresponding to different types of satellites in Figure 3, we can conclude that the prediction performance of IIR-M Rb is relatively better than that of the other four types of satellite.

###### 4.2.5. Affection of Fluctuation of Clock Bias Signal on Accuracy of Prediction Model

Among the 28 satellites, the clock bias signals of PRN 08 fluctuate in high frequency and large amplitude, followed by PRN 05, 21, and those of other satellites change smoothly. It can be seen from the results in Figure 4 that the character of fluctuation of the signal has a certain impact on the accuracy of the four models. The greater the amplitude or frequency of fluctuation is, the lower the accuracy is. In contrast, the improved GM (1,1) model is the least insensitive to the fluctuation of the clock bias signal. Take the satellite of PRN 08 as an example: the RMS of the QP model, ARIMA, and traditional GM (1,1) model are 4.88 ns, 4.13 ns, and 5.47, respectively, while that of the improved GM (1,1) model is only 2.78 ns.

###### 4.2.6. Accuracy of a Few Abnormal Satellites

When the forecast duration is 6 h, the RMS values of PRN 03, 08, 23, and 27 are much larger than those of other satellites, which is shown in Figure 4, where the clock bias signal of PRN 08 fluctuates up and down, the prediction accuracy of these four models is very low. For PRN03, 23, and 27 with smooth clock bias signals, the RMS values of both the traditional and improved grey model are much higher than that of the QP model, which indicates that the clock bias signal of these three satellites shows a quadratic polynomial change trend rather than exponential.

Through the above tables, figures, and qualitative and quantitative analysis, it is indicated that the accuracy of the improved GM (1,1) model is higher than that of conventional GM (1,1), QP, and ARIMA. Stability is another indicator to evaluate the prediction performance, which can be analysed qualitatively by the forecast error chart of multiple satellites and be analysed quantitatively with the width of the confidence interval. The forecast errors of the four models are shown in Figure 5. Because of the limited space for this article, only the prediction errors of the next 1 h and 6 h predicted by clock bias of the first day are displayed here and those of other forecast durations are not presented one by one. The values of the confidence interval and its width are listed in Tables 4 and 5.

By observing the signal curve in Figure 5, we can find that the error distribution of the traditional GM (1,1) model is relatively scattered both 1 h and 6 h forecast duration but that of the improved GM (1,1) model is the most compact. When the forecast duration is 1 h (the left side of Figure 5), the forecast error curves of almost all the satellites for the IGM model are distributed between minus 0.5 ns and plus 0.5 ns and that of QP and ARIMA are also the same with IGM except one or two curves. However, most error curves of the GM model are scattered, and many error values are beyond plus or minus 1 ns. We can find out easily from Figure 4(a) the corresponding satellite numbers with large error values. With the increase in forecast duration, the error curves of each model begin to diverge slowly, and there are a few curves deviating from the central axis seriously. Relatively speaking, the error curve of the IGM model diverges most slowly, and the maximum value of the error curve far away from the central position is also significantly smaller than that of other models. In addition, the change trend of the forecast error curve of the IGM and GM is almost the same, but the change range of IGM is far less than that of GM. In accordance with the above qualitative analysis, we can conclude that the improved GM (1,1) model has the highest stability among the four models.

The confidence interval and its width of the first day for the four models are listed in Table 4. Generally, the lower limit of the confidence interval is negative and the upper limit is positive, in clock bias prediction. However, when the prediction accuracy is very low, both the lower limit and upper limit are positive or negative such as that on March 10 and 11 here. Therefore, Table 5 only lists the average value of the 25 days’ confidence interval and its width. All the values in both tables are calculated at the confidence level of 0.95.

Comparing the data of the corresponding forecast duration in the two tables, we can see that the average confidence interval width of 25 days is slightly higher than that of the first day. Furthermore, it can be seen from the data in the single table that the models whose confidence interval width values change from large to small are GM, ARIMA, QP, and IGM when the forecast duration is within 3 h, but those of the models are ARIMA, GM, QP, and IGM while the forecast duration is greater than 3 h, except that the confidence interval widths of ARIMA for 15 mins forecast duration and QP for 12 h forecast duration are slightly lower than IGM. It is further proved by the quantitative data that the stability of IGM is the highest among the four models. Compared with GM, the average prediction stability of the IGM model of 25 days for the forecast duration of 15 mins, 30 m, 1 h, 3 h, 6 h, 12 h, and 24 h has been increased by approximately (56.99, 59.72, 56.23, 53.43, 45.63, 36.1, and 19.71) %, respectively.

#### 5. Conclusion

Aiming at the instability performance of the traditional GM (1,1) model in short-term clock bias prediction, an improved method is proposed by optimizing initial conditions. In order to verify the effectiveness of the improved GM (1,1) model, we apply it to predict the GPS rapid precision clock bias, then use RMS and the confidence interval width of forecast error to evaluate the accuracy and stability of the prediction model, and finally draw the following conclusions:(1)The improved GM (1,1) model can be used to predict SCB effectively, and we can obtain better prediction result even when the clock bias signal fluctuates seriously. In addition, the improved GM (1,1) model is not sensitive to the type of satellite clock, the fluctuation characteristics of clock bias signal, and its influence is much less than that of the QP model, ARIMA model, and traditional GM (1,1) model.(2)The prediction accuracy and stability of the improved GM (1,1) model are much higher than that of the QP model, ARIMA, and traditional GM (1,1) model, especially for the short forecast duration. However, with the increase in forecast duration, the improvement rate of accuracy and stability of IGM begin to decline, compared with other models. Therefore, in future research, we can combine the method mentioned in this article with other measures, such as preprocessing the initial sequence, reconstructing the background value, and optimizing the time response function to further enhance the accuracy and stability of the GM (1,1) model.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by National Outstanding Youth Science Found Project of National Natural Science Foundation of China (Grant nos. 41804076 and 61503404) and the Science and Technology Project of the Department of Education of Jiangxi Province (Grant no. GJJ211814). The authors also thank the IGS authorities for providing the data and products for this study.