Application of empirical mode decomposition to very low frequency signals for identification of seismic-ionospheric precursor phenomena

This study investigates the application of empirical mode decomposition to signals from very low frequency transmitters in Europe that were received in Thessaloniki, Greece, to provide a method for depicting seismic-ionospheric precursor phenomena that occur prior to an earthquake. The basis for ionosphere interactions with seismic phenomena has been well documented in past studies, and the depiction of disturbances applied from the earthionosphere waveguide on the received signals was the purpose of this study. Empirical mode decomposition is a method for processing of nonlinear and nonstationary signals, to decompose them into their functional components, known as intrinsic mode functions. This method can provide high pass filtering to signals, thus depicting a clearer image of any abnormal disturbances in the signals that are not part of the normal noise content. Observations of such precursor phenomena are presented and correlated to earthquakes, to demonstrate the effectiveness of this method.

It is generally accepted that there is a lithosphereatmosphere-ionosphere connection mechanism that works through one of two mechanisms: (i) By the generation of atmospheric gravity waves at acoustic frequencies during the earthquake preparation period.These begin in the area of the epicenter, and subsequently move upwards, to enrich the turbulent content of the ionosphere over the epicenter and to initiate gravity waves that propagate in the waveguide of the ionosphere [Molchanov et al. 2004]; or (ii) By ion exhalation and the subsequent variations in the electric field at the site of the earthquake preparation area, which produces variations in the ionosphere over the epicenter area [Pulinets and Ouzounov 2010].
It was observed a long time ago that very low frequency/low frequency (VLF/LF) electromagnetic signals that propagate in the Earth-ionosphere waveguide are greatly influenced in their amplitude and phase by the plasma condition of the lower ionosphere [Gokhberg 1989, Molchanov et al. 1998].These can be used as an integrated diagnostic for lower ionosphere (D layer) variations over the propagation path.For this reason, the follow-up of VLF/LF electromagnetic signal propagation is considered to be a reliable diagnostic means for pre-, co-and post-earthquake ionospheric perturbations for earthquakes that might occur along the propagation path of the signal.However, the ionosphere is a very complex nonlinear system, and there are a large number of factors that can contribute to its variability, and thus this does not always provide a clear link between cause and effect.To reveal uneven disturbances and to facilitate the detection of possible earthquake-connected perturbations, the use of the Hilbert-Huang transform [Huang andAttoh-Okine 1998, Huang et al. 2005] is proposed here, or rather the part of it that is described by the empirical mode decomposition (EMD) method.

Empirical mode decomposition
Traditional filtering methods are realized in the frequency domain only.However, filtering in the frequency domain is very difficult to effectively implement for nonstationary and nonlinear signals.This is mostly because frequency analysis of Special Issue: EARTHQUAKE PRECURSORS

Subject classification:
Processes and dynamics, Earthquake faults: properties and evolution, Seismic risk, Atmospheric electric field, LAI coupling.
nonlinear and nonstationary signals generates harmonics over a wide range, and as a result, any filtering in the frequency domain will eliminate some of the harmonics.This will eventually cause deformation of the waveforms of the fundamental modes if they lie outside the filtering range.
The Hilbert-Huang transform [Huang andAttoh-Okine 1998, Huang et al. 2005] is a data-driven signal-processing method.It consists of two parts; the first is EMD, where the signal is decomposed into a series of structural components, known as intrinsic mode functions (IMFs).An IMF is defined as any function that has the same number of zero-crossings and extrema, and also has symmetric envelopes that are defined by the local maxima and minima [Huang andAttoh-Okine 1998, Huang et al. 2005].IMFs allow well-behaved Hilbert transforms, and thus the second part of the method is just the Hilbert transform, which provides instantaneous frequencies as a function of time for each of the IMF components.
Compared to standard signal-processing methods, namely Fourier transform or fast Fourier transform, as well as more contemporary signal-processing methods, such as wavelets analysis, the Hilbert-Huang transform has two significant advantages.First, it is highly applicable to nonlinear and nonstationary signal processing, as it is based on the local characteristic time scale of the data.Secondly, it is totally adaptive and data driven, as there is no need for a-priori selection of a basis (i.e., mother wavelet) for the data analysis.
For a discrete time signal x(n), the EMD starts by defining the envelopes of its maxima and minima using cubic splines interpolation.Then, the mean of the two envelopes is calculated, as: . (1) The mean is then subtracted from the original signal: , and the residual is examined according to the IMF criteria.If it is an IMF, then the procedure stops and the new signal under examination is expressed as: .
(3) However, if no IMF is identified, the signal is processed using 'sifting', and the process is continued k times, until the first IMF is realized.The sifting process is continued until the last residual is either a monotonic function or a constant.It should be mentioned that as the sifting process evolves, the number of extrema from one residual to the next drops, thus guaranteeing that complete decomposition is achieved in a finite number of steps. . (4) Then the IMF is provided as: . (5) The final product is wavelet-like decomposition that goes from higher oscillations to lower oscillations, which means that the spectrum moves to lower frequencies as the order of the IMF increases.The large difference, however, with wavelet analysis is that while modes and residuals can intuitively be given a 'spectral' interpretation in the general case, their high versus low frequency discrimination applies only locally and corresponds in no way to predetermined sub-band filtering.Instead, the selection of modes corresponds to automatic and adaptive (signal-dependent) time-variant filtering [Flandrin et al. 2004a,b]; this can provide a high pass filter for a signal that can give a clear image of abnormal disturbances in the signal that are not part of the normal noise content

Data analysis
VLF signals transmitted by a number of European VLF transmitters (Table 1; Figure 1) were monitored for over a year in Thessaloniki (40.69˚N, 22.78˚E).The signals received were sampled and stored for off-line processing.The receiver was developed by Elettronika Srl, and it is part of the International Network for Frontier Research on Earthquake Precursors [INFREP; Biagi et al. 2011, http://beta.fisica.uniba.it/infrep/).
As examples of how the above process works and of the disturbances that can be detected, the September 6, 2009, earthquake that occurred in Albania (41.46˚N; 20.41˚E; M = 5.6) is examined.As there was no appreciable seismic activity prior to the main earthquake, this provides a better example of the detection of the precursor phenomena, and also considering that the Kp and Ap indices revealed no extraordinary disturbances, with the maximum kp index for that day of 1+ and a mean Ap of 3 (ftp://ftp.ngdc.noaa.gov/).The preparation zone for this earthquake was calculated as 255.85 km, according to the Dobrovolsky equation [Dobrovolsky et al. 1979, Pulinets et al. 2004]: The distance between the location of the event and the receiver in Thessaloniki, Greece is 244 km, and therefore Thessaloniki is considered to be in the preparation zone (Figure 2).Consequently, it is expected that there would be observable disturbances in all or most of the VLF signals received.
For the purpose of the present study, the signals received were sampled (sampling rate, 1 min) over a period of 4 days prior to (from September 2, 2009), to 2 days after (to September 8, 2009) this seismic event.During this period, there was no other notable ionospheric activity that could have provided disturbances of any type.

Results and discussion.
As indicated above, four VLF transmitters were monitored in Thessaloniki: ICE (Keflavik, Iceland), ICV (Tavolara, Italy), GBZ (Anthorn, UK) and DHO (Rhauderfehn, Germany).The signals received from the first three of these stations revealed disturbances on the days previous to the earthquake and a few hours prior to the earthquake, as compared to typical relevant signals (Figures 3,4,5;respectively).
As can be seen in the signals received from ICE (Figure 3), using the EMD method, disturbances were identified 1 day prior to the seismic event.For the signal received from ICV (Figure 4), it is possible to discern disturbances up to 2 days prior to the event, and for the signals received from GBZ (Figure 5), it is possible to identify disturbances a few hours prior to the event.On the other hand, in the signals received from DHO (Figure 6), there were no appreciable differences to the normal patterns of the signals received, although this can be attributed to the recurring disturbances that were already seen in the signals received from this station, which were due to other factors that might have masked any precursor phenomenon.
As observed from the figures referring to the periods in question, the processing of the signals with the EMD method and their decomposition into their IMFs can provide an accentuation of disturbances that can be attributed to precursor phenomena and a better view of the effects of the event in question.Additionally, it should be noted that the processing of the signals with this method clears the received signals of extraneous information that is not related to disturbances.As such, this method provides a common basis for analysis across the stations, with different characteristics thus providing an array of data that can be processed easier.

Conclusions
From the above-presented data, we can conclude that the processing of signals with the EMD method and their decomposition into their IMFs can provide an accentuation of disturbances that can be attributed to precursor phenomena.Moreover, the processing of the signals received with EMD clears them from extraneous information.Consequently, this can provide data arrays that can be used EMD, IDENTIFICATION OF SEISMICIONOSPHERIC PRECURSORS  as inputs in an automatic detection system for seismicionospheric precursor phenomena.

Figure 1 .
Figure1.Map of the transmitters, the receiver and the earthquake epicenter.Orange, defined earthquake preparation area.

Figure 2 .
Figure2.Map of the study area.Orange, defined earthquake preparation area.

Figure 3 .
Figure 3. Signals recorded from the ICE station (Keflavik, Iceland).IMF1 and IMF2 received from September 2, 2009, to September 8, 2009.Green line, the seismic event, Black arrow, points of interest in the signals.

Figure 4 .
Figure 4. Signals recorded from the ICV station (Tavolara, Italy) IMF1 and IMF2 received from September 2, 2009, to September 8, 2009.Green line, the seismic event, Black arrow, points of interest in the signals.

Figure 5 .
Figure 5. Signals recorded from the GBZ station (Anthorn, UK).IMF1 and IMF2 received from September 2, 2009, to September 8, 2009.Green line, the seismic event, Black arrow, points of interest in the signals.

Table 1 .
Details of the European VLF transmitters monitored in the present study.