Prediction of Earth Rotation Parameters with the Use of Rapid Products from IGS, Code and GFZ Data Centres Using Arima and Kriging – A Comparison

Earth Orientation Parameters (EOP), i.e., coordinates of the Earth’s pole (PMx, PMy), celestial pole offsets (dX, dY), and UT1-UTC describe the irregularities in the Earth’s rotation. An additional parameter that is the first negative time derivative from the UT1-UTC is the Length of Day (LOD) (Modiri et al., 2020). In this contribution, a subset of EOP, Earth Rotation Parameters (ERP), consisting of Earth’s pole coordinates (PMx, PMy) and UT1-UTC (LOD) is subject to prediction.

The knowledge of EOP is essential for transformation between celestial reference frame (CRF) and terrestrial reference frame (TRF) (Gambis and Luzum, 2011), which enables precise positioning in space and on the Earth’s surface, navigation, time-keeping and positional astronomy.

EOP are determined by space geodetic techniques such as VLBI (Very Long Baseline Interferometry), DORIS (Doppler Orbitography and Radiopositioning Integrated by Satellite), SLR (Satellite Laser Ranging) and GNSS (Global Navigation Satellite System), and are finally combined into a single solution (Petit and Luzum, 2010). This combination of techniques enables to determine polar coordinates with the accuracy of app. 50 μas or better and 10 μs in case of LOD (Dick and Thaller, 2020).

Complexity of measuring and data processing involved in EOP determination forces their prediction since it is impossible to obtain them in real-time mode.

The long EOP time series are published with latency by International Earth Rotation and Reference System Service (IERS) and many individual analysis centres or services (IERS Annual Report 2018). Therefore, in operational settings, this gap must be filled with rapid products, e.g., IGS (International GNSS Service), CODE (Center for Orbit Determination in Europe) or GFZ (German Research Centre for Geosciences). This may also be done by basing on rapid products published by IERS, e.g., Daily Rapid EOP Data.

A plenty of prediction methods have been applied to EOP prediction through the years; the following list is by no means exhaustive: least-squares collocation (Hozakowski, 1989), artificial neural networks (Schuh et al., 2002; Liao et al., 2012), autocovariance prediction (Kosek, 2002), autoregression (Kosek et al., 2007), combination of least-squares and multivariate stochastic methods (Niedzielski and Kosek, 2008), wavelets and fuzzy inference systems (Akyilmaz et al., 2011), extreme learning machine (Lei et al., 2015), Gaussian process regression (Lei et al., 2015), combination of least-squares extrapolation and auto-regressive modelling of EAM (Dill et al., 2019), singular spectrum analysis (Okhotnikov and Golyandina, 2019), one-step forecasting method based on the LS+AR (Wu et al., 2019), combination of singular spectrum analysis and Copula-based analysis (Modiri et al., 2020) and Kalman filter (Nastula et al., 2020), combination of Holt–Winters algorithm and angular momenta of global surficial geophysical fluids (Luo et al., 2022).

This study is intended in ultra-short-term (up to 15 days into the future) and short-term (up to 30 days into the future) ERP predictions using geostatistical method of ordinary kriging and autoregressive integrated moving average (ARIMA) model. Rapid GNSS products EOP 14 12h from IGS, CODE and GFZ and also IERS final products (IERS EOP 14 C04 12h(IAU2000A)) are used and the results are compared. For all the time series and data centres the first ERP prediction was made for 1 January 2017 (MJD 57754.5) and the last for 15 April 2022 (MJD 59684.5) with 1 day shift.

In this contribution, rapid EOP 14 12h (IGS, CODE, and GFZ) and final EOP 14 C04 12h (IERS) daily products were used to predict ERP. The obtained predictions cover the time interval from 1 January 2017 (MJD 57754.5) to 15 April 2022 (MJD 59684.5) for all evaluated ERP time series/data centres.

Prediction procedure is based on that one described in Michalczak and Ligas (2021) and Michalczak and Ligas (2022) set of parameters such as the number of the nearest neighbors (NN) for kriging, number of observations for structure function estimation, the time span (TS) necessary to estimate parameters of a linear trend and number of periodic components. The TS parameter for PMx and PMy is equal to 6.4 yrs and for LOD 9.2 yrs; NN for PMx and PMy equals to 10 immediate neighbors and for LOD it is equal to 25.

In case of ARIMA(p, d, q) model, the best set of parameters p and q was obtained by means of corrected Akaike Information Criterion (AICc) (Akaike, 1998). The loop with AICc criterion was run over the number of p and q varying from 1 to 5; the best set of parameters was selected and then a prediction was made based on them. For each 15-day and 30-day ERP prediction, the best set of p and q was selected separately. Parameter d determines a degree of differencing to be applied in order to transform a non-stationary time series into the stationary one in the mean sense. Stationarity of each ERP time series was checked using the KPSS test. Pole coordinates time series (after removing components mentioned in 2.4.) passed the stationary test for d = 2 and for LOD time series (after removing components mentioned in 2.2. and 2.4.) d = 1.

The prediction process starts with the selection of the first day of forecast. Then, the 30 observations preceding the first day of forecast were selected (rapid products from CODE, GFZ or IGS), while further back observations were taken from the final product (IERS). Afterward, the ERP predictions were carried out with ARIMA model and kriging for 15 and 30 days ahead. Obtained results were compared with each ERP IERS final time series. The whole process is then repeated, shifting the first day of prediction to next day and repeated until the last day of the forecast is reached, i.e., 15 April 2022. The process is illustrated in Figure 4. The prediction process based on final product (IERS) did not use any rapid products (this is to compare results of prediction with and without rapid product supplementation), which means that the 30 observations preceding the first day of prediction were selected from the final product – IERS EOP 14 C04 12h.

Figure 4.

For the IGS, GFZ, and IERS, the total number of 15-day predictions for each ERP and each method amounts to 1886 and the total number of 30-day predictions is equal to 1871. Due to gaps in the CODE time series, the number of each ERP forecasts is equal to 1822 and for 30-day predictions amounts to 1807.

The ERP prediction accuracy was assessed by means of the mean absolute prediction error (MAPE) given by 9 MAPEj=1n∑i=1n| ei,j |j≔prediction day number \[\begin{matrix} MAP{{E}_{j}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathop \sum }}\,|{{e}_{i,j}}| & j:=prediction~day~number \\ \end{matrix}\] where n is a number of predictions, e_i,j = O_i,j – P_i,j is the prediction error, O stands for observed and P for predicted ERP.

Figure 5 and 6 present the values of MAPEs for ultra-short (15-day into the future) PMx and PMy prediction using ARIMA and kriging for data time series from four analysis centres.

Figure 5.

Figure 6.

PMx prediction with ARIMA model offers better accuracy up to day 13; after this day the advantage is on the side of kriging. For the first 5 days, the advantage is around 0.3 mas in favour of the ARIMA model. As the forecast day increases, the difference between ARIMA model and kriging decreases. For the 15th day, the difference is around 0.15 mas in favour of kriging. MAPEs for PMx for each analysis centres are very similar. The average difference between MAPEs for PMx prediction with ARIMA model is around 0.01 mas. In case of kriging, prediction performance is similar, but only for first 10 days. Then, MAPEs for CODE time series increase and differ from the other centres by about 0.04 mas.

PMy prediction with ARIMA model offers better accuracy up to the last day. The differences between ARIMA and kriging are around 0.23 mas in favour of ARIMA. Similarly to PMx results, the average difference between MAPEs for PMy prediction with ARIMA for four centres is around 0.01 mas with maximum value 0.03 mas. MAPEs for various centres for kriging prediction up to 8th day of prediction are nearly the same (with difference around 0.01 mas). On the 9th day of prediction differences between various centres increase to 0.02 mas and then on the last day, they reach values around 0.04 mas.

For both PMx and PMy, prediction accuracy based on final data (IERS) does not differ significantly from accuracy obtained with rapid products (CODE, GFZ and IGS) supplementing. This insignificant difference between the forecasts based on either final or rapid data applies also to 30-day predictions of PMx/PMy and to both LOD forecast horizons.

Figure 7 presents the values of MAPEs for ultra-short (15-day into the future) LOD prediction using ARIMA and kriging for time series from four analysis centres. For the first 7 days, lower values of MAPEs are obtained, with advantage around 0.004 ms, for ARIMA. After this day, the advantage is on the side of kriging with maximum difference equal to about 0.016 ms. The accuracy of predictions within kriging and ARIMA for each time series is nearly identical.

Figure 7.

Figures 8 and 9 present the values of MAPEs for short (30-day into the future) PMx and PMy prediction using ARIMA and kriging.

Figure 8.

Figure 9.

Up to 11 day of PMx prediction, MAPEs for ARIMA are smaller than those for kriging with maximum advantage around 0.5 mas. After this day, kriging reaches better accuracy of forecast. The maximum differences between ARIMA and kriging are on the 30th day and are equal to around 1.91 mas. The accuracy of predictions based on various data centres is nearly identical; the maximum differences vary from around 0.01 mas to 0.05 mas.

For the first 18 days of PMy prediction, results indicate a better prediction performance of ARIMA model. From the 19th day, the advantage is in favour of kriging and reaches maximum value about 0.3 mas. The differences between predictions based on various data centres are small, around 0.01 – 0.03 mas.

Figure 10 shows a comparison of MAPEs for 30-day LOD prediction for ARIMA and kriging. Results indicate that for the first 6 days, prediction with ARIMA is more accurate. From the 7th day onward, kriging gives a more accurate forecast. The maximum difference between the forecasting methods is about 0.026 ms for CODE (on the last day of prediction). It is worth noting that MAPEs for IGS and IERS using kriging are most similar to each other. In majority, this also holds for ARIMA.

Figure 10.

Earth Rotation Parameters predictions based on ordinary kriging and ARIMA model have been presented in this contribution. Predictions have been performed on the basis of rapid (CODE, IGS and GFZ) and final (IERS) ERP products. The quality of polar motion forecast is visibly coordinate-dependent, i.e., PMx is worse predictable than PMy. The results indicate that the ARIMA-based prediction is better for ultra-short prediction; for PMx it turned out to be 11 -days of forecast, for PMy 18 -days, and for LOD 8 -days. The maximum differences for first few days of 15-day predictions between methods are around 0.32 mas (PMx), 0.23 mas (PMy) and 0.004 ms (LOD) in favour of ARIMA. On the last day, the differences reach 0.15 mas (PMx) and 0.016 ms (LOD) with advantage to kriging. For the first few days of 30-day prediction, ARIMA gives lower MAPEs of approximately 0.5 mas (PMx), 0.3 mas (PMy) and 0.007 ms (LOD). The maximum differences of MAPEs on the last days of 30-day predictions are 1.91 mas (PMx), 0.3 mas (PMy) and 0.026 ms (LOD) with advantage to kriging method. For all ERP, the differences of MAPEs for time series from various analysis centres are not significant and vary up to maximum value of around 0.05 mas (PMx), 0.04 mas (PMy) and 0.005 ms (LOD). These values are on the level of accuracy of corresponding ERP determination, i.e., polar coordinates with the accuracy of app. 50 μas or better and 10 μs in case of LOD (Dick and Thaller, 2020). It is also worth noting that the accuracy of prediction on the basis of rapid products (30-day gap filling) is on the same level as on final product only. The conducted comparisons did not give a clear answer as to which method, i.e., ARIMA or kriging, is better for a given ERP and for a given forecast length, with the exception of ultra-short-term (15 days) prediction for PMy where ARIMA systematically behaved better. The results of this contribution are supported by a considerable number of performed predictions; therefore, all intermediate results are available at the website mentioned in the acknowledgements section.