Improved Prediction of Polar Motions by Piecewise Parameterization

The linearized Liouville equations that describe the relations between the EOP and GAM are presented in Equations 1 and 2 1 ΔLODr∧0=χz(t) where ΔLOD_r is the variation of the length of day (ΔLOD) with zonal tidal components removed and Λ₀ is the nominal length of day of 86400 seconds. 2 p→(t)+iσ0dp→(t)dt=χ→(t) $$\vec p(t) + {i \over {{\sigma _0}}}{{d\vec p\left( t \right)} \over {dt}} = \vec \chi (t)$$

Given the coupling between x and y components of polar motion, the equation is written in complex notation, with χ→(t)=χx(t)+iχy(t),p→(t)=px(t)−ipy(t). Here, the negative sign in p→(t) accounts for p_y(t) being positive toward 90°W longitude, as the coordinates of polar motions are defined as a left-hand Cartesian coordinate system, whereas the angular momentum function components are defined in the TRF frame which is a right-hand Cartesian coordinate system.

The transformation of GAM and EOP in the z axis is straightforward; here we only give the discrete transformation formula for the x and y components for numerical use (the derivation of these two formula can be found in Barnes et al. (1983)) 3 χ¯(t)=i2σcwδtexp(−iπδtTcw)⋅{ p→(t+δt)+[ 1−exp(iσcwδt) ]p→(t)−exp(iσcwδt)p→(t−δt) } \[\bar{\chi }(t)=\frac{i}{2{{\sigma }_{cw}}\delta t}\exp (\frac{-i\pi \delta t}{{{T}_{cw}}})\cdot \{\vec{p}(t+\delta t)+[1-\exp (i{{\sigma }_{cw}}\delta t)]\vec{p}(t)-\exp (i{{\sigma }_{cw}}\delta t)\vec{p}(t-\delta t)\}\] 4 p→(t)=−iσcwδt2exp(iπδtTcw)⋅[ χ→(t)+χ→(t−δt) ]+exp(iσcwδt)p→(t−δt) \[\vec{p}(t)=\frac{-i{{\sigma }_{cw}}\delta t}{2}\exp (\frac{-i\pi \delta t}{{{T}_{cw}}})\cdot [\vec{\chi }(t)+\vec{\chi }(t-\delta t)]+\exp (i{{\sigma }_{cw}}\delta t)\vec{p}(t-\delta t)\] where δt =1 day is the time interval for the GAM and polar motion series, the complex Chandler frequency σ_cw = 2π(1 + i/2Q)/T_cw, with period T_CW = 433 days and dampling of Q = 179 (Gross et al. 2003).

Figure 2 shows the transformed x and y components of GAM, ESMGFZ EAM, and differences between GAM and EAM. It can be seen a remarkable correlation between GAM and EAM and seasonal oscillations in the Y component, whereas relatively lower correlation and seasonal oscillations are seen in the X component.

Figure 2.

In Figure 3, we present plots of GAM versus EAM to show such phenomenon. The magnitude of GAM and the correlation between GAM and EAM are higher for the Y component than that for the X component. Since the major continents of the Earth are more aligned in the Y axis that produce relatively larger excitation of PMY from atmospheric surface pressure loading, whereas atmospheric loading effect over the oceans is largely compensated by the inverted barometer response of the ocean surface (Boy et al. 2009).

Figure 3.

As shown in Figure 1, there are two steps of LS+AR. The first LS+AR is used to predict the 6 days GAM-EAM series and the second LS+AR is used to predict the 90 days GAM series. The strategy of LS is same for the two steps, we fit the x and y components together in complex form, 5 χ(t)={ a1+a2ta3+a4ta5+a6ta7+a8t+∑j=14[ bjcos(2πtTj+ϕj)+ibj*cos(2πtTj+ϕj*) ] \[\chi (t)=\{\begin{array}{*{35}{l}} {{a}_{1}}+{{a}_{2}}t \\ {{a}_{3}}+{{a}_{4}}t \\ {{a}_{5}}+{{a}_{6}}t \\ {{a}_{7}}+{{a}_{8}}t \\ \end{array}+\underset{j=1}{\overset{4}{\mathop \sum }}\,[{{b}_{j}}\cos (\frac{2\pi t}{{{T}_{j}}}+{{\phi }_{j}})+ib_{j}^{*}\cos (\frac{2\pi t}{{{T}_{j}}}+\phi _{j}^{*})]\] where a₁, a₂, …, a₈ are complex numbers with ai=(areal+iaimag*)i, b_j and bj* denotes the amplitude, ϕj* and ϕ_j denote the phase of the periodical componts, with T₁ = 1 years, T₂ = 1/2 years, T₃ = 1/3 years, and T₄ = 13.7 days. Same as Dill2019, the linear components are fitted every year, and the four periodic components are fitted with 4 years data. In total, we used 5 complexes, a₁, a₂, a₄, a₆, a₈, and 16 real numbers (b_j, bj*, ϕ_j, and ϕj*) in each step of the least squares fit. a₃, a₅ and a₇ are not free parameters, because of the piecewise continuity constraint.

Figure 4.

Up to now, the techniques we used are all identical to Dill2019. Here we present the major technical differences of our method.

For a stationary random sequence x(t)(t = 1, 2,…, N), the AR(p) model is expressed as 6 x(t)=c+∑i=1pϕix(t−i)+εt. where ϕ_i, ϕ₂, …, ϕ_p are autoregressive coefficients, that can be solved by the Yule-Walker equations by means of the Levinson-Durbin rercursion (Brockwell and Davis 1996), and p is the order of the AR model that is determined by Akaike’s Final Prediction Error (FPE) criterion (Akaike, 1971). ∈_t is the unmodeled random noise component.

Practically, p, the order of AR model, is a critical parameter which is very sensitive to the final accuracy of prediction. In Dill2019, the p is adopted as constants, with p = 20 days for 1-6 days prediction of the equatorial components of GAM-EAM, with p = 60 days for 1-6 days prediction of the axial component of GAM-EAM, with p = 2 days for 7-90 days prediction of the equatorial component of GAM and with p = 25 days for 7-90 days prediction of the axial components of GAM.

In our prediction, we aussmed there exists a best p_n for the nth-day prediction, which is not known at the beginning but can be searched by further evaluations. Frankly speaking, such trials are still on the testing stage, and we still have no physical explanation of why some p works well than others. Therefore, feasibility of such method still needs careful evaluation. Meanwhile, in the steps of AR, we used a modified AR formula that includes one more sampling parameter, lag, 7 χ(t)=c+∑i=1pϕiχ(t−i)sign(i)+εt. \[\chi (t)=c+\underset{i=1}{\overset{p}{\mathop \sum }}\,{{\phi }_{i}}\chi (t-i)sign(i)+{{\varepsilon }_{t}}.\] with i < p, and sign(i) is a signum function, 8 sign(i)={ 1i%lag=00i%lag≠0 $$[sign(i) = \{ \matrix{ 1 & {i\% \,lag\, = 0} \cr 0 & {i\% \,lag\, \ne 0} \cr } $$]

The effect of the new lag parameter is a kind of desampling smooth filter. When lag equals 1, the Eq. 7 will be same as Eq. 6, whereas when lag equals 2, which means setting the odd autoregressive coefficients, ϕ₁, ϕ₃, ϕ₅, …, in Eq. 7 as 0. Similar to the order p, the best lag is also determined by evaluation. In summary, in the process of AR, we have two parameter (p, lag), values of which are determined by evaluation.

Figure 5 shows exmaples of prediction for polar motions x and y, by adopting different (p, lag). It is noted that the AR is conducted for x and y components separately, whereas the LS is conducted in the complex domain for both x and y components simultaneously. It can be seen that the prediction of PMY is around 2 times better than the prediction of PMX. Given that the AR is prediction of noisy signals, and χ_x are much more noisy than χ_y, a better prediciton of PMY than PMX is explicable. Especially, within 6 days, the selection of (p, lag) is very sensitive for prediction of PMX.

Figure 5.

Asumming there exists best (p, lag) of AR, we make a large amount of trials and prediction of EOPs in 441 days from Sep. 2019 to Feb. 2021 to find the best AR parameter (p, lag). For the prediction of GAM-EAM within 6 days, we searched p within the range of [10,90] and searched the lag within the range of [1, 20]. For AR of GAM within 7 to 90 days, we searched p within the range of [2, 20], and lag within the range of [1,20]. In Tables 1 and 2, given are the best AR parameters (p, lag) we adopted for the polar motion x and y, respectively.

In Figure 6, shown are the prediction errors for all the 441 predictions from Sep. 2019 to Feb. 2021. Within 90 days, the maximum prediction error of PMX is 36 mas, the maximum prediction error of PMX is 16 mas; within 10 days, the maximum prediction errors of PMX and PMY are 4.5 and 2 mas. On an average, the prediction of PMY is around 2 times better than that of PMX. In addition, we found that the prediciton errors of PMX are much higher in 2020 than 2019. In 2019 September, the prediction errors of both PMX and PMY are lower than 10 mas, whereas in 2020 August and September, the 30-90 days prediction errors are around 2 or even 3 times higher than that of 2019. We suspect that such phenomenon might be due to the pandemic of coronavirus COVID-19, which result in a decrease in global airline by more than 60%. As the air-borne meteorological radar data are very important for numerical weather model, a decrease of flight number finally results in worse weather and AAM forecast, and then worse EOP prediction. A future study on evaluation of whether and/or how much the AAM data might be influenced by the lack of air-borne meteorological data can be helpful to verify or deny such hypothesis.

Figure 6.

To compare with IERS bulletin A prediction, in Figure 7, we show the MAE of both our prediction and IERS bulletin A prediction. Our prediction of PMY is better (~20%) than bulletin A prediction in all timescale. Especially, on the 5th day and 90th day, the MAE is reduced by 49.0% and 28.9%, respectively. For the prediciton of PMX, within 30 days, our predictions are slightly better (2% - 8%) than bulletin A, but become worse (−7%~ −19%) than bulletin A within 30-90 days. In Table 3, we present the prediction errors (MAE) at different future days.

Figure 7.