Spectro-temporal features applied to the automatic classification of volcanic seismic events

Accepted Manuscript Spectro-temporal features applied to the automatic classification of volcanic seismic events

Ricardo Soto, Fernando Huenupan, Pablo Meza, Millaray Curilem, Luis Franco PII: DOI: Reference:

S0377-0273(18)30030-1 doi:10.1016/j.jvolgeores.2018.04.025 VOLGEO 6368

To appear in:

Journal of Volcanology and Geothermal Research

Received date: Revised date: Accepted date:

18 January 2018 25 April 2018 25 April 2018

Please cite this article as: Ricardo Soto, Fernando Huenupan, Pablo Meza, Millaray Curilem, Luis Franco , Spectro-temporal features applied to the automatic classification of volcanic seismic events. The address for the corresponding author was captured as affiliation for all authors. Please check if appropriate. Volgeo(2018), doi:10.1016/ j.jvolgeores.2018.04.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Spectro-temporal features applied to the automatic classification of volcanic seismic events Ricardo Sotoa , Fernando Huenupana , Pablo Mezaa , Millaray Curilema , Luis Francob a Department

of Electrical Engineering, Universidad de La Frontera, Av. Francisco Salazar 01145. Temuco, Chile Vulcanol´ogico de los Andes Sur, Rudecindo Ortega 03850, Temuco, Chile

RI PT

b Observatorio

Abstract

MA

NU

SC

Feature extraction and selection are very relevant processes in the design of automatic classifiers. In the context of volcanic seismic signal classification, most of the features presented in the literature have been extracted separately from a single domain, such as time or frequency. However, spectrograms, which combine time and frequency information, are widely used by experts during the classification of manual seismic events. This paper proposes to evaluate the performance of classifiers trained with features extracted from the spectro-temporal domain, individually or combined with other conventional features. The parameters were extracted from the spectrogram, based on a curve which combines the high energy components and the frequency bandwidth information through the duration of the event. The tests were performed at the Llaima volcano and seismic events were classified into four classes: long-period, tremor, volcano-tectonic and tectonic, using a database of signals recorded between the years 2009 and 2017. The main achievements of this study were the reduction of more than 70% of the error and false positive rates and also a reduction of approximately 30% of the number of features, compared with a baseline established in previous studies. Thus, the inclusion of spectro-temporal information was considered relevant to complement the conventional features and to support classification.

PT E

D

Keywords: Volcano monitoring, Signal processing, Pattern recognition, Spectrogram, Spectro-temporal curves

1. Introduction

AC

CE

There are approximately 1,500 active volcanoes in the world producing, on average, 30 eruptions per year. Many of these massifs are close to populated areas, resulting in a permanent risk situation to the populations in their vicinity (Brown, 2017). Moreover, the strategic position of threatened cities, coupled with other consequences of eruptions, such as climatic change or water pollution, directly affect the economy of the region (Chester et al., 2000). Hence, monitoring volcanic activity, including ground deformation, volcanic gases, temperature variations, and seismic signals, are crucial in hazard assessment and in the development of contingency plans (Sparks et al., 2012). Volcano monitoring is based also on the assumption that the quantity and/or behavior of seismic volcano events change before an eruption. The growing need to monitor a high number of volcanoes has driven the development of automatic identification systems (AISs), either by volcanic events detectors or classifiers (Orozco-Alzate et al., 2012). In particular, the Preprint submitted to Journal of LATEX Templates

analysis of seismic signals based on pattern recognition (PR) techniques has been the principal method used to develop AISs. Designing a seismic events classifier is not a trivial task due to the large amount of volcanic earthquakes with similar features, hindering signal discrimination. Thus, to ensure a good performance of an AIS, two PR stages are required: feature extraction and classification. To the best of our knowledge, the most frequently used AIS classification techniques are hidden Markov model (HMM), support vector machine (SVM) and self organizing map (SOM). HMM has been widely used in volcanic seismic event identification, achieving a true positive accuracy of over 90% and an average accuracy of over 85% in classifying events, as shown in Bhatti et al. (2016) and Bicego et al. (2013), respectively. SVM has been less explored than HMM in seismic event identification. However, SVM has classification results similar to HMM, as in Masotti et al. (2006), which applied SVM to classify volcanic tremor data recorded during different states of activity at Etna volcano (Italy) achieving an accuracy over 90%. In May 3, 2018

ACCEPTED MANUSCRIPT

NU

SC

RI PT

(2017) and Unglert et al. (2016). Hibert et al. (2017) assessed the relevance of 60 features commonly exploited by human operators for volcanic event classification. The evaluation was made through the Random Forest algorithm, separating VTs from rockfall events, where two of the six best selected features were extracted from the spectrogram. In contrast, Unglert et al. (2016) developed a technique to automatically extract three relevant features to identify changes in volcanic activity by experts. These features are related to the variations in the spectral content of volcano seismicity from spectrograms. In this study, a novel approach for feature extraction is proposed, based on the spectro-temporal curves obtained from spectrograms. The spectrogram represents the evolution over time of the spectral content of a signal. It is a tool widely used by seismologists because it highlights patterns associated with the evolution of the bandwidth and energy behavior of the frequency over time. The proposed curves were designed using expert knowledge and describe the spectral content of seismic signals through the evolution of energy, weighted by a frequency bandwidth support. The proposed curves were developed using seismic signals from the Llaima volcano (between 2009 and 2017), one of the biggest volcanoes by volume in Chile and one of the most active volcanoes in South America. The seismic events considered in this study are from the classes LP, VT, tremor (TR) and tectonic (TC). The paper is structured as follows: A description of the Llaima volcano and the seismic events considered is given in section 2. Section 3 describes the proposed spectro-temporal curves and their feature extraction process. Section 4 analyzes the impact of the new features on classifier performance. Section 5 presents the experimental results and a discussion. Finally, conclusions and future research challenges are discussed.

AC

CE

PT E

D

MA

Lara-Cueva et al. (2016) and Curilem et al. (2016), two class and four class classifiers were used, respectively, achieving accuracies higher than 95%. Unlike other classification techniques, SOM is an unsupervised method widely applied in data exploratory analysis. Esposito et al. (2008) used this technique to associate physical features of the Stromboli volcano with waveform differences of VLP (Very Long period) events. Langer et al. (2009) contrasted SOM with SVM, Multilayer Perceptron (MLP) and Cluster Analysis (CA) at tremor data classification from Etna volcano. As an advantage of SOM the authors highlights the capability to identify transitional stages from a regime of volcanic activity to another one. Langer et al. (2011) developed a software based on SOM, in order to detect changes in the dynamic regime of the Etna volcano. This software was used for the monitoring of Etna volcano. Carniel et al. (2013) applied SOM in assessing low level seismic activity prior to small scale phreatic events at the Ruapehu volcano in New Zealand, distinguishing three regimes (pre-eruptive, eruptive and posteruptive) when classifying tremor spectral patterns. In the literature, most signal features are extracted from the time and frequency domains. Temporal features, used to describe signal shape and distribution, are computed from absolute signal values by applying statistical measures such as standard deviation, mean, kurtosis and skewness, to name a few (Curilem et al., 2009). On the other hand, frequency features are typically extracted from the Fourier spectrum, where characteristics related to spectral energy, frequency and spectral evolution are extracted from the spectrogram (Carniel, 2014). The wavelet transform (WT) is another tool widely used to analyze non-stationary time series as presented in Jones et al. (2012). It should be noted that the WT is commonly used to decompose signals into different frequency bands to distinguish volcanic event classes, like long-period (LP) and volcanotectonic (VT), by examining the wavelet energy levels per band (Lara-Cueva et al., 2016; Curilem et al., 2014). In general, feature extraction techniques have been inherited from other research fields where the input is a time series, as shown in Orozco-Alzate et al. (2012), where the authors described most of feature extraction techniques applied to recognition of volcano-seismic signals. In particular, the application of speech framework has been successfully proven, as presented in Beyreuther et al. (2008). In recent years, feature extraction techniques have been developed based on expert knowledge, extracting information from the visual assessment of spectrograms as was done in Hibert et al.

2. Llaima volcano seismicity The Llaima volcano is located in the Araucan´ıa Region of Chile on the western edge of the Andes (38◦ 41’ S – 71◦ 44’ W). It is part of the Andean volcanic arc in the segment of the so-called Southern Volcanic Zone (Stern, 2004), where volcanism results from the subduction between Nazca and South America plates. Llaima is a complex active stratovolcano with a mainly basalticandesite composition. Its base has an average altitude of 740 m a.s.l. and the volcanic cone rises to around 2,430 m, surpassing the surrounding summits by nearly 1,500 m. The main volcanic edifice consists of two peaks, the 2

ACCEPTED MANUSCRIPT

event is as follows: first the different stations around the volcano are analyzed to evaluate if the signal is a seismic event, discarding noise. This analysis also consists in estimate the amplitude and the arrival time of the signal; if the amplitude is greater and the arrival time (P wave) is earlier in the station(s) closer to the crater, the events is defined as volcanic, not tectonic. Second, the reference station of the volcano is used for event classification. For OVDAS the reference station is the most reliable sensor in the volcano related to location, low noise (< 0.2 μm/s) and stable data transmission. Under these criteria, the volcanic event classes VT, LP and TR are classified. VT events are earthquakes that typically have broad spectra extending up to 15 Hz. They are associated with shear fractures in the brittle rock of the volcanic edifice or in the crust beneath it, due to the pressure induced by magma, gas or sudden temperature changes. Their name is due to their close similarity to tectonic earthquakes (Chouet & Matoza, 2013). Their duration is, in general, shorter than other volcanic seismic events, often lasting less than 1 minute. In addition, the source mechanism of these earthquakes basically generates body waves that are longitudinal (P) and transversal (S) along with superficial waves (L and R). Their spectrogram is very characteristic: related to energy its begin is impulsive and shows an exponential decay towards the end of the event, being their more energetic frequency components over 5 Hz. In relation to frequency, the bandwidth (frequency range where most of the energy is concentrated) of VTs widens quickly, then it narrows progressively as shown in Fig. 2.g. According to Chouet & Matoza (2013), the excitation mechanism of LPs may be observed in different volcanic processes, such as formation of lava domes and solid extrusion processes in the craters of andesiticdacitic stratovolcanoes (Gil-Cruz & Chouet, 1997; Yamasato, 1998; Neuberg et al., 2000; Moran et al., 2008). Nevertheless, water is abundant in stratovolcanoes, where the steam exsolved from magma interact with the superficial aquifers (Branch, 1976; Wohletz & Heiken, 1992; Fournier, 1999; John et al., 2008). Also, seismicity of LPs is directly related to magmatic fragmentation, degasification and explosions (Arciniega-Ceballos et al., 2003, 2008; Molina et al., 2004; Chouet & Dawson, 2011), and may be related to a brittle failure of magma on volcanic conduits (Goto, 1999; Neuberg et al., 2006), or self-sustained oscillations into magma’s channels (Julian). A widely accepted mechanism for LP events is the resonance of fluid-filled structures (Aki et al., 1977; Chouet, 1985, 1986, 1996; Ferrazzini & Aki, 1987; Julian;

higher being the north one (3179 m a.s.l.), which is separated by 750 m from the south summit (2930 m a.s.l.). The tallest summit has an open crater about 300 m in diameter (Naranjo & Moreno, 2005).

MOT

ROC

LAV

PIC

CRU

RI PT

Llaima

38.7˚ S

AGU CON

SC

LAJ

LLA

PAI

0 km 38.9˚ S

10 km 71.8˚ W

71.7˚ W

MA

NU

38.8˚ S

71.6˚ W

D

Figure 1: Location of the Llaima volcano and its stations.

AC

CE

PT E

The Llaima volcano is one of the most active volcanoes in the southern Andes and it is considered the second most dangerous Chilean volcano (Lara et al., 2011). The last eruptive period of Llaima volcano began in May 2007 with a weak ash emission followed by two major Strombolian eruptions: January 1 st , 2008 and April 4th , 2009 (Franco et al., 2014). Currently, the volcano is monitored by the ten seismic stations shown in Fig. 1. The Observatorio Vulcanol´ogico de los Andes Sur (OVDAS) of the Servicio Nacional de Geolog´ıa y Miner´ıa (SERNAGEOMIN) is the state agency responsible for volcano monitoring in Chile. OVDAS uses the criteria defined in Latter (1981); Lahr et al. (1994); Chouet (1996); Chouet & Matoza (2013) to identify and label the seismic events recorded by the seismic stations of monitored volcanoes. An adequate installation of the sensors, in a robust instrumental network design, improves the discrimination of the volcanic seismicity from other kind of signals, such as noise or tectonic signals. Like in Rouland et al. (2009), OVDAS uses a SNR criterion to define a threshold detection for the seismic events. The procedure to detect and classify a volcanic 3

ACCEPTED MANUSCRIPT

Normalized Amplitude

1.0 0.5 0 -0.5

RI PT

15 10

SC

Frequency [Hz]

-1.0 20

5 0 10

20 30 Time [s]

20 10 0 -10

NU

0

30

MA

Figure 2: Time and spectral representations of the signals belonging to the four classes of events examined. From left to right: LP, TR, VT and TC. Their time representation is shown at top and their spectrograms at bottom. Signals were filtered between 1 and 10 Hz and spectrograms are in decibel scale, being calculated from signals measured in μm/s , windowed with overlap of 90% using a Hamming window of length 1 second and 1024 FFT samples.

accepted models to TR and LP events is the response of fluid-filled conduits to sustained pressure fluctuations (Chouet, 1988, 1992). In general terms TR is characterized by a continuous vibration that may last from minutes to years (e.g. Villarrica volcano). Its waveform could be described as an emergent increase in wave amplitude and its end is characterized by a slow attenuation. When the signal intensity decreases, the tremor amplitude is masked by the background noise, being difficult to determine the end of the event. OVDAS analysts establish the end of the tremor when the amplitude of the signal becomes similar to the background noise. The tremor in the Llaima volcano is not a permanent signal and it has become increasingly scarce in periods furthest from the last eruption. In the Llaima volcano, the tremor can be distinguished from other events by analyzing the energy variation in the frequency band 0.5–3.0 Hz and its duration characterized by continuous oscillations or sustained amplitude. Tremors have a characteristic spectrogram, similar to several LPs very close together resulting in a long duration and narrow low-frequency content, as shown in Fig. 2.f.

AC

CE

PT E

D

Neuberg et al., 2000; Jousset et al., 2003). These events generate oscillations that attenuate and produce resonance with dominant frequencies related to the geometry of the cavity formed by the volcanic duct and the physical and acoustic properties of the fluids (Chouet, 1992; McNutt, 1992; Julian). The LP spectral pattern is narrower than the spectral pattern of VT and it is bounded mainly within the range of 0.5 and 5.0 Hz in the Llaima volcano. Regarding the evolution of energy, the LP spectrogram often has a non-impulsive beginning followed by a slow decay at the end of the event, its maximum energetic components typically are below 5 Hz. LP’s bandwidth is characterized by a smooth increase and decay, as shown in Fig. 2.e. The origin of TR mechanism has been widely studied in Steinberg & Steinberg (1975); Aki et al. (1977); Chouet (1985, 1986); Crosson & Bame (1985); Gordeev (1993); Julian; Benoit & McNutt (1997); Ripepe et al. (2001). In the literature it is possible to find different source models for TR, such as related to the movement of gases along volcanic vents and self-oscillations Steinberg & Steinberg (1975), or associated to channels opening due to over-pressure fluids Aki et al. (1977); Aki & Koyanagi (1981); Chouet (1981), among others. However, one of the most widely

Other events that are not related to volcanic activity are also recorded by seismic stations. One of the most important are TC events. TC events are related to 4

ACCEPTED MANUSCRIPT

Table 1: Number of events of each class and the period in which they were recorded.

2.1. Seismic events database

Training / validation

Test

Total

Period

LP TR VT TC Total

497 324 202 940 1963

249 163 102 471 985

746 487 304 1411 2948

2011 - 2013 2009 - 2010 2011 - 2017 2011 - 2013 2009 - 2013

RI PT

Class

3.1. The spectrogram

NU

SC

The spectrogram is a visual representation of the time-dependent Fourier transform where the onedimensional sequence x(m) is converted into a twodimensional function, X(m, k), of the time variable m and the frequency variable k. X(m, k) constitutes the short-time Fourier transform of the signal where the frequency is sampled in K equally spaced 2π/K. The time dependent Fourier transform (TDFT) is defined as Oppenheim & Schafer (2009):

MA

the dynamics of the geological faults (Chouet, 1996) or subduction processes (e.g. Nazca and South American plates). Their waveform and spectral content depend on many factors, such as the distance from the epicenter and the site effect. In relation to the distance from the epicenter, a far TC has a lower frequency than a near TC. Meanwhile, site effect is related to the soil characteristics in which the station is located. The site effect imposes a spectral pattern on the signal, thereby affecting the measurement of all events. In TC, the frequencies present at the end of the event are related to this effect. Similar to VT, TC spectra are characterized by an impulsive beginning and a subsequent exponential decay, with more distinguishable P-S waves, as shown in Fig. 2.h. To distinguish TC from VT events the main criteria applied by OVDAS is the time difference of arrival (TDOA) between P and S waves (S − P). A sensor network composed by seismic stations near to the active crater (< 5 km) and other distant (> 10 km), allows to determine the origin of events. Thus, a minor TDOA in a sensor close to the crater suggests proximity to the volcanic edifice, otherwise the events are considered as TC.

X(m, k) =

N−1 X

x[m + n]w[n]e− j2πkn/K .

(1)

n=0

For a fixed m, X(m, k) is the discrete Fourier transform (DFT) of the sequence x[m + n]w[n], where w[n] is a window of length N < K, with w[n] = 0 outside the interval 0 ≤ n ≤ N − 1. The spectrogram S (m, k) is a representation of the squared modulus (in logarithmic scale) of X(m, k), and is given by (Lu & Zhang, 2009):   (2) S (m, k) = 10 log10 kX(m, k)k2 .

PT E

D

The database was generated by OVDAS analysts who selected the events from the Llaima volcano, specifically from the Laguna Verde (LAV) station. This station is located at 38.701◦ S 71.651◦ W, 7 km from the crater, as shown in Fig. 1. The data in the database were recorded between the years 2009 and 2017, where from 2009 to July 2010, the signals were recorded by the GeoSpace 1” sensor, from July 2010 to June 2011 by the Guralp 40TD 30” sensor and, currently, the Nanometrics 120” sensor. Unfortunately, the most active periods of the Llaima volcano had an instrumental network composed by sensors of short period. In June 2010, 30 seconds sensors were placed, but at that time, volcanic seismicity had drastically decreased (Franco, et al, 2014). The database was divided into two parts, training/validation and testing, where each part contained 2/3 and 1/3 of the total events, respectively, as presented in table 1.

CE

In this work, in order to avoid affecting the calculations due to event magnitude, the spectrogram is normalized as follows:

AC

S (m, k) − S min Sˆ (m, k) = , S max − S min

(3)

where S min and S max are the minimum and maximum spectrogram values, respectively. This procedure modified the original spectrogram range to values between 0 and 1, i.e. b S (m, k) ∈ [0, 1].

3.2. Proposed spectro-temporal curves 3. Spectro-temporal curves based on the spectrogram

In this proposal, two temporal curves were defined to characterize seismic signals, which were called the spectral evolution (SE) and the energy accumulation (ρ) ˆ curves. According to OVDAS analysts and seismologists, the spectrogram patterns of each event class are

This section presents the theoretical background required to define the proposed spectro-temporal curves. 5

ACCEPTED MANUSCRIPT

40 30 20 10

NU

SC

RI PT

0

MA

Figure 3: Sequence used to compute SE. (a) Spectrogram (S), (b) Filtered spectrogram (Sˆ μ ), (c) SE components, (d) SE curve.

calculated as:

essential when classifying an event. The spectrogram patterns correspond to the characteristic contours that result from its high energy components. Thus, the proposed curves are based on the energy and contour shape evolution given by the high energy components. The energy evolution in time is given by Eν (m), that represents the energy at m window considering only the high energy components. Eν (m) is computed as follows: K/2−1 X

Sˆ μ (m, k),

CE

Eν (m) =

PT E

D

νˆ (m) =

(4) ν(m) =

AC

1 Sˆ (m, k) ≥ μ 0 otherwise

B(m, k).

(8)

Then, the SE curve is obtained from the combination of the evolution of the high energy components over time Eν (m) and its contour shape νˆ (m), expressed as follows:

(5)

S E(m) = Eν (m) ∙ νˆ (m).

where B is a threshold matrix which determines the high ˆ energy components of S: B(m, k) =

K/2−1 X k=0

where Sˆ μ (m, k) represents the high energy components ˆ The high energy comof the normalized spectrogram S. ponents are the Sˆ elements above a threshold μ, defined by the average of Sˆ (Sˆ ). Then Sˆ μ (m, k) is calculated as:

(

(7)

where ν(m) is the bandwidth that concentrates most of the signal energy, νmax and νmin are the maximum and minimum ν values respectively. ν(m) quantifies the frequency bands that contains high energy components at m window, as follows:

k=0

Sˆ μ (m, k) = Sˆ (m, k) ∙ B(m, k),

ν(m) − νmin , νmax − νmin

(9)

Figure 3 shows the main steps for computing SE. First, the spectrogram is generated (Fig. 3.a) and filtered by hard thresholding (Fig. 3.b). Then, the SE components Eν and νˆ are calculated using Eq. (4) and Eq. (7), respectively (Fig. 3.c). Finally, the SE curve is computed by Eq. (9) (Fig. 3.d). Comparing figures 3.c and 3.d it can be observed that the curves Eν and SE are very

(6)

On the other hand, the contour shape was measured through the normalized high energy bandwidth νˆ (m), 6

RI PT

ACCEPTED MANUSCRIPT

1 ρmax

ρ(m) =

m 1 X

ρmax

S E(i),

i=0

(10)

MA

ρ(m) ˆ =

NU

similar, so Eν is able to characterize the events behavior. However, νˆ provides robustness, giving more weight to the energy components that have a wide bandwidth than the energy components that have a narrow bandwidth, which is a characteristic of the noise. The second curve used is the accumulation energy curve (ρ), ˆ which is calculated as follows:

SC

Figure 4: Characteristic behavior of the SE and ρˆ curves for each event class. The curves for LP (a), TR (b), VT (c) and TC (d) classes are shown from left to right, where the black line corresponds to the SE curve and the dashed brown line to the ρˆ curve.

where ρ(m) is the SE sum of accumulated value at m window and ρmax is the maximum ρ value. Hence, ρˆ ranges from 0 and 1 and the parameters extracted from it are limited in order to describe curve form. Figure 4 shows the characteristic behavior of each event class using the SE and ρˆ curves. The typical curves for LP (Fig. 4.a), TR (Fig. 4.b), VT (Fig. 4.c) and TC (Fig. 4.d) classes are shown from left to right, where the black line corresponds to the SE curve and the dashed brown line to the ρˆ curve. From both curves, it is possible to determine event class. However, ρˆ is smoother and less noisy than the SE curve. Thus, the parameters extracted from ρˆ are better for classifying as they are more robust. Figure 5 presents a contrast in terms of the shape of the characteristic ρˆ curves for each event class, where the different behavior of each class can be seen, including their slopes and inflection points, among other features useful to class discrimination.

LP TR VT TC

Normalized time (%)

Figure 5: Sample of the energy curve from each class.

AC

CE

PT E

D

ters of the function show an optimal fit. The adjustment parameters and the quality of the estimate are used as features for the classifying step. In this study, a first order and an exponential approximation function were selected due to their ease in fitting the energy curves. It should be noted that the exponential approximation provides a good representation of ρˆ for the LP and VT classes, as shown in Fig. 5. The features extracted from the first order approximation are the slope m s , the bias b and the approximation error (R2s ). For the exponential estimate, the curve is fitted by the function 1 − eδx where the parameter δ and the approximation error (R2δ ) are the extracted features. δ was calculated by solving the least-squares problem and the approximation errors, R2δ and R2s , were calculated using the squared correlation between ρˆ and its approximation (or other curve like a template), through Eq. (11). Let A be a curve fitted to a ρˆ curve. The squared correlation is defined as follow:

3.3. Features extracted from spectral temporal curves In this section, seven feature extraction techniques are presented, which were applied to obtain the characteristics of the energy accumulation curve (ρ). ˆ The features are extracted from two procedures: curve fitting and template similitude.

R2A =

3.3.1. Curve fitting features Curve fitting is the adjustment of the mathematical function to a series of data points where the parame-

cov(A, ρ) ˆ σA σρˆ

!2

,

(11)

where cov(A, ρ) ˆ is the covariance between A and ρ, ˆ σA and σρˆ are the standard deviation of A and ρ, ˆ respec7

ACCEPTED MANUSCRIPT

3.3.2. Template similitude features The template similitude is the correspondence measure between a template (ϕ) and an energy curve ρ. ˆ The similitude was measured using two techniques: dynamic time warping (DTW) and the squared correlation. DTW is a technique used to find an optimal alignment between two given (time-dependent) sequences under certain restrictions (M¨uller, 2007). In this study, the extracted feature is the distance between ϕ and ρˆ called DT Wϕ . The alignment was performed as in the study conducted by Turetsky & Ellis (2003), which is detailed in Appendix A. It provides the cost matrix D and the optimal path given by the sequences of length J and the p and q coordinates. Next, DT Wϕ is calculated as follows:

SC

RI PT

tively. Figure 6 shows an example of the two approximation functions. It is worth noting that the exponential approximation shows a better fit to ρˆ than the first order approximation. The representativeness of the functions is less relevant for the purpose of this article, in this paper the objective of curve fitting is to obtain effective parameters for discriminating between classes. Figure 7 shows the clustering of classes using curve fitting features. Figure 7.a shows the combination of two first order approximation features, which separate a considerable number of TRs, VTs and TCs. Figure 7.b shows the combination of two exponential approximation features. This combination allows all classes to be separated with low overlap.

DT Wϕ =

Normalized time (%)

NU

The squared correlation was defined in section 3.3.1 and calculated using Eq. (11), being the feature defined as R2ϕ . For each event, there are four template similitude features that have to be compared with the template of each class. Figure 8 shows the class clustering that results from the combination of the DT Wϕ and the R2ϕ features. The classes are separated with low overlap.

D

Figure 6: Energy curve and mathematical approximations.

CE

PT E

LP TR VT TC

Figure 8: Clustering of classes using similitude features.

AC

LP TR VT TC

(12)

MA

Energy curve Exponential approx. First order approx.

J 1X D(p( j), q( j)). J j=1

4. Experiments The performance and contribution of ρˆ features were measured by comparing them to the automatic classifiers proposed in Curilem et al. (2016). This work was established as the baseline for the proposed features. The event classes considered were LP, TR, VT and TC. The following subsections define the experimental set-up and the experiments.

Figure 7: Clustering of classes using curve fitting features. At the top (a) features extracted from the first order approximation and (b) from the exponential approximation are shown.

8

ACCEPTED MANUSCRIPT

(c)

x1 x2 x3

[ x1, x2, x7 ]

LP

x7 ]

TR

[ x 3,

[ x1, x2, x7 ]

xJ

[x1, x4, x5, x6]

RI PT

...

(d) 1 0

VT 0 1 TC

LP

Figure 9: Scheme of the implemented classifying structure one versus all.

SC

features combination for each SVM classifier are shown in 9.

NU

4.1.2. Pre-processing step Signals prior to 2011 were sampled at 50 Hz and then at 100 Hz. Hence, the signals were standardized to a 100 Hz sampling rate and the amplitude was measured in micrometers/s (μm/s). Next, all the signals were filtered by a 10th order Butterworth bandpass filter between 1 and 10 Hz as this is the frequency range of greatest interest for Llaima volcanic events. In order to generate SE, the spectrogram was set up using 1024 FFT samples and a hamming window with length N = 100 samples (1 second), which was overlapped in 90%. The window type was chosen because it provides a smooth truncation and it has given good results according to the literature (Oppenheim & Schafer, 2009). The window length and overlap were selected experimentally, where the selected values gave the best result. The window size was tested between 1 and 5 seconds, it was defined due to the lowest relevant frequency information of volcanic events begins between 0.3 and 1.0 Hz, according to analysts criteria/knowledge. The overlap was tested between 50% and 90% with step of 10%. The amount of DFT bins was defined according to the window’s length. It should be equal or higher to 128, thus 1024 bins were established, to give more resolution. Finally, ρˆ was considered from 1% to 99%, due to the fact that the start and end of the events are not uniform. The features based on template similitude were generated through four templates, one for each class. To define the templates, the ρˆ curve was generated for all the events of the training set. Then, one ρˆ curve that exhibited both a high correlation within the same class and low correlation with other classes was chosen for each class.

AC

CE

PT E

D

MA

4.1.1. Classifier design The classification was conducted using the ”one versus all” structure as recommended by Curilem et al. (2014, 2016). This structure requires at least N-1 classifiers to classify N classes whereas each classifier is designed to discriminate one class from the others. In this work the classifier structure was implemented with N SVM classifiers (one for each class), using the Radial Basis Function (RBF-kernel). The SVM is a binary classification algorithm that searches for an optimal hyperplane that separates the data into two classes. In this paper, the classifier was implemented using the Radial Basis Function (RBF-kernel). The kernel function performed the data projection from the original input space to the feature space where the two classes are linearly separable (Cortes & Vapnik, 1995). The scheme of the proposed classifying structure is presented in Fig. 9, where it can be observed that the SVM classifiers work simultaneously. The first step corresponds to feature extraction for all classifiers (J features), as shows in Fig. 9.a. The second step is feature selection process, where the characteristics are selected for each classifier, as shows in Fig. 9.b. The feature vectors can be different between classifiers, with unequal dimension and features selected. Then, at the same time each classifier evaluates its input (feature vector), assigning value 1 (true) if the input is localized into the region of the hyperplane assigned to interest class, in the other case the value assigned is 0 (false). As presented in Fig. 9.c, this may generate a conflict when more than one classifier has a positive output (true) or when all outputs are 0. In order to ensure a single output, a confidence stage was implemented as in Curilem et al. (2014). The confidence block selects the output of the most confident classifier, thus avoiding a conflict. For more details about this stage see Appendix C. To train the SVM classifiers, it is necessary to tune the regularization constant c = 2 x , parameter that controls the complexity of the SVM and the σ = 2y which control the standard deviation of the Gaussian kernel function. To choose these parameters, each feature combination was trained by sweeping a grid for the parameter x and y, where x took values from -5 to 15 and y from -15 to 9, with a step of 1, being independently adjusted for each classifier. The training was developed through a 4 folds cross-validation, with each fold containing approximately 490 events, the selected c and σ values and

(b)

(a)

Confidence block

4.1. Experimental set-up The experiments were developed using Matlab software version R2015. The classifier design and the preprocessing step are described as follows.

9

ACCEPTED MANUSCRIPT

Table 2: Overview of the extracted features used to classify each event class.

Description

Ax An Ad Fn E5 δ R2δ ms b R2s DT WLP DT WT R DT WVT DT WTC R2LP R2T R R2VT R2TC

Maximum value of the amplitude obtained from the time records of each event. Average value of the amplitude obtained from the time records of each event. Median value of the amplitude obtained from the time records of each event. Average of the five highest frequency peaks obtained from the Fourier transform. Percentage of energy in the 5th band [1.56–3.13] Hz, obtained from the wavelet transform of each event. Parameter extracted from curve fitting by the exponential function. Approximation error of curve fitting by the exponential function. Slope, parameter extracted from curve fitting by the first order function. Bias, parameter extracted from curve fitting by the first order function. Approximation error of curve fitting by the first order function. Similitude measure between ρˆ and a LP template using the DTW distance. Similitude measure between ρˆ and a TR template using the DTW distance. Similitude measure between ρˆ and a VT template using the DTW distance. Similitude measure between ρˆ and a TC template using the DTW distance. Similitude measure between ρˆ and a LP template using the squared correlation. Similitude measure between ρˆ and a TR template using the squared correlation. Similitude measure between ρˆ and a VT template using the squared correlation. Similitude measure between ρˆ and a TC template using the squared correlation.

MA

NU

SC

RI PT

Feature

D

4.2. Experiments description The following experiments were developed in order to measure the individual and combined performance of ρˆ features and the set proposed by (Curilem et al., 2016). Table 2 describes all the features used in the experiments where the first five correspond to those used in Curilem et al. (2016) and the others were extracted from the ρˆ curve.

tion considered only the features that had an EER average lower than 30%, according to Experiment 1. However, the number of features was limited to five to avoid augmenting the complexity of the classifier.

PT E

For these experiments, the performances were measured using statistical indices such as the exactitude, error, precision, false discovery rate (FDR), sensitivity and specificity. These indices are defined in Appendix B.

CE

• Experiment 1 - Individual feature evaluation: All features were individually tested by designing a SVM classifier for each class. The performance of the classifiers was measured through the equal error rate (EER) (Toh et al., 2008) average from the cross validation training.

5. Results and discussion This section presents and analyzes the results from the experiments detailed in section 4.2, which corresponds to an individual and a combined evaluation of both, new and conventional features. Table 3 shows the EER average for each feature and for each class. The features that have an EER average lower than 15% are highlighted. Thus, individually the features with a better performance in each class are δ for LP and VT, E5 for TR and R2δ for TC with 3.60%, 5.50%, 4.00% and 4.31%, respectively. As described in section 4.2, experiment 2 was developed to measure the contribution of ρˆ features compared to conventional features during classification. Thus, the tables 4, 5 and 6 show the contingency matrix of the

AC

• Experiment 2 - Combined features evaluation: In this experiment, three classifying structures were created. The first classifying structure was trained using the features of the previous work and its performance was established as the baseline for comparison to the overall performance. Then, two tests were performed to evaluate the contribution of the new parameters. In the first test, two ρˆ features were added to the baseline features. In the second test, a wrapper feature selection strategy (Guyon et al., 2006) was performed to find the best combination of all features (baseline and ρˆ features) for each class. The combina10

ACCEPTED MANUSCRIPT

Table 3: Average EER (%) of the best SVM classifier obtained from the individual evaluation of features.

35.99 34.76 23.82 30.15 4.00 33.82 27.10 35.53 34.04 28.01 10.77 11.63 10.09 33.56 32.15 32.84 28.42 34.42

25.33 23.77 44.16 36.53 12.14 12.64 20.36 13.43 5.50 22.35 10.83 15.54 37.34 25.43 10.26 10.45 15.32 11.02

31.82 35.10 39.88 12.85 31.25 11.58 29.04 35.88 21.17 4.31 13.13 19.98 39.39 24.77 7.36 8.11 7.99 17.32

D

Table 4: Contingency matrix of the baseline classifier structure.

Classifier LP TR

VT

TC

Total

LP TR VT TC Total

230 2 5 18 255

7 0 81 17 105

8 3 16 435 462

249 163 102 471 985

AC

CE

PT E

Expert

Table 5: Contingency matrix of the best classifier structure obtained from baseline combinations and additional ρˆ features.

Expert LP TR VT TC Total

Classifier LP TR

VT

TC

Total

241 1 4 1 247

3 0 98 4 105

4 1 0 466 471

249 163 102 471 985

1 161 0 0 162

Classifier LP TR

VT

TC

Total

LP TR VT TC Total

241 1 6 2 250

5 0 95 3 103

1 1 0 462 464

249 163 102 471 985

2 161 1 4 168

ture had an exactitude and precision of over 90% for each class, except for VTs which had 77.14% precision. This result was similar to the one obtained by Curilem et al. (2016) and it is difficult to improve it significantly. However, when adding less than three ρˆ features to the baseline features, the error and FDR were significantly reduced from 4.11 to 0.96% and from 10.39 to 2.69%, respectively. The performance of the best features combination gave similar results. Its performance was slightly lower than the previous, but with a smaller number of features. Table 8 shows a comparative statistic between the baseline classifying structure and the other structures in terms of error, FDR and the dimensionality of features. For the baseline + ρˆ features, the error and the FDR decreased around 75% for all classes. For the best features combination, the error and FDR were reduced on average by 67.90% and 61.60%, respectively, being the worst result for TRs and the best for TCs. For TRs, the error decreased by only 10% and the FDR increased by 35.83%. However, it is important to emphasize that the dimensionality of the features was reduced by 50%. In contrast, for the baseline + ρˆ features, the error and FDR decreased by 70% and 79.88%, respectively, and the dimensionality increased by 50%. For TCs, the results were similar for both experiments. The best features combination presented an error of 1.59% lower and a FDR above 10.78%. Thus, for all classes and especially for TRs, it was possible to improve the results significantly. Nevertheless, it is still necessary to increase the dimensionality. Table 9 summarizes the combination of features for each class and for each classifying structure, where the selected features are marked with an ”x”. In the second experiment (baseline + ρˆ features), most of the added features were related to the template similitude measures. For test 2 (best features combination), other features were selected, mostly extracted from ρˆ curves and combined with E5 or Fn baseline features. It is important to consider that the dimensionality of features was

baseline classifier, the combination of the baseline features with the ρˆ features and the combination of best features, respectively.

4 158 0 1 163

Expert

TC

RI PT

VT

SC

27.75 33.12 26.28 21.33 39.22 19.88 32.59 41.52 3.60 20.72 12.65 14.13 32.30 11.07 19.08 18.52 18.83 26.76

TR

MA

An Ad Ax Fn E5 b ms R2s δ R2δ DT WLP DT WTC DT WT R DT WVT R2LP R2TC R2T R R2VT

LP

NU

Feature

Table 6: Contingency matrix of the classifier structure obtained from the best features combination.

Table 7 shows the statistical indices obtained from the contingency tables. The baseline classifying struc11

ACCEPTED MANUSCRIPT

Table 7: Performance indices (%) for each classifier structure.

Baseline LP TR

Exactitude Error Precision FDR Sensitivity Specificity Features

95.53 4.47 90.20 9.80 92.37 96.60 5

98.98 1.02 96.93 3.07 96.93 99.39 4

VT

TC

X

95.43 4.57 77.14 22.86 79.41 97.28 4

93.60 6.40 94.16 5.84 92.36 94.75 5

95.89 4.11 89.61 10.39 90.27 97.01 4.50

Baseline + ρˆ features LP TR VT 98.58 1.42 97.57 2.43 96.79 99.18 7

99.70 0.30 99.38 0.62 98.77 99.88 6

98.88 1.12 93.33 6.67 96.08 99.21 5

TC

X

98.98 1.02 98.94 1.06 98.94 99.03 7

99.04 0.96 97.31 2.69 97.64 99.32 6.25

Best features combination LP TR VT TC 98.27 1.73 96.40 3.60 96.79 98.78 3

99.09 0.91 95.83 4.17 98.77 99.15 2

RI PT

Classifier Index

98.48 1.52 92.23 7.77 93.14 99.09 3

98.88 1.12 99.57 0.43 98.09 99.61 5

X 98.68 1.32 96.01 3.99 96.70 99.16 3.25

Table 8: Comparison between the baseline and the classifier structures with ρˆ features.

Error reduction FDR reduction ΔFeatures

68.18 75.22 40.00

70.00 79.88 50.00

75.56 70.83 25.00

TC

X

84.13 81.84 40.00

76.54 74.08 38.89

61.36 63.28 -40.00

10.00 -35.83 -50.00

66.67 66.02 -25.00

TC

X

82.54 92.62 0.00

67.90 61.60 -27.78

From the experiments related to the combination of features, it is possible to conclude that the features extracted from ρˆ complement the baseline features with the most relevant information related to frequency and the evolution of energy, such as E5, Fn, δ and DT WLP . Among the principal achievements of this study was the reduction of the error and the FDR by more than 70%. Furthermore, these features reduced the dimensionality of features by around 30%, decreasing the complexity of the classifying structures and speeding up the recognition process. Future research should seek to test these novel parameters for classifying seismic events that are detected and segmented automatically at the Llaima volcano and at other volcanos with similar physical features and activity.

MA

reduced around 30%. These tests prove that the features extracted from ρˆ complements the information provided by the conventional features, being that the most relevant information added was related to the frequency, the energy and the evolution of energy over time.

Best features combination LP TR VT

SC

Baseline + ρˆ features LP TR VT

NU

Classifier Variation (%)

D

6. Conclusions

CE

PT E

This paper proposed a new method for extracting features from the spectrograms of volcanic seismic signals, based on the knowledge of experts. First, a method was proposed for constructing the spectro-temporal curves (SE and ρ), ˆ which combine the information of the frequency and the energy of high energy components along the time duration of an event. The curves describe the behavior evolution of the energy, weighted by the bandwidth frequency behavior over time. The ρˆ curves showed a characteristic behavior for each class. The second step of the method was to extract features from the spectro-temporal curve. The features extracted from the ρˆ curves describe its shape and measure its similitude with templates or mathematical functions that represent the characteristic ρ-curve ˆ for each class. These features together with the features proposed by Curilem et al. (2016), considered as baseline features, were tested individually and then combined. Individually, most of the new features showed a good performance, highlighting features that described the exponential start and decay of events in terms of energy, such as δ and R2δ . Moreover, these two new individual features showed better performance than the baseline features, for the LP, VT and TC classes.

AC

Acknowledgments The authors would like to thank the project FONDEF IDeA IT15I10027 for financing this study. Also, the authors would like to thank OVDAS, which provided the data and geological knowledge that supported the analysis of the results.

Appendix A. Dynamic Time Warping (DTW) DTW is a well-known technique used to compare two (time-dependent) sequences X := (x1 , x2 , ..., xN ) of length N ∈ N and Y := (y1 , y2 , ..., y M ) of length M 12

ACCEPTED MANUSCRIPT

Table 9: Configuration of features and SVM-tune parameters for each classifier structure.

x

x x

x x x x x

x x

Baseline + ρˆ features LP TR VT TC x x x x x x

x x x

x x

x x x x x x

x x

x

Best features combination LP TR VT TC

x

x x

x x

x

x

x

x

x

x

x 2 29

3

2 215

3

2 211

3

2 213

1

2 23

5

2 213

3

2 29

1

2 23

x x

x

21 211

23 215

23 213

23 211

one class (positive class (P)) versus all other classes (negative class (N)). The contingency matrix contains true positive (TP), true negative (TN), false positive (FP) and false negative (FN) elements. TPs is the amount of P elements correctly classified as P, TNs is the amount of N elements correctly classified as N, FPs are the N elements wrongly classified as P, and FNs are the P elements wrongly classified as N. The interest indices extracted from the contingency table are (Fawcett, 2006):

PT E

D

MA

∈ N. These sequences may be discrete signals (timeseries) or, more generally, feature sequences sampled at equidistant points in time (M¨uller, 2007). To compare the sequences a similarity measure is calculated from a cost matrix C. The cost matrix is computed based on the two sequences to be synchronized (X and Y), where each element C(n,m) represents the distance between the elements xn and ym from each sequence. Each value in the N x M matrix is computed as follows (Foote, 1999):

x

x

x

1

σ c

x x x

TC

RI PT

x x x x x

VT

SC

An Ad Ax Fn E5 DT WLP DT WT R DT WVT R2LP R2TC b R2δ δ

Baseline LP TR

NU

SVM

Feature

Classifier Parameter

xn ∙ ym . (A.1) kxn kkym k Then, the cost-minimizing alignment path is determined from this matrix via dynamic programming (DP), tracing a path from the origin to (N - 1, M - 1), as described by Turetsky & Ellis (2003). DP consists of two stages, forward and traceback. In the forward step, the lowest-cost path is calculated for all the neighbour points, plus the cost to get from the neighbour to the point. Thus, the path cost matrix D is generated as described by the algorithm 1, between lines 2 and 14. In the traceback stage, the actual path is obtained itself by recursively looking up the point providing the best antecedent for each point on the path given by the sequences p and q coordinates.

Table B.10: Illustrative contingency matrix for two classes.

AC

CE

C(n, m) =

Expert

Classifier P

N

Total

P N Total

TP FP TP + FP

FN TN TN + FN

TP + FN TN + FP TP + FP + TN + FN

• Exactitude (Ex). Indicates the proportion of events correctly classified, calculated as:

Ex =

(T P + T N) ∙ 100 = 100 − Er. (B.1) T P + FP + T N + FN

• Error (Er). Measures the proportion of events wrongly classified, being calculated as:

Appendix B. Performance indices The performance indices are calculated from the contingency matrix, as shown in table B.10, considering

Er = 13

(FP + FN) ∙ 100 = 100 − Ex. (B.2) T P + FP + T N + FN

ACCEPTED MANUSCRIPT

• Precision or positive predictive value (PPV). Measures the proportion of events correctly classified as positive class, it is calculated as follows: PPV =

(B.3)

RI PT

• False discovery rate (FDR). Measures the proportion of events wrongly classified as positive class, being calculated as: FDR =

FP ∙ 100 = 100 − PPV. T P + FP

(B.4)

T PR =

T P ∙ 100 . T P + FN

(B.5)

NU

SC

• Sensitivity or true positive rate (TPR). Indicates the proportion of positive events correctly classified, it is calculated as follows:

• textbfSpecificity or true negative rate (TNR). Indicates the proportion of negative events correctly classified, being calculated as:

MA

T NR =

D

Appendix C. Bayes-based (BBCM)

T N ∙ 100 . T N + FP

(B.6)

confidence

measure

X

Feature selection stage

AC

CE

PT E

Bayes-based confidence measure (BBCM) (Huenup´an et al., 2008; Yoma et al., 2005) corresponds to the ordinary Bayes fusion weighted by the reliability of each individual classifier and provides a formal model for the reliability in a classification task. Has been used in speech recognition and combination of classifiers, the latter will be described in this appendix. CL1

CL2 . . .

Combination stage

Algorithm 1 Dynamic programming algorithm. 1: procedure DP(C) [N, M] = size(C) . Forward stage 2: D = zeros(N + 1, M + 1) . D Initialization 3: loc = zeros(N, M) 4: D(0, 1 : M) = In f 5: D(1 : N, 0) = In f 6: D(1 : N, 1 : M) = C 7: for i = 0 : N do 8: for j = 0 : M do 9: [val, loc(i, j)] = ... 10: min([D(i, j), D(i, j + 1), D(i + 1, j)]) 11: D(i + 1, j + 1) = D(i + 1, j + 1) + val 12: end for 13: end for 14: i= N−1 . Traceback from top left 15: j= M−1 16: p=i 17: q= j 18: while i > 0& j > 0 do 19: if loc(i, j) == 0 then 20: i=i−1 21: j= j−1 22: else if loc(i, j) == 1 then 23: i=i−1 24: else if loc(i, j) == 2 then 25: j= j−1 26: else 27: error 28: end if 29: p = [i, p] 30: q = [ j, q] 31: end while 32: return D, p, q 33: 34: end procedure

T P ∙ 100 = 100 − FDR. T P + FP

Final Decision

CLJ

Figure C.10: General scheme for combination of classifier using confidence measure.

Figure C.10 shows a typical scheme for combination of classifier with a feature selection stage, where X is the feature vector, X j is the input (selected characteristic from X) for the classifier j, CL j , with J the number 14

ACCEPTED MANUSCRIPT

of classifiers; S CL j is the score at the output of classifier j at input X j ; and dCL j is the local decision or classification due to classifier j. BBCMCL j is the confidence measureh estimation for the outputs of classifier j and i BBCM S CL j is the score of reliability. Finally, D(X) is the final decision or classification for the entire system, where the combination can be made by different schemes of merging the local decision as calculating the mean of all the BBCM, maximum, product, Majority Vote Rule (or weighted by the BBCM score), among others (Kittler et al., 1998). The BBCM associated with classifier CL j given a feature vector X j is specified by the following equation:

RI PT

where M is the total number of classes,  Pr S CL j (X j )|Cm is the a priori probability or likelihood for the score S CL j (X j ) given the class Cm and Pr (Cm ) is the a priori probability for class Cm .

NU

Figure C.11: Example of BBCM curves for TR (blue line) and VT (dashed red line) SVM classifiers.

Figure C.11 shows examples of BBCM curves, obtained with Eq. (C.2). A lower value of BBCM occurs when the classification score is close to the decision threshold. In the case of a SVM classifier the score is the distance between the sample and the hyperplane. If the sample is closer to the hyperplane, the decision of the classifier is less reliable, as a consequence a lower BBCM.

MA

where ”dCL j is ok” denotes that the local decision of   classifier j, CL j , is correct. Pr dCL j is ok|S CL j (X j ) means the probability of the local decision dCL j (ex. target or non-target for a SVM Classifier) is correct given the score S CL j (X j ). By applying the Bayes theorem, Eq. (C.1) can be written as:

SC

    BBCM S CL j (X j ) = Pr dCL j is ok | S CL j (X j ) , (C.1)

D

  BBCM S CL j (X j ) =     (C.2) Pr S CL j (X j )|dCL j is ok ∙ Pr dCL j is ok   . Pr S CL j (X j )   The BBCM S CL j (X j ) is a probability itself. More  over, the distribution Pr S CL j (X j )|dCL j is ok and the   probability Pr dCL j is ok provide information about the recognition engine’s performance. To compute the BBCM from Eq. (C.2), it is necessary three a priori p.d.f. and can be estimated in the training process:   • Pr S CL j (X j )|dCL j is ok it is calculated using the classifier score S CL j (X j ) of the samples that have been correctly classified or labeled by the classifier CL j .   • Pr dCL j is ok is the performance of the classifier, it is the a priori probability of a correct decision of the classifier CL j .   • Pr S CL j is the probability of the score S CL j (X j ) for the classifier CL j , and it is estimated as:

PT E

References

AC

CE

Aki, K., Fehler, M., & Das, S. (1977). Source mechanism of volcanic tremor: fluid-driven crack models and their application to the 1963 kilauea eruption. Journal of Volcanology and Geothermal Research, 2, 259 – 287. doi:10.1016/0377-0273(77)90003-8. Aki, K., & Koyanagi, R. (1981). Deep volcanic tremor and magma ascent mechanism under kilauea, hawaii. Journal of Geophysical Research: Solid Earth, 86, 7095–7109. doi:10.1029/JB086iB08p07095. Arciniega-Ceballos, A., Chouet, B., & Dawson, P. (2003). Longperiod events and tremor at popocatepetl volcano (1994–2000) and their broadband characteristics. Bulletin of Volcanology, 65, 124– 135. doi:10.1007/s00445-002-0248-8. Arciniega-Ceballos, A., Chouet, B., Dawson, P., & Asch, G. (2008). Broadband seismic measurements of degassing activity associated with lava effusion at popocat´epetl volcano, mexico. Journal of Volcanology and Geothermal Research, 170, 12 – 23. doi:10.1016/j.jvolgeores.2007.09.007. The 1994-Present eruption of Popocat´epetl: Background, current activity, and impacts. Benoit, J. P., & McNutt, S. R. (1997). New constraints on source processes of volcanic tremor at arenal volcano, costa rica, using broadband seismic data. Geophysical Research Letters, 24, 449– 452. doi:10.1029/97GL00179. Beyreuther, M., Carniel, R., & Wassermann, J. (2008). Continuous hidden markov models: Application to automatic earthquake detection and classification at las can˜adas caldera, tenerife. Journal of Volcanology and Geothermal Research, 176, 513 – 518. doi:10.1016/j.jvolgeores.2008.04.021.

M  X    Pr S CL j (X j )|Cm ∙ Pr (Cm ) , Pr S CL j (X j ) = m=1

(C.3) 15

ACCEPTED MANUSCRIPT

NU

SC

RI PT

doi:10.1029/JB090iB12p10237. Curilem, G., Vergara, J., Fuentealba, G., Acu˜na, G., & Chac´on, M. (2009). Classification of seismic signals at villarrica volcano (chile) using neural networks and genetic algorithms. Journal of Volcanology and Geothermal Research, 180, 1 – 8. doi:10.1016/j.jvolgeores.2008.12.002. Curilem, M., Huenupan, F., Beltr´an, D., Martin, C. S., Fuentealba, G., Franco, L., Cardona, C., Acu˜na, G., Chac´on, M., Khan, M. S., & Yoma, N. B. (2016). Pattern recognition applied to seismic signals of llaima volcano (chile): An evaluation of station-dependent classifiers. Journal of Volcanology and Geothermal Research, 315, 15 – 27. doi:10.1016/j.jvolgeores.2016.02.006. Curilem, M., Vergara, J., Martin, C. S., Fuentealba, G., Cardona, C., Huenupan, F., Chac´on, M., Khan, M. S., Hussein, W., & Yoma, N. B. (2014). Pattern recognition applied to seismic signals of the llaima volcano (chile): An analysis of the events’ features. Journal of Volcanology and Geothermal Research, 282, 134 – 147. doi:10.1016/j.jvolgeores.2014.06.004. Esposito, A., Giudicepietro, F., D’Auria, L., Scarpetta, S., Martini, M., Coltelli, M., & Marinaro, M. (2008). Unsupervised neural analysis of very-long-period events at stromboli volcano using the self-organizing maps. Bulletin of the Seismological Society of America, 98, 2449–2459. doi:10.1785/0120070110. Fawcett, T. (2006). An introduction to roc analysis. Pattern Recognition Letters, 27, 861 – 874. doi:10.1016/j.patrec.2005.10.010. Ferrazzini, V., & Aki, K. (1987). Slow waves trapped in a fluidfilled infinite crack: Implication for volcanic tremor. Journal of Geophysical Research: Solid Earth, 92, 9215–9223. doi:10.1029/JB092iB09p09215. Foote, J. (1999). Methods for the automatic analysis of music and audio. Technical Report FXPAL-TR-99-038 FX Palo Alto Laboratory, Inc. Fournier, R. O. (1999). Hydrothermal processes related to movement of fluid from plastic into brittle rock in the magmaticepithermal environment. Economic Geology, 94, 1193–1211. doi:10.2113/gsecongeo.94.8.1193. Franco, L., Palma, J., Gil-Cruz, F., & San Martin, J. (2014). La actividad sismica asociada a erupciones estrombolianas violentas en el volcan llaima, chile (2007-2010). Earth Sciences Research Journal, 18, 338–339. Gil-Cruz, F., & Chouet, B. A. (1997). Long-period events, the most characteristic seismicity accompanying the emplacement and extrusion of a lava dome in galeras volcano, colombia, in 1991. Journal of Volcanology and Geothermal Research, 77, 121–158. doi:10.1016/S0377-0273(96)00091-1. Gordeev, E. (1993). Modeling of volcanic tremor as explosive point sources in a singlelayered, elastic halfspace. Journal of Geophysical Research: Solid Earth, 98, 19687–19703. doi:10.1029/93JB00348. Goto, A. (1999). A new model for volcanic earthquake at unzen volcano: Melt rupture model. Geophysical Research Letters, 26, 2541–2544. doi:10.1029/1999GL900569. Guyon, I., Gunn, S., Nikravesh, M., & Zadeh, L. (2006). Feature Extraction: Foundations and Applications. Studies in Fuzziness and Soft Computing. Springer Berlin Heidelberg. doi:10.1007/978-3-540-35488-8. Hibert, C., Provost, F., Malet, J.-P., Maggi, A., Stumpf, A., & Ferrazzini, V. (2017). Automatic identification of rockfalls and volcano-tectonic earthquakes at the piton de la fournaise volcano using a random forest algorithm. Journal of Volcanology and Geothermal Research, 340, 130 – 142. doi:10.1016/j.jvolgeores.2017.04.015. Huenup´an, F., Yoma, N. B., Molina, C., & Garret´on, C. (2008). Confidence based multiple classifier fusion in speaker

AC

CE

PT E

D

MA

Bhatti, S. M., Khan, M. S., Wuth, J., Huenupan, F., Curilem, M., Franco, L., & Yoma, N. B. (2016). Automatic detection of volcano-seismic events by modeling state and event duration in hidden markov models. Journal of Volcanology and Geothermal Research, 324, 134 – 143. doi:10.1016/j.jvolgeores.2016.05.015. Bicego, M., Acosta-Mu˜noz, C., & Orozco-Alzate, M. (2013). Classification of seismic volcanic signals using hiddenmarkov-model-based generative embeddings. IEEE Transactions on Geoscience and Remote Sensing, 51, 3400–3409. doi:10.1109/TGRS.2012.2220370. Branch, C. (1976). Development of porphyry copper and stratiform volcanogenic ore bodies during the life cycle of andesitic stratovolcanoes. Volcanism in Australasia. Elsevier, Amsterdam, (pp. 337–342). Brown, G. E. (2017). Volcanic eruptions and what triggers them. Elements, 13, 3–4. doi:10.2113/gselements.13.1.3. Carniel, R. (2014). Characterization of volcanic regimes and identification of significant transitions using geophysical data: a review. Bulletin of Volcanology, 76, 848. doi:10.1007/s00445-014-0848-0. Carniel, R., Jolly, A. D., & Barbui, L. (2013). Analysis of phreatic events at ruapehu volcano, new zealand using a new {SOM} approach. Journal of Volcanology and Geothermal Research, 254, 69 – 79. doi:10.1016/j.jvolgeores.2012.12.026. Chester, D., Degg, M., Duncan, A., & Guest, J. (2000). The increasing exposure of cities to the effects of volcanic eruptions: a global survey. Global Environmental Change Part B: Environmental Hazards, 2, 89 – 103. doi:10.1016/S1464-2867(01)00004-3. Chouet, B. (1981). Ground motion in the near field of a fluiddriven crack and its interpretation in the study of shallow volcanic tremor. Journal of Geophysical Research: Solid Earth, 86, 5985–6016. doi:10.1029/JB086iB07p05985. Chouet, B. (1985). Excitation of a buried magmatic pipe: A seismic source model for volcanic tremor. Journal of Geophysical Research: Solid Earth, 90, 1881–1893. doi:10.1029/JB090iB02p01881. Chouet, B. (1986). Dynamics of a fluiddriven crack in three dimensions by the finite difference method. Journal of Geophysical Research: Solid Earth, 91, 13967–13992. doi:10.1029/JB091iB14p13967. Chouet, B. (1988). Resonance of a fluid-driven crack: Radiation properties and implications for the source of long-period events and harmonic tremor. Journal of Geophysical Research: Solid Earth, 93, 4375–4400. Chouet, B. (1992). A seismic model for the source of long-period events and harmonic tremor. In P. Gasparini, R. Scarpa, & K. Aki (Eds.), Volcanic Seismology (pp. 133–156). Berlin, Heidelberg: Springer Berlin Heidelberg. Chouet, B., & Dawson, P. (2011). Shallow conduit system at kilauea volcano, hawaii, revealed by seismic signals associated with degassing bursts. Journal of Geophysical Research: Solid Earth, 116. doi:10.1029/2011JB008677. Chouet, B. A. (1996). Long-period volcano seismicity: its source and use in eruption forecasting. Nature, 380, 309. doi:10.1038/380309a0. Chouet, B. A., & Matoza, R. S. (2013). A multi-decadal view of seismic methods for detecting precursors of magma movement and eruption. Journal of Volcanology and Geothermal Research, 252, 108 – 175. doi:10.1016/j.jvolgeores.2012.11.013. Cortes, C., & Vapnik, V. (1995). Support-vector networks. Machine learning, 20, 273–297. Crosson, R. S., & Bame, D. A. (1985). A spherical source model for low frequency volcanic earthquakes. Journal of Geophysical Research: Solid Earth, 90, 10237–10247.

16

ACCEPTED MANUSCRIPT

NU

SC

RI PT

2004–2006 (pp. 27–60). US Geological Survey Professional Paper Virginia volume 1750. M¨uller, M. (2007). Information retrieval for music and motion. Springer Science & Business Media. doi:10.1007/978-3-540-74048-3_4. Naranjo, J., & Moreno, H. (2005). Geolog´ıa del volc´an llaima, regi´on de la araucan´ıa, escala: 1:50.000. Carta Geol´ogica de Chile, Serie Geolog´ıa B´asica, 88, 33. Neuberg, J., Luckett, R., Baptie, B., & Olsen, K. (2000). Models of tremor and low-frequency earthquake swarms on montserrat. Journal of Volcanology and Geothermal Research, 101, 83 – 104. doi:10.1016/S0377-0273(00)00169-4. Neuberg, J. W., Tuffen, H., Collier, L., Green, D., Powell, T., & Dingwell, D. (2006). The trigger mechanism of low-frequency earthquakes on montserrat. Journal of Volcanology and Geothermal Research, 153, 37–50. doi:10.1016/j.jvolgeores.2005.08.008. Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing. (3rd ed.). Upper Saddle River, NJ, USA: Prentice Hall Press. Orozco-Alzate, M., Londo˜no-Bonilla, J. M., & Acosta-Mu˜noz, C. (2012). The automated identification of volcanic earthquakes: Concepts, applications and challenges. In S. D’Amico (Ed.), Earthquake Research and Analysis chapter 19. InTech. doi:10.5772/27508. Ripepe, M., Coltelli, M., Privitera, E., Gresta, S., Moretti, M., & Piccinini, D. (2001). Seismic and infrasonic evidences for an impulsive source of the shallow volcanic tremor at mt. etna, italy. Geophysical Research Letters, 28, 1071–1074. doi:10.1029/2000GL011391. Rouland, D., Legrand, D., Zhizhin, M., & Vergniolle, S. (2009). Automatic detection and discrimination of volcanic tremors and tectonic earthquakes: An application to ambrym volcano, vanuatu. Journal of Volcanology and Geothermal Research, 181, 196 – 206. doi:https://doi.org/10.1016/j.jvolgeores.2009.01.019. Sparks, R. S. J., Biggs, J., & Neuberg, J. W. (2012). Monitoring volcanoes. Science, 335, 1310–1311. doi:10.1126/science.1219485. Steinberg, G. S., & Steinberg, A. S. (1975). On possible causes of volcanic tremor. Journal of Geophysical Research, 80, 1600–1604. doi:10.1029/JB080i011p01600. Stern, C. R. (2004). Active andean volcanism: its geologic and tectonic setting. Revista geol´ogica de Chile, 31, 161–206. doi:10.4067/S0716-02082004000200001. Toh, K.-A., Kim, J., & Lee, S. (2008). Biometric scores fusion based on total error rate minimization. Pattern Recognition, 41, 1066 – 1082. doi:10.1016/j.patcog.2007.07.020. Turetsky, R. J., & Ellis, D. P. W. (2003). Ground-truth transcriptions of real music from force-aligned midi syntheses. In ISMIR. Unglert, K., Radi´c, V., & Jellinek, A. (2016). Principal component analysis vs. self-organizing maps combined with hierarchical clustering for pattern recognition in volcano seismic spectra. Journal of Volcanology and Geothermal Research, 320, 58 – 74. doi:10.1016/j.jvolgeores.2016.04.014. Wohletz, K., & Heiken, G. (1992). Volcanology and geothermal energy volume 432. University of California Press Berkeley. Yamasato, H. (1998). Nature of infrasonic pulse accompanying low frequency earthquake at unzen volcano, japan. Bulletin of the Volcanological Society of Japan, 43, 1–13. doi:10.18940/kazan.43.1_1. Yoma, N. B., Carrasco, J., & Molina, C. (2005). Bayes-based confidence measure in speech recognition. IEEE Signal Processing Letters, 12, 745–748. doi:10.1109/LSP.2005.856888.

AC

CE

PT E

D

MA

verification. Pattern Recognition Letters, 29, 957 – 966. doi:10.1016/j.patrec.2008.01.015. John, D. A., Sisson, T. W., Breit, G. N., Rye, R. O., & Vallance, J. W. (2008). Characteristics, extent and origin of hydrothermal alteration at mount rainier volcano, cascades arc, usa: Implications for debris-flow hazards and mineral deposits. Journal of Volcanology and Geothermal Research, 175, 289 – 314. doi:10.1016/j.jvolgeores.2008.04.004. Jones, J., Carniel, R., & Malone, S. (2012). Subband decomposition and reconstruction of continuous volcanic tremor. Journal of Volcanology and Geothermal Research, 213-214, 98 – 115. doi:10.1016/j.jvolgeores.2011.07.006. Jousset, P., Neuberg, J., & Sturton, S. (2003). Modelling the timedependent frequency content of low-frequency volcanic earthquakes. Journal of Volcanology and Geothermal Research, 128, 201 – 223. doi:10.1016/S0377-0273(03)00255-5. Julian, B. R. (). Volcanic tremor: Nonlinear excitation by fluid flow. Journal of Geophysical Research: Solid Earth, 99, 11859–11877. doi:10.1029/93JB03129. Kittler, J., Hatef, M., Duin, R. P. W., & Matas, J. (1998). On combining classifiers. IEEE Transactions on Pattern Analysis and Machine Intelligence, 20, 226–239. doi:10.1109/34.667881. Lahr, J., Chouet, B., Stephens, C., Power, J., & Page, R. (1994). Earthquake classification, location, and error analysis in a volcanic environment: implications for the magmatic system of the 1989–1990 eruptions at redoubt volcano, alaska. Journal of Volcanology and Geothermal Research, 62, 137 – 151. doi:10.1016/0377-0273(94)90031-0. Langer, H., Falsaperla, S., Masotti, M., Campanini, R., Spampinato, S., & Messina, A. (2009). Synopsis of supervised and unsupervised pattern classification techniques applied to volcanic tremor data at mt etna, italy. Geophysical Journal International, 178, 1132–1144. doi:10.1111/j.1365-246X.2009.04179.x. Langer, H., Falsaperla, S., Messina, A., Spampinato, S., & Behncke, B. (2011). Detecting imminent eruptive activity at mt etna, italy, in 2007–2008 through pattern classification of volcanic tremor data. Journal of Volcanology and Geothermal Research, 200, 1 – 17. doi:10.1016/j.jvolgeores.2010.11.019. Lara, L., Orozco, G., Amigo, A., & Silva, C. (2011). Peligros volc´anicos de chile. Carta Geol´ogica de Chile, Serie Geolog´ıa Ambiental, 13, 34. Lara-Cueva, R. A., Ben´ıtez, D. S., Carrera, E. V., Ruiz, M., & Rojo´ Alvarez, J. L. (2016). Automatic recognition of long period events from volcano tectonic earthquakes at cotopaxi volcano. IEEE Transactions on Geoscience and Remote Sensing, 54, 5247–5257. doi:10.1109/TGRS.2016.2559440. Latter, J. H. (1981). Tsunamis of volcanic origin: Summary of causes, with particular reference to krakatoa, 1883. Bulletin Volcanologique, 44, 467–490. doi:10.1007/BF02600578. Lu, W.-k., & Zhang, Q. (2009). Deconvolutive short-time fourier transform spectrogram. IEEE Signal Processing Letters, 16, 576– 579. doi:10.1109/LSP.2009.2020887. Masotti, M., Falsaperla, S., Langer, H., Spampinato, S., & Campanini, R. (2006). Application of support vector machine to the classification of volcanic tremor at etna, italy. Geophysical Research Letters, 33. doi:10.1029/2006GL027441. McNutt, S. R. (1992). Volcanic tremor. Encyclopedia of earth system science, 4, 417–425. Molina, I., Kumagai, H., & Yepes, H. (2004). Resonances of a volcanic conduit triggered by repetitive injections of an ash-laden gas. Geophysical research letters, 31. doi:10.1029/2003GL018934. Moran, S. C., Malone, S. D., Qamar, A. I., Thelen, W. A., Wright, A. K., & Caplan-Auerbach, J. (2008). Seismicity associated with renewed dome building at mount st. helens, 2004–2005. In A Volcano Rekindled: The Renewed Eruption of Mount St. Helens,

17

ACCEPTED MANUSCRIPT

Highlights • Pattern recognition was applied to identify events of Llaima volcano, Chile. • Development of spectro-temporal curves based on expert knowledge.

RI PT

• Development of spectro-temporal features extracted from the spectrogram, based on high energy components and frequency bandwidth information.

AC

CE

PT E

D

MA

NU

SC

• Improvement at classification of seismic events including spectro-temporal information.

18

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9