1. Introduction
The automatic or semi-automatic techniques for the detection of seismic events in active volcanic areas are challenging issues for a more detailed understanding of the volcano dynamics. Signal detection methodologies have considerably increased in number and improved in performance in recent years. This is due to the availability of huge amount of data that are not directly manageable by analysts through a visual inspection of the signal waveform. Therefore, many efforts have been made during the last decades to implement very sensitive detection algorithms. Methodologies based on spectral energy in a specific frequency band and STA/LTA filter [
1], wavelet analysis [
2], signal cross-correlation techniques [
3] and neural networks [
4,
5] are the most notable cases. The development of new signal detection procedures can be of great help for monitoring and surveillance purposes especially when they are designed to identify small-amplitude events that occur in areas characterized by high background noise. This is the case in densely populated volcanic areas where the identification and the classification of low-amplitude seismic events becomes a very difficult task. High levels of background seismic noise and low-amplitude seismic events are, in fact, the typical impediments that reduce the ability to detect and analyze seismicity in urbanized volcanic areas with classical methodologies [
6]. Active volcanoes are characterized by a variety of seismic events that produce signals which differ in amplitude, waveform and spectral content. The most common are volcanic tremor, low-frequency (LF) earthquakes, long-period (LP) events, hybrid events, and volcano tectonic (VT) earthquakes [
7,
8,
9,
10]. In this work, we refine and test an automatic detection procedure to identify LF earthquakes in the continuous seismic data recorded at Mt. Vesuvius by the seismic monitoring network managed by INGV-Osservatorio Vesuviano.
The seismicity of Mt. Vesuvius is characterized by VT earthquakes of small magnitude (M
D < 3), mostly located at depths between 0 and 2 km below sea level (b.s.l.) and sometimes occurring as low-energy swarms. During the last decade, the average rate of this seismicity has been of about 900 VT earthquakes per year [
11]. LF events characterized by low amplitude and deeper localization than VT earthquakes (between 5 and 7 km depth) have been identified in the seismic signals, sometimes as single sporadic events and in other cases as volcanic tremor consisting of low-frequency swarms [
12]. The source of this LF seismicity was studied in another work [
13], where the authors hypothesized the involvement of fluids in the genesis of these deep earthquakes. Over the past 10 years, the use of high-dynamic instrumentation, broadband sensors, and seismic arrays for surveillance, monitoring and research purposes has permitted the collection of high-quality data appropriate for automatic analysis. This facilitated the development of effective machine-learning techniques for analysis and classification of seismic signals recorded in active volcanic areas. For instance, an unsupervised technique based on a neural network approach for the detection and classification of seismic signals recorded at Mt. Vesuvius is provided by Esposito et al. [
14]. Previously, Scarpetta et al. [
15] provided automatic classification of seismic signals for the same volcano through neural networks. They introduced a new method for feature extraction, based on the combination of a waveform parametrization in the time domain, and an LPC (linear prediction coding) algorithm that provides a compressed and robust representation of the data in the frequency domain. More recently, the authors of the present work have used the statistical moments in the frequency domain to characterize the seismicity of Vesuvius [
16].
Following this latter work, the detection procedure we provide in the present paper is based on the coherence and statistical moments of seismic signals and turns out to be more empirical than many modern approaches that are quite articulated and based almost entirely on neural networks. Our main objective is the detection of deep LF events occurring beneath Mt. Vesuvius at depths greater than the predominant VT seismicity [
12]. The aim of improving the analysis of recorded seismic signals is to identify and classify in a semi-automatic manner the occurrence of LF events, possible signals due to landslides and other phenomena such as thunderstorms, in order to distinguish them from seismic signals associated with events arising outside the volcano, such as regional earthquakes and transient signals of artificial origin. To this end, the proposed approach focuses on the analysis of the spectral content of seismic signals using two techniques jointly, one based on the spectral parameters analysis through the computation of statistical moments, and another based on the coherence analysis. The effective detection of low-amplitude seismic signals such as LF events and volcanic tremor is the first step toward the detailed analysis of their source and driving mechanisms, which is of fundamental importance for assessing the eruptive potential of volcanoes [
17]. In the work of Galluzzo et al. [
16], spectral parameters were used for earthquake classification. In the present paper, the innovative aspect lies in the application of statistical moments for detection purposes and the additional use of coherence as a multichannel technique. The newly designed procedure was tested on continuous seismic signals recorded at Mt. Vesuvius from 2019 to 2021. After a brief outline of used data and methodologies, the main results are illustrated and discussed in the final part of the paper.
2. Data and Methods
The seismic monitoring network operating at Mt. Vesuvius comprises a permanent network, whose signals are transmitted in real time to the monitoring center for surveillance purposes [
18], and an offline mobile network [
19], whose data are characterized by a better continuity of signals and are mainly used for monitoring and scientific research. The dataset used in this work consists of continuous signals recorded in the years 2019, 2020 and 2021 by stand-alone seismic stations. The latter consisted of three-channel MarsLite Lennartz digitizers equipped with Guralp CMG40T broadband seismometers (60s). For our analysis, we selected those with the lowest seismic background noise (
Figure 1), i.e., those installed at high altitude (BKSG, BKWG) and in less-populated areas (SVAG). Once a possible seismic event has been identified by our automatic procedure, we also take into account signals recorded by the other seismic stations available in the area (BKE, VBKN, etc.;
Figure 1) and catalogs of local and regional earthquakes. In the present work, two analysis methodologies operating in the frequency domain have been used jointly: (a) the spectral analysis obtained with the statistical moments and (b) the coherence analysis. A brief description of the two adopted techniques is given in the following two sections.
2.1. Spectral Parameters Analysis
The spectrum of a seismic signal can be parametrized by using the statistical moments of the power spectrum [
20]. The statistical moments
λn of the seismic power spectrum
G were evaluated using the relation, Equation (1):
where the index n corresponds to the n-th moment of the seismic power spectrum
G. The integration is evaluated in frequency domain in the range
f1–
f2. For a single spectrum, the first three statistical moments
λ0,
λ1 e
λ2 were used to evaluate the two spectral parameters, i.e., central frequency
Ω and shape factor
ẟ, following the relationships (Equations (2) and (3) respectively):
The central frequency
Ω is a measure of the frequency where the signal power is higher, while the shape factor
ẟ indicates the dispersion of the power spectral density around the central frequency [
21]. The latter takes values between 0 and 1, with higher values corresponding to larger bandwidths. The power spectrum
G and then the spectral parameters were calculated on adjacent 30 s signal windows for the continuous seismic signals. For any seismic station, the power spectrum parameters in Equations (2) and (3) were evaluated, respectively, for each component of ground motion, and then the average central frequency and the average shape factor among the three components were computed. In this way, for a given time window, we obtained a pair of values
Ω and
ẟ for each seismic station. A duration of 30 s corresponds approximately to an earthquake of M
D = 1.3. We considered this window size a good compromise between the minimum detectable earthquake and the calculation cost of the procedure used. In particular, statistical moments were evaluated through Equation (1) in the frequency band 1–40 Hz. The lower integration limit was set at 1 Hz to reduce the effects of low-frequency seismic noise, which is often quite high at a frequency below 1–2 Hz in broadband signals, while the upper integration limit is given by the upper frequency limit of the flat bandwidth of the analyzed data. The calculation of the statistical moments was carried out on the continuous signals recorded at the BKSG, BKWG and SVAG sites for the years 2019, 2020 and part of 2021. In order to highlight potentially interesting signals for the search of LF events, signal windows were selected according to the results obtained by Galluzzo et al. [
16]. Specifically, signal windows characterized by a central frequency
Ω and a shape factor
δ in the ranges 3 Hz
≤ Ω ≤ 6 Hz and 0.45 ≤
δ ≤ 0.65, respectively, were considered in the present study. With this procedure, we obtain independent values of
Ω and
δ for any analyzed seismic stations.
2.2. Coherence Analysis
The coherence
γij between two signals
i and
j provides a measure of the similarity of their respective spectra in a specific frequency band, Equation (4):
where
Fi(ω) and
Fj(ω) are the spectra of the signals i and j, respectively, and F
i*(ω) is the conjugate of F
i(ω). A smoothing on a chosen bandwidth is necessary to obtain values between 0 and 1, where higher values correspond to signals with very similar spectral characteristics [
22]. The coherence defined in Equation (4) can be extended to more than two signals through the computation of the cross-spectral matrix, and thus, it is a multi-station analysis method. In our work, we used the signals recorded by the three stations BKSG, BKWG, and SVAG, which were utilized for the computation of spectral moments. As described above, we adopted adjacent 30 s signal windows to apply the coherence analysis to continuous data. The calculation was performed with a smoothing parameter of 20 and focused on a frequency of 4 Hz, in agreement with the spectral features of low-frequency events observed at Vesuvius by La Rocca and Galluzzo [
12]. The result of coherence analysis provides one value that is the coherence of the seismic signals among the three sites. Since low-amplitude LF events are almost always characterized by emergent onset and unclear waveforms in the time domain, the coherence provides more robust information than other mathematical tools in the time domain (such as cross-correlation) because it is based on the spectral characteristics of the analyzed signals. The coherence was computed on the vertical component of ground motion, and more than 28 months of continuous data have been analyzed. The results show that most of signals, at the selected frequency, have a value of coherence less than 0.4. Moreover, we have observed that the signals with coherence ≥ 0.45 are characterized by spectral amplitude higher than the background noise. On the basis of this observation, we set as a coherence threshold the value of 0.45.
3. Results
The methodological approach followed in this work consisted of three steps: (a) first, the two methods (statistical moments and coherence) described in
Section 2 were applied separately; (b) the obtained results were inspected and the signal windows that simultaneously exceeded the predetermined thresholds for both methods were considered; (c) the seismic signals from the selected time windows were visually inspected and their timing was searched in the seismic catalogs available online [
23]; “
www.ov.ingv.it; www.emsc-csem.org/ (accessed on day December 21, 2022)” looking for local and/or regional earthquakes. A synthetic flow diagram of the proposed procedure is shown in
Figure 2. Following the calculation of the spectral parameters, all signal windows characterized by shape factor
δ and central frequency
Ω values typical of LF events (see
Section 2.1) were selected in order to identify such possible seismic events.
The analysis of the continuous signal showed that most of the analyzed windows are characterized by spectral parameters and coherence values typical of random seismic noise. In particular, with regard to spectral parameters, most of the signal windows show central frequency
Ω around 10 Hz and shape factors
δ greater than or equal to 0.6, as found for noise records by Galluzzo et al. [
16]. As an example, we consider the spectral parameters and coherence analysis performed for May 2, 2020 (
Figure 3 and
Figure 4).
Figure 3 illustrates the results obtained for the three stations, BKWG, BKSG and SVAG. Plots in the top row of the figure show the results of shape factor vs. central frequency (blue points) for the three sites. Most of the dots are characterized by a central frequency greater than 6 Hz (histograms in the second row plots) and by a shape factor mostly greater than 0.5 (histograms in the bottom row plots). It is noteworthy that the three sites have different spectral features from each other, particularly with regard to the bimodal distribution of the shape factor at BSKG. Different results were expected because the spectral features of the background seismic signal are locally modified by site effects. In order to identify the windows that potentially contain signals associated with LF events, we need to inspect the results in the time domain.
Figure 4 shows the two spectral parameters (central frequency and shape factor) for the BKWG station, the coherence and the average spectral amplitude vs. time evaluated for 2 May 2020. In this case, the three larger yellow symbols identify the windows that simultaneously have the features of LF events at the three stations and coherence values ≥ 0.45.
Figure 5,
Figure 6 and
Figure 7 show the seismic signals within the selected windows. Three local earthquakes are recognized in the seismograms of the first window (
Figure 5). On the three components of ground motion, the direct P and S waves are identifiable only after appropriate filtering, and the time difference Ts − Tp = 1.3 s can be estimated for the BKSG and BKWG stations. The spectral content features are consistent with LF spectral parametrization. In this case, the three events have been classified as deep LF earthquakes. The second selected window contains signals with spectral features similar to those of local LF events, but in this case, they are generated by a regional earthquake (20200502 12:51:05 UTC Mw6.6 Crete, Greece) (
Figure 6). Visual inspection and examination of seismic catalogues revealed that the selected signal did not correspond to a local LF event. The third window selected from the automatic detection is summarized in
Figure 7. Seismograms of some stations of the monitoring network are shown in
Figure 7a. There appears to be no signal with a signal-to-noise ratio clearly observable in the seismic records. However, after filtering the seismograms in the frequency band of interest (2–8 Hz), a low-amplitude seismic event is visible at all the stations of the monitoring network. This is the case of a very low amplitude event, with spectral characteristics typical of LF events and recorded by the whole seismic network. It is not possible to identify any particular phases in these seismograms due to the emergent onset and the low amplitude of the signals and, for these reasons, the source of this event cannot be localized. Considering the spectral and shape characteristics, this event has been classified as an LF volcanic tremor.
Further examples of LF events detected by our procedure are shown in
Figure 8,
Figure 9 and
Figure 10. The two events shown in
Figure 8 and
Figure 9 are recognizable by visual inspection of the seismograms, but those of
Figure 9 are not easily identifiable and classifiable through a simple visualization because of their very low signal-to-noise ratio. In
Figure 10, we show the spectrogram of LF event 202009071932. Most of the frequency content is included in the 2–6 Hz band for both the small initial event and the subsequent larger one. In both cases, the spectral and coherence analysis provides parameters fully compatible with LF earthquakes. Specifically, our analysis applied to the selected dataset (years 2019, 2020 and part of 2021) highlighted the presence of 80 seismic signals potentially classifiable as LF events. Visual inspection of the seismograms and catalogue examination allowed us to recognize at least 12 LF earthquakes, a few LF tremors and 62 events not clearly classifiable in the LF category because they have extremely small amplitudes, almost undistinguishable from the background noise even in the filtered signals. The 12 LF events discovered by our analysis and shown in
Table 1 do not appear in the catalog of Vesuvius seismicity.
4. Discussion and Conclusions
To automate the detection and classification of seismic events, methodologies based on neural networks or even the use of time domain templates are often used to identify seismic events in the recorded signals. In the present work, we proposed a detection methodology based on the parametrization of spectral shape and coherence of the seismic signals among selected stations. The procedure, tuned for low-frequency earthquakes, was applied to continuous signals recorded at Mt. Vesuvius over a 28-month time span. The results highlighted 80 signal windows with spectral features compatible with LF events. For about 12 of these events (
Table 1), those with the higher amplitude and signal-to-noise ratio, the filtered seismograms allowed the recognition of P-S wave pairs that permit a rough location and unambiguous classification of these events as deep LF earthquakes, very similar to those reported by La Rocca and Galluzzo [
12]. In contrast, very low amplitude events are visible with a low signal-to-noise ratio only in filtered signals recorded at the lowest noise sites. Although they have spectral characteristics compatible with LF events, we could not provide a reliable classification because the lack of identifiable P-S wave pairs prevents any possible location of the source. These events constitute the majority of the detected signals (more than 80%). The hypothesis of an endogenous source is very reasonable because these events are visible only at the sites in the central area of Vesuvius and not at the seismic stations around the volcano, at a higher distance from the crater. A better investigation of these seismic events would be feasible with higher quality data, such as dense array recordings [
24] or deep borehole recordings.
In summary, the results of our analysis show that:
- (1)
The joint use of spectral parameters and coherence applied to a few selected stations of the seismic monitoring network operating at Mt. Vesuvius, which are characterized by low background noise, allowed for the detection of signals classifiable as LF earthquakes or LF tremors.
- (2)
The detected LF events reveal an amount of deep LF seismicity at Vesuvius much larger than that reported in the catalog.
- (3)
The applied procedure also detects a large number of low-amplitude events hidden in the background noise that have the characteristics of LF events, but their effective classification requires further investigation and possibly higher quality seismic data.
- (4)
The overall outcomes suggest that the proposed procedure could be appropriate for real-time applications aimed at LF signal detection, but its suitability must be further investigated.
The analysis procedure applied in this work is strongly focused on the detection of LF events characterized by spectral features in narrow ranges, as observed for events detected in the past years at Vesuvius. It must be pointed out that the use of narrow ranges for tuning parameters can lead to failure to detect seismic events with spectral characteristics different from those observed previously. The application of the proposed automatic detection procedure to other volcanic areas is feasible, but it could require a revision of the search parameters to be effective.