Analysis and prediction of polar motion using MSSA method

Polar motion is the movement of the Earth's rotational axis relative to its crust, reflecting the influence of the material exchange and mass redistribution of each layer of the Earth on the Earth's rotation axis. To better analyze the temporally varying characteristics of polar motion, multi-channel singular spectrum analysis (MSSA) was used to analyze the EOP 14 C04 series released by the International Earth Rotation and Reference System Service (IERS) from 1962 to 2020, and the amplitude of the Chandler wobbles were found to fluctuate between 20 and 200 mas and decrease significantly over the last 20 years. The amplitude of annual oscillation fluctuated between 60 and 120 mas, and the long-term trend was 3.72 mas/year, moving towards N56.79 °W. To improve prediction of polar motion, the MSSA method combining linear model and autoregressive moving average model was used to predict polar motion with ahead 1 year, repeatedly. Comparing to predictions of IERS Bulletin A, the results show that the proposed method can effectively predict polar motion, and the improvement rates of polar motion prediction for 365 days into the future were approximately 50% on average.


Introduction
The rotation of the Earth not only characterizes the overall motion state of the Earth, but also reflects the coupling processes between the solid Earth and the atmosphere, oceans, mantle, and core on spatial and temporal scales (Nastula and Ponte 1999;Zotov and Bizouard 2015). Polar motion (PM) is an objective representation of the Earth's rotation, so its changes also indicate variations in various geophysical factors. Therefore, the study of PM changes is of great value and significance for understanding geophysical processes, such as the movement of mass in the crust (mantle convection, topographic coupling, electromagnetic coupling, and plate movement), the flow of Earth's surface materials (ocean circulation and tides, atmospheric motion, postglacial rebound, and glacier melting), and the external effects on the Earth (the influence of the Moon, Sun, and other celestial bodies) (Chen and Wilson 2005;Lambeck 2005;Dill and Dobslaw 2010;Jin et al. 2013;Wang et al. 2016;Adhikari and Ivins 2016).
PM is the movement of the Earth's rotational axis relative to the crust, which is usually expressed by using X and Y coordinates, and these coordinate axes point to 0° and 90°W, respectively. Long-term drift, Chandler wobbles and the annual term are important components of PM, which is of great significance in studying the temporally varying characteristics of PM (Nastula et al. 2011). The Chandler wobbles were proposed by Chandler in 1891, which is subject of ongoing research since its discovery (Zotov and Bizouard 2015;Wang et al. 2016). As time changes, the period of the Chandler wobbles fluctuates from 1.13 to 1.20 years (Schuh et al. 2001), and its amplitude varies between 100 and 200 mas (Gross 2000; Wang et al. 2016). Although current studies have shown that the Chandler wobbles are related to interactions among the atmosphere, ocean and other Earth fluids, its excitation mechanism remains to be further studied (Chen and Wilson 2005;Malkin and Miller 2010;Nastula et al. 2011). The period and amplitude of the annual term remain basically stable, and the long-term drift is approximately 3.3-3.5 mas/year (Gross 2000;Schuh et al. 2001).
To better extract the periodic terms of PM and study their time-varying characteristics, bandpass filters have typically been used. For instance, Popinski et al. (1995) used Fourier transform bandpass filter to analyze PM variations and extracted the Chandler wobbles and annual term from IERS93C01 and IERS90C04. Wavelet transform was applied to determine and analyze the amplitudes of the Chandler wobbles and annual term (Schuh et al. 2001;Liu et al. 2007). Since 2010, singular spectrum analysis (SSA) has been applied for the analysis of PM. Malkin and Miller (2010) used SSA to extract the Chandler wobbles and annual term of PM, and the results showed that this method can obtain a more detailed PM series structure than other digital filters. Wang et al. (2016) proposed a Fourier basis pursuit bandpass filtering that can effectively suppress the edge effect in the process of extracting prograde and retrograde annual wobbles. However, these methods can only analyze the time series in the X-or Y-direction individually, whereas the PM is actually an entire movement, which is only expressed by X-and Y-coordinates. Therefore, the correlation of changes in PM between the two directions should be considered in PM analysis.
Multi-channel singular spectrum analysis (MSSA) is widely applied in the fields of climatology (Luzano 2020), seismology (Oropeza and Sacchi 2011) and physical geodesy (Jin et al. 2021). This analysis method extracts the spectral characteristics of a multi-channel time series and extends the analysis of the observed dynamic characteristics from a limited single channel to multiple channels, taking into account the correlation among different measurement signals. Compared with SSA and principal component analysis (Rangelova et al. 2007), MSSA can better extract useful information (such as periodic term and long-term variations) from signals (Zhou et al. 2018). Therefore, in order to better extract the principal components such as the trend, Chandler wobbles and annual term in the X-and Y-directions of PM, MSSA is used in this study to analyze PM variations.
With the development of observation technology, modern space geodetic technology has become the main method for solving precise PM parameters. However, owing to the complexity of the data collection and processing process, these parameters are usually delayed by hours to days. To satisfy the real-time requirements of spacecraft tracking, deep space detection and satellite positioning (Schuh and Behrend 2012;Dill et al. 2019), it is significant to research PM prediction.
Moreover, precise transformations between the international celestial and terrestrial reference frames, as well as the autonomous navigation and positioning of artificial Earth satellites have an increasing demand for precise PM prediction parameters with long term (Kalarus et al. 2010). To improve prediction of PM, many methods have been proposed. For instance, Kosek et al. (2007) used the combined least-square (LS) and autoregressive (AR) algorithm (as the official PM prediction algorithm of International Earth Rotation and Reference System Service (IERS)) to predict PM for one year into the future from 2003 to 2012. The mean absolute errors (MAEs) of 10-day-lead predictions with this method in the X-and Y-directions are 3.1 and 1.9 mas, respectively, and the MAEs of 180-day-lead predictions are 23.5 and 26.9 mas, respectively. Fuzzy-wavelet was used to predict PM during the period from 2000 to 2005. The root mean squares of 10 days and 180 days ahead prediction in the X-and Y-directions are 4.86 and 3.28 mas, 33.02 and 34.29 mas, respectively (Akyilmaz et al. 2011 applied artificial neural network (ANN) to predict PM. In the ultra-short-term forecast, its precision in the X-and Y-directions are only 21.84 and 16.87 mas, respectively, and the MAEs of 180-day-lead predictions are 22.74 and 27.67 mas, respectively. The Earth Orientation Parameters Prediction Comparison Campaign (EOP PCC) has compared these methods and shows that the LS + AR model is one of the most efficacious methods for predicting PM (Kalarus et al. 2010). Many scholars have improved the LS + AR model (Xu et al. 2012;Xu and Zhou 2015;Wu et al. 2018). For instance, Xu et al. (2012) added the Kalman filter method to the LS + AR model and achieved good results in short-term prediction. The MAEs of 10-day-lead predictions in the X-and Y-directions are 2.88 and 2.14 mas, respectively, and the precision of 50-day-lead predictions are 9.27 and 8.93 mas, respectively. And no long-term prediction was made.
However, these methods do well for short-term prediction and do not perform satisfactorily in mid-term and long-term predictions, because they cannot extract the variable amplitude of Chandler wobbles and annual terms (which have amplitudes up to several mas, or even dozens of mas) (Guo and Han, 2009;Su et al. 2014;Modiri et al. 2018). The aforementioned MSSA can effectively extract the temporally varying Chandler wobbles and annual term. Therefore, MSSA will be applied to predict PM for the first time in this paper, and a new method combining the linear model (LM) and autoregressive moving average model (ARMA) is proposed.
The structure of this paper is as follows: Sect. "Analysis of PM time series" introduces the MSSA method and analyzes the principal components of PM, such as the Chandler wobbles and annual term extracted by this method. Section "Prediction of PM parameters" explains how to use MSSA method to predict PM parameters for 1-year into the future and compares these predictions with the predictions of IERS Bulletin A. Section "Conclusion" contains the conclusion.

Data sources
The EOP 14 C04 PM time series provided by IERS used in this paper were downloaded from https:// www. iers. org/ IERS/ EN/ DataP roduc ts/ Earth Orien tatio nData/ eop. html. The data covers January 1, 1962, to May 12, 2020, with an interval of 1 day, which is calculated by solutions of various spatial measurement technologies (Very Long Baseline Interferometry (VLBI), Global Navigation Satellite Systems (GNSS), Satellite Laser Ranging (SLR) and Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS)), which is consistent with the conventional reference frame ITRF 2014, with a precision of 0.02 mas (Bizouard et al. 2019).

Multi-channel singular spectrum analysis
MSSA is a non-parametric analysis technique of nonlinear time series that can extract valid information from time series with unknown prior information. It has been widely used in the fields of climatology (Luzano 2020) and seismology (Oropeza and Sacchi 2011). As an extension of SSA, MSSA can better consider the correlations among various variables (Zhou et al. 2018). Therefore, MSSA was selected herein to analyze PM variations.
PM is divided into two components in the X-and Y-directions, so two channels are selected for analysis. First, the linear trends in the X-and Y-directions of PM are removed. It is assumed that the de-trended PM time series in the X-and Y-directions are x 1,t and x 2,t (t = 1, 2, · · · , N ) , respectively, where t is time and N is the length of the series. L represents the window size, which usually meets 1 < L < N /2 . To better extract the periodic terms, the window size is generally a multiple of the period (Golyandina and Zhigljavsky 2013). The Chandler wobbles and annual term are a reasonable choice as the main PM period (Zotov 2010;Shen et al. 2018), so 2190 is chosen herein as the window size. The total trajectory matrix H constructed by PM in the X-and Y-directions is expressed as follows: The trajectory matrix H is an L × 2 K matrix and is transformed into the sum of biorthogonal matrices of (1) rank 1, S = HH T , that is decomposed by singular value decomposition (Golyandina and Stepanov 2005). The singular value decomposition of trajectory matrix H is expressed as follows: where are the eigenvalues of the matrix S , d represents the rank of H , U 1 , . . . , U d are the corresponding eigenvectors, and Considering that H i is associated with the i-th largest singular value of H , the proper eigenvalues and corresponding eigenvectors are selected to reconstruct H: where r is the number of reconstruction, r < d.
Using anti-diagonal averaging on each block of H i , the reconstruction component (RC), a new N-length time series, is obtained. Therefore, the mth reconstruction component h m is expressed as: Through the above steps, 2190 RCs in the X-and Y-directions can be obtained, respectively. The variance contribution rates of the first 12 RCs were statistically analyzed, as shown in Fig. 1.
As shown in Fig. 1, RC1 and RC2 have the highest contribution rates, at 35.08% and 34.37%, respectively. RC3 and RC4 have the second-highest contribution rates, at 13.68% and 13.60%, respectively. The contribution rates of RC5, RC6 and RC7 are 1.06%, 0.53% and 0.48%, respectively. The contribution rates of the other RCs are relatively low and gradually decrease with the increasing order. The total variance contribution rate of the first 7 RCs is 98.81%. Under the statistical null hypothesis, F statistics can easily test the significance of any mode (Smith et al. 2007;Karegar et al. 2015), and the results show that the first 7 modes are dominant. Therefore, the first 7 RCs in the X-and Y-directions of PM were selected for analysis.

Periodic analysis of PM
Weighted correlation (w-correlation) analysis (Hassani 2007) can effectively reflect the correlations among various variables. To better determine the relationships among the RCs of PM, first, the reconstructed components in the X-and Y-directions were analyzed using the w-correlation analysis method. The w-correlation between any two reconstructed RCs of the de-trended time series of PM ( X i ) can be expressed as: w k x i k x j k , w k is the weight coefficient, and its value is min (k, M, N − k) . The closer the absolute value of ρ w i,j is to 1, the greater the correlation between the two signals is, and the higher the degree of similarity between the two signals is. It is necessary to group them together when extracting the periodic term.
The w-correlation analysis of the first 12 RCs is shown in Fig. 2. RC1 and RC2 have high correlation, and RC1 + RC2 has good independence relative to other components, so they can be divided into a group of signals. Similarly, RC3 and RC4 also have high correlation and independence, combing into a group of signals. There is a strong correlation between RC6 and RC7, but the two RCs are still weakly correlated with RC9 and RC10. Then, a fast Fourier transform (FFT) was used to perform a power spectrum analysis on the combined signals.
The de-trended PM original time series, merged results of the same periodical RC terms in the X-and Y-directions of PM and the corresponding power spectrum analysis are shown in Figs. 3 and 4, respectively. The periodic terms in the X-and Y-directions are relatively consistent, and the first principal component (RC1 + RC2) period is all 1.18 years, which represents the Chandler wobbles of the PM. The second principal component (RC3 + RC4) period is 1 year, which represents the annual term. The third principal components (RC5) are all nonlinear trends that represent the long-term nonlinear trend changes of As shown in Figs. 3 and 4, the periodic terms of PM reveal nonlinear variations by using MSSA. The amplitude of the periodic term changes with time, which is more in line with the actual variation in PM. Traditional least-square spectrum analysis can only provide a constant amplitude of the period. The amplitude of RC1 + RC2 (Chandler wobbles) varies from 20 to 200 mas, and has decreased significantly over the last 20 years. Since the length of PM series is only a hundred years, the last time a similar Chandler amplitude minima occurred was in the 1940s (Malkin and Miller 2010). The physical mechanism on this phenomenon is currently unclear, and there is still no relevant excitation mechanism that can explain it. RC3 + RC4 (the annual variation in PM) is relatively stable, and its amplitude varies between 60 and 120 mas. This change also reflects that the seasonal redistribution of the mass of Earth has good stability. For the periodic term of approximately 1.34 years, no corresponding excitation mechanism is currently known. This term may be a combination of multiple other periodic terms and thus only represent a mathematical expression of MSSA, and further research and analysis are needed.

Long-term trend analysis of PM
The trajectory and long-term trend of PM from January 1, 1962, to May 12, 2020, are shown in Fig. 5. Assuming that the PM series are X t (t = 1, 2 . . . N ) , the linear trend is extracted by using the least square (LS) method, expressed as follows: where ω 0 is a constant and ω 1 is the change rate of PM. The red dashed line represents the nonlinear trend of PM, ω 0 + ω 1 t + RC5 t (t = 1, 2 . . . N ) , which is formed by the superposition of the PM linear trend and RC5 term obtained by MSSA. The change rates of nonlinear trend of PM are also directly calculated by LS method. Table 1 (6) X t = ω 0 + ω 1 t, Fig. 3 RCs in the X-direction of PM and its perform power spectrum analysis shows the linear and nonlinear long-term trend directions and the movement rate of PM, and the results are relatively consistent. The change rates of PM in the Xand Y-directions are approximately 2.04 ± 0.02 mas/year and 3.11 ± 0.02 mas/year, respectively, and the long-term trend is 3.72 mas/year. As shown in Table 1, the long-term trends of PM obtained from different methods and different data periods have relatively large differences. During the period used in this paper, the North Pole moves to 56.79 °W in the longitudinal direction with respect to the crust, pointing to Greenland. As shown in Fig. 5, after July 1990, Earth's North Pole began to move to the west, which is also consistent with the results of Schuh et al. (2001). At the beginning of the twenty-first century, the North Pole of the earth began to move slowly eastward. In recent decades, with the increase in global average temperatures, the glaciers in Greenland have been experiencing accelerated melting (Chen et al. 2013;Adhikari and Ivins, 2016), which may be the main cause for this trend.

Principle of prediction
To meet the real-time requirements of satellite navigation for PM parameters (Jayles et al. 2016), the MSSA method combining LM and ARMA model was used herein to predict PM parameters. The prediction flowchart is shown in Fig. 6. First, the LS is used to extract the linear trend of PM series, and a linear model is established to obtain linear trend predictions. Then, MSSA is used to decompose and analyze the de-trended PM series, selecting appropriate principal components of PM for predictions. The ARMA model is used to predict the remaining components. Finally, the resulting PM predictions can be obtained by adding predictions of the three components.

Method of prediction
The prediction process includes three steps (Fig. 6): a) A linear model is used for prediction. Equation (6) shows that the linear trend of PM is X t = ω 0 + ω 1 t . The LM model is built for the linear trend predictions. b) The MSSA model is used for prediction (Golyandina and Stepanov 2005). According to step (a), the result-ing PM time series after removing the linear trends in the X-and Y-directions are Y (2) t to obtain the last L-1 values of the reconstructed signal, denoted by F as follows: where the vector R L = (a 1 , a 2 , . . . a L−1 ) can be expressed as: where U ∇ i is the vector containing the first L-1 values of the eigenvectors U i , π i is the last value of the eigenvectors U i , r is the number of reconstruction, and ν = r i=1 π 2 i .

The recurrent prediction formula for M predictions is
where M is the number of predictions for future, the coefficients a j , j = 1, 2, . . . , L − 1 derived in formula (9). The series f predictions in the X-direction of PM. Similarly, the Y-direction predictions of PM are obtained in the same way. c) The ARMA model is used for prediction. According to step (b), the remaining components series are ; this components are modeled by ARMA, which is expressed as follows: where p is the AR model's order, q is the MA model's order, c is the constant term, φ i is the autoregressive coefficient at time t − i , ε t is the error term at time t , and θ j is the moving average coefficient at time t − j . The extended autocorrelation function (EACF) is selected to determine the orders p and q of the AR model and the MA model, respectively, which can be used to build the ARMA model for the prediction of the remaining components.

Results and analysis
The PM series from January 2000 to December 2019 was selected for the prediction experiment, and the predictions were compared with the Bulletin A released by IERS to evaluate the accuracy of the LM + MSSA + ARMA model. The experiments conducted in this paper are all based on historical data, and the training period data were used for modeling to predict PM for 1-year into the future. The PM from January 1, 2000, to the predicted starting time was used as the training data. The predicted starting times are January 7, , February 4, 2016, March 3, 2016, …, October 11, 2018, November 8, 2018, and December 6, 2018. Each predicted starting time is moved forward by 4 weeks (28 days) for the next prediction. From 2016 to 2018, 13 predictions were made per year, for a total of 39 predictions. First, the PM time series from January 1, 2000 to January 7, 2016 were analyzed using MSSA to determine the principal components of PM. The reconstructed series with the first 6 and 7 RCs in the X-and Y-directions of PM and their corresponding de-trended series are shown in Fig. 7, and the correlation coefficients are shown in Table 2, showing that the reconstructed series with the first 6 and 7 RCs can represent the original series well. However, it can be seen from the w-correlation analysis results described in Sect. 2 that the combined components of RC6 and RC7 contain some other weak periodic signals. To better predict PM, this paper uses the first 6 and 7 RCs as the principal components for prediction, respectively. As shown in Fig. 7, the ARMA model is used to predict the differences (remaining components) between the de-trended PM series and reconstructed principal components by using MSSA. Through EACF analysis, it was found that for the first 6 and 7 RCs reconstructions, the ARMA (2, 9) has good applicability for modeling the remaining components in the X-and Y directions of PM, respectively.
Second, the mean absolute error (MAE) was used to evaluate the prediction accuracy (Kalarus et al. 2010), that is: where P i is the predicted value of the i-th period, O i is the observed value in the corresponding period, and n is the total number of predictions. The MAE comparison between the LM + MSSA + ARMA model and IERS Bulletin A predictions are shown in Fig. 8, and the statistics are shown in Table 3.
As shown in Fig. 8, compared to the predictions of Bulletin A, the MAEs of the predictions produced by the proposed method were smaller in most periods. Our prediction accuracies were excellent for predictive period (0-180 days and more than 250 day) for PMX and predictive period (110-320 days) for PMY. This may be because MSSA can effectively model and predict the variable amplitude of Chandler wobbles and annual terms. In addition, for the LM + MSSA + ARMA model, the prediction precision of PMX was higher than that of PMY.
As shown in Table 3, for the ultra-short-term and short-term prediction (Kalarus et al. 2010), the MAEs of 10-day-lead and 30-day-lead predictions with MSSA(6) model for PMX were 2.43 and 5.73 mas, respectively,  , the prediction precision of the proposed method was equivalent to that of Bulletin A. In summary, the proposed method can effectively and accurately predict PM parameters. Figure 8 can only reflect the overall prediction accuracy of the model. To further discuss the actual correction effect of the model, Bulletin A was used as the benchmark, and the improvement rate was used to evaluate the accuracy of the two models. When the absolute values of the measured value minus the predictions of Bulletin A is greater than that of the measured value minus the predictions of the proposed method, the predictions of Bulletin A were considered to be improved, otherwise, the improvement failed. The improvement rate formula is: The improvement analysis is shown in Fig. 9 and Table 4, yellow indicates improvement, and green indicates failure. The improvement in the prediction by the proposed method was approximately 50% on average. Overall, the proposed model with 6 RCs was slightly better than the model with 7 RCs. The highest (13) Improvement rate = Number of improvements in PM prediction Total predicted number of PM × 100%.

Conclusions
Variations in polar motion (PM) are accompanied by the physical processes of the Earth. To study the temporally varying characteristics of PM, multi-channel singular spectrum analysis (MSSA) was used to analyze approximately 58-year PM series from 1962 to 2020. Principal components such as the trend, annual term and Chandler wobbles were extracted effectively. The analysis shows that the amplitude of Chandler wobbles fluctuates between 20 and 200 mas, and its amplitude has decreased significantly over the last 20 years. The physical mechanism behind this phenomenon is currently unclear, and further research is needed. The annual variation in PM is relatively stable, and its amplitude varies between 60 and 120 mas. This change also reflects that the seasonal redistribution of Earth materials has good stability. The trends of PM in the X-and Y-directions are about 2.04 ± 0.02 mas/year and 3.11 ± 0.02 mas/year, respectively, and the long-term trend is 3.72 mas/year, moving towards N56.79 °W.
To improve prediction of PM, the MSSA method combining linear model (LM) and autoregressive moving average model (ARMA) was used to predict PM parameters for 1-year into the future, repeatedly. Comparing to predictions of IERS Bulletin A, the results show that the proposed method can effectively predict PM parameters, and the improvement rates of PM prediction for 365 days into the future are approximately 50% on average.