Layout 6 Application of empirical mode decomposition to very low frequency signals for identification of seismic-ionospheric precursor phenomena Christos Skeberis1, Dimitrios T. Xenos1, Thomas D. Xenos1, Michael E. Contadakis2, Dimitrios Arabelos2, Georgia Chatzopoulou1 1 Aristotle University of Thessaloniki, Department of Electrical and Computer Engineering, Thessaloniki, Greece 2 Aristotle University of Thessaloniki, Department of Survey Engineering, Thessaloniki, Greece ANNALS OF GEOPHYSICS, 55, 1, 2012; doi: 10.4401/ag-5312 ABSTRACT 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 earth- ionosphere 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. 1. Introduction A large amount of research over the last 20 years has indicated the existence of pre-, co- and post-earthquake ionospheric perturbations at all levels of the ionosphere (i.e. E, D, F layers). These have included ground-based experiments [Molchanov et al. 2004, 2005, Liperovsky et al. 2005, Shvets et al. 2004, Rozhnoi et al. 2004, 2009, Contadakis et al. 2007, 2008, Biagi et al., 2009, Tsolis and Xenos 2009, 2010], space-born studies [Hayakawa et al. 2000, Parrot 2006, He et al. 2011], and combined space-born and ground-based studies [Rozhnoi et al. 2007, Muto et al. 2008]. It is generally accepted that there is a lithosphere- atmosphere-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 and Attoh-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. 2. 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 199 Article history Received July 14, 2011; accepted January 19, 2012. 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 and Attoh-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 and Attoh-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: , (2) 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 3. 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]: . (6) 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. ( ) ( ) ( ) m n E n E n 2 max min 1 = + ( ) ( ) ( )h n x n m n1 1= - ( ) ( ) ( )x n x n h n1 1= - ( ) ( ) ( )h n h n m n11 1 11= - ( ) ( ) ( )IMF h n h n m nk k k1 1 1 1= = -- 10 kmr ( . )M0 43= ) 200 SKEBERIS ET AL. 201 4. 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. EMD, IDENTIFICATION OF SEISMICIONOSPHERIC PRECURSORS Table 1. Details of the European VLF transmitters monitored in the present study. Station code Location Latitude (˚N) Longitude (˚E) Frequency (Hz) GBZ Anthorn, UK 54.912 –3.277 19.580 ICV Tavolara, Italy 40.923 9.731 20.270 DHO Rhauderfehn, Germany 53.082 7.616 23.400 ICE Keflavik, Iceland 63.959 –22.542 37.500 Figure 1. Map of the transmitters, the receiver and the earthquake epicenter. Orange, defined earthquake preparation area. Figure 2. Map of the study area. Orange, defined earthquake preparation area. SKEBERIS ET AL. 202 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. 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. 203 5. 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 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. Figure 6. Signals recorded from the DHO station (Rhauderfehn, Germany). IMF1 and IMF2 received from September 2, 2009, to September 8, 2009. Green line, the seismic event. as inputs in an automatic detection system for seismic- ionospheric precursor phenomena. References Biagi, P.F., L. Castellana, T. Maggipinto, D. Loiacono, L. Sciavulli, T. Ligonzo, M. Fiore, E. Suciu and A. Ermini, (2009). A pre-seismic radio anomaly revealed in the area where the Abruzzo earthquake (M=6.3) occurred on 6 April 2009, Nat. Hazard Earth Sys., 9, 1551-1556. Biagi, P.F., T. Maggipinto, F. Righetti, D. Loiacono, L. Schiavulli, T. Ligonzo, A. Ermini, I.A. Moldovan, A. S. Moldovan, A. Buyuksarac, H.G. Silva, M. Bezzeghoud and M.E. Contadakis (2011). The European VLF/LF radio network to search for earthquake precursors: set- ting up and natural/ man-made disturbances, Nat. Ha- zard Earth Sys., 11, 333-341. Contadakis, M.E., D. Arabelos, G. Asteriadis, S. Spatalas and C. Pikridas (2007). TEC variations over the Medi- terranean during the seismic activity of 20th October, in the area of eastern Aegean, General Assembly of EGU (European Geosciences Union), Vienna, Austria, April 2007, Geophys. Res. Abstr., vol. 9,. Contadakis, M.E., D.N. Arabelos, G. Asteriadis, S.D. Spa- talas and C. Pikridas (2008). TEC variations over the Mediterranean during the seismic activity period of the last quarter of 2005 in the area of Greece, Nat. Hazard Earth Sys., 8, 1267-1276. Dobrovolsky, I.P., S.I. Zubkov and V.I. Miachkin (1979). Estimation of the size of earthquake preparation zones, Pure Appl. Geophys., 117, 1025-1044. Flandrin, P., G. Rilling and P. Goncalves (2004a). Empirical Mode Decomposition as a filter bank, IEEE Signal Proc. Lett., 11 (2), 112-114. Flandrin, P., P. Goncalves and G. Rilling (2004b). Detren- ding and denoising with Empirical Mode Decomposi- tions, Eusipco, 12th European Signal Processing Conference, Vienna, Austria, 6-10 September 2004. Gokhberg, M.B., I.L. Gufeld, A.A. Rozhnoi, V.F. Marenko, V.S. Yampolshy and E.A. Ponomarev (1989). Study of sismic influence on the ionosphere by superlong wave probing of the Earth-ionosphere waveguide, Phys. Earth Planet. Int., 57, 64-67. Hayakawa, M., O.A. Molchanov, T. Kodama, V.V. Afonin and O.A. Akentieva (2000). Plasma density variations observed on a satellite possibly related to seismicity, Adv. Space Res. Lab., 26 (8), 1277-1280. He, Y., D. Yang, J. Qian and M. Parrot (2011). Response of the ionospheric electron density to different types of seismic events, NHESS, 11, 2173-2180. Huang, N.E., Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung and H.H. Liu (1998). The empirical mode decomposition and Hilbert spectrum for nonlinear and nonstationary time series analysis, P. Roy. Soc. Lond. A Mat., 454, 903-995. Huang, N.E. and N.O. Attoh-Okine (2005). The Hilbert- Huang transform in engineering, Taylor & Francis. Liperovsky, V.A., C.-V. Meister, E.V. Liperovskaya, N.E. Va- sil’eva and O. Alimov (2005). On Es-spread effects in the ionosphere before earthquakes, Nat. Hazard Earth Sys., 5 (1), 59-62. Molchanov, O.A., M. Hayakawa, T. Ondoh and E. Kawai, (1998). Precursory effects in the subionospheric VLF si- gnals for the Kobe earhquake, Phys. Earth Planet Int., 105, 239-248,1998 Molchanov, O., P.F. Biagi, M. Hayakawa, A. Lutikov, S. Yunga, D. Iudin, S. Andreevsky, A. Rozhnoi, V. Surkov, V. Chebrov, E. Gordeev, A. Schekotov and E. Fedorov (2004). Lithosphere-atmosphere-ionosphere coupling as governing mechanism for preseismic short-term events in atmosphere and ionosphere, Nat. Hazard Earth Sys., 4 (5/6), 757-767. Molchanov, O., A. Schekotov, M. Solovieva, E. Fedorov, V. Gladyshev, E. Gordeev, V. Chebrov, D. Saltykov, V.I. Si- nitsin, K. Hattori and M. Hayakawa (2005). Near sei- smic effects in ULF fields and seismo-acoustic emission: statistics and explanation, Nat. Hazard Earth Syst., 5, 1- 10. Muto, M., T. Yoshida, M. Horie, M. Hayakawa, M. Parrot, and O.A. Molchanov (2008). Detection of ionospheric perturbations associated with Japanese earthquakes on the basis of reception of LF transmitter signals on the satellite DEMETER, Nat. Hazards Earth Sys., 8, 135- 141. Parrot, M., J.J. Berthelier, J.P. Leberton, J.A. Sauvaud, O. Santolik and J. Blecki (2006). Examples of unusual io- nospheric observations made by the DEMETER satel- lite over seismic regions, Phys. Chem. Earth, 31, 486-495. Pulinets, S.A., T.B. Gaivoronska, A. Leyva Contreras and L. Ciraolo (2004). Correlation analysis technique re- vealing ionospheric precursors of earthquakes, Nat. Hazards Earth Sys., 4, 697-702. Pulinets, S. and D.D. Ouzounov (2010). Lithosphere-At- mosphere-Ionosphere Coupling (LAIC) model. A uni- fied concept for earthquake precursors validation. J. Asian Earth Sci.; doi: 10.1016/j.jseaes.2010.03.05. Rozhnoi, A., M.S. Solovieva, O.A. Molchanov and M. Ha- yakawa (2004). Middle latitude LF (40kH) phase varia- tions associated with earthquakes for quiet and disturbed geomagnetic conditions, Phys. Chem. Earth, 29, 589-598. Rozhnoi, A., O. Molchanov, M. Solovieva, V. Gladyshev, O. Akentieva, J.J., Berthelier, M. Parrot, F. Lefeuvre, P.F. Biagi, L. Castellana and M. Hayakawa (2007). Possible seismo-ionosphere perturbations revealed by VLF si- SKEBERIS ET AL. 204 205 gnals collected on ground and on a satellite, Nat. Ha- zards Earth Sys., 7, 617-624. Rozhnoi, A., M. Solovieva, O. Molchanov, K. Schwingen- schuh, M.Y. Boudjada, P.F. Biagi, T. Maggipinto, L. Ca- stellana and M. Hayakawa (2009). Anomalies in VLF radio signals prior the Abruzzo earthquake(M=6.3) on 6 April 2009, Nat. Hazards Earth Sys., 9, 1727-1732. Shvets, A.V., M. Hayakawa, O.A. Molchanov and Y. Ando (2004). A study of ionospheric response to regional sei- smic activity by VLF radio sounding, Phys. Chem. Earth, 29, 627-637. Tsolis, G.S. and T.D. Xenos (2009). Seismo-ionospheric coupling correlation analysis of earthquakes in Greece, using empirical mode decomposition, Nonlin. Proces- ses Geophys., 16, 123-130. Tsolis, G.S. and T.D. Xenos (2010). A qualitative study of the seismo-ionospheric precursors prior to the 6 April 2009 earthquake in L’Aquila, Italy, Nat. Hazard Earth Sys., 10, 133-137. *Corresponding author: Christos Skeberis, Aristotle University of Thessaloniki, Department of Electrical and Computer Engineering, Thessaloniki, Greece; email: cskeberis@auth.gr. © 2012 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. EMD, IDENTIFICATION OF SEISMICIONOSPHERIC PRECURSORS