Search for the 531-day-period wobble signal in the polar motion based on EEMD

In this study, we use a nonlinear and nonstationary time series analysis method, the ensemble empirical mode decomposition method (EEMD), to analyze the polar motion (PM) time series (EOP C04 series from 1962 to 2013) to find a 531-day-period wobble (531 dW) signal. The 531 dW signal has been found in the early PM series (1962–1977), but cannot be found in the recent PM series (1978–2013) using conventional analysis approaches. By virtue of the demodulation feature of EEMD, the 531 dW can be confirmed to be present in PM based on the differences of the amplitudes and phases between different intrinsic mode functions. Results from three sub-series divided from the EOP C04 series show that the period of the 531 dW is subject to variations, in the range of 530.9–524 days, and its amplitude is also time-dependent (about 2–11 mas). Synthetic tests are carried out to explain why the 531 dW can only be observed in recent 30-year PM time series after using EEMD. The 531 dW is also detected in the two longest available superconducting gravimeter (SG) records, which further confirms the presence of the 531 dW. The confirmation of the 531 dW existence could be significant in establishing a more reasonable Earth rotation model and may effectively contribute to the prediction of the PM and its mechanism interpretation.


Introduction
It is recognized that the polar motion (PM) contains two dominant components: the annual wobble (AW) with a 12 month period and the Chandler wobble (CW) with a 14 month period.Some researchers suggested that the CW is highly variable with respect to its amplitude (e.g., Carter, 1981Carter, , 1982;;Höpfner, 2003;Chen et al., 2009), some considered it having double or multiple frequencies (e.g., Chao, 1983;Pan, 2012), and some considered its frequency being invariant (e.g., Okubo, 1982;Vicente and Wilson, 1997;Gross et al., 2003;Seitz and Schmidt, 2005).If the CW is frequency modulated as Carter (1981 and1982)  by the magnitude, it will creates an infinite number of sidebands, arranged symmetrically about the carrier and spaced at integer multiples of the modulating frequency (Carter, 1981).The first upper and lower sidebands could be at 1 cpy (cycles per year) and 0.69 cpy respectively when the beat frequency is 0.157 cpy, where the beat frequency or beat period is the time period required for the Earth's pole to complete a cycle of the combined AW (12 months) and CW (14 months).Because the first upper sideband was contaminated in the spectrum of the AW, it was concluded that only the first lower sideband, which is located at 0.686 cpy, could most likely be detected.Based on a 16-year time series of International Polar Motion Data (spanning from 1962 to 1977), a 0.686 cpy component with its amplitude being around 10 to 17 mas (milliarcsecond) was very weakly detected in Catrer (1982), but its signal-to-noise ratio (SNR) is very low.Based on the PM series (spanning from 1974 to 1981) obtained by the Doppler satellite tracking and lunar laser ranging, Morgan et al. (1982) identified two spectral peaks at 532±10.8 days and 537±15.2days, with their amplitudes around 8.6±2.0 mas and 7.4 ± 2.2 mas, respectively.However, after then (1980s), this signal was seldom announced in the PM time series with higher SNR (especially for the records after the mid-1990s).Recently, Na et al. (2011) found a 500-day period component in the PM data with an average amplitude of 20 mas, and argued that if assuming the existence of this component, the RMS error of the prediction of the PM could be reduced by 50 %; and suggested that this phenomenon should be caused by resonance of unidentified oscillating mode of the Earth (possibly Earth's inner core wobble).In addition, this wobble (or referred to as an 18 month wobble) was found in the analysis of the atmospheric angular momentum data by Wahr (1983) and Chen et al. (2010).Furthermore, a signal with a period about 1.5 year could be also found from the variation of length of day (LOD) by using the wavelet analysis as suggested by Chao et al. (2014).prediction of the PM, the main purpose of this study is to detect the ∼ 530-day-period wobble (simply 531 dW hereafter for convenience, the reason of which could be found in Sect.3) component and provide a possible explanation why in general cases it could not be detected from observation series (e.g.EOP C04 series) based on the Fourier anaylsis.
Generally, the traditional Fourier analysis method cannot observe this 531 dW signal, hence, in this study a nonlinear and non-stationary time series analysis method, the ensemble empirical mode decomposition (EEMD) (Huang and Wu, 2008;Wu and Huang, 2009) is applied to a PM time series as a dyadic filter bank to detect the 531 dW signal.
Here, the PM time series, the EOP C04 series spanning from 1962 to 2013 with one-day sampling interval from the International Earth Rotation and Reference System Service (IERS) (http://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/), is used for observing the signal of interest.

Method
EEMD was proposed to overcome the disadvantages existing in the empirical mode decomposition (EMD) (Huang et al., 1998;Huang and Wu, 2008), such as the modemixing problem and the end effect (see as Fig. 2 in Shen and Ding, 2014), although EMD has demonstrated its applicability in a wide range of geoscience studies over the last 15-year (e.g., Pee and McMahon, 2006;Thomas et al., 2009;Franzke, 2009;Jackson and Mound, 2010;Jeng and Chen, 2011;Lee et al., 2011;Shen and Ding, 2014;Chambers, 2015).The details of EMD and EEMD can be found in many relevant literatures (e.g., Huang et al., 1998;Huang and Wu, 2008;Wu and Huang, 2009;Shen and Ding, 2014).
A given time series x(t) can be decomposed into a number of IMFs using the following steps (Huang and Wu, 2008;Shen and Ding, 2014): Introduction

Conclusions References
Tables Figures

Back Close
Full 1. Identify all local maxima and local minima of x(t), where its upper (lower) envelope can be formed using a cubic spline line to connect all the local maxima (minima).
2. Use h 1 (t) = x(t) − m 1 (t) to get a new series, where m 1 (t) is the mean of the upper and lower envelopes of x(t).Steps 1 and 2 are called a sifting procedure.
3. Generally, h 1 (t) is not an IMF, hence, treat h 1 (t) as the data given just as x(t), and repeat the sifting procedure (step 1 and 2) k times until h k (t) is an IMF, namely, , which is designated as the first IMF.
5. Repeat the above steps, and then the j -th target IMF, c j (t), can be obtained.
The above process is the EMD.In these steps, two different iterative loops exist, and the criteria to stop them can be found in Huang et al. (1998).
Given that EMD has the mode-mixing problem and end effects, EEMD is developed with the following aspects (Huang and Wu, 2008;Shen and Ding, 2014): 1. Add a white noise series to the targeted time series x(t).
2. Decompose the series with added white noise into IMFs.
3. Repeat step 1 and step 2 iteratively, but with different white noise series each time.Introduction

Conclusions References
Tables Figures

Back Close
Full By using EEMD, the chosen time series can be decomposed into a finite number of simple IMFs, and different-scales signals in the series will be re-combined by proper IMFs based on the fact that different IMFs have different frequency bands.Hence, for the chosen EOP C04 series, f (t), and its IMFs satisfy the following linear equation where n is the number of the IMFs and r(t) is a residue term.In addition, EEMD can be used to demodulate a frequency-modulated time series (Huang et al., 1998;Huang and Wu, 2008).That means EEMD can be used not only as a dyadic filter bank but also as a demodulator simultaneously.Since EEMD has these two advantages, we use it to detect the 531 dW in the PM.
A frequency modulation signal can be expressed as follows (Cater, 1981 and1982), where e t (x, y) is the expected value of the x or y component, f m the frequency of the modulating signal, here we set it to be equal to 0.157 cpy (the beat frequency between the CW and AW); f c and C c are the frequency and amplitude of the CW, and M is the modulation index, defined as M = ∆f /f m , where ∆f is the maximum variation of f m ; φ 0 is the initial phase, and it is simply set as zero (same as the AW and the 531 dW).

Results from the PM series with/without using EEMD
The waveforms of the two components, the x and y components of the chosen PM series are shown in Fig. 1a, and their corresponding Fourier spectra are shown in Fig. 1b and c (see the black bold curves, we call them the original spectra).The vertical dashed lines (Fig. 1b and c) locate at the target freuquency, 0.6875 cpy, and there is Introduction

Conclusions References
Tables Figures

Back Close
Full no peak that can be identified for the 531 dW based on original PM series.Meanwhile, if the EEMD is used to these two components, 11 intrinsic mode functions (IMFs) can be obtained, respectively.In the Fourier spectra of IMF5 and IMF6, the 531 dW can be clearly detetcted (see the blue and red curevs in Fig. 1b and c).The spectra for other IMFs (sum) of x and y components are denoted by the green curves (we call them the residual spectra).No peak can be identified for the 531 dW in the residual spectra.The amplitudes of the 531 dW signals in IMF5 and IMF6 are denoted by the green and red arrows respectively.Comparing with the original spectra, the amplitudes in IMF5 and IMF6 are outstanding (Fig. 1b and c), but the phases of them are almost opposite (Fig. 1d and e).According to the characteristics of EEMD, the mean amplitudes of 531 dW in x and y components achieve 7.1 and 7.2 mas, respectively; while the corresponding noise level is about 4 mas.This might be the reason why we cannot identify the 531 dW signal directly from the original PM time series spanning from 1962 to 2013.However, previous studies show that this signal has been found in the 1962-1977 PM series.One possible inference is that the amplitude of the 531 dW may be varing with time.
Based on the conventional Fourier anaylsis approach, the results as shown in Fig. 2 clearly indicate that only the target peak in the spectra of the 1962-1977 series (see Fig. 2b) is over their corresponding background noise level (see Fig. 2d; note that a 0.57 cpy peak in Fig. 2f is also over the noise level, but it is not the interesting signal of this study), which is consistent with previous studies (Carter, 1981 and1982;Morgan et al., 1982).Without using EEMD, our estimates for the target 531 dW, CW, and AW (the annual wobble) are tabulated in Table 1.For the x and y components of the 531 dW signal from the 1962-1977 series, the corresponding amplitudes are 11.3 and 14.6 mas, while the estimates of previous studies are about 8 mas (Carter, 1981 and1982;Morgan et al., 1982).However, this wobble cannot be found in the 1978-1994 and 1995-2013  After apppling EEMD to the three sub-series, only IMF5 and IMF6 contain the 531 dW signal, the corresponding spectra being shown in Fig. 3 (for x component) and Fig. 4 (for y component).As for other IMFs, the corresponding peaks at 531 dW frequency can neither exceed the background noise levels of their spectra just as in the whole 1962-2013 PM series they cannot (see the green curves in Fig. 1).Note that we do not rule out the possibility that some part (energy) of the interesting signals may also present in the adjacent IMFs, but they must be very weak (see Fig. 1) so that they can be neglected.Hence, in the present case, we will only concern the IMF5 and IMF6.
As shown in Figs. 3 and 4, the phases for CW (and AW) in IMF5 and IMF6 are same (except for Figs.5f and 6f, where there is no peak for AW), whereas the phases for the 531 dW in IMF5 and IMF6 are opposite with each other.The corresponding estimates are listed in Table 1.Taking into account the opposite phases as shown in Figs. 5 and 6, one can find that the amplitudes from the IMFs for CW and AW are consistent with the results without using EEMD (see Table 1).For example, for the x component of the 1962-1978 series, the amplitude of CW after using EEMD is about 128.7 mas (103.6 + 25.1 mas), whereas its amplitude without using EEMD is 129.2 mas.As for the amplitude of the 531 dW, the results from the x and y components of the 1962-1978 series after using EEMD are respectively (44.1 mas − 33.2 mas =) 10.9 mas and (44.9 mas − 32.7 mas =) 12.2 mas (note that they have opposite phases), whereas the corresponding results without using EEMD are respectively 11.3 and 14.6 mas; they are consistent with each other very well.Moreover, the results from the x and y components of the 1978-1994 series after using EEMD are (50.0-45.3mas =) 4.7 and (50.5-44.1 mas =) 6.4 mas, whereas the corresponding noise levels as shown in Fig. 2d are 3.63 and 3.43 mas; and for those of the 1995-2013 series after using EEMD, the results are (33.3-31.2mas =) 2.1 and (33.3-32.8mas =) 0.5 mas, whereas the corresponding noise levels as shown in Fig. 2f are 5.63 and 5.56 mas.Considering the estimation errors, the results using EEMD also indicate that the 531 dW cannot be found in the spectra by using the Fourier analysis (as shown in Fig. 2e and f).Introduction

Conclusions References
Tables Figures

Back Close
Full Wahr (1983) implied that there exists a 18 month polar wobble (which can be considered as coinciding with the 531 dW signal mentioned in the present paper), which could be caused by oceanic excitation.Chen et al. (2010) further demonstrated that the atmospheric excitation may give rise to the 531 dW signal with an amplitude 15 mas using the atmospheric angular momentum (AAM) dataset spanning from 1948 to 2006.But from the above results, those excited 531 dW signals almost disappear from the 1978-1994, and 1995-2013 series.This might be due to the fact that the 531 dW signal has varying amplitudes and phases with time.However, after using EEMD, the 531 dW signals can be detected in two IMFs with opposite phases.Hence, we may carefully suggest that this phenomenon can be explained by the demodulation feature of EEMD.And this feature of EEMD gives us a chance to confirm that the 531 dW is present in PM due to the differences of the phases between different IMFs.

Synthetic tests for the frequency modulation and demolulation
By carefully examining Table 1 we can find that the amplitudes of the 531 dW in IMF6 clearly have some proportional relation with their corresponding amplitudes of CW, whereas the amplitudes of the 531 dW and CW in IMF5 have no obvious relationship.According to previous studies, here we assume that the amplitude of the 531 dW in IMF6 is caused by the frequency modulation of CW.
The expression of the frequency modulation signal can be found in Eq. (2) in Sect. 2. From Table 1, we can get the amplitudes of the x and y components of the three sub-series (without using EEMD), hence, we use those parameters to generate three synthetic noise-free time series based on Eq. ( 2), and they have the same sampling intervals and lengths as the 1962-1977, 1978-1994, and 1995-2013 series, respectively.As previous studies (Carter, 1981(Carter, , 1982) suggested, we first test the modulation index M of CW which equals 0.23 and 0.38 respectively, and set different amplitudes for the CW and AW to ensure that the synthetic spectra coincide with the spectra of the actual observations (here we use the x components of the three sub-series as examples), and the frequencies of the CW and AW are set to 0.8437 and 1.00 cpy.The Introduction

Conclusions References
Tables Figures

Back Close
Full set amplitude parameters (input) and their corresponding synthetic output are listed in Table 2.The corresponding spectral results are shown in Fig. 5. From Table 1 and Fig. 5, one can clearly find that if M = 0.23, for the 1962-1977 series, the amplitude (x-component) of the sideband that is caused by the frequency modulation of the CW is 15.17 mas, and our estimates from the 1962-1977 series are about 11.3 and 14.6 mas for x and y components, respectively.Obviously, our synthetic result is consistent with our observed results very well.If we only consider this fact, we may conclude that the observed target wobble from the 1962-1977 series is very likely to be originated from the frequency modulation of CW.But here we do not want to infer whether the CW is frequency modulated, but obviously, the synthetic results for the 1978-1994 and 1995-2013 series based respectively on M = 0.23 and 0.38 clearly show the appearance of the target wobble, whereas there is no significant peak for the target wobble in the corresponding actually observed spectra.Namely, the modulation index M = 0.23 or 0.38 of CW cannot explain the observed results from the later two sub-series.Although we cannot ensure that M is a constant, we try to find an M that is suitable for all the three sub-series.We find that M = 0.5 is a good choice.Now we set M = 0.5, and construct three synthetic noise-free time series based on Eq. ( 2) for the 1962-1977, 1978-1994, and 1995-2013 series, and still set different amplitudes for the CW and AW to guarantee that the theoretical spectra coincide with the observational spectra (the input parameters are as same as the parameters for the three series without considering frequency modulation of CW, which are listed in Table 3; note here only CW and AW are concerned).The corresponding results are shown in Fig. 6  531 dW in IMF6 is really caused by the frequency modulation of CW with a modulation index M = 0.5, the 531 dW in IMF5 may be excited by some geophysical processes, such as the atmospheric/oceanic excitation (Wahr, 1983;Chen et al., 2010), then we can appropriately explain the observations after using EEMD and the disappearance of the 531 dW in the recent PM series.Here, we would like to further demonstrate that EEMD can directly demodulate a frequency-modulated time series.We generate six synthetic noise-free time series.The length of the synthetic series I and II is equal to that of the 1962-1977 series, the length of the synthetic series III and IV is equal to that of the 1978-1994 series, and the length of the synthetic series V and VI is equal to that of the 1995-2013 series, with one-day sampling interval.
The synthetic series I (III/V) contains three sinusoidal components without frequency modulation, and the frequencies and amplitudes of CW and AW (see Table 3) are set to make the results for CW and AW after considering the frequency modulation of CW be as same as with the observations from the x component of the 1962-1977 series (1978-1994 and 1995-2013; see Table 1 and Fig. 7b-f).The third signal is the 531 dW which is listed in Table 3.As for the synthetic series II (IV/VI), its parameters are as same as the synthetic series I (III/V), but a mechanism of frequency modulation of CW is considered with the modulation index M = 0.5.The spectra of the six synthetic series are shown in Fig. 7. Obviously, after considering the frequency modulation of CW, the amplitude of the 531 dW in the synthetic series II is reduced to only 11.05 mas, which is almost equal to the observed result obtained from the x component of the 1962-1977 series without using EEMD (11.3 mas, see Fig. 2 and Table 1).The corresponding spectra of the IMF5 and IMF6 of the synthetic series II obtained after using EEMD are shown in Fig. 7b1 and b3.Similar to the real data, the phases for CW (or AW) in IMF5 and IMF6 are same, whereas the phases for the 531 dW in IMF5 and IMF6 are opposite with each other (Fig. 7b2).Comparing those synthetic results with the results from the x component of the 1962-1977 series as listed in Table 1, the amplitudes of the 531 dW in IMF5 (44.1 mas) and IMF6 (33.2 mas) of the real data are almost equal to those of the synthetic results.The results as shown in Fig. 7c-f  nature with the results as shown in Fig. 7a and b.Namely, the input 531 dW signal is in fact demodulated into IMF5, whereas the 531 dW signal caused by the frequency modulation of CW is demodulated into IMF6.
The frequency modulation mechanism of CW is a pending question.If we accept this mechanism, the results obtained from the PM series after using EEMD can be appropriately explained.However, the excitation sources of the 531 dW in IMF5 still cannot be confirmed.We do not intend to figure out this question in this study, but we will try to find the 531 dW in the gravity records in the following section to further confirm its presence.

Results from the SG records
If the 531 dW can be excited by the atmospheric/oceanic angular momentum, or it is the normal mode of the Earth, it may be found in the gravity records.Hence, we choose the SG records from the GGP network to further confirm the existence of this signal.The GGP network has been operating about 25 yr (from 1997), but the longest available continuous records (without very large gaps) are only about 16 yr long, namely the mb (1995-2011) and mc (1998-2014) records.After removing the tidal effects (solid and ocean), the local atmospheric effect (without 531 dW signal in it), and the pole tides, the spectra of the residual gravity records are shown in Fig. 8 (see the black curves).Two clearly spectral peaks for the 531 dW can be found in Fig. 8a (for mb) and b (for mc).We have confirmed that the local atmospheric effect and the pole tides have no contributions to the 531 dW.Since the pole tides have been removed, namely the CW, AW and the possible 531 dW signals caused by the polar motion have been removed from the gravity records, the residual 531 dW must be casued by some other sources.
Here we further consider the loading (global) effect except the atmospheric loading (global) and non-tidal ocean loading effects.The two latter effects are the same possible sources which excite the 531 dW in PM series.However, Fig. 8 clearly shows those three effects cannot explain the residual 531 dW, only the residual AW can be appropriately explained.This finding seems being not consistent with previous studies, but they

Conclusions References
Tables Figures

Back Close
Full can be explained.As Wahr (1983) and Chen et al. (2010) suggested, the 531 dW can be excited by the atmospheric/oceanic effects, and the excited signal by atmospheric or oceanic effect can be found in the PM series only after convoluting with the CW term.However, the excited 531 dW signals by atmospheric and oceanic effect seem have different phases, hence the 531 dW cannot be found in the PM series (Chen et al., 2010).
If we accept this, the 531 dW signal in gravity records comes from two different sources: indirectly from polar motion which are affected by the atmospheric/oceanic effects, and directly from the atmospheric/oceanic effects.The excited 531 dW from PM has been removed by taking away the pole tide effects.The excited 531 dW directly from the atmospheric/oceanic effects has just the same nature as the excitation in the PM; if no convolution is processed with the CW, the 531 dW signal cannot be found in the PM or gravity records.Hence, our findings are actually consistent with previous studies.
No matter how, our results clearly show the demodulation feature of EEMD is helpful for detecting the 531 dW signal in the PM series, and we confirm that the amplitude and frequeny of the 531 dW are varying with time.

Discussions
After applying EEMD to the 1962-2013 PM time series (EOP C04), a 531 dW is clearly found with a mean amplitude about 7 mas (with much larger amplitudes in IMF5 and IMF6 series respectively), but in the spectra without using EEMD, this signal cannot be found.The 531 dW has been found by previous studies in the 1962-1977 PM series with a lower SNR.To confirm previous observations (Carter, 1981 and1982) and further study this signal, we divide the whole PM series into three sub-series, 1962-1977, 1978-1994 and 1995-2013  the decomposed IM5 and IMF6 series), which (taking x-components as example) are respectively about (44-33 mas =) 11 mas in 1962-1977, (50-45 mas =) 5 mas in 1978-1994 and (33-31 mas =) 2 mas in 1995-2013, while the corresponding noise levels of the three sub-series are about 2.7, 3.5 and 5.5 mas recpectively.That is reason why the 531 dW signal can only be directly found in the 1962-1977 PM series without using EEMD.
Although the frequency modulation mechanism of CW is a pending question, we find that if the modulation index M of CW equals 0.5, the results obtained from the PM series after using EEMD can be appropriately explained.Furthermore, using synthetic tests we confirmed the demodulation feature of EEMD which can help us find the 531 dW signal in PM series, but we cannot explain the possible exctiation sources of the 531 dW.Given that the 531 dW signal can be directly found in the SG records without using EEMD, it should be considered in the study of the long period effects in some relevant geophycial datasets.However, owing to the frequency and amplitude of the 531 dW signal being time-varying, it becomes quite difficult to explore its excitations.It might be caused by the core dynamics, or even by some random process (see as Chao et al., 2014), but this is only a conjecture.Undoubtedly, whether the 531 dW belongs to a class of normal modes of the Earth is also a pending question; more studies are needed to further understand this signal and its excitation mechanisms.Introduction

Conclusions References
Tables Figures

Figure 1 .Figure 2 .
Figure 1.The x (black) and y components (red) of the 1962-2013 PM series (a) and their corresponding Fourier amplitude spectra (original spectra, denoted by the black curves in b and c).(b) x component; (c) y component.The spectra for IMF5 and IMF6 after applying EEMD to x and y conponments are indicated by blue and red curevs; and their corresponding phases are shown in (d and e), respectively.The spectra (residual spectra) for other IMFs (sum) of x and y components are denoted by the green curves.