Noise Reduction, Atmospheric Pressure Admittance Estimation and Long-Period Component Extraction in Time-Varying Gravity Signals Using Ensemble Empirical Mode Decomposition

Time-varying gravity signals, with their nonlinear, non-stationary and multi-scale characteristics, record the physical responses of various geodynamic processes and consist of a blend of signals with various periods and amplitudes, corresponding to numerous phenomena. Superconducting gravimeter (SG) records are processed in this study using a multi-scale analytical method and corrected for known effects to reduce noise, to study geodynamic phenomena using their gravimetric signatures. Continuous SG (GWR-C032) gravity and barometric data are decomposed into a series of intrinsic mode functions (IMFs) using the ensemble empirical mode decomposition (EEMD) method, which is proposed to alleviate some unresolved issues (the mode mixing problem and the end effect) of the empirical mode decomposition (EMD). Further analysis of the variously scaled signals is based on a dyadic filter bank of the IMFs. The results indicate that removing the high-frequency IMFs can reduce the natural and man-made noise in the data, which are caused by electronic device noise, Earth background noise and the residual effects of pre-processing. The atmospheric admittances based on frequency changes are estimated from the gravity and the atmospheric pressure IMFs in various frequency bands. These timeand frequency-dependent admittance values can be used effectively to improve the atmospheric correction. Using the EEMD method as a filter, the long-period IMFs are extracted from the SG time-varying gravity signals spanning 7 years. The resulting gravity residuals are well correlated with the gravity effect caused by the Earth’s polar motion after correcting for atmospheric effects.


INTRODUCTION
Continuous observations of superconducting gravimeters (SGs) deployed by the Global Geodynamics Project (GGP) provide valuable information on the environmental changes and geodynamic processes. Due to its high sensitivity, stability and relatively low drift, the SG is capable of detecting gravity variations as small as 0.01 nm s -2 (Hinderer and Crossley 2000;Jentzsch et al. 2004). The continuous SG record has a broad spectrum, with periods varying from a few seconds to one year or longer. A global network of SGs compiles significant data for a range of important studies spanning a number of disciplines concerned with the Earth's gravity, tides, environment, and geodetics (Crossley et al. 1999). As a result of the fact that the time-varying gravity signals involve nonlinear and non-stationary processes -and not only the solid-Earth and ocean tides -certain other gravity effects, caused by the non-tidal ocean loading, by the atmospheric, terrestrial water, as well as by the Earth polar motions, are recorded. Also recorded are gravity changes due to earthquakes, free Earth oscillations, crustal movement, instrument noise and other phenomena (Crossley and Hinderer 2008). SGs are influenced by all of the effects noted above due to their high precision and broad recording spectrum. This makes SGs extremely useful, but requires much care during processing, in order to be able to separate the numerous contributions. The systematic study of one or several signal components in observed gravity time series requires the subtraction of contributions that can be calculated using the theoretical model. Subjecting the residual signals to more subsequent detailed refinement or interpretation is required to provide better analysis of the various contributions in terms of amplitude (nm s -2 ) versus period. The multi-scale analysis method is helpful in distinguishing the separate gravimetric signals of various periods.
The Fourier Transform (FT) and Wavelet Transform (WT) are two of the most widely used signal processing methods, applied broadly to time-varying gravity signal analysis and processing (e.g., Vauterin 1998;Benciolini 1994;Chao and Naito 1995;Venedikov et al. 2003;Hu et al. 2005Hu et al. , 2007. The FT shows the global characteristics of signals in the frequency domain and requires many additional harmonic components to simulate non-stationary data that are non-uniform globally (Huang et al. 1998). The FT has difficulties when confronted with certain nonlinear or nonstationary time-varying signals. As a result, the resulting FT will have little physical sense for many natural phenomena that cannot be approximated by linear and temporally stationary signals. By comparison, the WT works in both the time and frequency domains to identify the variations in features over time. This method is a localised time-frequency analysis that provides for uniform evaluation of different time scales. However, the WT is limited by the size of the basic wavelet function, its disadvantage being that uniform resolution results in uniformly poor resolution (Huang et al. 1998;Yang and Tavner 2009). Meanwhile, no general guidelines have been proposed for properly selecting the wavelet basis function (Qiu et al. 2006). Little attention has been paid to inherent deficiencies in the WT, such as border distortion, energy leakage, etc. (Peng et al. 2005).
The empirical mode decomposition (EMD) is a new multi-scale analysis method for signal processing that is intuitive, direct, a posteriori and adaptive. EMD is based on and derived from the data (Huang et al. 1998). EMD is used to decompose a signal into locally narrow band components, termed intrinsic mode functions (IMFs), which do not require any prior known basis (Huang et al. 1998). This method has unique advantages in processing non-stationary and nonlinear data. It was introduced for the processing of nonlinear geophysical signals and the extraction of anomalies (Battista et al. 2007;Huang and Wu 2008;Jeng and Chen 2011). The ensemble empirical mode decomposition (EEMD) is an improved method based on EMD, which identifies the IMFs components using an ensemble of trials, each involving the signal plus a white noise of finite amplitude . EEMD largely eliminates the mode mixing (scale-mixing) problem and the end effect caused by intermittency signal noise in the EMD method (Huang and Wu 2008), although EMD has demonstrated its applicability in a wide range of geoscience studies over the last decade (Vasudevan and Cook 2000;Huang and Wu 2008;Jackson and Mound 2010;Chen et al. 2012). When EEMD analysis is used to decompose nonlinear and non-stationary data the added white noise is averaged out over a sufficient number of trials. The persistent part that survives the averaging process is the original signal component, which is then treated as the true and more physically meaningful result. This new approach allows obtaining a uniform reference frame in the time-frequency plane and significantly reducing the problem caused by mode mixing. A number of applications for geoscience studies using EEMD (Breaker and Ruzmaikin 2011;Ehrhardt et al. 2012;Shen and Ding 2014) and algorithm research based on EEMD parameter selection have recently been developed (Niazy et al. 2009;Yeh et al. 2010;Zhang et al. 2010).
In this study an introduction to the EMD and EEMD methods is given first, with the IMF properties explained. The EEMD method is then applied as a dyadic filter bank to the gravity and atmospheric pressure records from the SG (GWR-C032) station in Wuhan, China. The original data after pre-processing are decomposed into IMFs which have a distinct central frequency. Further analysis of the relationship between the IMFs and their mean period provides a better way to perform certain applications, including those for reducing high-frequency noise, estimating atmospheric admittance based on various frequencies and extracting long-period gravimetric signals.

METHODS
EMD is an adaptive time-frequency data analysis method based on signal frequency decomposition characteristics, and thus the results may be used to retain the various characteristics of the original signal. EMD is a sifting process in which the signal is decomposed into IMFs from nonlinear and non-stationary signals based on the idea that any time series is composed of various IMFs and each IMF can be characterised by a well-defined frequency. The identification of IMFs and the principles of the EMD method have been discussed in various studies (e.g., Huang et al. 1998;Huang and Wu 2008;Franzke 2009;De Michelis et al. 2012). Therefore, in this paper we will not describe these fundamentals in detail. EEMD takes full advantage of the statistical characteristics of white noise to perturb the signal to bring it in its true solution neighbourhood. The white noise is then cancelled out after serving its purpose (Wu and Huang 2004;Flandrin et al. 2004). This method represents a substantial improvement over the original EMD method and a truly noise-assisted data analysis method (Huang and Wu 2008). EEMD is an algorithm consisting of the following steps: (1) add a white noise series to the targeted data; (2) decompose the data with the added white noise into IMFs; (3) repeat steps (1) and (2) several times, but with a different white noise series each time; and (4) obtain the (ensemble) mean values of the corresponding IMFs of the decompositions as the final result.
EEMD empirically satisfies all of the usual primary mathematical requirements of a time-series decomposition method, including convergence, completeness, orthogonality and uniqueness. It is important, particularly in filtering applications, to identify the correspondence between the IMFs and specific frequencies. When using an ideal digital filter, the cut-off frequency is usually set to limit the signal to a specific frequency band. In contrast, EEMD is based on recent white noise statistical property studies and is derived from the data. Thus, when applied to time-series data, EEMD is an effective adaptive dyadic filter bank (Wu and Huang 2004). This dyadic filter bank is a collection of bandpass filters that have a constant band-pass shape but with neighbouring filters covering half or double the frequency range of any single filter in the bank. The frequency ranges of the filters can be made to overlap; for example, a simple dyadic filter bank can include filters spanning such frequency bands as 50 -120, 100 -240, and 200 -480 Hz. However, because the actual signal is complex, the frequency distribution of an IMF is often not uniform and may lack a certain frequency band in the dyadic filter bank. When the frequency range of a certain IMF is intermittent, the next level of IMF central frequency (mean period) is then determined using the decomposition of the actual frequency band, which further demonstrates that EEMD is an adaptive process.

DATA
The gravimetric data used in this study are the longterm continuous SG observations recorded at the Wuhan stations. The data provided by the GGP database (http://ggp. gfz-potsdam.de/) include the current station's location, the instrument type, the sampling rate, the scale value of gravimetric and barometric pressure and other parameters. We start with the uncorrected data in the database, which are given in time intervals of minutes. Before using the SG data series in the multi-scale analysis study using EEMD, their spikes, gaps, steps, and large-amplitude variations caused by earthquakes are corrected using the graphical, interactive software T-soft (Vauterin 1998). After this correction, we use the EEMD method to evaluate one month of the data, including the gravimetric residuals and the barometric pressure data. The gravimetric residuals are obtained by subtracting the synthetic earth tides, which are calculated using the Tamura potential (Tamura 1987), and tidal parameters calculated using the tidal analysis program VAV (Venedikov et al. 2003). Later, the noise is reduced in the corrected gravimetric data by removing the high-frequency IMFs, and the manner in which atmospheric admittance may be effectively estimated using the relationships between the gravimetric residuals and barometric pressures in the IMFs of various frequencies are discussed.
We used 7 years (from 2002 -2008) of the continuous SG observation data to analyse the long-period gravity variations caused by perturbations in the Earth's rotation and in other geodynamic surface loads (e.g., Sun and Luo 1998). The corrected minute-interval data series are resampled by using a low-pass filter with a cut-off frequency of 12 cycle per day (cpd) to yield an hourly data series, and are decomposed into IMFs of various periods. In this study, we focus primarily on the gravimetric effects of the long-period polar motion, comparing theoretical gravity variations induced by polar motion and the long-period gravity IMFs after atmospheric IMF corrections. Theoretically, for a rigid Earth, the gravimetric variations induced by polar motion can be predicted by the relationship (e.g., Wahr 1985): where R (6371 km) is the mean radius of the Earth, X (7.292115 × 10 -5 rad s -1 ) is the mean rate of the Earth's rotation, θ and m are the latitude and west longitude of the observation site, respectively, and x(t) and y(t) are the instantaneous polar coordinates of the Celestial Ephemeris Pole relative to the International Reference Pole. We use the EOPC04 polar motion data set, provided by the International Earth Rotation Service (IERS), and Eq. (1) to determine the theoretical polar tides at the Wuhan SG stations.

Reductions of Noise in Gravimetric Data
We need to pre-process the raw data before analysing the SG time series with high precision (0.01 nm s -2 ) and high sensitivity. The purpose of the pre-processing is to eliminate bad records such as spikes, steps, and large-amplitude vibrations caused by earthquakes. Still, the results cannot be fully scrubbed of artificial operational effects. The traditional method used for removing noise is through Fourier analysis-based filtering in the frequency domain. However, the traditional method has shortcomings due to the underlying stationarity and linearity assumptions. Using the developed adaptive time domain of the EEMD method as a filter, we are able to reduce the noise variations involving nonstationary and nonlinear processes. The result of EEMD time space filtering preserves the full nonlinearity and nonstationarity in the physical space ). For example, the result of low-pass filtered signal having IMFs components can be expressed as: The high-pass filtering result can be expressed as: and the band-pass filtering result can be expressed as: We perform the EEMD on SG data corresponding to one month (sample rate is 1 minute). The resulting IMF components are shown in Fig. 1. In the first two IMFs of high-frequency components, one can clearly see the residual error signals, which are caused by electronic device noise, Earth background noise and the residual effects of pre-processing. The subsequent 3 rd -10 th IMF long period components belong to the tidal frequency band. Thus, based on the characteristics of each component of low to high frequencies, we reconstruct the IMFs so that they reach the low-pass effect.
To examine the filtering effect in detail, we plot, as shown in Fig. 2a, the Fourier spectra corresponding respectively to the observed data, the theoretical tides, the EEMDfiltered data (the gravity signals without IMF1 and IMF2), and the removed noise (IMF1 and IMF2). The spectral magnitudes of the theoretical tides, of the original records, and of the filtered data are almost identical in the low-frequency band (i.e., these datasets contain the primary tidal-wave fluc-tuations, such as diurnal, semidiurnal, and third-day waves). All of the substantial differences are in the higher frequency range. They consist of the first two IMFs (IMF1 + IMF2; the curve has been shifted by a factor of 100 for better visibility in Fig. 2a). Specifically, the spectra after the filter indicate a drastic decrease in energy at frequencies higher than 100 cpd, whereas the spectral density exhibits a negligible impact in the low-frequency range.
The difference between the filtered and unfiltered SG data is shown in Fig. 2b. This comparison indicates that the high frequency component (IMF1 + IMF2, after removing the theoretical tides from the observed data) is dominated by the Earth's background noise. This noise is either in the form of spikes or intermittent noise, with a peak-to-peak range of 10 nm s -2 . Usually earthquake related disturbance signals lasting hours are removed first as noise in the data pre-processing. However, the large magnitude of the removed noise appears near the middle of the dataset and may be due primarily to the effects of earthquakes for a few days. These processing results also indicate that using the EEMD filter is an adaptive process, determined a posteriori based on criteria derived from the data, in contrast with the selection of basic wavelets and cut-off frequencies required in wavelet and Fourier analyses. Using the mean IMF periods we are able to recognise the underlying processes of a given time-varying gravimetric signal and identify the relationships between the Fig. 1. The intrinsic mode function (IMF) components of the SG (GWR-C032) data obtained using the ensemble empirical mode decomposition (EEMD) method. The central frequencies of IMF1 and IMF2 higher than 100 cpd and the residual represent the overall trend of the signal.
IMFs and the various scales in the information.

Estimates of Atmospheric Admittance
The variable atmospheric gravity load is the second largest effect (up to 10% of Earth tides) in the SG time series. It spans a wide frequency range from minutes to seasonal periods. Merriam (1992) found that about 90% of the atmospheric gravity comes from the local region of the gravimeter (i.e., a cylinder of 50 km radius with the same pressure as the station). For a good approximation, a scalar admittance relates observed gravity to local pressure variations. Earth tides are efficiently subtracted from SG records using synthetic local tides and the gravity residuals are strongly negatively correlated with the barometric pressure data (see input data in Fig. 3). However, removing the atmospheric pressure effects precisely from gravity measurements is not easy (Hu et al. 2005). The atmospheric admittance estimation is based on fitting the barometric pressure data to the gravity residuals. However, the barometric signal also exhibits nonlinear and non-stationary characteristics, and the signals in various frequency bands correspond to various energy levels. Thus, a single atmospheric admittance cannot be used to properly correct the gravity effects of atmospheric pressure. Warburton and Goodkind (1977) first suggested a frequency-dependent admittance might be a better approximation. The idea was revived by Neumeyer (1995) and Crossley et al. (1995) to improve the accuracy of the atmospheric correction in the sense of reducing the amplitude of the corrected signal.
In this study we obtained the atmospheric pressure effects on the gravity variations (IMFs) in various frequency bands to estimate the time-and frequency-dependent admittance based on the reconstruction of the IMFs in various frequency bands. The central frequency of the IMFs acts essentially as a dyadic filter bank resembling those involved in wavelet decompositions. The decomposition results (the 1 st -14 th IMFs) are divided into various frequency ranges: the first three IMFs, for periods shorter than 30 minutes, are part of the so-called seismic frequency band (400 -50 cpd); the subsequent 4 th -6 th IMF components, for periods shorter than 4 hours, belong to the long-period seismic frequency band (50 -6 cpd), which is of interest taking into account the long-period of Earth's free oscillations; the IMF components from 7 th -9 th for the periods in tidal band (6 -0.5 cpd); IMF components from 10 th -11 th are for periods equal or longer than 2 days, which belong to the so-called the longperiod band (0.5 -0.125 cpd); finally, the last three IMFs and residual (R) represent the overall trend of gravity and atmospheric pressure (less than 0.125 cpd) (see Fig. 3).
The correlation coefficients between the gravimetric and the atmospheric pressure, as well as the atmospheric admittance in each band are listed in Table 1. In the seismic frequency band the IMFs (IMF1 -IMF3) reflect the highfrequency gravimetric signals from ocean noise, seismic background noise and instrument noise. The correlation between the gravity residuals and the atmospheric signal is very low. These results show that the small atmospheric loads have little effect on the gravity in the high frequency band. The degree of similarity between the gravimetric and atmospheric signals is then seen to begin to increase in the long-period seismic frequency band (IMF4 -IMF6). The results indicate that the local atmospheric pressure starts to influence the gravimetric variations in this band: the overall (a) (b) Fig. 2. Comparison between the Fourier spectra and amplitudes of observed data, theoretical tides, EEMD-filtered data and the removed noise. (a) The Fourier spectra of the observed data (blue), the EEMD-filtered data (red), and the theoretical tides (black) show a good agreement in the relative low frequency (less than 100 cpd). The first two IMFs (IMF1 + IMF2, the black curve) have been shifted by a factor of 100 for better visibility.
(b) The comparison of EEMD filter results, including the differences between the theoretical and the observed data, and the substantive differences between the observed and the EEMD-filtered data, are all in the high-frequency range, with a peak-to-peak range of 10 nm s -2 .
correlation between the gravity and pressure is -0.6165, whereas the admittance is -3.612 nm s -2 hPa -1 . The longperiod seismic frequency band is of interest, as the longperiod free-Earth oscillations and atmospheric fluctuations play an important role in the gravity residual variations, which provide indications of the long-period normal mode and spectral splitting. In the tidal band (IMF7 -IMF9), the gravity residual variations are caused primarily by the global effects of atmospheric tides S 3 , S 2 , and S 1 . The correlation associated with this tidal band is clearly higher than that of the long-period seismic frequency band and the admittance is -2.79 nm s -2 hPa -1 . Because of very weak tidal Fig. 3. Comparison of reconstructed IMFs based on gravity and atmospheric pressure data from various bands. components and strong local atmospheric pressure signals in the long-period band (IMF10 -IMF11) and in the secular trend (IMF12 -R), the gravity variations highly correlate with local atmospheric pressure fluctuations and the admittance is -2.116 and -8.069 nm s -2 hPa -1 , respectively. Especially in the secular trend band there is a very high overall correlation of -0.971 between pressure and gravity, and the pressure fluctuations fit the gravity variations very well in the time domain (see Fig. 3, IMF12 -IMF14 and residual). We can clearly see that there are strong gravity variations on longer time scales than a few days or even one month. Thus, the atmospheric pressure is corrected in the low frequency (IMF12 -R) IMFs of gravity variations by multiplying the admittance value of -8.069 nm s -2 hPa -1 for this band. Needless to say, there are possible variations in admittances on time scales from minutes to days caused by the local atmosphere. They are clearly variable on these short time-scales. Local weather systems can also move rapidly over a gravity observation station within a few hours (Müller and Zürn 1983). Furthermore, a single admittance factor fails to account for the entire load effect that arises from the total pressure distribution (e.g., local, regional, and global contributions). The residual gravity using a simple admittance will be unreliable at the diurnal harmonics or the longwavelength contributions. This is because the solar harmonics for a long-period exist at a global scale, although they can still be correlated to the local pressure. The amplitude of the gravity residuals is reduced from 100 -20 nm s -2 (the red line in Fig. 4) after subtracting the various atmospheric gravity effects (pressure multiplied by admittance) from the gravity signal in the various bands. These results demonstrate that the effectiveness of the time-and frequencydependent admittance. The EEMD method can effectively help with providing a correction for the atmospheric gravity effect. However, Kroner and Jentzsch (1999) pointed out that the smaller amplitude value was not necessarily better, due to the probable contamination of the gravity residuals by other signals such as hydrology.

Analysis of Long-Period Signal
The work described above has demonstrated EEMD capability in terms of filtering the high-frequency and tidal bands. In analysing the long-period gravity signals caused by polar motion and surface loading (e.g., atmospheric and hydrological effects), we considered 7 years (from 2002 -2008) high-quality gravity and barometric records from the Wuhan SG station to confirm that the EEMD method is suitable for extracting the long-period gravity and pressure components of SG signals.
The main purpose of this study is to verify the advantages of multi-scale analysis using EEMD on broad spectrum, non-stationary and nonlinear gravity data. As the impact of environmental factors such as hydrological effects and non-tidal ocean loading are not considered, the gravity residuals after atmospheric correction still contain these effects. Fifteen IMFs were obtained using the EEMD method to decompose the data. Only the power spectra of the IMF11 and IMF12 are shown in Fig. 5, because the frequencies of IMF1 -IMF11 are more than 0.004 cpd (periods are less than 256 d), and we only consider the frequency band range from 0.004 -0.001 cpd as being of interest (periods of polar motion). For instance, the low-frequency components (long-period IMFs, IMF12 -IMF15) may be thought of as constituting long-period, non-tidal and other overall patterns, whereas the high-frequency components (shortperiod IMFs, IMF1 -IMF11) may be thought of as noise and as gravity effects of earthquakes and tides (see Fig. 6). The data related to atmospheric pressure are also filtered in a similar manner to estimate the atmospheric effects in the long-period band.
The gravity residuals were obtained by subtracting the atmospheric pressure load effects from the long-period IMFs. The gravity variations induced by polar motions are theoretically predicted using Eq. (1) and the EOPC04 polar motion data (see Fig. 7). The polar motion consists of two primary frequency components, including the Chandler Wobble and the forced annual wobble, which recur at periods of approximately 432 and 365 d (Harnisch and Harnisch 2006;Hu et al. 2007), respectively. The peak-topeak magnitude of the gravity fluctuations induced by these Fig. 4. The results after correcting for the atmospheric gravity effect. Fig. 5. The comparison of Fourier spectra between the IMF11 (black) and IMF12 (red). polar motions is approximately 100 nm s -2 . Figure 7 shows a comparison between the gravity residuals and the gravity effects of this polar motion. The results show high correlation between observed data and theory predicts. The high-precision SG may be used to clearly identify the gravitational field perturbations caused by the Earth's polar motion. The polar motion gravity effects cannot be ignored when using the time-varying gravity signals in long-period or long-term geodynamic phenomena studies (e.g., seasonal variations in terrestrial water and crustal movements).

CONCLUSIONS
Gravimetric information from various time scales can be identified and isolated owing to the advantages of EEMD multi-scale decomposition, significantly increasing the credibility of gravity-effect corrections and interpretation during any post-processing. EEMD was used in this study to identify and analyse surface SG gravimetric and atmospheric pressure time series data. By decomposing gravity and pressure data into dyadic frequency bands using EEMD band-pass filters, we obtained a series of IMFs in which the central frequency displayed binary variations from high to low. Combining IMF features we carried out three application aspects. First, we applied a high-frequency filter to one month of SG data in order to reduce the continuous gravimetric observation data noise and the uncertainty caused by pre-processing, greatly reducing level of noise in the data. Second, EEMD was used to take into account both the time and frequency information of the atmospheric pressure effects. The atmospheric admittance estimates reduced the atmospheric effects from the gravimetric records in the search for weak geodynamic signals. Finally, using a low-pass filter to obtain the long-period gravimetric signal effectively retained the polar-motion information. The results indicate a high degree of correlation between the gravity residuals after correcting for the atmospheric effects and polar motion gravimetric effects.
However, the environmental effects (e.g., the global non-ocean loading, the influence of terrestrial water) in the tidal and pole tidal band could not be removed using the EEMD method. To achieve further improvement these effects should be removed by introducing more auxiliary data and environmental mathematical models. In our view, EEMD offers a potentially viable method for time-varying gravity data analysis. In addition to considering the influence factors of some basic principles (e.g., the selection of the spline, the extension of the end), the physical meanings of each IMF component and the best IMF selection and  scale of multi-scale decomposition are the most important aspects of scientific research in this area. These gravity research aspects may constitute new research possibilities for areas of future study.