Application of Subspace-Based Detection Algorithm to Infrasound Signals in Volcanic Areas

Abstract

Infrasonic signals investigation plays a fundamental role for both volcano monitoring purpose and the study of the explosion dynamics. Proper and reliable detection of weak signals is a critical issue in active volcano monitoring. In particular, in volcanic acoustics, it has direct consequences in pinpointing the real number of generated events (amplitude transients), especially when they exhibit low amplitude, are close in time to each other, and/or multiple sources exist. To accomplish this task, several algorithms have been proposed in literature; in particular, to overcome limitations of classical approaches such as short-time average/long-time average and cross-correlation detector, in this paper a subspace-based detection technique has been implemented. Results obtained by applying subspace detector on real infrasound data highlight that this method allows sensitive detection of lower energy events. This method is based on a projection of a sliding window of signal buffer onto a signal subspace that spans a collection of reference signals, representing similar waveforms from a particular infrasound source. A critical point is related to subspace design. Here, an empirical procedure has been applied to build the signal subspace from a set of reference waveforms (templates). In addition, in order to determine detectors parameters, such as subspace dimension and detection threshold, even in presence of overlapped noise such as infrasonic tremor, a statistical analysis of noise has been carried out. Finally, the subspace detector reliability and performance, have been assessed by performing a comparison among subspace approach, cross-correlation detector and short-time average/long-time average detector. The obtained confusion matrix and extrapolated performance indices have demonstrated the potentiality, the advantages and drawbacks of the subspace method in tracking volcanic activity producing infrasound events. This method revealed to be a good compromise in detecting low-energy and very close in time events recorded during Strombolian activity.

Introduction

Amplitude transient detection plays a fundamental role in volcano monitoring, allowing counting amplitude transients, identifying amplitude and occurrence rate variations. Besides, it is an essential step to localize seismic sources and their possible migration, which could be related to changes in volcano state and dynamics.

Classical methods for signal detection in seismology and volcano-seismology are grouped into two main categories: energy and correlation detectors. The former, called incoherent energy detectors, include algorithms searching for signals which are not or poorly known, such as STA/LTA (short-time average/long-time average) algorithm (Allen, 1978; Trnkoczy, 2012). These techniques, routinely used in volcano seismology, need no data pre-processing, beside the filtering applied to identify the desired signals, such as volcano-tectonic (VT) earthquakes or long period (LP) events. This approach suffers from high rate of false alarms or even of missed detections, due mainly to the background noise strongly affecting the reliability of this technique. This is especially true on active volcanoes with continuous volcanic tremor, which could dramatically reduce the signal (intended as the amplitude transients) to noise ratio. The latter group, correlation detectors, consists of algorithms based on the cross-correlation between a known waveform and the continuously recorded signal. These algorithms are very sensitive, give low false alarm rate but have the disadvantage of being able to detect only signals which are very similar to the template waveform, which in turn needs to be well known (Withers et al., 1999; Gibbons et al., 2007). In volcano acoustics, similar techniques (e.g., Cannata et al., 2013a; Cannata et al.,2013b; Thompson, 2015; Hotovec-Ellis and Jeffrines, 2016; Matoza et al., 2019a; Senobari et al., 2019), or methods making use of advanced signal processing techniques (Bueno et al., 2019), are implemented to identify and extract amplitude transients from the real-time streaming of signals, that characterize explosive or degassing activity. In particular, energy detectors, such as STA/LTA, are efficient algorithms when multiple infrasound sources are active (as at multi-vent volcanoes) and exhibit space-time variations, while correlation detectors are a powerful tool when we want to identify amplitude transients produced by a single and/or stable infrasound source in order to study its physical properties (Montalto et al., 2010; Sciotto et al., 2013; Cannata et al., 2013a; Hotovec-Ellis and Jeffrines, 2016; Yokoo et al., 2019).

Subspace-based detectors overcome the aforementioned limitation, in that they operate a comparison between the continuous signal and a set of reference waveforms hereafter called templates (Harris, 2006). One of the strong points of this method is the assumption on noise statistical features: it is supposed to be uncorrelated zero-mean gaussian noise. Signals acquired on active volcanoes generally are affected by band overlapped noise (e.g., volcanic tremor and volcanic infrasound tremor; e.g., Cannata et al., 2013b; Cannavò et al., 2019). In the light of it, without loss of generality, well known sources of noise, like infrasound tremor, can be preventively filtered. While in correlation detectors the waveform is a single template or a stacked waveform (Gibbons and Ringdal, 2006), in the subspace approach the designed set of templates is built by means of the Singular Value Decomposition (SVD) of a matrix whose columns are a variable number of templates. Subspace methods have been carried out mainly in seismology, where they have been applied for earthquakes tracking, especially in case of aftershock sequences (Harris and Dodge, 2011; McMahon et al., 2017), as well as to identify low-frequency earthquakes in non-volcanic tremor (Maceira et al., 2010).

Volcano acoustic plays a fundamental role for both monitoring purpose and the study of the explosion dynamics and revealed to be a reliable tool to characterize eruptive activity and shed a light into the shallow plumbing structure system at Etna (Cannata et al., 2013a; Sciotto et al., 2013, Spina et al., 2015; Cannavò et al., 2019; Sciotto et al., 2019). Proper detection of signal of interest is a crucial, and at the same time critical, issue in volcano seismology, in that it allows extracting and collecting amplitude transient waveforms (events), which are therefore analyzed to provide information about spectral content and source location. In particular, in volcano monitoring, events detection has direct consequences in pinpointing the real number of generated events and identifying amplitude and occurrence rate variations. This information can be of support to follow the explosive activity and to improve the assessment of volcanic activity. This is particularly true on Etna, where multiple open-conduit vents exist, whose activity often consists of persistent Strombolian explosions, producing low amplitude and very close in time infrasound events. In order to accomplish the detection task, several algorithms have been proposed in literature; in particular, to overcome limitations of classical approaches such as short-time average/long-time average and cross-correlation detector, in this paper a subspace-based detection technique has been implemented.

In this paper, we attempt to clear the way to the application of subspace detection method in volcano-acoustics, previously investigated in Sciotto et al., (2011) in a preliminary study, comparing its performance with correlation and STA/LTA detectors. In particular, we test this technique on signals recorded by the infrasound permanent network deployed at Mt. Etna, which represents an ideal dataset to lead tests on this matter. Indeed, infrasound activity at Mt. Etna is almost continuous, and is also characterized by both discrete amplitude transients and continuous tremor, produced by several summit craters and eruptive fractures often opening on the flanks of the summit cones. Moreover, the infrasound signals are generated by different source mechanisms related to explosive activity, such as Strombolian activity and lava fountaining, as well as to degassing phenomena (Cannata et al., 2013a; Sciotto et al., 2013; Spina et al., 2015; Sciotto et al., 2019). Therefore, in a multi-vent and open-conduit volcano such as Mt. Etna, where volcanic activity is almost persistent and prone to eruptive fracture opening, infrasound signal can consist of a superposition of signals from different time-variant and stationary infrasound sources. If on the one hand each infrasound source is repetitive, on the other hand it can undergo modifications in time. In these cases, correlation detector may fail in detection of the variation in infrasound waveforms caused by these factors. Subspace detector is supposed to accomplish these two tasks: high sensibility and high flexibility.

Data and Method

For the purpose of subspace-based detection implementation, theory of detection problem is first introduced. Successively subspace approach is explained and an empirical procedure used to build and design signal subspace described. Other two sub-sections are dedicated to discuss the statistical analysis of noise related to parameter estimation, both for subspace, correlation and STA/LTA detectors.

Subspace Detection Problem

Detectors usually implement a binary hypothesis test on the presence or absence of the signal of interest in a data observation window (Van Trees, 1968). In particular, the test is aimed to choose between the null hypothesis H₀, where noise only is present, and the alternative hypothesis H₁, where both the signal of interest and noise are present.where x[n] is the n-long window of continuous data, s is the signal of interest and η is the background noise assumed zero-mean Gaussian and temporally and spatially uncorrelated. In general, considering a multi-channel data acquisition, if N_t is the number of samples of observation window and N_c is the total number of data channel streams, the total number of samples N of a multiplexed data stream vector x[n] is:

In our framework, infrasound sensor acquires only one channel, so in Eq. 3.

Signal s in Eq. 2 is assumed deterministic and dependent on a vector of an unknown parameters a and expressed by unknown linear combination of a basis waveform:where U is a matrix of d unknown signals that represent the subspace bases. The subspace dimension d takes value from 1 to length of vector : . Without loss of generality, U can be made orthonormal:where I is matrix.

Under these assumptions, the probability densities function (pdf) for the recorded data under the null hypothesis (no events present) is:while, under the null hypothesis of (events present), pdf can be expressed as:

As formulated by Harris (2006), the detection rule is a likelihood ratio test comparing the probability that the observed data are due to signal and noise to the probability that they are due to noise alone:using Eqs. 6, 7, the likelihood ratio test expressed in Eq. 8 can be rewritten as a Generalized Likelihood Ratio Test (GLRT; Van Trees, 1968):Using natural logarithm, Eq. 9 can be rewritten as:

where when the pdfs are in the exponential family, is the least-squares estimate of the signal in the detection window:and , known as the subspace detection statistics, represents the ratio of the energy projected into the signal subspace U to the energy in the original data, and is given by:

The generalized likelihood ratio test (Eq. 9) detects an event of interest if the generalized log likelihood ratio (Eq. 10) exceeds a certain threshold :

Considering the subspace detection statistics , an event is detected if:where is the threshold for the subspace and needs to be defined. In order to apply subspace detector based on Eq. 14, the first step is the construction of the signal subspace U, starting from the template matrix.

This matrix has peculiar characteristics, which are described in the Template Matrix, and consists of templates, representing previously observed events of interest, and is a fundamental tool for building the subspace. The number and type of templates needed to build the matrix depends on detector design. Signal subspace is the core of the algorithm, since it is the vector subspace used to represent the reference templates to be found into the signal. In order to extract orthonormal bases, Singular Value Decomposition (SVD) is applied to the template matrix, and then the dimension of the subspace is chosen. The dimension determines the amount of energy that the subspace is able to capture. Once the SVD is applied, d singular values are used to build the subspace (Eqs. 4, 11). A few approaches have been implemented in literature to set this parameter, aiming to gain a compromise between detecting weak and less represented events (characterized by waveforms quite different from the reference templates) and having low false alarm or loss of significant events. In this paper, following Harris (2006), we selected the dimension parameter by means of an empirical approach making use of two different graphs as explained in Subspace Design.

Regarding the definition of the threshold , Harris (2006) studied the distribution of statistics and derived the threshold using the Neyman-Pearson criterion (Van Trees, 1968). Under this criterion, the subspace dimension d is firstly determined by maximizing the probability of detection for a fixed false alarm rate using the following equations:where is evaluated from the cumulative central F distribution with d and N-d degrees of freedom under the null hypothesis and is expressed in terms of the cumulative doubly non-central F distribution (Mudholkar et al., 1976) with the same degrees of freedom, is the non-centrality parameter for the numerator, is the non-centrality parameter for the denominator, f^¯c is the average fraction of energy for all design set events, N is the embedding space dimension, and (N-d) is the dimension of the orthogonal complement of the signal subspace; finally SNR is the signal-to-noize ratio in the detection window.

Dataset

In order to design a dataset for the analysis, we chose a 1-h-long time interval (13:30–14:30 of May 30, 2019) of infrasound continuous signal recorded at EMFO station. This station belongs to the Infrasound Permanent Network run by Istituto Nazionale di Geofisica e Vulcanologia (INGV), is equipped with a GRAS 40AN microphone with a flat response at a sensitivity of 50 mV/Pa in the frequency range of 0.3–20,000 Hz and sampling rate of 50 Hz, and is located about 8 km far from Etna summit craters and about seven from the eruptive fracture (Figure 1). This station, together to ESLN (which is deployed at about the same distance from the summit area), was the only station able to record the explosive activity, is one of the less noisy stations among the permanent network, and, if compared with summit stations (located at higher altitude) is less affected by wind noise that can hide weak amplitude transients.

FIGURE 1

Tarquini et al., 2007) with the location of the infrasonic station used in this work (green triangle EMFO), summit crater acronyms (VOR, Voragine; BN, Bocca Nuova; NEC, North-East Crater; SEC, South-East Crater; NSEC, New South-East Crater), infrasonic tremor source locations (red circles) and infrasonic event source locations (blue circles).

Signal buffer is characterized by infrasound events generated by an intense Strombolian activity that occurred at an eruptive fracture opened southeast of New Southeast Crater on the firsts hours of May 30, 2019 (here after NSEC, INGV-OE Internal Report, 2019). Here lava flows, ash emission, Strombolian and spattering activity took place. Explosive activity produced infrasound amplitude transients characterized by most of energy in the band 2.5–10 Hz, which are identifiable in the spectrogram from about 07:00 UTC, and with amplitude varying in a wide range (Figure 2). From about 14:00 UTC infrasonic activity at this fracture became more energetic, explosion generated infrasound events were more energetic and very close in time (Figure 2). Higher amplitude infrasonic events were detected and located by the real-time automatic system in force at INGV-OE (INGV-OE Internal Report, 2019) in correspondence of the eruptive fracture, as shown by blue circles in Figure 1.

FIGURE 2

In addition to infrasound events from the eruptive fracture, an overload continuous low frequency infrasonic tremor (∼0.6 Hz, Figure 2), whose source was located at Bocca Nuova crater (BN; red circles in Figure 1), was recorded. These characteristics make the dataset particularly useful to be used as test for an automatic detection algorithm. In particular, the chosen signal is suitable for verifying the subspace capability to detect the maximum number of infrasound events, especially of low amplitude ones, and to compare its performance with other detection algorithms. Furthermore, it allows us to verify this triggering technique in presence of noise, which is represented by the overlying low frequency infrasonic tremor.

Subspace Algorithm Implementation

Three buffers of infrasound signal were selected for subspace method application (Figure 2): 1) a 1 h-long time interval of signal consisting of background noise, and with no infrasound events, recorded during the same day of the dataset of analysis, to carry out statistical parameter estimation (00:00–01:00 of May 30, 2019; all times are in GMT); 2) a 3 h-long time interval of signal characterized by infrasound events of interest, for waveform templates selection (12:00–15:00 of May 30, 2019); and 3) 1 h-long time interval to test the subspace detector and searching for events of interest (13:30–14:30 of May 30, 2019).

We performed a first test by detrending and filtering each signal window between 1 and 10 Hz, to get rid of noise such as wind and low frequency tremor generated by a second infrasound source (Figures 1, 2), and set to zero mean and unit variance. A second test was performed by filtering signal between 0.5 and 10 Hz, in order to include the low frequency infrasonic tremor, and verify its influence on the detection.

Template Matrix

Building the template matrix is the preparatory step for subspace design. The template matrix is thought to consist of events of interest we are searching into the continuous signal. These can be manually selected, or, for a more robust procedure, waveforms can be automatically detected by a trigger algorithm (Maceira et al., 2010; Song et al., 2014). We made use of this last approach, and first triggered the events by means of STA/LTA energy detector. Secondly, we applied waveform cross-correlation, choosing an appropriate threshold, and selected the first event of each family, related to the infrasound source of interest (Figure 3A) in which the events were grouped. Once extracted, waveforms were aligned (Figure 3B); the algorithm has been designed to allow the operator to choose the alignment method. Waveforms can be aligned based on maximum or minimum amplitude value, or by means of manual picking. Successively, they were placed as columns in the template matrix.

FIGURE 3

Subspace Design

Signal subspace (U in Eq. 4) is the vector subspace used to represent the reference templates we want to find into continuous signal. The SVD provides the singular values allowing to build the subspace of the signal, that is a low-dimension representation of signal. Meaning of the dimension of subspace relies in the amount of energy that it is able to capture, and hence in the degree of waveform variation the algorithm is capable to detect. A few approaches have been implemented in literature (e.g., Harris, 2006; Song et al., 2014) to set this parameter, aiming to gain a compromise between detecting weak and less represented events (that is events having waveform quite different from reference template) and having low false alarm rate and possible loss of significant events. In this paper, following Harris (2006), we selected the dimension parameter by means of an empirical approach, making use of two different graphs. First, the fractional energy captured for each event is calculated:where is the fraction of energy captured by the ith template and is the ith unknown parameter of the coefficient matrix (Eq. 4, see Song et al., 2014 for further details). Figure 4A shows f_c of each template in function of the dimension of representation. The second plot is built by calculating the difference between the average captured energy in function of the dimension (Figure 4C):where is the average fraction of energy captured by each of the d templates. Values of and determined for our dataset of analysis are plotted in Figure 4. In particular, in Figure 4A the dimension d, instead of maximizing P_D value (Eq. 16), is graphically determined as the lower value so that all curves (or the average curve) lay above the assumed percentage of energy we want to capture (e.g., 80 or 90%). As an alternative, in Figure 4C an adequate dimension is the value beyond which the amount of energy increase is negligible. Once the dimension is selected, the subspace is built.

FIGURE 4

The sufficient statistic for the subspace can be now calculated by implementing Eqs. 11, 12, by means of windows of signal subspace sliding against the continuous signal, and compared against the threshold (Eqs. 13, 14) to declare if an event is present.

In subspace detector algorithm, signal in a detection window is projected into a subspace spanned by the d columns of the subspace representation. The statistic is therefore the ratio of the squared norm of the projected vector to the squared norm of the original data vector (Eq. 12). It ranges between 0 and 1 and is a measure of the linear dependence between the signal and the orthonormal bases constituting the signal subspace. Every time the sufficient statistic exceeds the given threshold () a detection is declared.

Threshold Setting

The choice of the threshold is always a compromise between an aggressive value, with the highest number of detections, even of less energetic events and leading to a high number of false detections, and a conservative value, when we want a minimum number of false detections at the cost of less true detections. Usually, threshold choice is based on the operator background experience and signal characteristics, that makes its value pretty subjective. In this paper, we tried to derive an empirical threshold value based on the data statistics and on the nature of the problem, e.g., waveform of events to be identified into incoming signals. In order to objectively compare the performance of subspace detector, we determined threshold of the three different applied triggering algorithms (subspace detector, correlation detector and STA/LTA) by means of the same approach.

Following Harris (2006) and Song et al. (2014), we implemented the Neyman-Pearson decision criterion (Van Trees, 1968). In the criterion, the threshold γ is derived from the false alarm rate with Eq. 15. In order to obtain the detection threshold, a few parameters need to be determined: 1) the false alarm probability (Eq. 15), 2) subspace dimension and 3) N. Regarding this latter, we should discuss about noise. Indeed, noise in the detection windows is assumed to be statistically uncorrelated. As Wiechecki-Vergara et al. (2001) point out, the effective dimension of the embedding space can be significantly lower than N if the data are filtered prior to detection. Noise could be correlated and could reduce the effective dimension of the embedding space even if data are not filtered. As Song et al. (2014) suggested, we applied the correction for the influence of the correlated noise, and estimated the effective embedding space dimension of detection windows used in subspace/correlation, and in STA/LTA detectors. According to Wiechecki-Vergara et al. (2001), the effective dimension of the embedding space is related to the variance of the sample correlation coefficient between noise data and event signal .

In particular, once calculated the cross-correlation values using the specific window length N, the variance is obtained and the effective embedding space of the respective detection window (subspace, correlation, STA/LTA) are calculated by means of:in the light of it, Eq. 15 can be rewritten as:

Hence, simple correlation detector can be written as:where is the master event data, is data to be detected. A comparison between Eq. 12 and Eq. 22 shows that the correlation coefficient is equivalent to the square root of the subspace detection statistics with a signal subspace dimension of d = 1. Here cross-correlation threshold is and both detectors have a false alarm rate:

After obtaining for subspace, next step is related to estimation of false alarm probability P_F by means of Eq. 23.

In order to accomplish this task, threshold is obtained by cross-correlating each template with noise data by means of Eq. 19 (see Song et al., 2014 for further details). Cross-correlation value distribution is then plotted and the detection threshold is set aiming to obtain the minimum number of false detection (Figure 5). In this paper, we empirically estimated the correlator threshold as the value corresponding to a chosen percentile of the distribution.

FIGURE 5

Once the cross-correlation threshold and the effective embedding space of detection window are known, false alarm probability can be estimated by inverting Eq. 23.

At this stage of the processing, we own all parameters needed to derive γ from Eq. 21, that is: 1) the false alarm probability P_F (Eq. 23), 2) the effective dimension of the embedding space , and 3) the subspace dimension (whose method of derivation is exposed in following sub-section).

With the aim of evaluating the advantages/effectiveness of the subspace-based detector, we make a comparison with the performance of STA/LTA and the simple correlator trigger algorithms. As regards the former, the detection statistic is calculated:where are data to be scanned, and the detection problem is formalized as:The STA/LTA threshold was determined by means of γ_r the same approach implemented for subspace, once computed STA, LTA (obtained by Eqs. 19, 20) and P_F (obtained by Eq. 23), solving the following equation:

Regarding the STA/LTA detector, based on event spectral content, we empirically set STA and LTA window length equal to 3 and 25 times the dominant period (Withers et al., 1998) respectively, corresponding with the lower frequency characterizing infrasound events.

Correlation detector was instead implemented as a subspace detector in which d = 1, by means of Eq. 22. Detection window length chosen for subspace and correlator scan was set to 4 times the dominant period to include the entire waveform (Figure 3B).

Results

Subspace-based algorithm, as well as correlation detector and STA/LTA, were applied to the dataset of 30/May/2019 (13:30–14:30), which consists of the infrasound signal recorded by EMFO station, and is characterized by infrasound events located at the eruptive fracture and infrasonic tremor located at BN (Figures 1, 2). Concerning the subspace method, we used, as event templates, waveforms extracted from signal recorded at the same station (Figure 3A), by means of the approach described in Template Matrix, setting a cross-correlation threshold equal to 0.6. Once template matrix has been designed, waveforms were cut in 62 points-long windows, filtered, normalized and aligned by their positive peaks (Figure 3B). A first test was carried out by filtering signal in the frequency band 1–10 Hz. This frequency band was chosen with the aim of filtering out the correlated low frequency tremor.

We performed two computations, by using 99.9 and 99.99, as percentile for statistic threshold estimation (Eqs. 22, 23, Figure 6), and built the subspace by using a dimension of representation equal to 4 (Figures 4A,C). The estimated thresholds and other setting parameters are reported in Table 1.

FIGURE 6

One subspace detection statistic c[n] value is calculated in each detection window, with a sliding step of three points, we obtain N_tot/3 c[n], where N_tot is the signal buffer length. In Figure 6, the recorded signal and the sufficient statistics of subspace detector (c[n]) above the estimated threshold are shown.

With the aim of avoiding more values of detection statistics for each triggered event, we extrapolated one detection in a 1 s long window. For the purposes of the comparison among subspace, correlator and STA/LTA performance, the detections computed by the aforementioned algorithms are overlapped to the continuous signal in Figure 7A. Histograms, obtained by counting detections in 2 min long windows for all algorithms, are shown in Figure 7B.

FIGURE 7

Plots of the occurrence rates show that subspace succeeds in detecting a higher number of amplitude transients, especially if compared with the STA/LTA triggered events (Figure 7B).

Results reveal the capability of the subspace-based algorithm to detect infrasound events of lower amplitude, while STA/LTA algorithm is able only to detect high amplitude transients. Figures 8A,B, which reports a zoom of continuous signal and the detected event positions, shows that subspace-based algorithm detected even more transients than correlator.

FIGURE 8

Successively, with the aim of highlighting the influence and the importance of parameter setting in this kind of approach, we run the detectors by setting the percentile equal to 99.99 instead of 99.9, used to obtain the respective detection threshold. The estimated parameters are reported in Table 1. By using a 99.99 percentile, detection threshold of the three detectors raises and, as a consequence, we observe a decrease of detection number as also demonstrated by the occurrence rates (Figure 9 and Table 1). Furthermore, subspace and correlator exhibit a similar trend of the event occurrence rate.

FIGURE 9

In the first case (percentile 99.9), the three detectors trigger more events with respect to the second one (percentile 99.99) (Figures 8, 9). In particular, the subspace succeeds in the detection of very low amplitude transients. Nevertheless, a lower threshold can imply the identification of a higher number of false detections. By an inspection of signal buffer reported in Figures 8, almost all the detections are real and not false positives.

Reliable estimation of the detection capability when correlated noise, due to low frequency infrasonic tremor, is overlapped, was tested using the three described methods. In the light of it, signal was previously filtered in the band 0.5–10 Hz. In such a way, signal to be scanned is characterized by both infrasound events, exhibiting a frequency content in the band 2.5–10 Hz, and the infrasonic tremor, whose spectral peak is at ∼0.6 Hz (Figure 2).

The resulting cross-correlation value distribution is shown in Figure 5B. A value of 99.9 as percentile was chosen, and the subspace was built by means of five SVD (Figures 4B,D). By introducing a correlated continuous tremor on the signal, the subspace and correlator thresholds result higher, due to the increase of variance, as expected (Table 1).

In this case, comparison with STA/LTA detector is inconsistent due to a not well defined statistics (Eqs. 19, 20, 26) used in the threshold γ_r computation. Despite this, simple correlation detector and subspace detector, from a practical point of view, exhibit robustness in the detection of waveforms transient as reported in Figure 10.

FIGURE 10

Comparison between Figures 7B, 10B, whose results were obtained by using similar false alarm probability (Table 1), reveals that, if the low frequency infrasonic tremor is overlapped to the signal and not filtered out, the number of detections by means of both subspace and correlation based methods is lower. Furthermore, even in this case, subspace method succeeds in event detection with respect to the correlator detector.

Discussion and Conclusion

In the present work, we applied a subspace-based trigger algorithm for the automatic detection of infrasound amplitude transients in volcanic area. A 1 h-long buffer of continuous infrasound signal characterized by amplitude transients related to Strombolian activity, taking place at an eruptive fracture, was analyzed (Figure 1). In order to test the feasibility and performance of this technique, we made a comparison among results obtained by implementing STA/LTA, correlation and subspace detectors. Several computations were carried out by using different setting parameters (Table 1).

Results highlight that subspace detector succeeds in detection of explosive activity related infrasound events. Indeed, Figures 7B and 9B, show that subspace detector turns out to detect a higher number of amplitude transients with respect to both correlation and STA/LTA trigger algorithms. This is particularly true for a lower false alarm probability (that is choosing 99.99 as percentile of cross-correlation distribution values in Eq. 22, Figure 5A). Indeed, the ratio between the total number of detections of subspace with respect to correlation and STA/LTA is higher than for a higher false alarm probability (that is choosing 99.9) (Table 1).

It is worth noting that subspace detector is able to detect even infrasound events of low amplitude (Figure 8), which, especially in presence of noise, are difficult to be triggered by means of an energy detector. Several runs were performed by using variable setting parameter values. These tests demonstrate that the choice of setting parameter plays a fundamental role in the outcomes of elaborations. In particular, an in-depth statistical analysis of noise needs to be carried out, in that distribution value of cross-correlation between template waveform and noise determines the threshold of the detector.

Succeeding in the detection of all amplitude transients allows to monitor the time variation of occurrence rate, and thus to follow the evolution of explosive activity. Indeed, energy-based trigger algorithm (as STA/LTA) often fails in detection of low amplitude transients, events with low signal to noise ratio and events too close in time to each other.

With the aim of performing a quantitative estimate of the subspace effectiveness, especially in terms of false alarms and missed detections, and in general, to validate the results, a visual inspection of the 1 h-long signal buffer (Figure 6) and its comparison with the results obtained by means of the three algorithms were carried out. In particular, concerning the performance assessment goal, we chose to analyze the results of the second test presented into the work (Figures 8C,D, 9), that is the one with percentile equal to 99.99 (which is the most conservative test), and inspected the true and false positive, and true and false negative events. We collected a dataset of 669 true events belonging to the same family (same source). Successively, by comparing waveforms with the detections (Figure 9A), a confusion matrix was calculated, for each of the three algorithms (Figure 11). The performance indices were then extrapolated from each confusion matrix (Table 2).

FIGURE 11

The validation results show that the subspace detector is generally characterized by the best indices among those derived from the confusion matrix (Figure 11; Table 2). In particular, it performs the best error rate (number of all incorrect predictions divided by the total number of dataset), accuracy (correct predictions divided by the total number of the dataset), and sensitivity, also called true positive rate, which is the number of correct positive predictions divided by the total number of positive, compared to the other two methods. In particular, STA/LTA gives the poorest results, as expected (Table 2). Nevertheless, it has to be highlighted that the correlator and STA/LTA show the best precision, which quantifies the number of correct positive predictions divided by the total number of the positive prediction, specificity (representing the true negative rate) and false positive rate, due to the lower false positive values and higher true negative values (Figure 11; Table 2). Finally, F-Score, providing a single score balancing both precision and sensitivity, shows the highest value in case of the subspace detector, due to the high rate of the actual triggered waveforms (Figure 11; Table 2). In the light of this, the subspace detector method can be considered a good compromise in recognizing low-amplitude waveforms, particularly useful in tracking volcanic activity producing low-energy and very close in time events.

As regards the better performance of subspace over correlation detector, in terms of error rate, accuracy, and sensitivity, this can be ascribed to the events waveform variability. Indeed, subspace method allows making a comparison between the continuous signal and a set of waveforms of interest (templates) and the linear combination among them, instead of a single one or multiple in few cases (Hotovec-Ellis and Jeffrines, 2016; Senobari et al., 2019), as the correlation detector does. This makes the subspace detector particularly attracting and suited to detect infrasound events undergoing slight modifications in waveforms due for example to geometrical characteristic variations (e.g., vent/crater enlargement) (e.g., Cannata et al., 2011; Fee et al., 2017; Matoza et al., 2019b).

Automatic detection of seismo-volcanic events (seismic and infrasound) in monitoring framework and/or analyses of huge dataset is successfully currently carried out by means of energy detectors, correlation-based algorithms or methods making use of advanced signal processing techniques (e.g., Cannata et al., 2013a, Cannata et al., 2013b; Thompson, 2015; Hotovec-Ellis and Jeffrines 2016; Bueno et al., 2019; Matoza et al., 2019b; Senobari et al., 2019). Among more recent developed trigger algorithms, for example, VINEDA, designed for infrasound event detection (Bueno et al., 2019), parses the original signal into a characteristic function, whose amplitude is proportional to the sharpness of the original explosion onset. This latter method is very useful to increase signal to noise ratio, much better than the STA/LTA approach, and is able to detect even low-amplitude transients. Another recent method, REDPy (Hotovec-Ellis and Jeffrines, 2016), designed for earthquake detection, consist in performing a STA/LTA, storing the triggered events, and then cross-correlating waveforms to find similar events. REDPy is an effective tool in order to discover events belonging to the same family and offers the advantages and the disadvantages of both STA/LTA and cross-correlation. In the framework of our goal, that is to identify and extracting waveforms similar to each other (hence, belonging to the same source), subspace detector exhibits the advantages of being able to detect low-amplitude transients associated to a specific source, as well as transients showing slight waveform modifications. The latter feature is related to the use of SVD and the subspace method, which consists of the scanning of the signal buffer with all the linear combinations among basis vectors derived from the template matrix (previously designed). An additional feature of the subspace method lies in its reliability in the detection of events that are close in time to each other.

Furthermore, in the proposed subspace-based detection method, algorithm parameters are automatically tuned by implementing a statistical analysis of the background noise.

In general, as reported in Song et al., 2014, the limitation of the subspace detector is the complexity and relatively large computation cost in building the signal subspace. Nevertheless, the aim of our work focused on the statistical analysis of noise in order to optimize the quality of the trigger, in terms of low-amplitude event number.

In future developments, the subspace detector method has the potentiality to be adapted to multi-station framework, considering multiple instances or by multiplexing the data (e.g., Song et al., 2014; McMahon et al., 2017). In the perspective of a real-time implementation, a detector should scan a 2-min-long buffer. Actually, using simple Matlab scripts, it takes 0.4 s on average to process 2 min of infrasound data, seeming suitable for real-time processing. Future work will be dedicated to optimization and code compilation to build a module working in a real-time framework.

Findings of this work reveal the potentiality of subspace-based method in infrasound event detection. Advantages of this technique are particularly interesting in infrasound recorded in open-conduit volcanoes such Mt. Etna, where activity at summit craters often consists of persistent Strombolian explosions, producing infrasound events very close in time, that can take place at several vents, thus giving rise to multiple and time varying infrasound sources.

This first application of a subspace-based detection algorithm to infrasound signal proves that this is an efficient technique for identification and triggering of events in volcanic area.

References

• 1

  AllenR. V. (1978). Automatic earthquake recognition and timing from single traces. Bull. Seismol. Soc. Am.68 (5), 1521–1532.

  • Google Scholar
  • View reference in article

• 2

  BuenoA.Diaz-MorenoA.ÁlvarezI.De la TorreA.LambO. D.ZuccarelloL.et al (2019). VINEDA—volcanic INfrasound explosions detector algorithm. Front. Earth Sci.7, 335. 10.3389/feart.2019.00335

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  CannataA.Di GraziaAliottaG. M.CassisiC.MontaltoP.PatanèD. (2013a). Monitoring seismo-volcanic and infrasonic signals at volcanoes: Mt. Etna case study. Pure Appl. Geophys.170 (11), 1751–1771. 10.1007/s00024-012-0634-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  CannataA.MontaltoP.PatanèD. (2013b). Joint analysis of infrasound and seismic signals by cross wavelet transform: detection of Mt. Etna explosive activity. Nat. Hazards Earth Syst. Sci.13, 1669–1677. 10.5194/nhess-13-1669-2013

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CannataA.SciottoM.SpampinatoL.SpinaL. (2011). Insights into explosive activity at eruptive ssure closely-spaced vents by infrasound signals: example of Mt. Etna 2008 eruption. J. Volcanol. Geoth. Res.208, 1–11.

  • Google Scholar
  • View reference in article

• 6

  CannavòF.SciottoM.CannataA.Di GraziaG. (2019). An integrated geophysical approach to track magma intrusion: the 2018 christmas eve eruption at mount Etna. Geophys. Res. Lett.46, 8009. 10.1029/2019GL083120

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  FeeD.IzbekovP.KimK.YokooA.LopezT.PrataF.et al (2017). Eruption mass estimation using infrasound waveform inversion and ash and gas measurements: evaluation at Sakurajima Volcano, Japan. Earth Planet Sci. Lett.480. 10.1016/j.epsl.2017.09.043

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  GibbonsS. J.RingdalF. (2006). The detection of low magnitude seismic events using array-based waveform correlation. Geophys. J. Int.165, 149–166. 10.1111/j.1365-246X.2006.02865.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  GibbonsS. J.SørensenM. B.HarrisD. B.RingdalF. (2007). The detection and location of low magnitude earthquakes in northern Norway using multi‐channel waveform correlation at regional distances. Phys. Earth Planet. In.160, 285–309. 10.1016/j.pepi.2006.11.008

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  HarrisD. B. (2006). Subspace detectors: theory. Livermore, CA: Lawrence Livermore National Laboratory Internal Report UCRL-TR-222758

  • Google Scholar
  • View reference in article

• 11

  HarrisD. B.DodgeD. A. (2011). An autonomous system for grouping events in a developing aftershock sequence. Bull. Seismol. Soc. Am.101, 763–774. 10.1785/0120100103

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  Hotovec-EllisA. J.JeffriesC. (2016). Near real-time detection, clustering, and analysis of repeating earthquakes: application to mount st. Helens and redoubt volcanoes. Reno, Nevada: Seismological Society of America Annual Meeting

  • Google Scholar
  • View reference in article

• 13

  INGV-OE Internal Report (2019). 23/2019. BollettinoEtna20190604. Available at: www.ct.ingv.it.

  • Google Scholar
  • View reference in article

• 14

  MaceiraM.RoweC. A.BerozaG.AndersonD. (2010). Identification of low-frequency earthquakes in non-volcanic tremor using the subspace detector method. Geophys. Res. Lett.37, L06303. 10.1029/2009GL041876

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  MatozaR. S.Arciniega-CeballosA.SandersonR. W.Mendo-PerezG.Rosaldo-FuentesA.ChouetB. (2019a). High-Broadband seismoacoustic signature of vulcanian explosion at Popocatépetl volcano, Mexico. Geophys. Res. Lett.46 (1), 148–157. 10.1029/2018GL080802

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  MatozaR. S.FeeD.GreenD.MialleP. (2019b). Volcano infrasound and the international monitoring system: challenges in middle atmosphere dynamics and societal benefits Book: infrasound monitoring for atmospheric studies. Switzerland, Europe: Springer Nature. 10.1007/978-3-319-75140-5_33

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  McMahonN. D.AsterR. C.YeckW. L.McNamaraD. E.BenzH. M. (2017). Spatiotemporal evolution of the 2011 Prague, Oklahoma, aftershock sequence revealed using subspace detection and relocation. Geophys. Res. Lett.44, 7149–7158. 10.1002/2017GL072944

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  MontaltoP.CannataA.PriviteraE.GrestaS.NunnariG.PatanèD. (2010). Towards an automatic monitoring system of infrasonic events at Mt. Etna: strategies for source location and modeling. Pure Appl. Geophys.167, 1215–1231. 10.1007/s00024-010-0051-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  MudholkarG. S.ChaubeyY. P.Ching-ChuongL. (1976). Approximations for the doubly noncentral-F distribution. Commun. Stat. Theor. Methods5, 49–63. 10.1080/03610927608827331

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  SciottoM.CannataA.GrestaS.PriviteraE.SpinaL. (2013). Seismic and infrasound signals at Mt. Etna: modelling of North-east Crater conduit and its relation with the feeding system of the 2008-2009 eruption. J. Volcanol. Geoth. Res.254, 53–68. 10.1016/j.jvolgeores.2012.12.024

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  SciottoM.CannataA.PrestifilippoM.ScolloS.FeeD.PriviteraE. (2019). Unravelling the links between seismo-acoustic signals and eruptive parameters: Etna lava fountain case study. Sci. Rep.9, 16417. 10.1038/s41598-019-52576-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  SciottoM.RoweC. A.CannataA.ArrowsmithS.PriviteraE.GrestaS. (2011). Investigation of volcanic seismo-acoustic signals: applying subspace detection to lava fountain activity at Etna volcano. AGU–American Geophysical Union, Fall Meeting. San Francisco, CA.

  • Google Scholar
  • View reference in article

• 23

  SenobariN. S.FunningG. J.KeoghE.ZhuY.YehC. M.ZimmermanZ.et al (2019). Super‐efficient cross‐correlation (SEC‐C): a fast matched filtering code suitable for desktop computers. Seismol Res. Lett.90, 322–334. 10.1785/0220180122

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  SongF.WarpinskiN. R.Nafi ToksőzM.Sadi KuleliH. (2014). Full-waveform based microseismic event detection and signal enhancement: an application of the subspace approach. Geophys. Prospect.62, 1406–1431. 10.1111/1365-2478.12126

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  SpinaL.CannataA.PriviteraE.VergniolleS.FerlitoC.GrestaS.et al (2015). Insights into Mt. Etna’s shallow plumbing system from the analysis of infrasound signals. Pure Appl. Geophys.172, 473–490. 10.1007/s00024-014-0884-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  TarquiniS.IsolaI.FavalliM.BattistiniA. (2007). TINITALY, a digital elevation model of Italy with a 10 m-cell size. Ercolano, Italy: Istituto Nazionale di Geofisica e Vulcanologia (INGV).

  • Google Scholar
  • View reference in article

• 27

  ThompsonG. (2015). “Seismic monitoring of volcanoes,” in Encyclopedia of earthquake engineering. Editors BeerM.KougioumtzoglouI. A.PatelliE.AuS. K. (Berlin, Heidelberg: Springer). 10.1007/978-3-642-35344-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  TrnkoczyA. (2012). IASPEI New manual of seismological observatory practice 2 (NMSOP-2). Deutsches GeoForschungsZentrum GFZ, Understanding and parameter setting of STA/LTA trigger algorithm. Berlin, Germany: Helmholtz-Zentrum. 10.2312/GFZ.NMSOP-2_IS_8.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  Van TreesH. L. (1968). Detection, estimation and modulation theory. Hoboken, NJ: John Wiley & Sons.

  • Google Scholar
  • View reference in article

• 30

  Wiechecki-VergaraS.GrayH. L.WoodwardW. A. (2001). Tech. Rep. DTRA-TR-00–22, Statistical development in support of CTBT monitoring. Dallas, TX: Southern Methodist University.

  • Google Scholar
  • View reference in article

• 31

  WithersM.AsterR.YoungC. (1999). An automated local and regional seismic event detection and location system using waveform correlation. Bull. Seismol. Soc. Am.89, 657–669.

  • Google Scholar
  • View reference in article

• 32

  WithersM.AsterR.YoungC.BeirigerJ.HarrisM.MooreS.et al (1998). A comparison of select trigger algorithms for automated global seismic phase and event detection. Bull. Seismol. Soc. Am.88 (1), 95–106.

  • Google Scholar
  • View reference in article

• 33

  YokooA.IshiiK.OhkuraT.KimK. (2019). Monochromatic infrasound waves observed during the 2014–2015 eruption of Aso volcano, Japan. Earth Planets Space71, 1210.1186/s40623-019-0993-y

  • CrossRef
  • Google Scholar
  • View reference in article