Natural Hazards and Earth System Sciences Stochastic analysis of geo-electric field singularities as seismically correlated candidates

The study of the Earth’s electromagnetic field prior to the occurrence of strong seismic events has repeatedly revealed cases where transient electric potential anomalies, often deemed as possible earthquake precursors, were observed on electromagnetic field recordings. In an attempt to understand the nature of such signals, several models have been proposed based upon the exhibited characteristics of the observed anomalies, often supported by different mathematical models simulating possible generation mechanisms. This paper discusses a candidate Electric Earthquake Precursor (EEP) signal, accompanying the Kythira Mw=6.9 earthquake in Greece (occurred on 8 January 2006). Neuro-Fuzzy along with stochastic models are currently incorporated for the modelling and analysis of the recorded Earth’s electric field. The results of the study indicate that the Neuro-Fuzzy model treats the observed possible EEP signal as an external additive component to the recorded Earth’s electric field, while the stochastic TARMA models accurately represent the recorded electric signals in both the time and the frequency domains. The complementary findings of both methodologies might potentially contribute to the future development of a more accurate and generalized framework for the efficient recognition and characterization of possible EEP’s.


Introduction
Earthquake prediction remains one of the biggest challenges in modern earth science, as it aims to predict origin time, hypocenter (or epicenter) and magnitude of an impending earthquake.Intermediate and long-term predictions are mainly based on past and current records of seismic activity Correspondence to: F. Vallianatos fvallian@chania.teicrete.gr(earthquake forecasting).Short-term earthquake prediction, in addition, relies on the detection of precursors (seismoassociated phenomena).Only short-term and intermediateterm time scales can be considered as true earthquake prediction deterministic methods.During the last couple of decades, a variety of types of earthquake precursors has been increasingly reported by many researchers, with predominant category that of pre-seismic and/or co-seismic variations of the Earth's electric, magnetic and electromagnetic fields Park et al. (1993).The latter seems to cover almost the entire electromagnetic spectrum, but the most promising lay on the ULF-ELF-VLF bands Hayakawa (2004).
Ground-based measurements reveal subsurface transient electrotelluric and/or magnetic signals Varotsos (2005); Console (2001), whilst radio and ionospheric soundings along with space observations, manifest VLF/LF subionospheric signals, ionospheric perturbations and electromagnetic phenomena in the coupled system lithosphere-atmosphereionosphere associated with large seismic events Park et al. (1993).Furthermore, scientists have tried to verify and explain most of the observations mentioned above by conducting laboratory experiments on rock samples Hayakawa, 2004;Varotsos, 2005;Vallianatos et al., 2008Vallianatos et al., , 2004;;Stavrakas et al., 2004.Several physical generation mechanisms have been proposed such as electrification due to microfracturing, fluid diffusion, electrokinetic effects, charged particle generation and motion, moving charged dislocations, electrical conductivity changes, piezo-stimulated currents, piezomagnetic and piezoelectric effects etc., with secondary effects in the atmosphere-ionosphere coupled system such as total electron content changes, ionospheric motions, joule heating etc., but no definitive conclusions have been drawn yet Tzanis and Vallianatos (2002).
A. Konstantaras et al.: Stochastic analysis of seismic related geo-electric singular Despite the richness of observed earthquake precursory phenomena, significant scientific efforts are still required to ensure qualitative and quantitative rigorous definitions of earthquake precursors that could lead in developing robust and reliable earthquake prediction methods and this is mainly due to lack of thorough understanding the physics of pre-seismic processes, the physics governing the generation, propagation and observation of the various precursory phenomena and the physics of the earthquake source, as well as due to large recurrence times involved Varotsos, 2005;Console, 2001;Colangelo et al., 2000Colangelo et al., , 2001. .In the present study, the observed characteristics and features of a recorded candidate Electric Earthquake Precursor (EEP) signal of a strong earthquake in Greece, are investigated.Previous studies Konstantaras et al., 2007Konstantaras et al., , 2002Konstantaras et al., , 2004Konstantaras et al., , 2006a,b ,b provide information regarding the nature of EEP signals.These results indicate that EEP signals are transient electric potential anomalies external to the natural (of ionospheric origin) electromagnetic field of the Earth.
Despite these similarities, there is an open issue as it concerns this particular candidate EEP signal accompanying the Kythira earthquake, namely the fact that it was only observed upon the recordings of the electric field, whilst there is no indication of the latter upon the simultaneous recordings of the magnetic field.To enhance and strengthen the detection and identification of the observed transient electric field anomaly as a possible EEP candidate signal associated with the Kythira earthquake, two approaches are currently addressed: ( The aforementioned approaches are mutually complementary, as non-linear Neuro-Fuzzy systems focus on the signal's lower frequency band and any kind of non-linearities (i.e.transients, etc.), while Time-dependent ARMA (TARMA) models aim at capturing the signal's higher frequency components and any kind of non-stationarities.The rest of this paper, is organized as follows: the recording station along with the recorded data are presented in Sect. 2. Neuro-Fuzzy and stochastic modelling and analysis results are discussed in Sects.3 and 4, respectively, and the concluding remarks of the study are finally summarized in Sect. 5.

Data description
The present study focuses on providing evidence on the dynamical structure of the recorded signal aiming to reinforce its identification as a possible EEP associated with the M w =6.9 Kythira earthquake occurred on the 8 January 2006, located 36.21 • N (latitude) and 23.41 • E (longitude), with an epicentre depth of approximately 70 km.The signal under investigation was recorded by the MVC-2DS recording station, which was designed and developed by the Institute of Terrestrial Magnetism, Ionosphere and Radio-wave Propagation, Saint Petersburg Filial (SPbF IZMIRAN), Russian Academy of Sciences, in the frame of the INTAS-99-1102 project entitled: "Study of the ULF electromagnetic phenomena related to earthquakes (SUPRE)", Kopytenko et al., 2001a,b.The MVC-2DS recording station Kopytenko et al. (2001b); Vallianatos et al. (2002), measures in GPS stamped time the horizontal components of the electric field and the three (two horizontal plus vertical) components of the magnetic field of the Earth, with a variable sampling frequency, currently set to 5 Hz.In particular, magnetic field recordings are obtained using a highly sensitive torsion photoelectric magnetometer, whilst electric field recordings are obtained using two electric dipoles that have length of about 100 m and are deployed in the NS and EW direction, respectively.The electric dipoles use Pb electrodes, buried at a depth of about 1 m.The station has been in operation since 2001, which has enabled us to distinguish various distortions from random and/or artificial noise sources.Indeed, it is characterized by self-consistency as shown in Figs. 1 and 2 (Fig. 1 shows the recorded Earth's electromagnetic fields during June 2006, while magnetotelluric signatures are evident in Fig. 2).
The location of the recording station with respect to the main earthquake and the associated seismic sequence are presented in Fig. 3, see also Vallianatos et al. (2006).The recorded signals under investigation, last from approximately 9 p.m. on the 29 December 2005, till 1 p.m. on the 11 January 2006, Fig. 4. It is interesting to note, that while a candidate EEP is evident in both the electric field recordings (Fig. 4a  and b), there is no obvious indication of it upon the simultaneous magnetic field recordings (Fig. 4c and d, the square in all subplots of Fig. 4, indicates the occurrence time of the main seismic event).
Figure 5a demonstrates the seismic sequence in parallel with the observed electric transient potential anomaly, Fig. 5b.No foreshocks have been observed before the main earthquake whilst there were only a few aftershocks.It is interesting, though, to note that the main swarm of the aftershocks practically ends with the return of the electric field recordings to the background level, although the main event occurred almost in the middle of the signal's "duty cycle".It is also worth noting the large duration of the recoded anomaly and the fact that it outlasts the main seismic event.

Neuro-fuzzy modelling and analysis
Neuro-fuzzy models are neural networks with intrinsic fuzzy logic abilities where each layer of the network emulates the input Membership Functions (MFs), rules, output MFs, and defuzziffication function of a fuzzy inference system, respectively Jang (1993).The use of Hybrid Adaptive Neuro-Fuzzy Inference Systems aims for the recovery of the possible EEP signature from the electric field background Varotsos (2005); Konstantaras et al. (2007), which enables significant information to be extracted regarding its nature and possible association with the accompanying main seismic event.
To investigate whether this EEP candidate signal is indeed an external addition to the natural electric field of the Earth, a pattern recognition experiment has been devised with the incorporation of soft computing technology.A neuro-fuzzy model, i.e. a neural network with intrinsic fuzzy logic abilities Jang (1993) has been developed and trained to identify the recorded electric field signal using the data recorded prior to the occurrence of the possible electric earthquake precursor.Then, propagating through the electric field recordings, the neuro-fuzzy model is used to forecast the next sample of the recorded signal based upon a number of previously recorded data.The purpose of the experiment is to identify whether the neuro-fuzzy model follows the detected EEP signal as if it was part of the natural electric field of the Earth; or rejects it as an external additive component by considerably suppressing the EEP, aiming for the actual value of the natural electric field alone.

Training and evaluation procedure
To train and evaluate the reaction of the neuro-fuzzy model, 4096 data samples of electric field recordings have been selected Konstantaras et al. (2004), corresponding approximately to the time-period starting at 9 p.m. on the 29 December 2005 and ending at 1 p.m. on 11 January 2006, which include the possible electric earthquake precursor.Although the initial sampling frequency of the recorded data is f s =5 Hz, the overall data set has been decimated by a factor of 1280 as it is very costly in terms of processing time Kosko (1991) to train a neural network with such a heavy workload.A sliding window consisting of four previous inputs, at n−12, n−24, n−36 and n−48 Addison and Wermter (2002), propagating through the time-series, determines the input vectors fed to the neuro-fuzzy model.The first half of the input vectors, i.e. the first 2048 samples in the time-series, is used to train the neuro-fuzzy model to predict the next sample at time (n+1) in the time-series Jang et al. (1997), whilst the second half remains unseen during training.An initial neuro-fuzzy model is obtained by applying grid partitioning Konstantaras et al. (2002) on the first half of the input data set.This initial model is subjected to training with a hybrid algorithm Jang (1993), a combination of the least squares method and the back-propagation algorithm.During a forward pass, an input vector is fed to the neuro-fuzzy model and the least squares estimator is used to adapt its consequent parameters, which define the rules and output membership functions (MF's) of the model.A training error is computed by subtracting the output of the neuro-fuzzy model, for the current set of parameters, from the required output (the actual value of the electric field signal at sample n+1).The training error is deployed during the backward pass through the neuro-fuzzy model by the back-propagation algorithm to adapt its premise parameters, which determine the shape and dimensions of the input MF's.After every training epoch the neuro-fuzzy model is tested with the first 500 samples of the unseen data (also recorded before the occurrence of the possible electric earthquake precursor) to prevent overtraining Jang et al. (1997).The final neuro-fuzzy model holds the set of parameters that minimized the checking error.

Neuro-fuzzy model development and operation analysis
To generate an initial fuzzy inference system, grid partitioning was applied on the input data of the input/output data set.For this particular application the two-dimensional input space of every input was partitioned into four overlapping fuzzy regions, each of which is governed by a fuzzy ifthen rule.The structure of the neuro-fuzzy model depicted in Fig. 6, depends on the number of inputs and input membership functions per input.The developed neuro-fuzzy model has four inputs (layer 1), each of which has assigned to it two input membership functions (layer 2), and is guided by   sixteen rules (layer 3).The contribution of each rule to the output of the neuro-fuzzy model is determined by the output MF (layer 4) allocated to it, whilst the bias neuron (dashed line in Fig. 6) sets a weighting factor to each rule.The neuron in layer 5 defuzzifies the normalized weighted outputs of all rules to produce a crisp output (layer 6).The layer-by-layer operation of the developed neuro-fuzzy model is described as follows Jang (1993): -Layer 1: The present and three previous samples of the recorded electric signal are used as inputs (A to D) to the network.
-Layer 2: Every node i in this layer is an adaptive node with a node function: where x (or y, or z, or k) is the input to node i and A i (or B i , or C i , or D i ) is the equivalent membership function denoted by µ(•).The type of MF's A, B, C and D is that of the generalized bell function: where {a i , b i , c i } are the premise parameters of the network which determine the shape and size of the MF, Jang (1993).
-Layer 3: Every node in this layer is a fixed node calculating the normalized firing strength of either rule: where: -Layer 4: Every node i in this layer is an adaptive node using an output membership function to compute the weighed output of the equivalent rule, according to the following node function: where {p i , q i , m i , n i , r i } are the consequent parameters of the network that specify the rules of the fuzzy inference system, Jang (1993).
-Layer 5: The single node in this layer is a fixed node, which converts the weighted fuzzy outputs of all rules in the system into a single crisp output, as described by the following node function: -Layer 6: The node describes the actual output of the neuro-fuzzy model for a given input data set.

Neural-fuzzy pattern recognition results
The proposed neuro-fuzzy model was initially trained upon the first 2548 data samples of the recorded electric field signal (E x ), aiming to identify the main characteristics of the electric field variations and thus predict the next sample in the time-series.The initial 2048 data samples were used for training the neuro-fuzzy model whilst the next 500 data samples (2049 to 2548) remained unseen for validation purposes.
The outcome of the neuro-fuzzy pattern recognition application on the full electric field signal after training the neurofuzzy model is shown in Fig. 7.The comparison of Fig. 7a and b, reveals a significant suppression of the apparent magnitude of the possible recorded EEP signal.In detail, the output of the neuro-fuzzy model closely follows the recorded electric field signal until the moment of the occurrence of the possible EEP at approximately sample 2752.The rapid rise in magnitude of the recorded signal over the next few samples "confuses" the neuro-fuzzy model and for a short time it becomes unstable, hence the large spikes observed between samples 2780 and Nat.Hazards Earth Syst.Sci., 8, 1451Sci., 8, -1462Sci., 8, , 2008 www.nat-hazards-earth-syst-sci.net/8/1451/2008/

Stochastic FS-TARMA modelling and analysis
Non-stationary stochastic signals, that is signals with timedependent (evolutionary) characteristics, are frequently encountered in engineering, and have been receiving increasing attention in recent years Newland (1993); Kitagawa and Gersch (1996); Fouskitakis andFassois (2001, 2002).Typical application areas include seismic and structural systems, mechanical vibrations, speech analysis, reliability analysis, automotive and aircraft systems, rotating machinery, and so on.From a physical standpoint, non-stationarity stems from either time-dependent dynamics and transient phenomena and/or inherently non-linear behavior, Poulimenos and Fassois (2006).Functional Series (FS) methods, are physically motivated for many engineering systems where the evolution of the signal characteristics cannot be regarded as random.The imposition of maximum (deterministic) structure on parameter evolution (through properly selected functional spaces) allows Functional Series methods to achieve: (i) High parsimony, (ii) tracking of "fast" and "slow" evolution, and (iii) high accuracy and resolution.
Functional-Series (FS) Time-dependent (T) AutoRegressive (AR) Moving-Average (MA) (FS-TARMA) models, provide an important generalization of their stationary ARMA counterparts Pandit and Wu (1983), in that their parameters are explicit functions of time, belonging to functional spaces spanned by sets of preselected deterministic functions of time, referred to as the basis functions.A FS-TARMA(na, nc) p model, with (na, nc) indicating its Au-toRegressive (AR) and Moving-Average (MA) orders, respectively, and p the corresponding functional base dimensionality, is of the form (referred to as the shifted form): Hence, the time-dependent AR and MA parameters are expressed in terms of the mutually orthogonal basis functions as: with a i,j , c i,j representing the corresponding coefficients of projection.Functional spaces comprised by shifted polynomials (such as Chebyshev, Legendre, Laguerre, and so on Abramowitz and Stegun, 1970), trigonometric (sine/cosine) functions or discontinuous Haar or Welch functions Fouskitakis and Fassois (2002), may be considered In the case where no MA part is included in the model, the model of Eq. ( 7) drops to the FS-TAR case, in which it may be easily shown that the model is linear in its pa- rameters, and thus its coefficients of projection may be estimated via ordinary least-squares Poulimenos and Fassois (2006).In the full TARMA case, its coefficients of projection may be estimated via the Generalized Polynomial-Algebraic (GPA) method Fouskitakis and Fassois (2001), or the Two-Stage Least-Squares (2-SLS) method Poulimenos and Fassois (2006).Model selection may be based upon the Bayessian Information Criterion (BIC) Fouskitakis and Fassois (2001), which strikes a compromise between achievable accuracy and model complexity.Final model acceptance (validation) is based upon formal assessment of the uncorrelatedness hypothesis of the model's residuals.

FS-TARMA model estimation and identification results
The basic aim of the stochastic TARMA modelling and analysis procedure, is to reveal and accurately represent the recorded electric field's inherent dynamical structure -in both the time and the frequency domains -during (as indicated by the dashed frames shown in Fig. 8a and c) the possible EEP (when the recoded electric field is "high").The rejection of any kind of homogeneous non-stationarities (such as trends, transients, mean value non-stationarity, cyclostationarities, etc.) from the recorded electric field signals, is achieved by taking the signal's (raw data) first-order differences Brockwell and Davis (1996).The differenced signals are subsequently low-pass filtered in order to reduce noise effects and sample size (and thus computational complexity), and finally: f s =0.0278 Hz, T s =36 s, N=4160.
The resulting E x , E y signals shown in Fig. 8b and d, are thus modelled by full FS-TARMA models (estimated by the 2-SLS estimation algorithm, Poulimenos and Fassois (2006)).The application of several TARMA models of various orders (na, nc), various functional spaces (spanned by trigonometric, discontinuous Haar and shifted orthogonal polynomials) and various subspace dimensionalities (p), leads (based upon the BIC criterion) to a TARMA(40, 20) 4 model for both E x , E y , the functional subspace of which is comprised by the first 4 Haar functions in both cases.The selected TARMA models constitute very accurate representations -in the time domain -of the signals under study, as they achieve very low percentage RSS/SSS ratios [RSS: Residual (1-step-ahead prediction error) Sum of Squares; SSS: Series Sum of Squares], RSS/SSS=1.6381%for the E x case; RSS/SSS=1.6761%for the E y case, see also Fig. 9a  and c, respectively.They are also formally validated, as the normalized sample autocorrelation function of their residuals, falls within the limits of statistical insignificance at the α=0.05 level, Fig. 9b and d.
The "frozen" spectrum Fouskitakis and Fassois (2002), of the selected TARMA models are now presented in Fig. 10, and compared with their non-parametrically estimated (by the Short-Time Fourier Transform (SFFT), Fouskitakis and Fassois, 2002), counterparts.It is interesting to note, that in both cases the TARMA model-based spectrum is in very good agreement with its non-parametric counterpart.In addition, high frequency components are more evident present after sample 3000.

Concluding remarks
The observed anomaly accompanying the Kythira M w =6.9 earthquake, demonstrated some special characteristics and features with respect to earlier observations of possible electric earthquake precursors Stavrakas et al. (2004), namely its large duration, the fact that it outlasts the main earthquake and the lack of any signature on the simultaneous magnetic field recordings, thus necessitating the establishment of a reliable signal recognition and identification procedure.
In the present study, two mutually complementary techniques, e.g.non-linear Neuro-Fuzzy and non-stationary stochastic models, were utilized for the accurate representation of the dynamical structure of the recorded Earth's electric field.The two model classes presented, offered accurate representation of the signal's special characteristics in the lower and the higher frequency bands, respectively.Both of them are characterized by high representation accuracy, with the Neuro-Fuzzy model classifying the observed anomaly as an extrinsic additional component to the Earth's natural electric field recordings.In addition, TARMA models provided very accurate signal representations indicating a significant change in the signal's frequency content.
The complementary findings attained by both methodologies, might contribute to the future development of a more accurate and generalized framework for the efficient recognition and characterization of potential EEP's including the aforementioned special characteristics and features, especially when they will not be distinct upon the Earth's recorded electric field.
a) Non-linear Hybrid Adaptive Neuro-Fuzzy Inference Systems, and (b) Non-Stationary Stochastic systems based upon Functional-Series (FS) Time-dependent (T) AutoRegressive (AR) Moving-Average (MA) (FS-TARMA) models.The aim of this paper is twofold: (a) To access the possibility of efficiently representing the Earth's recorded electric field data by both system classes, and (b) To efficiently detect and identify a possible electric field anomaly on the Earth's recorded electric field.

Fig. 1 .Fig. 2 .
Fig. 1.Recorded electric and magnetic field components from the 7 June 2006, till the 21 June 2006, f s =5 Hz: (a) Electric field component E x ; (b) electric field component E y ; (c) magnetic field component H x ; and (d) magnetic field component H y .

Fig. 4 .
Fig. 4. Recorded electric and magnetic field components from approximately 9 p.m. on the 29 December 2005, till 1 p.m. on the 11 January 2006, f s =0.0039 Hz, T s =256 s: (a) Electric field component E x ; (b) electric field component E y ; (c) magnetic field component H x ; and (d) magnetic field component H y (the squares indicate the occurrence time of the main seismic event).

Fig. 5 .
Fig. 5. (a) Main seismic event ( ) and aftershocks ( ) sequence accompanying the Kythira earthquake; (b) the observed transient electric anomaly in parallel with the seismic sequence.

Fig. 7 .
Fig. 7. (a) Recorded Earth's electric field (the square indicates the occurrence time of the main seismic event); (b) neuro-fuzzy model output indicating rejection of the possible EEP signal as an external addition upon the electric field recordings; (c) neuro-fuzzy model-based prediction error (the error signal highlights the close proximity of the neuro-fuzzy model's output signal to the recorded electric field before and after the occurrence of the possible recorded EEP signal as well as the continuous incremental rejection of the latter at the time of its occurrence).

Fig. 8 .
Fig. 8. (a) Recorded electric field component E x ; (b) first-order differenced E x ; (c) recorded electric field component E y ; (d) first-order differenced E yf s =0.0278 Hz, T s =36 s, the squares " " in subplots (a, c) and the vertical dashed lines in subplots (b, d) indicate the occurrence time of the main seismic event, the dashed frames is subplots (a, c) indicate the corresponding segments of E x , E y signals that are subsequently differenced and then shown in subplots (b, d), respectively.
2, . . ., N), y[t] the stochastic non-stationary signal being modelled, e[t] an innovations (uncorrelated) sequence with zero mean and variance σ 2 e and N the sample size of y[t].a i [t], c i [t] designate the model's timedependent AR and MA parameters, respectively (notice that a 0 [t]≡c 0 [t]≡1), which belong to a functional space of dimensionality p of the form: