Extraction of periodic signals in Global Navigation Satellite System (GNSS) vertical coordinate time series using the adaptive ensemble empirical modal decomposition method

Empirical modal decomposition (EMD) is an efficient tool for extracting a signal from stationary or non-stationary time series and is enhanced in stability and robustness by ensemble empirical mode decomposition (EEMD). Adaptive EEMD further improves computational efficiency through adaptability in the white noise amplitude and set average number. However, its effectiveness in the periodic signal extraction in Global Navigation Satellite System (GNSS) coordinate time series regarding the inevitable missing data and offset issues has not been comprehensively validated. In order to thoroughly investigate their impacts, we simulated 5 years of daily time series data with different missing data percentages or a different number of offsets and conducted them 300 times for each simulation. The results show that high accuracy could reach the overall random missing rate below 15 % and avoid consecutive misses exceeding 30 d. Meanwhile, offsets should be corrected in advance regardless of their magnitudes. The analysis of the vertical components of 13 stations within the Australian Global Sea Level Observing System (GLOSS) monitoring network demonstrates the advantage of adaptive EEMD in revealing the time-varying characteristics of periodic signals. From the perspectives of correlation coefficients (CCs), root mean square error (RMSE), power spectral density indices (κ) and signal-to-noise ratio (SNR), the means for adaptive EEMD are 0.36, 0.81, −0.18 and 0.48, respectively, while for least squares (LS), they are 0.27, 0.86, −0.50 and 0.23. Meanwhile, a significance test of the residuals further substantiates the effectiveness in periodic signal extraction, which shows that there is no annual signal remaining. Also, the longer the series, the higher the accuracy of the reasonable extracted periodic signal concluded via the significance test. Moreover, driving factors are more effectively facilitated by the time-varying periodic characteristics compared with the constant periodic signal derived by LS. Overall, the application of adaptive EEMD could achieve high accuracy in analyzing GNSS time series, but it should be based on properly dealing with missing data and offsets.

1 Introduction

With the development of the GNSS (Global Navigation Satellite System) and the International GNSS Service (IGS), many countries have established GNSS networks, such as the worldwide IGS, SAPOS (the satellite positioning service of the German National Survey) in Germany, GEONET (the GNSS Earth Observation Network system) in Japan and CMONOC (the Crustal Movement Observation Network of China) in China. Also, there is a GNSS data center named SONEL (Système d'Observation du Niveau des Eaux Littorales), which serves the Global Sea Level Observing System (GLOSS). These GNSS networks are widely applied in geodesic and geodynamic studies (Gülal et al., 2013; Scaramuzza et al., 2017). GNSS coordinate time series are a crucial type of geophysical data. They can be used not only to obtain accurate positions and velocities of stations (Johnson et al., 2021; Zhou et al., 2022), but also to better study Earth's dynamic processes such as crustal movement (Calais et al., 2023), post-ice rebound and sea level change (Din et al., 2019), inversion mass change (Willen et al., 2022) and other applications (Dong et al., 2002; Munekane, 2021; Costantino et al., 2023). These coordinate time series are generally comprised of a periodic signal, trends and noise. The periodic signals are caused by a complex interplay of the Earth's own physical properties and environmental factors, which is particularly obvious in the elevation direction (Dong et al., 2002; Ray et al., 2008; Wang et al., 2018). The trend terms are composed of long-term trends and transient trends. For the transient trend, this mainly includes the effects of earthquakes and tides. The 26 December 2004 Sumatra earthquake (M_w 9.2–9.3) caused an after-slip of larger than 10 mm within minutes as far away as India (Blewitt et al., 2006). It is usually modeled with logarithmic or exponential functions (Hetland and Hager, 2006; Nishimura, 2014; Tobita, 2016). Moreover, modulations in plate motion are compatible with the solar year, the period of the lunar perigee and the lunar nodes, which clearly the influence of lunar and solar tidal forces on the plate motion (Zaccagnino et al., 2020). Also, this motion by oscillating tides can be shown by the deep tremors (Ide and Tanaka, 2014). Since our focus is the long-term trend, we mainly investigate the impact of periodic signals on it. Therefore, accurately extracting periodic signals from GNSS coordinate time series is crucial (Van Dam et al., 2001; Sorin et al., 2021). However, the periodic signals estimated by the traditional least squares (LS) method did not agree well with the actual periodic variations (Bennett, 2008; Li, 2020). Therefore, some scholars aimed at more suitable periodic signal extraction methods. The Kalman filtering, assuming signals randomly, is contrary to the real GNSS time series containing power-law noise (Davis et al., 2012). Also, the high sensitivity of F-test detection to noise may lead to the incorrect determination of the presence of periodic signals (Li and Shen, 2014). The principal component analysis (PCA) method found that most GNSS stations around the world have a nonlinear periodic pattern, with annual and semi-annual terms dominating (Shen et al., 2014). These studies verified the existence of time-varying periodic signals in GNSS time series and separated the signals according to different methods.

Moreover, empirical mode decomposition (EMD) introduced by Huang et al. (1998) initially focused on the analysis of signals that are a mixture of harmonic oscillations with different periods. EMD requires the construction of the upper and lower envelopes based on the local extrema, by which approximating curves are drawn using cubic splines (Kopsinis and McLaughlin, 2008). Because of this, it works unstably once local extrema are poorly defined, e.g., when the data contain relatively long plateaus of constant values or short-lived waves of small amplitudes and short periods (Huang et al., 2009). Adding white noise to the original data could create auxiliary local extrema, which regularizes the determination of the intrinsic mode function (IMF) sequence, together with the subsequent averaging, sharply reduces the influence of white noise due to its independence (Wu and Huang, 2009). This is a feature of ensemble empirical mode decomposition (EEMD). However, the selection of empirical white noise magnitude and set average number may vary from person to person. Due to the adaptivity in the two variables, this paper employs an adaptive ensemble empirical mode decomposition (EEMD) method in the GNSS coordinate time series analysis. The adaptive EEMD has significant advantages in terms of noise immunity, avoidance of local extremes, stability of IMFs adaptivity and consistency (Liu et al., 2023). It is suitable for signal components of different scales and frequencies, which improves the stability and reliability of decomposition. Nevertheless, it has high requirements for data completeness (Agnieszka and Dawid, 2022). In addition, offsets which usually occur in GNSS coordinate time series have not yet been comprehensively included in the application of adaptive EEMD.

To obtain the absolute sea level from the tide gauges, their vertical trend should be removed, which can be measured by GNSS stations nearby the tide gauges. The GNSS stations belong to the SONEL aim of providing high-precision continuous monitoring. However, the trend is affected by the periodic signals contained in the GNSS time series (Bos et al., 2010; Klos et al., 2018). Therefore, it is of paramount significance to get a periodic signal of high accuracy. This paper applies adaptive EEMD to the periodicity assessment of 13 Australian stations in SONEL. Section 2 describes the principles in detail, Sect. 3 verifies their effectiveness in synthetic time series, and 300 simulations were conducted for each of the experiments in the focused missing data and offset issues. The practical analysis for the selected 13 stations located in Australia is carried out in Sect. 4, and the conclusions are summarized in the last one.

2 Principle and method

Due to the uncertain quantity of IMFs, there may be issues of frequency overlap and mode mixing among different IMFs, which can affect the decomposition and analysis of EMD signals (Huang et al., 1999). The employed EEMD method effectively mitigates the impact of mode mixing phenomena, enhancing the precision and robustness of signal decomposition (Agnieszka and Dawid, 2022). This method consists in adding white noise to the original signal and uses the uniform distribution of Gaussian noise to alter the distribution of extreme points, ensuring signal continuity within each frequency band. Subsequently, the components of the multiple decomposed IMFs are averaged to mitigate the impact of random noise (Huang et al., 1998). The specific steps are as follows (Peng et al., 2018; Liu et al., 2023).

• Step 1: add the same length Gaussian white noise sequence to the data x(t) multiple times.

  $$\begin{array}{}\text{(1)}& X_i\left(t\right)=x\left(t\right)+u\cdot a_i\left(t\right)\end{array}$$

  u is the magnitude of the added white noise, and a_i(t) is the white noise added for the ith time, where $t=\mathrm{1},\mathrm{2},\mathrm{\dots },T$, $i=\mathrm{1},\mathrm{2},\mathrm{\dots },m$ and m is the set average number.

• Step 2: decompose the signal X_i(t) with noise to obtain the IMF components by EMD.

  $$\begin{array}{}\text{(2)}& X_i\left(t\right)=\sum _{j=\mathrm{1}}^n{c}_{i,j}\left(t\right)+r_{i,n}\left(t\right)\end{array}$$

  c_i,j is the jth IMF of EMD with white noise added for the ith time, and r_i,n is the remainder after decomposition.

• Step 3: calculate the mean value of the corresponding IMF component obtained from each decomposition to get the final IMF component of EEMD.

  $$\begin{array}{}\text{(3)}& {\mathrm{IMF}}_j\left(t\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{i=\mathrm{1}}^m{c}_{i,j}\left(t\right)\end{array}$$

  IMF_j(t) is the jth IMF component obtained after EEMD of the decomposed signal.

We can see that those two variables (added white noise magnitude u and set average number m) in Step 1 are crucial and directly affect the EEMD. For too large white noise, this can mask the real signal. For tiny white noise, its modal effect is not obvious. The larger the set average number, the smaller the decomposition error. However, this leads to a reduction in computational efficiency. Therefore, it is key to balancing both white noise amplitude and set average number, which represents the adaptivity in EEMD (Liu et al., 2023). The adaptive EEMD method can avoid setting variables manually and greatly improves decomposition efficiency and accuracy. However, the adaptive EEMD still suffers from boundary effects, which refers to the problem of unstable or distorted decomposition due to the lack of extreme points at the end. In this paper, it is mitigated by the extreme value extension method (Wu and Qu, 2008).

While preprocessing is also essential for adaptive EEMD application, which involves the selection of good continuity and sufficient data, outlier removal, offset detection and missing data filling. Here, the robust interquartile range (IQR) method is used for outlier removal (He et al., 2017), and iterative empirical mode decomposition (IEMD) is utilized for missing data filling (Qiu et al., 2022). The specific IEMD algorithm is seen in Fig. 1. All the missing data are firstly filled with the mean of the observed data to construct an initial complete time series x^′(t). Subsequently, traditional EMD decomposes x^′(t) into a series of IMF components with frequency from high to low and one residual term. Generally, the boundary IMF index k is determined by the correlation coefficient between the IMF component and the time series x^′(t) (Wang et al., 2023), at which it reaches the minimum value for the first time. Then the signals ${\widehat{x}}^{\prime }\left(t\right)$ are reconstructed by summing the IMF components after the boundary IMF component and the residual term. The missing data are updated by the corresponding part of the reconstructed signals ${\widehat{x}}^{\prime }\left(t\right)$. When the difference in the missing data |Δx^′(t)_miss| between two iterations is smaller than ε, the iterative process is terminated, resulting in the final complete time series. Otherwise, the iteration continues. ε represents a threshold (0.005 in this paper) to terminate the iteration process.

Figure 1 Flowchart of IEMD for filling.

Figure 2 Data filling of the (a) random missing 15 % and (b) consecutive missing 30 d.

With the complete data, the periodic signals can be extracted using the adaptive EEMD method. Firstly, the time series are decomposed into multiple IMFs and corresponding residuals. Secondly, the high-frequency and low-frequency IMFs are determined by the energy value method (Cheng et al., 2006). Thirdly, Lomb spectral analysis is applied to the low-frequency IMFs to identify the periodic signals (Ray et al., 2008). Finally, the IMFs containing periodic signals are selected and combined. Meanwhile, classic LS is employed for comparison. Generally, the function model includes annual, semi-annual, the first to sixth draconitic-year periodic signals (Ray et al., 2008), offsets and after-slips. The after-slip term is typically expressed in the form of logarithmic or exponential functions (Hetland and Hager, 2006; Nishimura, 2014; Tobita, 2016), and in this paper, we adopt an exponential decay function model. Specifically, the functional model can be expressed as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}x\left(t\right)& =a+bt+\sum _{i=\mathrm{1}}^{\mathrm{8}}(c_i\cdot \sin \left(\mathrm{2}\mathit{\pi }\cdot f_i\cdot t\right)\\ & +d_i\cdot \cos \left(\mathrm{2}\mathit{\pi }\cdot f_i\cdot t\right))+\sum _j{e}_i\cdot H\left(t-t_j\right)\\ & +\sum _z{g}_z\cdot \left(\mathrm{1}-e^{-\frac{(t-t_z)}{{\mathit{\tau }}_z}}\right)+n\left(t\right).\end{array}\end{array}$$

t is the observation time, a is the initial position constant, b is the linear trend, c_i and d_i are the coefficients of the periodic signal (c₁ and d₁ represent the annual periodic coefficients, while c₂ and d₂ represent the semi-annual periodic coefficients and others are draconitic-year periodic coefficients), f_i is the frequency, e_i is the offset magnitude, t_j is the moment of offset occurrence, H is the Heaviside function, g_z is the after-slip magnitude, t_z is the seismic event occurrence time, τ_z is the relaxation time and n(t) is the noise term.

Last but not least, the assessment of the signal reliability is crucial, which can be measured from different perspectives. Generally, it is measured by correlation coefficients (CCs), root mean square error (RMSE), mean absolute error (MAE), power spectral density index (κ) and signal-to-noise (SNR) (Bos et al., 2008; Chen et al., 2020; Ran et al., 2022). Their specific calculation formulas are as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{CC}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(X_i-\overline{X})(Y_i-\overline{Y})}{\sqrt{{\sum }_{i=\mathrm{1}}^N(X_i-\overline{X}{)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^N(Y_i-\overline{Y}{)}^{\mathrm{2}}}}},\text{(6)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(X_i-Y_i{)}^{\mathrm{2}}}{N}}},\text{(7)}& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left|X_i-Y_i\right|}{N}},\text{(8)}& {\displaystyle }P\left(f\right)=P_{\mathrm{0}}{\left({\displaystyle \frac{f}{f_{\mathrm{0}}}}\right)}^{\mathit{\kappa }},\text{(9)}& {\displaystyle }\mathrm{SNR}=\mathrm{10}\cdot {\log }_{\mathrm{10}}\left({\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{X}_i^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N(X_i-Y_i{)}^{\mathrm{2}}}}\right).\end{array}$$

X_i and Y_i are the true signal and prediction signal at the ith sampling point, respectively. $\overline{X}$ and $\overline{Y}$ are the means of the two signals. N is the length of the signal. P(f) is the power spectral density of the signal, f is the frequency, P₀ and f₀ are constants representing the amplitude and reference frequency and κ is the power spectral density index obtained through linear fitting. However, in the practical analysis of stations, there is no true signal present. Therefore, X_i and Y_i represent the values at the ith sampling point for the post-processed vertical components (outliers removed, missing data filled and offset corrected) and the extracted periodic signal, respectively. P(f) denotes the power spectral density of the residual signal (X_i−Y_i). Due to the true signal known in the simulated experiments, the extracted signal was evaluated using the RMSE and MAE. For the real data analysis, the CC, RMSE, κ and SNR were employed. Specifically, the higher the CC of the extracted periodic signal and the original time series, the more accurate the extraction. A smaller RMSE indicates a more accurate extraction of the periodic signal. The κ of the residuals is used to evaluate the noise characteristics and indirectly reflects the accuracy of the extracted signal, of which the lower absolute values imply more accurate extraction. Similarly, the larger the SNR, the better the periodic model established.

Figure 3 Decomposition of filled series with adaptive EEMD: (a) random missing 15 % and (b) consecutive missing 30 d.

Figure 4 The obtained periodic signal and the real periodic signal of (a) randomly missing and (b) consecutively missing.

3 Synthetic time series analysis

To test and validate the performance of the adaptive EEMD method, simulation experiments are carried out. In a simplified scenario, only annual and semi-annual periods are simulated as follows.

$$\begin{array}{}\text{(10)}& \begin{array}{rl}s\left(t\right)& =a+bt+\sum _{i=\mathrm{1}}^{\mathrm{2}}(c_i\cdot \sin \left(\mathrm{2}\mathit{\pi }\cdot f_i\cdot t\right)\\ & +d_i\cdot \cos \left(\mathrm{2}\mathit{\pi }\cdot f_i\cdot t\right))\end{array}\end{array}$$

Specifically, the 5-year data x(t) are simulated and are generated with the parameters in Table 1. To be more realistic, not only white noise but also power-law noises are included (Eq. 11). The magnitudes of white noise p(t) and power-law noise q(t) are 0.6 mm and 2 mm^−0.6. This means that the power-law spectral index here is −0.60.

$$\begin{array}{}\text{(11)}& x\left(t\right)=s\left(t\right)+p\left(t\right)+q\left(t\right)\end{array}$$

Since the missing data and offset are the general cases, these are mainly investigated in this section.

Table 1 The simulated parameters in synthetic time series (mm).

Figure 5 RMSE and MAE of (a) randomly missing and (b) consecutively missing.

3.1 Impact of missing data on adaptive EEMD

GNSS provides critical daily position time series for geodetic and geophysical studies. However, due to unforeseeable factors such as receiver malfunctions, human errors or deteriorating environmental conditions, the occurrence of missing data is inevitable (Bao et al., 2021). The IEMD method described in Sect. 2 is used in data filling due to its advantages of being efficient and accurate. However, its performance could be influenced by the percentage of missing data. Since there are no true data in the missing epochs, it is impossible to evaluate the performance. Therefore, the missing data in the simulations are artificial deletion from the original data. In practice, data missing may occur randomly or consecutively (Shen et al., 2014). To be more realistic, randomly missing and consecutively missing cases are simulated by deleting data from the original data. For randomly missing, we removed 5 %, 10 %, 15 % and 20 % from the original time series data, while for the simulations of consecutive missing we randomly deleted data for 15, 30, 45 and 60 consecutive days. Here, we present an example with a missing rate of 15 % and a consecutive 30 d missing (Fig. 2). We can see that data marked in red are filled by IEMD, which preserves the same variation of the original time series. The decomposition of the filled data using adaptive EEMD is shown in Fig. 3. From the IMFs, it is easy to identify the simulated annual (IMF6) and semi-annual (IMF5) signals. Figure 4 demonstrates the obtained periodic signal with the real periodic signal. As expected, the results indicate that, with the increase in the percentage of missing data, the extracted periodic signal gradually deviates from the true values. In order to obtain more reliable statistical results, each simulation is carried out 300 times. The mean RMSEs and MAEs of each simulation are presented in Fig. 5. As depicted in the figure, it can be observed that, with the increase in the missing data, both RMSE and MAE gradually increase, which indicates that the stability of the filling results is decreasing and may lead to larger errors. Especially when the random missing data rate reaches 15 %, RMSE and MAE show rapid increases at 0.11 and 0.05, respectively, when consecutively missing data points are between 30 and 45 and the increase in RMSE and MAE is more pronounced, at 0.32 and 0.06, respectively. Therefore, to ensure the accuracy and reliability of adaptive EEMD, the overall random missing rate should be less than 15 %, and the consecutive missing epochs should be less than 30 d.

3.2 Impact of offset on adaptive EEMD

The offsets caused by satellite orbit changes, receiver hardware or software problems, data processing algorithm improvements, and other unknown factors are common in GNSS time series (Gazeaux et al., 2013). For the known offsets, in data preprocessing, the magnitudes of offsets are estimated according to the epochs of earthquake or instrument replacement, and then the known offsets are corrected. However, for the unknown factors, it is difficult to detect the offset time due to the outliers and noises included. While the percentage due to unknown factors can be up to 20 %, its influence cannot be ignored (Griffiths and Ray, 2016). Therefore, it is urgent to identify and correct the offset. Although there are many methods to detect the offset (Wang and Herring, 2019; Tehranchi et al., 2020), neither is 100 % effective. To verify whether the offset correction is essential or not, tests (with and without offset correction) are carried out. Traditionally, manual visual inspection is combined with the STARS method for offset detection, and the LS method is used for offset correction (Rodionov, 2004; Cai, 2020). Since the number and magnitude are the two factors of offsets, these are variables in the simulations. Specifically, the offset number ranges from 2 to 8, with an interval of 2. The simulated offset magnitudes are 0.5σ, 1.0σ, 1.5σ, 2.0σ and 2.5σ mm (Su et al., 2023), with a standard deviation (σ) of the simulated time series (σ=12.70 mm). Also, the simulations are repeated 300 times. Since our focus is the extraction of periodic signals, the effect of simulated offsets on the periodic signals is analyzed. Accordingly, the RMSE and MAE between the periodic signals obtained using the adaptive EEMD method with and without offset correction and the real periodic signal are shown in Fig. 6.

Figure 6 Periodic signal RMSE and MAE with and without offset correction of different number and magnitude offsets: (a) two offsets, (b) four offsets, (c) six offsets and (d) eight offsets.

It is found that, with no correction beforehand, RMSE and MAE rapidly increase as the magnitude and the number of offsets increase. With the offset corrected, both RMSE and MAE reveal a substantial decrease. This indicates that, in adaptive EEMD application, the presence of offsets significantly leads to the reduction in the accuracy of periodic signal extraction. Further, the improvement ratio of the offset correction to without correction is shown in Table 2. We can see that, even when the offset magnitude is 0.5σ and the offset number is 2, the ratios for RMSE and MAE are 68.15 % and 67.39 %, respectively. Meanwhile, as the offset number increases, the improvement ratios for RMSE and MAE gradually rise. This indicates the necessity for offset correction even under conditions of small magnitudes. However, these findings are based on all the offsets detected as accurate, although, in reality, this is hard to realize. Thus, we recommend performing offset correction for the adaptive EEMD application in the GNSS coordinate time series analysis.

Table 2 The improvement ratios of RMSE and MAE with an offset magnitude of 0.5σ.

Table 3 GNSS station information used in this study.

Table 4 Offsets of the 13 investigated GNSS stations.

4 Results and analysis

4.1 Data sources

Thirteen GNSS stations located in the Australian region are selected for the study, and their locations are shown in Fig. 7. The vertical components for these selected stations can be downloaded from the SONEL website (https://www.sonel.org/spip.php?page=cgps&lang=en, last access: May 2023). Table 3 provides essential station information, including longitude and latitude coordinates, the time span, the raw missing rate (%) and the missing rate after gross error elimination (%), among which the longest time span is the MAC1 station and the shortest one is the BRO1 station. Additionally, the XMIS station exhibits the highest missing rate after outlier removal, which reaches 13.60 %. The known offset epochs and magnitudes available from the SONEL website, together with the offset that occurred on 12 August 2011 at the NORF station detected with the STARS method along with its magnitude estimated by LS, are shown in Table 4.

Figure 7 Location of the 13 selected GNSS stations in Australia.

Figure 8 Signals derived by the adaptive EEMD and LS of the vertical components at 13 stations.

4.2 Data processing and analysis

In this section, we carry out a comparative analysis of the periodic signal in the vertical component of the 13 selected GNSS stations derived by the adaptive EEMD method and the LS method. Figure 8 presents the post-processed vertical components (outliers removed, missing data filled and the offset corrected) and the extracted signal. It is apparent that the signal extracted is significantly different, from which the adaptive EEMD method shows its advantage in time-varying signal extraction. For example, the signals of XMIS and ESPA demonstrate an apparent peak signal in 2016 and 2015. To further evaluate its performance, comprehensive assessment indicators are displayed in Fig. 9. It is observed that, the higher the CC (Fig. 9a), the lower the RMSE (Fig. 9b), the lower the absolute of κ (Fig. 9c) and the higher the SNR (Fig. 9d) of the adaptive EEMD, which shows its outstanding advantage over LS.

Figure 9 Reliability indicators for periodic model assessment: (a) CC, (b) RMSE, (c) κ and (d) SNR.

Figure 10 Stacked comparison of power spectrum analysis.

In addition, power spectral stacking analysis is performed on residuals obtained from 13 stations. From Fig. 10 we can see that residuals derived by adaptive EEMD have the lower power of the frequency band in periodic signals than LS. However, it is also noted that the residual signal using the adaptive EEMD method still shows an annual peak. To confirm the significance of this signal within the residual signals, a power spectral significance test is conducted (Horne and Baliunas, 1986; Press et al., 2010), as shown in Fig. 11. The false alarm probability, which is used to estimate the significance of peaks in the power spectrum, indicates the probability of being a result of true signals as opposed to a random noise distribution. Generally, lower values (< 0.01) indicate higher probability. The results indicate that the 1.00 cpy signal is not significant in the residuals. In contrast, the residual signal using the LS method still shows a significant presence of the 1.00 cpy signal. This further highlights the superiority of the adaptive EEMD method in extracting periodic signals.

Moreover, the length of the time series in EEMD application is an important factor to consider, which is essential for investigating the impact on periodic signal extraction. Thus, the aforementioned 13 stations are categorized into two groups based on data length: stations with time series exceeding 10 years (mean time series length of 14.13 years) and stations with time series shorter than 10 years (mean time series length of 5.21 years). Also, significance tests are conducted on stacked residual signals for these two groups of stations. As illustrated in Figs. 12–13, the results demonstrate that, the longer the time series, the better the performance of both methods in extracting periodic signals. Furthermore, the adaptive EEMD method outperforms in both scenarios, which is no matter for data longer or shorter than 10 years.

Figure 11 Significance test of the residual signal derived from (a) adaptive EEMD and (b) LS.

Figure 12 Significance test of the residual signal derived from (a) adaptive EEMD and (b) LS for stations with data lower than 10 years.

Figure 13 Significance test of the residual signal derived from (a) adaptive EEMD and (b) LS for stations with data more than 10 years.

4.3 Discussion

Figure 14 displays all the periodic signals of the 13 stations in the vertical direction. It is evident that all the stations exhibit significant periodic patterns but with varying amplitudes. This may be significantly influenced by geophysical phenomena (Abraha et al., 2017). It is noteworthy that an anomalous peak is observed in the periodic signal of the XMIS station on 2 May 2016 (approximately 2016.35 years) characterized by an approximately twice-larger amplitude than its mean amplitude (which is approximately 0.5 cm). Similar phenomena are also observed at the ESPA station on 10 February 2015 (approximately 2015.11 years). To find a rational explanation, detailed analysis is conducted and utilized hydrological loading data from the EOST Loading Service (http://loading.u-strasbg.fr/listdata.php?dirn=_dicf, last access: July 2023), which possesses a spatial resolution of ${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$ and a time resolution of 3 h. Additionally, the daily rainfall time series provided by the Australian Bureau of Meteorology (http://www.bom.gov.au/climate/, last access: July 2023) is adopted. Rainfall anomalies are defined as the difference between the daily rainfall and its mean value. The accumulated rainfall anomalies are obtained by summing the daily anomalies with the preceding days' anomalies (Singh et al., 2021). All the possible factors of the two example stations are demonstrated in Fig. 15.

Figure 14 Periodic signals derived by adaptive EEMD of vertical components of 13 stations.

We can see that the extremely low rainfall in the second half of 2015 led to the minimum accumulated rainfall anomaly at the XMIS station on 7 May 2016. The increasing hydrological loading responded to the low rainfall but with a time lag reaching its maximum. According to the reports from the Australian Bureau of Meteorology, the presence of a strong El Niño phenomenon in 2016 (Kennedy et al., 2017) caused the spread of moist water vapor across the Australian continent. The hydrological loading decreased, which was affected by the rainfall directly. Meanwhile, the hydrological loading shows the visible annual features. The CC of 0.50 between hydrological loading data (detrended) and periodic signals for the year 2016 indicates that hydrological loading is one of the primary factors contributing to periodic changes. In addition, the ESPA station also exhibits a similar phenomenon around February 2015, with extremely low rainfall during the first 3 months of the year. This leads to a gradual decrease in accumulated rainfall anomalies, reaching its minimum on 29 March 2015. Hydrological loading during the same period reaches its current peak. The CC between the hydrological loading data (detrended) and the 2015 periodic signal is 0.57. This further emphasizes the significance of hydrological loading variations in relation to periodic changes. In summary, the primary characteristics of periodic variations are induced by hydrological loading, and other factors may also affect the periodic signals to some extent, which provides a direction for future research to explore in depth.

Figure 15 Abnormal periodic signal analysis at (a) the XMIS station (2015–2017) and (b) the ESPA station (2014–2016).

5 Summary and conclusions

This paper employed the adaptive EEMD method for more efficient and accurate extraction of periodic signals. The simulations are conducted focusing on the inevitable missing data and offset issues. The results revealed that, with the missing data filled with the IEMD method, the adaptive EEMD method exhibits efficiency in signal extraction. As expected, when the percentage of missing data increases, the extracted periodic signals may gradually deviate from the true values. Specifically, for the missing random data, when the missing rate reaches 15 %, RMSE and MAE show rapid increases at 0.11 and 0.05, respectively. For the missing consecutive data, when missing lasts for 1 month, the increases in RMSE and MAE also reach their maximum values at 0.32 and 0.06, respectively. Therefore, it is recommended that high accuracy could reach the overall random missing rate below 15 % and avoid consecutive missing epochs exceeding 30 d when dealing with GNSS coordinate time series. Regarding the offset, it is apparent that correction beforehand is essential since RMSE and MAE are markedly reduced. This also means that the offset detection and correction cannot be ignored for the adaptive EEMD application. For the processing vertical component from 13 GNSS stations in Australia, adaptive EEMD shows its advantage over LS. The mean values for CC, RMSE, κ and SNR for the adaptive EEMD are 0.36, 0.81, −0.18 and 0.48, respectively, while the mean values for LS are 0.27, 0.86, −0.50 and 0.23. Moreover, significance tests are conducted on the residuals, which further confirms that annual signals still exist in the residuals for LS but not for the EEMD. The periodic signals of the 13 investigated GNSS stations reflect the time-varying characteristics. It is found that there are significant variations of the XMIS in 2016 and ESPA in 2015. Their CC with hydrological loading are 0.50 and 0.57, respectively, indicating that hydrological loading is one primary factor contributing to the periodic variations.

In conclusion, this research highlights the advantages of the adaptive EEMD method for periodic signal extraction no matter the synthetic or real GNSS coordinate time series analysis. The handling of the missing data and offset is firstly comprehensively investigated, which provides a helpful solution to preprocessing EEMD. Furthermore, it reveals the time-varying characteristics of the periodic signal, which is favorable in the subsequent trend analysis of the GNSS stations in the SONEL network. Without doubt, the adaptive EEMD can also be applied not only to the GNSS coordinate time series, but also other geophysical data.