Application of the Huang-Hilbert transform and natural time to the analysis of Seismic Electric Signal activities

TL;DR

This paper combines the Huang-Hilbert transform and natural time analysis to improve the identification of Seismic Electric Signals by decomposing signals into intrinsic modes and distinguishing them from noise.

Contribution

It introduces a novel approach integrating Huang-Hilbert transform with natural time analysis for better SES characterization and noise discrimination.

Findings

IMFs' instantaneous amplitudes help distinguish SES from noise

Natural time analysis enhances SES detection accuracy

Huang-Hilbert transform effectively decomposes SES signals

Abstract

The Huang-Hilbert transform is applied to Seismic Electric Signal (SES) activities in order to decompose them into a number of Intrinsic Mode Functions (IMFs) and study which of these functions better represent the SES. The results are compared to those obtained from the analysis in a new time domain termed natural time after having subtracted the magnetotelluric background from the original signal. It is shown that the instantaneous amplitudes of the IMFs can be used for the distinction of SES from artificial noises when combined with the natural time analysis.

Tables

Table 1. Table 1: Analysis in natural time of IMF sums.

Sums of IMF:	${\kappa }_1$	$S$	$S_{-}$
1	$0.077\pm 0.008$	$0.087\pm 0.012$	$0.086\pm 0.018$
1-3	$0.074\pm 0.003$	$0.073\pm 0.005$	$0.097\pm 0.003$
1-5	$0.081\pm 0.004$	$0.104\pm 0.006$	$0.077\pm 0.005$
1-9	$0.060\pm 0.012$	$0.075\pm 0.014$	$0.052\pm 0.017$
2-4	$0.080\pm 0.003$	$0.080\pm 0.005$	$0.098\pm 0.002$
2-9	$0.062\pm 0.014$	$0.078\pm 0.012$	$0.056\pm 0.016$
3-9	$0.056\pm 0.006$	$0.069\pm 0.004$	$0.050\pm 0.010$

Table 2. Table 2: Analysis in natural time of the instantaneous amplitudes of IMF sums.

Sums of IMF:	${\kappa }_1$	$S$	$S_{-}$
1	$0.078\pm 0.005$	$0.064\pm 0.008$	$0.064\pm 0.009$
1-3	$0.067\pm 0.002$	$0.068\pm 0.003$	$0.071\pm 0.004$
1-9	$0.062\pm 0.006$	$0.083\pm 0.008$	$0.060\pm 0.006$
2-4	$0.0601\pm 0.0009$	$0.0400\pm 0.0008$	$0.0418\pm 0.0005$
2-9	$0.067\pm 0.007$	$0.081\pm 0.006$	$0.063\pm 0.007$
3-9	$0.082\pm 0.005$	$0.083\pm 0.005$	$0.083\pm 0.007$

Table 3. Table 3: Summary of the results of the analysis in natural time.

Analysis in natural time of:	${\kappa }_1$	$S$	$S_{-}$
Signal after eliminating the background (dichotomous)	$0.062\pm 0.006$	$0.080\pm 0.007$	$0.057\pm 0.007$
Sum of IMF2-9 (dichotomous)	$0.062\pm 0.014$	$0.078\pm 0.012$	$0.056\pm 0.016$
Signal after eliminating the background (using $P$)	$0.0626\pm 0.0002$	$0.0789\pm 0.0002$	$0.0593\pm 0.0003$
A of sum of IMF1-9 (using $A^2$)	$0.062\pm 0.006$	$0.083\pm 0.008$	$0.060\pm 0.006$

Equations

$$\Phi (\omega )=\sum _{k=1}^N{p}_k\exp (i\omega \frac{k}{N})$$

$$\Pi (\omega )=1-0.07{\omega }^2+\dots$$

$${\kappa }_1=<{\chi }^2>-<\chi {>}^2=0.07$$

$$S=<\chi \ln \chi >-<\chi >\ln <\chi >$$

$$S,S_{-}<S_u$$

$$H[x(t)]\equiv \hat{y}(t)=\frac{1}{\pi }PV{\int }_{-\infty }^{+\infty }\frac{x(\tau )}{t-\tau }d\tau$$

$$z(t)=x(t)+i\hat{y}(t)=A(t)e^{i\theta (t)}$$

$$A(t)=(x^2+{\hat{y}}^2{)}^{\frac{1}{2}}$$

$$X(t)=\sum _{i=1}^n{c}_i+trend$$

$$x_4=\frac{x_3{e}^{(-\frac{1}{20})}+x_2{e}^{(-\frac{2}{20})}+x_1{e}^{(-\frac{3}{20})}}{e^{(-\frac{1}{20})}+e^{(-\frac{2}{20})}+e^{(-\frac{3}{20})}}$$

$$data=\frac{data-min(data)}{max(data)-min(data)}$$

$$P={A(t)}^2-{\hat{y}}^2$$

Full text

Keywords: seismic electric signal, natural time, Huang-Hilbert transform, instantaneous amplitude

Seismic Electric Signal (SES) activities are low frequency electric signals precursory to earthquakes that have been found to exhibit critical dynamics. Their distinction from noise which is usually characterized by different dynam- ics is an important task. Here, the Huang decomposition method and the Hilbert transform, used for the analysis of non-stationary time series, have been applied for the first time to time-series of SES activities of short duration in order to study their features. By decomposing the SES activity using the Huang method into some functions, called Intrinsic Mode Functions (IMFs), it is found that the first of the IMFs, hereafter IMF1, does not contain enough information to secure a classification of the signal (i.e., SES or noise), but it might contain some information about the distance of the source emitting the electric sig- nal, from the recording station. Furthermore, the Hilbert transform is applied to the IMFs in order to determine their instantaneous amplitudes $A(t)$. Earlier studies revealed that the classification of the signal as an SES can be achieved upon employing natural time analysis and the instantaneous power $P$ of the time-series as described in Ref. Varotsos et al., 2009, provided that the usual sinusoidal background has been first determined and then eliminated. Here however, it is shown that the squared instantaneous amplitude $A(t)^{2}$ of the sum of all the IMFs (i.e., the $A^{2}$ of the original time-series, including the background) can be used directly for the analysis in natural time in order to classify the signal as an SES.

I Introduction

Seismic Electric Signals (SES)Varotsos and Alexopoulos (1984a, b); Varotsos (2005) are low frequency $(\leq 1Hz)$ changes of the electric field of the earth, that have been found to precede earthquakes with lead times ranging from several hours to a few monthsUyeda (1996); Varotsos and Lazaridou (1991). They are emitted when the gradually increasing stress before an earthquake reaches a critical value Varotsos et al. (1993) in which the electric dipoles formed due to point defectsVarotsos and Miliotis (1974); Lazaridou et al. (1985); Kostopoulos et al. (1975) in the future focal area that affect the usual dielectric propertiesVarotsos and Alexopoulos (1982) exhibit cooperative orientation, thus leading to an emission of a transient electric signal. Upon analysing these signals in a new time domain termed natural timeVarotsos et al. (2001, 2002), they can be distinguishedVarotsos et al. (2003a, b) from other electric signals of different origin, e.g., signals emitted from man-made electrical sources. The latter are hereafter called “artificial” signals.

In general, the analysis in natural time may uncover properties hidden in complex time-seriesAbe et al. (2005). In particular, for the case of SES, the natural time analysis of the small earthquakes subsequent to a series of SES, termed SES activityVarotsos and Lazaridou (1991), may identify the occurrence time of the impending mainshockSarlis et al. (2008).

For a time-series consisting of N events, the natural time can be definedVarotsos et al. (2011a) as $\chi_{k}=k/N$ and it serves as an index for the occurrence of the k-th event. We studyVarotsos et al. (2002, 2006, 2005); Varotsos et al. (2003a); Varotsos et al. (2004); Varotsos et al. (2011b) the evolution of the pair ($\chi_{k},Q_{k}$) where $Q_{k}$ is a quantity proportional to the energy released in the k-th event. For a signal of dichotomous nature, $Q_{k}$ is proportional to the duration of the k-th event; otherwise, the quantity $Q_{k}$ is determinedVarotsos et al. (2009) by means of the instantaneous power $P$.

The normalized power spectrum is introduced asVarotsos et al. (2001) $\Pi(\omega)=|\Phi(\omega)|^{2}$ where

[TABLE]

with $p_{k}=Q_{k}/\sum_{n=1}^{N}Q_{n}$. As $\omega\rightarrow 0$ Eq. (1) gives some useful statistical properties. For example, it has been shownVarotsos et al. (2001) that

[TABLE]

which for an SES activity leads toVarotsos et al. (2001) the variance of $\chi$ hereafter labelled $\kappa_{1}$:

[TABLE]

Other parameters we take under consideration are the entropy S in natural time given byVarotsos et al. (2003b)

[TABLE]

where $<f(\chi)>\equiv\sum_{k=1}^{N}p_{k}f(\chi_{k})$, and the entropy under time reversal $S_{-}\equiv\hat{T}S$, where the effect of the time-reversal operator $\hat{T}$ on $Q_{k}$ is given byVarotsos et al. (2011a) $\hat{T}Q_{k}=Q_{N-k+1}$ which positions the first pulse (k=1) as last in the new time reversed time-series etc. It has been shownVarotsos et al. (2003b) that an SES activity obeys the inequalities:

[TABLE]

where $S_{u}$ the entropy of the uniform distribution $S_{u}=0.0966$. The relations (3) and (5), which involve the parameters $\kappa_{1}$, $S$ and $S_{-}$, are used in order to classify as SES an electric signal of short duration after eliminating its sinusoidal background, termed magnetotelluricVarotsos et al. (2003c); Sarlis and Varotsos (2002); Varotsos (2005) background. This background comprises electrical variations induced by frequent small variations of the Earth’s magnetic field originated from extraterrestrial sources. An example of an electric signal analysis following the procedure described above will be given below in Section III.

The Huang-Hilbert transformHuang et al. (1998, 2003); Yang et al. (2004); Bowman and Lees (2013); Xie and Wang (2006) shortly reviewed in Section II is a novel method used in the analysis of non-stationary time-series by computing instantaneous amplitudes and frequencies. One of its components is the Empirical Mode Decomposition (EMD) which decomposes a signal into some functions called Intrinsic Mode Functions (IMF) that, in combination with the Hilbert transform, can give information about the spectrum of the signal. It is the scope of the present paper to study the application of the Huang-Hilbert transformHuang et al. (1998, 2003); Yang et al. (2004); Bowman and Lees (2013); Xie and Wang (2006) to SES activities of short duration by considering also natural time analysis. In particular, the latter analysis is used for the classification of an SES activity in Section III, while the IMFs of the SES activity are extracted in Section IV. The interconnection of these results is investigated in Section V. Specifically, in order to determine which of the IMFs represent the SES, we show in Section V, that the squared instantaneous amplitudes as derived after the application of the Hilbert transform can also be used instead of the instantaneous power for the classification of the electric signal as SES. The reason we choose to study the instantaneous amplitudes is because, as shown in Section V (in particular see Fig. 7c that will be discussed later), the instantaneous amplitude of the sum of all the IMFs extracted from the signal can approximate the original time-series.

II The Huang-Hilbert Transform

The Hilbert transform of a function $x(t)$ isHuang et al. (1998):

[TABLE]

where $PV$ indicates the Cauchy principal value. The original time-series and its Hilbert transform are the real and imaginary part respectively of an analytic function $z(t)$

[TABLE]

where

[TABLE]

stand for the instantaneous amplitudes.

An IMF should satisfy two conditionsHuang et al. (1998): First, the number of extrema and the number of zero crossings must be equal or differ at most by one. Second, the mean value of the envelopes defined by the local maxima and minima is zero. If we apply the Hilbert transform to the IMFs we derive the instantaneous amplitudes $A(t)$ which will be used later for the analysis of the signal in natural time.

In order to begin the EMD (the process which extracts the IMFs) the signal has to have at least two extremaHuang et al. (1998). If there are not any extrema then the signal can be differentiated to reveal them. Finally, the characteristic time scale is defined by the time lapse between the extrema and not the zero crossings. The decomposition, i.e. the method to extract the IMFs, is designated as the sifting process and is described belowHuang et al. (2003); Yang et al. (2004); Bowman and Lees (2013); Xie and Wang (2006):

a) The maxima are determined and connected with a cubic spline, thus creating an upper envelope.

b) The minima are determined and connected with a cubic spline so that a lower envelope is created. These two envelopes should contain all the data.

c) The mean value of the envelopes m=\big{(}(upper-envelope)+(lower-envelope)\big{)}/2 is calculated and subtracted from the signal $X(t)$ deriving a component $h$, $X(t)-m=h(t)$.

d) The sifting process starts anew treating the component $h$ as the data and a new component $h^{\prime}$ is derived etc.

This process will stop when one of the $h$ components satisfies the definition of an IMF (zero mean value, the number of extrema and zero crossings is equal or differ at most by one). The first IMF will be the last $h$ component extracted. We set $h=c$.

e) The IMF c is subtracted from the original signal and the process begins anew (steps a-d) treating the residual $res=X(t)-c$ as the data until the second IMF is extracred etc.

The original signal can be recomposed using the following relation

[TABLE]

where trend stands for the last IMF. When the final residual is a monotonic function the process stops.

III Natural time analysis of an SES activity of short duration

Figure 1a depicts the electric signal that is used as a typical example in the present analysis. It was recorded at a measuring station located close to Volos city, Greece, on November 21, 2012 and preceded the $M_{b}$=4.6 earthquake with an epicentre at 38.26oN, 22.22oE that occurred on April 28, 2013.

In order to eliminate the magnetotelluric background we work as follows:

First, we plot the histogram of the signal along with its probability density function (pdf) see Fig. 2. We need to determine a level which will separate the high-level states (i.e., those having the largest deflections of the electric field amplitude) from the low-level states designated as the backgroundVarotsos et al. (2003a). For this purpose, the minimum value of the pdf, i.e., 22mV in this case, (Fig. 2) is selected as background level. Thus, two time series of subsequently high-level and low-level states are constructed. The values greater than 22mV (the high-level states) will represent the signal while the values below this level (i.e. the low-level states) correspond to the magnetotelluric background.

The time-series of the background appears to have some discontinuities (“missing points” in Fig. 3). We re-calculate these “missing points” following the process below:

The mean duration of an SES pulse lies in the rangeVarotsos et al. (2011a) 11-14sec. In every “missing point” of Fig. 3 a new value is assigned. We assume that this value is a past average with weight $e^{(-\frac{\Delta t}{20})}$. Thus, we employ the following procedure in order to estimate how many of these “past values” should be used: we apply an exponential fit moving from the “present” towards the “past”. For example, if the first missing point is the $x_{4}$ then after the fitting it will take the value

[TABLE]

Applying the exponential fit to the entire signal of Fig. 1a we plot in Fig. 1b the signal alongside the resulting background, while in Fig. 1c we plot the final signal after subtracting the resulting background from the original time-series.

Next, we compute the instantaneous power P and use it for the analysis of the electric signal in natural time by following the procedure described in Ref.Varotsos et al. (2009). We then find $\kappa_{1}=0.0626\pm 0.0002$, $S=0.0789\pm 0.0002$ and $S_{-}=0.0593\pm 0.0003$.

An alternative way to analyse a short duration electric signal in natural time and evaluate the parameters $\kappa_{1}$, $S$ and $S_{-}$, is by approximating the signal under consideration with another one which is of dichotomous nature. Along this line we do the following: First, we rescale the signal of Fig. 1c so as its values are in the region [0,1] according to the relation:

[TABLE]

We then determine a threshold $V_{thres}=0.5$ for the vertical axis designating the measured values of the potential difference $V$ and, starting from the beginning of the signal, we compare each value of $V$ with $V_{thres}$. All the points for which $V$ exceeds $V_{thres}$ correspond to the high-level states and the value 1 is assigned; to all the other points we assign the value 0. Thus, we have a new, dichotomous signal consisting of subsequent aces and zeros. The duration $Q_{k}$ of each pulse will now be the sum of the consecutive aces. The dichotomous signal along with the original signal are depicted in Fig. 4. The analysis in natural time gives $\kappa_{1}=0.062\pm 0.006$, $S=0.080\pm 0.007$ and $S_{-}=0.057\pm 0.007$. These values, after considering that the minimum $\kappa_{1}$ value measuredVarotsos et al. (2011a) to date for an SES activity is $\kappa_{1}=0.063\pm 0.003$, indicate that they obey Eqs. (3) and (5), thus this signal can be classified as an SES activity. This conclusion is consistent with the one obtained above after using the instantaneous power $P$.

IV Application of the Huang-Hilbert transform to an ses of short duration

In Fig. 5 we depict the IMFs that were extracted by applying the procedure described in Section II to the signal depicted in Fig. 1a.

The first IMF is actually the oscillation with the highest frequencies and we want to determine if it contains information about the signal. Thus, we subtract it from the original time-series (Fig. 1a) and analyse the difference in natural time after making this new signal dichotomous using the method described in Section III. We find $\kappa_{1}=0.061\pm 0.006$, $S=0.079\pm 0.005$ and $S_{-}=0.056\pm 0.009$ while during the analysis in natural time of the same signal after eliminating the magnetotelluric background using the exponential fit, we found as mentioned in Section III $\kappa_{1}=0.062\pm 0.006$, $S=0.080\pm 0.007$ and $S_{-}=0.057\pm 0.007$. Thus, we conclude that by subtracting the first IMF we can still classify the signal under consideration as SES, thus probably indicating that the first IMF does not contribute to the signal.

In order to examine which of the IMFs better represent the SES, we analyse in natural time various sums of IMFs using the corresponding dichotomous signals. The results are shown in Table  1. Comparing the results in Table  1 with those of Section III we can see that the better representation of the SES comes after summing IMFs 2-9, i.e. by subtracting only IMF1 from the original signal.

In Fig. 6 we depict the original time-series (blue dotted line) alongside the IMF sums (red line). Their comparison shows that complete concurrence exists only in Fig. 6c where all the IMFs are summed up. The original time-series already had a magnetotelluric background so the sum of all the IMFs, which in reality represents the original time-series, contains a background as well. The two plots better coincide in Fig. 6e regarding the pulses of the signal. In this case the sifting process produces some erroneous estimations in the values of the extrema with that near the 750sec mark being the most characteristic. This problem might be due to the use of the cubic spline used for the extraction of the IMFs as described in Section II. The cubic spline can add large swings near the ends of the signal which can spread inwards and corrupt the data especially the low-frequency components. This problem though, still allows for the sifting process to extract the essential scales from the dataHuang et al. (1998) and does not affect the basic characteristics of the signal.

V Analysis of the instantaneous amplitudes of the Hilbert transform

The squared instantaneous amplitude ${A(t)}^{2}$ of the sum of all the IMFs can be considered proportional to the instantaneous power $P$ of the original time-series, i.e. $P=|x(t)|^{2}={x(t)}^{2}$. From Eq. (8), we have $x^{2}={A(t)}^{2}-{\hat{y}}^{2}$. Thus, we derive the following interrelation between $P$ and $A(t)$:

[TABLE]

Our scope is to investigate whether the instantaneous amplitudes can be used instead of the instantaneous power for the analysis in natural time. The process is the following: We sum the IMFs and by applying the Hilbert transform we find the instantaneous amplitude of each sum. We then square the instantaneous amplitude and begin the analysis in natural time by applying the method of the instantaneous power as described in Ref.Varotsos et al. (2009), treating the squared instantaneous amplitude as the instantaneous power. We test various IMF sums in order to examine which of the amplitudes better represent the SES. The results are shown in Table 2 while in Fig. 7 we depict the original time-series alongside the instantaneous amplitudes of the IMF sums.

An inspection of Fig. 7 reveals that the instantaneous amplitude of the IMF1-9 sum, i.e., the amplitude of the Hilbert transform of the original time series, better approximates the original signal.

We remind that in Section III for the same signal using the instantaneous power method we had found $\kappa_{1}=0.0626\pm 0.0002$, $S=0.0789\pm 0.0002$ and $S_{-}=0.0593\pm 0.0003$. Comparing these results with the ones of Table 2 we find that we have better results by adding to the computations the highest-frequency component (i.e. IMF1). Thus, the $A^{2}$ of the original time-series (i.e. the $A^{2}$ of the sum of all the IMFs) can be used for the analysis in natural time instead of the instantaneous power.

VI Summary and concluding remarks

Table 3 summarizes the results of the analysis in natural time after the application of the Huang-Hilbert transform to an SES activity and compares them to those obtained from the natural time analysis after eliminating the background from the original time-series. First, we showed that IMF1 does not affect the durations of the pulses of the signal since, by subtracting it from the original time-series, the natural time parameters remain the same. This might indicate that IMF1 contains information about the distance of the source emitting the electric signal, from the recording station. Second, we showed that the $A^{2}$ of the sum of all the IMFs (i.e. the $A^{2}$ of the original time-series) can be used instead of the instantaneous power of the time-series (after having subtracted the background) for the analysis in natural time in order to classify the signal as an SES.