In this section, we explain the technical side of the method A6-DFMSD in accordance with the flow diagram in Fig. 2. Later in the discussion section we expand the advantages and the capabilities of the method. A user manual for the detector A6-DFMSD is also provided in the supplementary file to explain more details for potential users [note: this manual will be open and provided after acceptance of the manuscript]. The main technical ideas of A6-DFMSD program are explained in the following subsections:
- Input Parameters
- Configuration of the Detection Model (CDM-step)
- Detection Field preparation (DF preparation)
- Single Station Detection (SSD-step)
- Multi Station Detection (MSD-step).
2.1 Input Parameters
A6-DFMSD is designed to detect the seismic signals which originate from a pre-defined target zone (e.g., a volcanic field) and are recorded by a local seismological network. For this reason, in addition to the three continuous seismic velocity records (Z: vertical, N: north-south, E: east-west), the following information is needed to establish a detection model:
(a) coordinates and codes of the seismological stations,
(b) P- and S-wave seismic velocity models,
(c) a center location, a radius length of the detection area, an upper and a lower depth range for determining the event locations. These data allow to define a cylindrical geometry of the target zone.
2.2 Configuration of the Detection Model (CDM-step)
In this step, the program A6-DFMSD uses the input information to establish a model and calculates some outputs which later are used in DF preparation and MSD steps. The CDM-step consists of the following sub-steps:
(a) Defining the target seismic zone: synthetic source positions are distributed on the circles that represent the upper and lower boundaries of the cylindric target zone. Figure 3A provides an exemplary top view of distributed synthetic source positions which are the same as in the application of this study (see section 3). It presents a circle with sources inside a radius of 25 km across the EEVF.
(b) Calculating the P- and S-wave travel times (Tp and Ts, respectively) for each individual combination of the synthetic sources and the recording stations: for this calculation the CDM-step considers the local seismic P- and S-wave velocity models, the synthetic source positions, the positions of the seismological stations and Snell’s law for the ray path approximation.
(c) Obtaining the time intervals based on the potential travel times between the recording stations and possible source locations: these time limits are used in the following MSD-step when the seismic signals, which originate from the target zone, are discriminated from locations outside the target zone. In this part, the program calculates the relative difference in arrival time of the seismic phases between each two stations for each synthetic source. After this calculation, for each pair stations, the maximum and minimum values are taken. They are used for the MSD-step to limit the search for the arriving seismic phases for each pair of stations in a certain time interval. Figs. 3 B and C, visualize these time limits for the recording stations DEP02 and ABH, respectively.
(d) Determination of the four closest stations relative to each other station: for each station, the four closest stations are determined and saved. This information is used for the local coherency search in the MSD-step. Fig. 4 A shows how the stations are grouped for the test example in this study.
(e) Estimation of an optimum time window (called ”minStatisDur”) for the seismic trace windowing: minStatisDur is obtained by a statistical approach to estimate the typical half duration of the micro-earthquake signals which are supposed to occur inside the target zone and are recorded at the seismological network. This value is the geometric mean value minus the geometric standard deviation value of all the calculated Ts-Tp values, for which all the synthetic sources of the target model in the shallow depth and all the stations are considered. In the following step (section 2.3) minStatisDur is used for windowing the so-called energy mono bands (see section 2.3).
2.3 Detection Field Preparation (DF preparation)
In this part, the routine reads in the three components of the continuous ground motion velocity records and produces a new set of traces which are proportional to the recorded energy in one-Hertz wide frequency bands (Fig. 5). Each component of the original velocity records is decomposed into 29 frequency mono band traces (with a width of one Hertz) covering the frequency range 1–30 Hz (Figs. 5B, C and D). This frequency range is adjusted to our case of volcanic and tectonic earthquakes in the EEVF and may be changed for other purposes. The mono band traces are converted to energy-proportional (E*) mono bands by squaring the velocity (v) proportional data points to achieve a quantity which is proportional to the kinetic energy EK of the ground motion (EK = m· v2/2) ignoring the mass m and the coefficient. For each sample i of each mono band j the total energy* (E*j,i) is obtained by the summation:
$${E}_{j,i}^{*}={Z}_{j,i}^{2}+{N}_{j,i}^{2}+{E}_{j,i}^{2}$$
1
with Z, N and E the velocity-proportional seismic recording components (Fig. 5E). The 1-Hertz energy mono bands are integrated measures of the ground motion which increase coherent signals between the three recording components and which are especially helpful to search for frequency-dependent signals, e.g., such as LF magmatic tremor.
In the next step, the energy mono bands (E*j,i ) are windowed by a window length (routine minStatisDur) which is equal to the half duration of the typical micro-earthquake signals recorded at the seismic network (Fig. 5F). ”minStatisDur” is obtained statistically during the CDM-step. As results of this procedure we obtain 29 windowed energy mono bands called the Detection Field (DF). The DF of each mono band is indexed by j (DFj).
2.4 Single Station Detection (SSD-step)
In this step, signal anomalies are detected independently in different frequency bands called signal classes. In our example 11 signal classes are chosen. Each signal class is defined by the two values denoting the upper and lower corner frequency of the signal class in Hz. For each signal class a short-term amplitude (STA) to long-term average amplitude (LTA) ratio is determined with all DFs, which are within the upper and lower corner frequencies of the signal class. This leads to the detection criterion. A signal anomaly is detected, if all DFs within a signal class reach to the sample that reached amplitude values higher than the mean value of the three previous sample amplitudes plus a multiplication of 0.7 (an empirical constant) to their standard deviation. As example, in Fig. 5G a signal anomaly detection in signal class n = 6 is presented. The signal class n = 6 includes the detection fields of DF8, DF9, DF10, DF11, DF12, DF13, DF14, DF15, and DF16. In the example in Fig. 5G the red dot (at sample i) is found as a signal anomaly, because at all DFs the amplitude values at sample i have a higher value than the mean plus a multiplication with 0.7 to the median absolute deviation of the reference samples (shown as three green dots right before the sample i in the DFs). For each signal anomaly, which is detected in a signal class, 6 unique values are saved as labels of the signal anomaly in a file called ”SSD-file”. These six values are:
- the station code (implicitly, the latitude and longitude of the station location)
- the sample number of the detection point (corresponding to the detection time)
- the signal class
- a value called ”energy* smoothness factor” which represents the energy* distribution status at the point of detection in the corresponding signal class. This value is the standard deviation of the DF amplitudes in the signal class at the sample of detection divided by their mean value.
- a value called ”sigClasPower” which represents the amplitude level of energy* in the signal class at the detection point sample relative to the 3 previous samples just ahead the detection point sample. This value is equal to the result of the division of two values: the numerator is the stack value of the DF amplitudes of the signal class at the detection point minus the mean value of the three samples right before the samples at the detection point (in Fig. 5H this is equal to the length of the red line), and the denominator is the standard deviation of the values of the three samples right before the detection point sample at the stacked trace (in Fig. 5H this is equal to the length of the blue line).
The ”energy* smoothness factor (Λ)” and the ”sigClasPower (Γ)” in signal class (n) with the upper border (j2) and lower border (j1) at the detection sample (i) are obtained according to the equations (2) and (4), respectively:
$${{\Lambda }}_{n,i}=\frac{\sqrt{\frac{1}{{j}_{2}-{j}_{1}+1}\sum _{j={j}_{1}}^{{j}_{2}}{({DF}_{j,i}-{\mu }_{n,i})}^{2}}}{{\mu }_{n,i}}$$
2
\({\mu }_{n,i}=\frac{1}{{j}_{2}-{j}_{1}+1}\sum _{j={j}_{1}}^{{j}_{2}}{{DF}_{j,i}}^{ }\)
|
(3)
|
\({{\Gamma }}_{n,i}=\frac{{\mu }_{n,i}-{M}_{n,i}}{\frac{1}{{M}_{n,i}}\sqrt{\frac{1}{3}{\sum }_{i=i-3}^{i-1}{({\mu }_{n,i}-{M}_{n,i})}^{2}}}\)
|
(4)
|
\({M}_{n,i}=\frac{1}{3}\sum _{i=i-3}^{i-1}{{\mu }_{n,i}}^{ }\)
|
(5)
|
where: µn,i is the mean value of the DF amplitudes in the class (n) at the detection sample (i) and Mn,i is the mean value of the 3 samples right before the sample of detection (i) at the stacked trace in the signal class (n)
2.5 Multi Station Detection (MSD-step)
After the SSD-step is finished for all network stations, in the MSD-step the search for the coherent signal anomalies begins among the labelled signal anomalies at all recording stations. This search is done by applying a 6-dimensional floating-search frame (6-DFF). 6-DFF uses labels to bound its searching ranges and limits. This means that depending on what is saved as the label of each signal anomaly (the 6 unique values for each signal anomaly), 6-DFF sets the searching ranges accordingly. The searching ranges, which are used in this part regarding the place and time, are determined station-wise at the CDM-step. 6-DFF uses these ranges to limit the search of the coherent labels locally. In this way, each labelled signal anomaly is taken as a ”reference label” and a 6-DFF is established around it. Based on the reference label content:
- the coherency search in place is limited to the four closest stations (called ”orbit stations”) to the station with the reference label. This supports the principle that, regarding the local station network configuration and the direct wave propagation of the local events, at stations closer to the source the waves arrive earlier than at more distant ones. As a result of applying the local search, the number of false detections was remarkably reduced in our tests.
- the coherency search in time is limited to the certain time window around the time sample in the reference label (the time window borders are unique and differ from station to station depending on the relative position of the reference and orbit stations in relation with the position of the target zone). The time limits, which are applied here, are already determined in the CDM-step.
- the coherency search in a signal class is limited to the corresponding signal class of the reference label.
- the coherency search in sigClasPower (Eq. 4) is limited to the labels with a sigClasPower value more than 2.
- the coherency search is limited to the borders around the energy* smoothness factor of the reference label (Λ) which is called ”smoFacVici” and obtained according to the empirical Eq. (6):
$$smoFacVici\left(\varLambda \right)={1.6}^{-({\Lambda }+iniSmoFacVici)}$$
6
The value for iniSmoFacVici is provided in the control file of the program. Depending on the upper depth limit of the target region it can reach values larger than 2.7. In Eq. (6), when Λ has a low value, the search around it has a wider border limit and vice versa. This allows the detector to conduct a robust coherency search around the dominant frequency of the signal anomalies. This approach is supported by the fact that signal anomalies around their dominant frequencies have a lower value in their energy* smoothness factor relative to the remaining part of their spectra. Although this equation is obtained empirically (simply by varying the involved values and parameters in the equation and then observe their effect on the number of the true and false detections), a simple rule can be assigned for allocating a value to iniSmoFacVici. This rule is applied when the upper depth boundary of the target zone is considered to be at a deeper position (deeper than 3 km). In such a case, we can assume a weaker ray path effect on the frequency content of the recording signals, and consequently expecting a narrower coherency search vicinity for the energy* smoothness factor. Based on this assumption, we allocated 2.7, 2.9 and 3.1 to the iniSmoFacVici value when the upper boundary values of the target zones are respectively set to the 3 km, 30 km or 45 km depth. Figure 4B shows how the smoFacVici value changes according to the energy* smoothness factor of a reference label while detecting events in different layers of defined target zones.
- for each label, that is considered as a reference label, if two coherencies are found (three coherencies including the reference label) through the 6-DFF step, the routine accepts the event belonging to the reference label as a detected event which originated from the target zone. It is then written into the detection list. The detection list of the program includes the following information regarding each detection:
- the date and arrival time of the event signal in UTC
- the code of the stations with coherent labels for each detection
- the signal class of the reference label.
The two latter information helps users while manually checking the detection outputs to select the events based on their dominant frequency content and/or based on the stations in which the signal of the event appears.