Empirical Wavelet Transform-based Detection of Anomalies in ULF Geomagnetic Signals Associated to Seismic Events with a Fuzzy Logic-based System for Automatic Diagnosis

Abstract

Owing to the relevance and severity of damages caused by earthquakes (EQs), the development and application of new methods for seismic activity detection that offer an efficient and reliable diagnosis in terms of processing and performance are still demanding tasks. In this work, the application of the Empirical Wavelet Transform (EWT) for seismic detection in ultra-low-frequency (ULF) geomagnetic signals is presented. For this, several ULF signals associated to seismic activities and random calm periods are analysed. These signals have been obtained through a tri-axial fluxgate magnetometer at the Juriquilla station localized in Queretaro, Mexico, longitude -100.45° N and latitude 20.70°E. In order to show the advantages of the proposal, a comparison with the discrete wavelet transform (DWT) is presented. The results shown a better detection capability of seismic signals before, during, and after the main shock than the ones obtained by the DWT, which makes the proposal a more suitable and reliable tool for this task. Finally, a fuzzy logic (FL)-based system for automatic diagnosis using the variance of the EWT outputs for the tri-axial fluxgate magnetometer signals is also proposed.

1. Introduction

Among the natural disasters, the EQs have attracted the interest of many researchers around the world due to the huge amount of human, material, and economic losses [1]. They have been focused on finding pre-seismic precursors for prediction and forecasting tasks [2-6]. In this regard, different electromagnetic phenomena (EP) encompassing a large frequency range, being the ultra-low-frequency (ULF) range one of the most promising, have been associated with EQs because they typically occur during, but sometimes prior to seismic activity [7-14] Although many detection methods have been proposed, the relevance and severity of EQs damages demand still more efficient and reliable diagnosis methods.

Techniques and methods such as polarization or spectral density ratio analysis [15-16], transfer function analysis [17], fractal analysis [18-22], singular value decomposition [23], principal component analysis [16, 21], and the discrete wavelet transform (DWT) [2, 24], among others have been proposed to analyse the ULF geomagnetic signals associated to EQs. Yet, despite showing promising results, the inherent characteristics of the ULF geomagnetic signals such as high noise levels, weak amplitude, and interferences from other sources due to the distance between the epicenter and sensor, among others may compromise the performance and reliability of the analysis. From this point of view, the development and application of new detection methods to make a more efficient and reliable diagnosis in terms of processing and performance are still interesting research fields.

In this chapter, the application of the empirical wavelet transform (EWT) to ULF signals for detecting seismic precursor anomalies is presented. Besides, a comparison with the DWT is carried out in order to show the EWT advantages. Moreover, a fuzzy logic (FL) system for automatic diagnosis using the variance of the EWT results is proposed. For this, three ULF signals associated to seismic activities and random calm periods are analysed. The obtained results show the usefulness and effectiveness of the proposed methodology, making it a suitable and reliable tool to detect ULF anomalies.

2. ULF geomagnetic data

In order to investigate the relationships between ULF geomagnetic signals and pre-seismic anomalies, ULF geomagnetic data from Juriquilla seismic station, located in Queretaro, Mexico, with geographic coordinates: longitude -100.45° N and latitude 20.70° E, are used. The ULF geomagnetic signals are monitored by means of a fluxgate magnetometer. It allows monitoring three mutually orthogonal components of the magnetic field, two horizontal (Mx: North-South and My: East-West) components and a vertical component (Mz). The three geomagnetic components are measured using a sampling frequency of 1 Hz to obtain 65,000 samples during a time window of 18 hrs, which comprise 9 hours before the main shock and 9 hrs after it. In this research, three recent seismic events with magnitude greater than 6.0 are analysed. Further, for comparison purposes, random analyses during periods of seismic calm are used. Table 1 summarize the characteristics of the studied EQs.

To discriminate the geomagnetic activity of the magnetosphere because of the solar activity and cultural noise, the analysed EQs data are compared with the geomagnetic activity expressed by Dst index (http://wdc.kugi.kyoto-u.ac.jp/dstdir/), where those indices apparently had no correlation with EQs variation.

Event	Year	Month	Day	Hour	Min	Longitude	Latitude	Magnitude, M	Depth, km	Distance, km	ρ, km	Distance/ρ
1	2009	8	3	13	0	-112.24	28.48	6.9	10	1473	2292	0.64
2	2009	9	24	2	16	-107.43	17.72	6.2	21	8052	1109	0.73
3	2010	6	30	2	22	-98.03	16.22	6.0	8	563	684	0.82

Table 1.

Characteristics of the Earthquakes occurred in Mexico during 2009–2010. Their magnitudes are presented in bold (Catalogue of National Seismological Service, Mexico)

Note: The Year / month / day / hour / minute are the exact time of the EQ (Local Time); Latitude and Longitude are the geographic coordinates of the epicenter, Magnitude and Depth are the EQ measures, Distance is the distance between the epicenter and Juriquilla station, and p is the radius of the EQ preparation.

3. Wavelet transform

This section presents the theoretical background of the Discrete Wavelet transform and the Empirical Wavelet Transform used for the analysis of ULF signals.

3.1. Discrete wavelet transform

Discrete Wavelet Transform (DWT) is a useful method for analysis of non-stationary, no linear and transient signals because it decomposes the time series signal into multiple time-frequency levels retaining the characteristics of the analysed signal [2]. DWT is defined by Eq. (1), where x(n) and h(n) denote the discrete signal and the wavelet basis function, respectively, of the total number N of samples contained in the signal x(n). j and k represent the time scaling, and the shifting of the discrete wavelet function, respectively.

$DW{T}_{j,k}={\displaystyle \sum }_{N}x(n)\overline{h_{j,k}(n)}$E1

The DWT is calculated using a set of low- and high-pass filters bank called approximations (AC_L) and details (DC_L) into desired levels L (Mallat algorithm), respectively, as shown in Figure 1. Based on the Mallat algorithm, the approximation obtained from the first level is split into a new decomposition and approximation and this process is repeated [25]. Once the discrete time signal x(n) has been decomposed into the desired levels, the signal is reconstructed by applying the decomposition process in an inverse way, which is known as the inverse discrete wavelet transform (IDWT).

Figure 1.

DWT basis construction.

According to the Mallat algorithm, the frequency band for the approximations $A{C}_L$ and decompositions $D{C}_L$ DC are given by Eq. (2) and Eq. (3), respectively, where $F_s$ represents the sampling frequency of the signal.

$A{C}_L\Rightarrow \left[0,\frac{F_s}{2^{L+1}}\right]$E2

$D{C}_L\Rightarrow \left[\frac{F_s}{2^{L+1}},\frac{F_s}{2^L}\right]$E3

Different types of wavelet mother function have been proposed to analyse ULF signals in order to find anomalies related with EQs such as Daubechies, Haar, Morlet, Symlets, Coiflets, and Meyer (Figure 2). However, it has been demonstrated that the most effective to analyse ULF signals is the Daubechies mother function [2, 24, 26]. For this reason, Daubechies as mother wavelet is used in this work.

Figure 2.

Wavelets used in ULF signals.

3.2. Empirical wavelet transform

EWT is a new adaptive wavelet transform capable of decomposing a time series signal x(t) into adaptive time-frequency sub-bands according to its contained frequency information [27]. This advantage allows generating narrow time-frequency sub-bands, unlike the DWT where the calculated time-frequency sub-bands depend on sampling frequency of the time signal. To provide an adaptive wavelet with respect to the analysed signal, the segmentation of the signal is an important step in the EWT. It can be performed either manually or by means of the Fourier spectrum. If the signal is segmented manually, the boundaries of the wavelet filters can be user selected as contiguous segments. On the other hand, using the Fourier spectrum, first, the local maxima of the Fourier spectrum $x(\omega )$ x are calculated. Next, the boundaries ${\omega }_i$ ω of each segment are defined as the center between two consecutive maxima. Thus, the Fourier support [0, π] is segmented N into contiguous sets or frequency bands. In both segmentations, each frequency band is indicated by ${\text{Λ}}_{\text{n}}=\left[{\text{ω}}_{\text{n}-1},{\text{ω}}_{\text{n}}\right]$ Λn=[ω and satisfy, as shown in Figure. 3. A transition phase of width $2{\tau }_n$ 2τ is defined to obtain a tight frame around each ${\text{ω}}_{\text{n}}$ ω. A more detailed selection of ${\tau }_n$ τ is presented in [27]. Observing Figure. 3, the empirical wavelets are defined by one low-pass represented by LPF ${\varphi }_n(\omega )$ and $N-1$ N band-pass ${\psi }_n(\omega )$ filters represented by BPF corresponding to the approximation and details components, respectively, on each ${\text{Λ}}_{\text{n}}$ Λn.

Figure 3.

EWT basis construction.

Following the idea used in deriving the Meyer’s wavelet, [27] defines the empirical scaling function to estimate the low-pass wavelet filter coefficients according to following Equation:

${\varphi }_n\left(\omega \right)=\{\begin{array}{c}1 if\left|\omega \right|\le {\omega }_n-{\tau }_n\\ cos\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_n}\left(\left|\omega \right|-{\omega }_n+{\tau }_n\right)\right)\right] if {\omega }_n-{\tau }_n\le \left|\omega \right|\le {\omega }_n+{\tau }_n\\ 0 otherwise\end{array}$E4

And an empirical wavelet function to build the N-1 band-pass filters as:

${\psi }_n\left(\omega \right)=\{\begin{array}{c}1 if {\omega }_n+{\tau }_n\le \left|\omega \right|\le {\omega }_{n+1}-{\tau }_{n+1}\\ cos\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_{n+1}}\left(\left|\omega \right|-{\omega }_{n+1}+{\tau }_{n+1}\right)\right)if{\omega }_{n+1}-{\tau }_{n+1}\le \left|\omega \right|\le {\omega }_{n+1}+{\tau }_{n+1}\right]\\ sin\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_n}\left(\left|\omega \right|-{\omega }_n+{\tau }_n\right)\right)if{\omega }_n+{\tau }_n\le \left|\omega \right|\le {\omega }_n+{\tau }_n\right]\\ 0 otherwise\end{array}$E5

where β(x) is a polynomial function taken as $\beta \left(x\right)=x^4(35-85x+70x^2-20x^3)$

After having built the wavelet filters, the signal x(t) is decomposed into different frequency bands through empirical wavelet transform defined by

$W_f^{ϵ}\left(n,t\right)=F^{-1}\left(x(\omega ){\psi }_n\left(\omega \right)\right)$E6

$W_f^{ϵ}\left(0,t\right)=F^{-1}\left(x(\omega ){\varphi }_n\left(\omega \right)\right)$E7

where the details $W_f^{ϵ}(n,t)$ and approximation $W_f^{ϵ}(0,t)$ coefficients are obtained by the inner products of the signal with the empirical wavelets low-pass and band-pass filters, respectively, andF^-1is the inverse Fourier transform.

4. EWT and fuzzy logic system

This section presents the proposed methodology. It consists of the ULF geomagnetic signals analysis through the EWT, then a statistical parameter based on the variance is applied and, finally, an automatic diagnosis by means of a FL system is computed as shown in Figure. 4.

Figure 4.

Proposed methodology.

4.1. EWT analysis

Generally, the pre-seismic anomalies are too much weak to be detected by the Fourier transform (Chavez et al., 2010); hence, the frequency bands are selected manually using the EWT. Several experimental using both algorithms are carried out to estimate the best frequency band of the ULF geomagnetic signal to detect anomalies associated to the EQs. After the experimental runs, it is found that for the EWT algorithm the frequency band from 0.0470 to 0.0781 Hz and for the DWT algorithm the frequency band or third decomposition from 0.0625 to 0.125 Hz with a Daubechies wavelet of order 5 generate the best results, enhancing correlation with associated seismic anomalies events. Figure. 5 presents the obtained results for the EWT and the DWT, where both the seismic calm signal (left-side plots) and the seismic activity signal (right-side plots), which corresponds to event 1 (Table 1), can be observed. Figure.5(a) shows that in both analysis, the seismic calm period do not present significant spikes over time, indicating the absence of seismic activity. On the other hand, observing the results shown in Figure. 5(b), both time-frequency analysis can detect the occurrence of peaks prior to seismic (Pre-seismic event zone) and another peaks after the main shock (Post-seismic event zone). These magnetic perturbations occur about 8hrs before the main shock and about 2 hrs after it. Open circles remark perturbations in the signal.. But, it is noticeable that, using the EWT method, it is better noticeable of the pre-seismic and post-seismic anomalies.

Figure 5.

Comparison between the EWT time-frequency analysis and DWT time-frequency analysis; for (a) seismic calm, and (b) with seismic activity.

In order to evaluate the significance of these results, a complementary statistical analysis based on the variance of the EWT and DWT results is computed to measure the fluctuations between seismic activity and seismic calm period as follows:

$V=\frac{1}{N}{\displaystyle \sum }_{n=1}^N\left[x\left(n\right)-\overline{x}\right]$E8

where V is the variance of each data window at the frequency band, N is the total of data analysed for the region of interest, x(n) is the input sequence, and $\overline{x}$ is the mean value of x(n). Figures. 6(a) 6(b) show the variance (V) results obtained by the EWT and the DWT method, respectively, for seismic activity and seismic calm period. The results correspond to running data windows each 1000 samples. It is observed that the EWT method presents a better performance to detect seismic anomalies than DWT method. Hence, it can be established that the EWT analysis allows the observation of ULF signal perturbations with low amplitude and embedded in high level of noise.

Figure 6.

Variance (V) of the seismic activity and seismic calm period using: (a) EWT and (b) DWT.

4.2. Study cases

After showing in the previous section that EWT improves the correlation between the seismic event and the ULF electromagnetic signal, the proposed EWT time-frequency analysis is applied to seismic calm and three seismic events with different geographical location. Figure. 7 shows the EWT time-frequency analysis for seismic calm period and the 3 seismic events. The three components of the magnetic field, Mx, My, and Mz, are analysed as shown in Figure. 7. The main shock position is indicated by an arrow and two circles. It shows the pre-seismic and post-seismic anomalies associated with the EQs. The EWT results for seismic calm period do not present significant spikes over time, indicating the absence of seismic activity, as shown in Figure. 7(a-c); unlike the signals with seismic activity where several spikes appear before and after the main shock. All the analysed signals comprise 9 hrs before and 9 hrs after each seismic event, considering the time zero as the specific time of the occurrence of the EQ.

Figure 7.

EWT analysis of the three geomagnetic components: Mx, My and Mz, for a seismic calm (a-c), and with seismic activity: (d) to (f) EQ1, (g) to (i) EQ2, and (j) to (l) EQ3.

Similar to previous section, to evaluate the significance of the results, a complementary statistical analysis based on the variance of the EWT results is computed to measure the fluctuations between seismic activity and seismic calm period. Figure. 8 shows the variance V results obtained for three analysed seismic events as well as for their three components (Mx, My, and Mz). The results correspond to running data windows each 1000 samples. According to the obtained results in Figure.8(a-c), there are important variations of V for the three components that could be associated with occurrence of the EQs. As observed, V increases before, during and, after the seismic event on the three components: Mx (a), My (b), and Mz (c). Observing Figure. 8, the Mx geomagnetic component presents an important variance in three different ranges: from 8 to 5 hrs before seismic event, between 2hr before and 2hr after the main shock, and from 3 to 7 hrs after the main EQs (post-seismic zone). The My component shows variance in different ranges, from 8 to 6 hrs before seismic event, between 2 hrs before and 2 hrs after main shock, and from 3 to 8 hrs after the main EQs. Finally, the Mz geomagnetic component also presents variance in three different ranges: from 8 to 6 hrs before seismic event, 2 hrs before and 2 hrs after main shock, and from 4 to 8 hrs after seismic event. In summary, these results show that the EWT is adequate to find electro-magnetic seismic precursors related to the variance magnitude.

Figure 8.

Variance of the EWT for the three geomagnetic components: (a) Mx, (b) My, and (c) Mz.

4.3. Fuzzy logic-based system

Once the variance of the EWT signals is computed, a FL-based system is used for automatically diagnosing the severity of the ULF geomagnetic variations associated to seismic events. A FL system represents a group of rules for reasoning under uncertainty in an imprecise or fuzzy manner. It is usually used when a mathematical model of a process does not exist or does exist but is too difficult to encode and too complex to be evaluated fast enough for real time operation. Besides, it can use several sources of information in order to take a decision according to a particular objective.

The designed and implemented FL system to perform the diagnosis process is a Mamdani-type fuzzy inference system with two inputs, one output, and 16 rules. The system uses Max–Min composition, and the centroid of area method for defuzzification. The inputs are the variance of EWT results for the signals Mx and My, the Mz signal is not considered since it presents a low difference between seismic activity and calm period. These inputs are partitioned into four trapezoidal membership function sets, as shown in Figures. 9(a) and 9(b). They are labelled as NV (normal variance), LV (low variance), MV (medium variance), and HV (high variance). The output is also divided into four trapezoidal membership functions as shown in Figure. 9(c), their labels are NF (normal fluctuations), LF (low fluctuations), MF (medium fluctuations), and HF (high fluctuations). The crisp output of the Mamdani FL system can assume values between 0.5 and 4.5, where normal variations = 1, low variations = 2, medium variations = 3, and high variations = 4. The parameters of membership functions are determined according to the interpretation of the variance results by the authors. The set of rules that classifies the inputs variance is show in Table 2; there, one rule can be read as follows if (variance Mx is NV and variance My is NV) (light gray) then the geomagnetic fluctuations magnitude is NF (dark gray), and so on.

Geomagnetic Fluctuations Magnitude		Variance (My)
NV	LV	MV	HV
Variance (Mx)	NV	NF	NF	LF	MF
LV	LF	LF	MF	HF
MV	LF	MF	HF	HF
HV	MF	HF	HF	HF

Table 2.

Rules table for the proposed FL

Figure 9.

Membership functions: (a) Input Mx variance, (b) Input My variance, and (c) Output diagnosis.

The FL output for a calm period signal is shown in Figure. 10(a), where most of the results are Normal and only a few data indicate Low ULF geomagnetic variations. On the other hand, the outputs for the three geomagnetic signals associated to EQs indicate many Medium and High variations as shown in Figure. 10(b-d); therefore, if these results are obtained in future data they could be associated to seismic activities.

Figure 10.

FL-based diagnosis for the analysed cases: (a) calm period and (b-d) seismic activities associated to ULF geomagnetic variations.

5. Conclusions

In this chapter, a new time-frequency study based on the EWT for analysing ULF geomagnetic signals, at Juriquilla station in Queretaro Mexico, is presented. In order to prove the effectiveness of the EWT algorithm; three real data of EQs are analysed. The results demonstrate that the proposal has a greater detectability than DWT for detecting anomalies before, during, and after the main shock since EWT allows selecting narrow time-frequency sub-bands, unlike the DWT where the calculated time-frequency sub-bands depend on sampling frequency of the time signal. Further, the variance, a statistical complementary analysis, shows that a seismic event can be detected from 8 to 5 hrs before it occurs. It also indicates that relevant information can be obtained from 563 to 1473 km (epicenter distance) to the testing station. Therefore, the proposed time-frequency analysis can extract the abnormal signals in the ULF range of the EP related to different stages of the EQ preparation. Finally, the proposed FL system can automatically classify the magnitude variations of the EP into Normal, Low, Medium, and High variations using both Mx and My signals, where is found that Medium and High variations could be associated to seismic activities.

In a future work, the overall methodology will be implemented into a digital signal processor (DSP) for online and continuous monitoring of ULF geomagnetic variations and possibly used for obtain a implicit correlation between the seismic magnitude and the ULF geomagnetic signals.