# Multiband array detection and location of seismic sources recorded by dense seismic networks

Natalia Poiata, Claudio Satriano, Jean-Pierre Vilotte, Pascal Bernard, Kazushige Obara

► To cite this version:

Natalia Poiata, Claudio Satriano, Jean-Pierre Vilotte, Pascal Bernard, Kazushige Obara. Multiband array detection and location of seismic sources recorded by dense seismic networks. Geophysical Journal International, 2016, 205, pp.1548-1573. 10.1093/gji/ggw071 . insu-03581658

HAL Id: insu-03581658

https://insu.hal.science/insu-03581658

Submitted on 20 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Multiband array detection and location of seismic sources recorded by dense seismic networks

Natalia Poiata,<fnref n="1" /> Claudio Satriano,<fnref n="1" /> Jean-Pierre Vilotte,<fnref n="1" /> Pascal Bernard<fnref n="1" /> and Kazushige Obara<fnref n="2" />

Accepted 2016 February 17. Received 2016 February 17; in original form 2015 June 25

## SUMMARY

We present a new methodology for detection and space–time location of seismic sources based on multiscale, frequency-selective coherence of the wave field recorded by dense large-scale seismic networks and local antennas. The method is designed to enhance coherence of the signal statistical features across the array of sensors and consists of three steps: signal processing, space–time imaging, and detection and location. The first step provides, for each station, a simplified representation of seismic signal by extracting multiscale non-stationary statistical characteristics, through multiband higher-order statistics or envelopes. This signal processing scheme is designed to account for a priori unknown transients, potentially associated with a variety of sources (e.g. earthquakes, tremors), and to prepare data for a better performance in posterior steps. Following space–time imaging is carried through 3-D spatial mapping and summation of station-pair time-delay estimate functions. This step produces time-series of 3-D spatial images representing the likelihood that each pixel makes part of a source. Detection and location is performed in the final step by extracting the local maxima from the 3-D spatial images. We demonstrate the efficiency of the method in detecting and locating seismic sources associated with low signal-to-noise ratio on an example of the aftershock earthquake records from local stations of International Maule Aftershock Deployment in Central Chile. The performance and potential of the method to detect, locate and characterize the energy release associated with possibly mixed seismic radiation from earthquakes and low-frequency tectonic tremors is further tested on continuous data from southwestern Japan.

**Key words:** Time-series analysis; Earthquake source observations; Seismic monitoring and test-ban treaty verification; Subduction zone processes.

## 1 INTRODUCTION

A key challenge in both fundamental and applied earthquake seismology is to detect and characterize the seismic signature of Earth's deformation, and to identify and classify the associated seismic sources.

The rapid evolution of seismological instrumentation has led, in the last decades, to the progressive densification of observation and monitoring systems at all the scales, from micro-earthquake studies to large plate interface observatories. Modern seismic networks can today be regarded as arrays, or antennas, of increasing density and spatial coverage. These systems record time-continuous, high-sensitivity signals, progressively enlarging the seismic observation window, and making it possible to explore time and frequency scales at which seismic sources were essentially unknown. As a result, spectacular discoveries unveiled the complexity and variety of seismic sources associated with active deformations in different

tectonic and industrial exploitation contexts over a wide range of spatial and temporal scales. Among these breakthroughs is the detection of repeating micro-earthquakes on active fault surfaces (e.g. Nadeau & McEvilly 1999); the discovery of episodic tremor and slip (e.g. Obara 2002; Rogers & Dragert 2003; Nadeau & Dolenc 2005) and of low and very low frequency earthquakes (Ide et al. 2007; Ito et al. 2007; Beroza & Ide 2011); the observation of long-period, long-duration seismic events during hydraulic stimulation in shale reservoirs (Das & Zoback 2013a,b).

Today, seismological monitoring systems generate increasing volumes of data, challenging the current source detection and extraction methodologies. Thus, development and implementation of advanced automated and scalable methods continues to be a central issue. The main focus being on developing innovative schemes that can fully exploit multiscale coherence of the seismic signal—over wide range of wavelengths and frequencies—across dense arrays of stations and improve source detection levels and resolution in the

<footnote n="1">¹Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Univ Paris Diderot, UMR 7154 CNRS, F-75005 Paris, France. E-mail: poiata@ipgp.fr</footnote>

<footnote n="2">²Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan</footnote>
case of multiple sources exhibiting wide dynamic ranges and large spatial extents.

The most common approach to automated earthquake detection and location is based on a three-step procedure: (1) reduction of individual waveform records to series of arrival time picks, that is, phase picking; (2) association of time picks with specific seismic phases and a potential common source; (3) location of earthquake sources through traveltime inversion methods. Earthquake location methodologies based on traveltime inversions have been pioneered by Geiger (1912). These methods rely on local or global minimization schemes of a cost function generally expressed as a term linked to the data, that is, misfit between the observed and theoretical arrival times of multiple seismic phases, and a regularization term incorporating a-priori knowledge of the source (see Lomax et al. 2009, and references therein). A significant limitation of such detection and location approaches is that automated picking is often done on each seismogram individually, making little or no use of the coherency information between stations. In the case of multiple and possibly overlapping seismic events, like for example during foreshock and aftershock sequences or seismic swarms, phase identification and event association can quickly become a difficult task, leading to missed detections and/or reduced location resolution.

To address these shortcomings, alternative signal-based strategies for automated detection and location of seismic sources, which bypass the picking and association steps, have emerged in recent years (see Cesca & Grigoli 2015). The general idea of these methods is to convert seismic signals recorded at different stations into time-series of 3-D spatial images, through reverse propagation, delay-and-sum or migration techniques. Depending on the approach, each 3-D image can either represent a snapshot of the full seismic wave field, an energy function or a likelihood of source location. The locations of the seismic sources are then extracted as points in space and time where the reconstructed wave field refocuses, or as the local maxima of the energy or likelihood function (using peak search algorithms or Bayesian approaches).

A first class of methods belonging to these alternative strategies makes use of simulated reverse space-time propagation of the recorded signals to produce time-series of 3-D spatial images of the full seismic wave field (e.g. McMechan et al. 1985; Gajewski & Tessmer 2005; Larmat et al. 2006; Xuan & Sava 2010). These methods are computationally demanding and require accurate knowledge of the adjoint Green functions and the structure of the medium (e.g. Baker et al. 2005), which, in case of most passive earthquake monitoring systems, is rather problematic, and becomes the principal factor limiting the accuracy and robustness of the methods.

A second class of signal-based methods performs, as preliminary step, a signal transformation at the level of each individual record (by constructing characteristic functions (CFs) based on a particular defined transform) to enhance signal features for further spatial conversion and source extraction processing stages. The transformed signal is then converted into time-series of 3-D spatial images representing an energy or likelihood function, through delay-and-sum (or migration) techniques, using arrival time prediction combined with a grid search strategy. Early examples of this second class of methods involve CFs based on short-term over long-term average (STA/LTA) in the context of earthquake location at regional and local scales (Withers et al. 1999), and the Source Scanning Algorithm, initially developed by Kao & Shan (2004, 2007) in the context of tectonic tremors location, extended by Liao et al. (2012) to earthquake aftershocks.

Many automated phase pickers or signal-based detection and location methods follow the approach of Allen (1978, 1982) in defin-

ing the nonlinear transform of the seismic signal, called characteristic function (CF). Signal CFs typically rely on the identification of changes in energy, frequency content, polarization and/or other characteristics of the seismic signal, with respect to the background noise (see Withers et al. 1998 for a detailed review). The most frequently used class of CF is based on the ratio between the signal STA to its LTA (Freinberger 1962; Allen 1978, 1982). Recent works introduced and investigated various alternative CFs, pointing out that higher-order statistics (HOS), such as skewness and kurtosis, provide significantly better results than energy-based STA/LTA in the conditions of low signal-to-noise ratio (SNR), and intense seismic activity (e.g. Saragiotis et al. 2002; Gentili & Michelini 2006; Kuperkoch et al. 2010; Baillard et al. 2014). Different types of CFs, based on STA/LTA or second and fourth order statistical moments, have been recently developed and incorporated into the detection and location schemes applied to data sets from regions with intense seismic activity, such as induced micro-seismicity (Grigoli et al. 2013), and active volcanic areas (Drew et al. 2013; Langet et al. 2014). Most of these studies are based on an a-priori assumption of known seismic environment with particular types of the seismic sources and predominant frequencies of the signals (e.g. Kao & Shan 2004; Drew et al. 2013; Grigoli et al. 2014; Langet et al. 2014). Such a-priori information remains problematic in the case of complex seismic environments, with multiple sources spanning wide spatial scales and high dynamic ranges.

In this paper we present a scalable strategy and several methodological contributions that are a first step towards the development of an advanced automated multiscale framework for extraction (i.e. detection, location and classification) of multiple sources with a wide range of dynamic and spatial scales (e.g. earthquakes, tremors, low-frequency events), minimizing the amount of necessary a-priori information. The proposed strategy belongs to the class of signal-based methods and consists in a number of adaptive processing steps, designed to enhance coherent statistical features of the wave field recorded across the array of sensors. Through the sequence of these steps the raw signals are transformed into time-series of 3-D spatial images, representing the likelihood of each pixel to be the part of a seismic source.

We evaluate the performance of the proposed method on two test cases: earthquake records from the International Maule Aftershock Deployment (IMAD) in Central Chile after the 2010 Maule earthquake ($M_w$ 8.7; Beck et al. 2014); and several earthquakes and low-frequency tectonic tremors recorded in southwestern Japan by the dense high-sensitive Hi-net monitoring system (Okada et al. 2004; Obara et al. 2005) operated by the National Institute for Earth Science and Disaster Prevention (NIED).

## 2 METHODOLOGY

We introduce here the main steps of the proposed adaptive strategy for extraction of multiple seismic sources exhibiting wide dynamic ranges and spatial scales (Fig. 1a).

In the first step (Fig. 1b), the signal recorded at each station is transformed into CF, designed to extract and enhance multi-scale non-stationary statistical features, through multiband HOS or envelope-based structural components. This stage is designed to account for a priori unknown transient signals potentially associated with a variety of sources, and to prepare data for a better performance in posterior steps.

In the second step, station-pair time-delay estimate (TDE) functions are evaluated by calculating local cross-correlations (LCCs) of
<figure><img src="image_3.png" /><figcaption>Figure 1. Diagram of the detection and location scheme: (a) Overall flow of the method. (b) Step-by-step flow of the signal processing scheme transforming an input seismic record into broad-band characteristic function (CF).</figcaption></figure>

the transformed signal. The imaging functions are then constructed through combination of the station-pair TDE functions and their spatial mapping, according to theoretical phase time difference of arrivals (TDOA). This step produces time-series of 3-D spatial images representing the likelihood of each pixel being part of a source.

In the third step, space-time location of the source is evaluated using extraction algorithms, applied to the time-series of 3-D images.

## 2.1 Signal transformation and CF

Seismic signals are time sequences regularly sampled and for the most part non-stationary in time. Non-stationarity may be expressed in various ways: varying powers from one time instance to the next, changing correlation between neighbouring time samples, or even non-stationary statistics, that is, higher moments, of the time samples. When dealing with such signals, the most widely accepted

approach is to consider them as the realization of some stochastic process.

The initial step of our detection and location scheme is statistical signal analysis, aiming to construct a CF that extracts and enhances the main properties and signatures of the signal non-stationarity. We are going to deal with two main families of signals: short transients, possibly associated with impulsive, earthquake-like sources and emergent signals, characterized by slow variation of energy with time, such as those related to tremor-like sources. We therefore use two different signal transformations: higher-order statistics (HOS), for impulsive signals, and energy envelope, for emergent signals.

HOS of a probability distribution of a random variable X are given by the standardized central moments, which can be expressed as

$$S_n = \frac{E[(X - \mu)^n]}{(E[(X - \mu)^2])^{n/2}} = \frac{m_n}{m_2^{n/2}}, \quad (1)$$
where $E[\cdot]$ is the expectation operator, $\mu$ is the mean defined by $\mu = E[X]$, $m_2 = \sigma^2$ is the variance (second moment), and $m_n$ is the nth central moment. Most widely used for the detection of transients in non-stationary random signals are the third and fourth normalized central moment HOS, known as skewness and kurtosis. Nonetheless, some of the recent applications of HOS in signal processing algorithms (e.g. Lokajíček & Klíma 2008) demonstrated that HOS of the fifth and sixth order can provide robust results in detecting signals corresponding to acoustic emissions in low SNR conditions. Here we will only focus on the fourth central moment HOS (kurtosis), leaving the potential extensions to other statistic moments as subject of future research.

Kurtosis is defined as a fourth moment about the mean normalized by the square of the second moment or variance. Following eq. (1) it can be written as:

$$S_4 = \frac{E[(X - \mu)^4]}{(E[(X - \mu)^2])^2} = \frac{m_4}{m_2^2}, \qquad (2)$$

where $m_4$ is the fourth central moment. The above expression defines a population kurtosis. For a discrete signal $u(t) = \{u(t_1), u(t_2), \dots, u(t_M)\}$, we define a HOS CF as a sample kurtosis, expressed as:

$$\mathrm{CF}_{\mathrm{HOS}}(t_i) \equiv \hat{S}_4(t_i) = \frac{\frac{1}{M} \sum_{i=1}^{M} (u(t_i) - \hat{\mu})^4}{\left(\frac{1}{M} \sum_{i=1}^{M} (u(t_i) - \hat{\mu})^2\right)^2}, \qquad (3)$$

where $M$ is the number of samples in the chosen window, and $\hat{\mu}$ is the estimate of the mean for $M$ samples.

As discussed in Saragiotis et al. (2002), presence of transient changes in the non-stationary statistics of the signal, like phase arrivals associated with earthquake source, temporally generates rapidly increasing values of kurtosis function. This makes kurtosis HOS an efficient CF for extracting the information on the onset of the short transient signals. However, kurtosis function by itself cannot be used as the detector of the seismic phases, since the large values of kurtosis can be associated with any transient providing a significant change in signal statistics. A following step, including coherency analysis across the network, is necessary for differentiating a kurtosis maxima corresponding to an accidental, non-correlated spike/noise from the arrival of the phase that is coherent over the array's stations.

On the other hand, the signal associated with tremor-like sources is characterized by slow energy transients, lasting several tens of seconds and not having clear separate impulsive phases, embedded within even slower variations of energy with time, of the order of minutes to hours (e.g. Obara 2002). For this kind of signals a more appropriate characterization is one based on the amplitude or energy modulation, that is, the energy envelope (e.g. Obara 2002; McCausland et al. 2005; Wech & Creager 2008; Maeda & Obara 2009; Suda et al., 2009; Ide 2010), which is sensitive, at the same time, to the slow-varying background energy level and to the time-localized energy transients. We shall see, in the next sections, that the relative timing of localized energy transients can be used to detect and locate a tremor source.

We thus define an envelope CF based on root mean square (RMS) envelope, expressed as:

$$\mathrm{CF}_{\mathrm{env}}(t_i) \equiv \mathrm{RMS}(t_i) = \sqrt{\frac{\sum_{i=1}^{M} u(t_i)^2}{M - 1}}, \qquad (4)$$

where $M$ is the number of samples in the analysed discrete signal.

### 2.1.1 Recursive scheme for CF computation

To improve computation efficiency and to make the approach suitable for possible real-time implementations, we estimate the HOS and envelope CFs using a recursive scheme, similar to the one described in Langet et al. (2014), which provides an efficient way for accumulating time-averaged statistics. An estimate of the recursive exponential average $\hat{\mu}(t)$ for a signal $u(t)$ is expressed as:

$$\hat{\mu}(t_i) = C \cdot u(t_i) + (1 - C) \cdot \hat{\mu}(t_{i-1}), \qquad (5)$$

where $u(t_i)$ is the ith sample of the signal, and $C: 0 \le 1 - C \le 1$ is the decay constant.

By directly extending the above definition to the estimates of the second and fourth central moments defined in eqs (2) and (3) we obtain following expressions:

$$\begin{aligned} \hat{m}_2(t_i) &= C \cdot (u(t_i) - \hat{\mu}(t_{i-1}))^2 + (1 - C) \cdot \hat{m}_2(t_{i-1}) \\ \hat{m}_4(t_i) &= C \cdot (u(t_i) - \hat{\mu}(t_{i-1}))^4 + (1 - C) \cdot \hat{m}_4(t_{i-1}). \end{aligned} \qquad (6)$$

The final formulation of recursive HOS CF for each time step $t_i$ can then be defined by the ratio of two quantities in eq. (6) according to the kurtosis definition provided in eq. (2):

$$\begin{aligned} \mathrm{CF}_{\mathrm{HOS}}(t_i) &\equiv \hat{S}_4(t_i) = \hat{m}_4(t_i) / (\hat{m}_2(t_i))^2 \\ &= \frac{C \cdot (u(t_i) - \hat{\mu}(t_{i-1}))^4 + (1 - C) \cdot \hat{m}_4(t_{i-1})}{(C \cdot (u(t_i) - \hat{\mu}(t_{i-1}))^2 + (1 - C) \cdot \hat{m}_2(t_{i-1}))^2}. \end{aligned} \qquad (7)$$

Following the same idea, we define the recursive envelope CF as a single component recursive RMS envelope of the recorded signal $u(t)$, using the following expression:

$$\mathrm{CF}_{\mathrm{env}}(t_i) \equiv \widehat{\mathrm{RMS}}(t_i) = \sqrt{C \cdot u^2(t_i) + (1 - C) \cdot (\widehat{\mathrm{RMS}}(t_{i-1}))^2} \qquad (8)$$

where $C$ is a decay constant, similarly to eqs (5) and (6).

In the above formulations, the decay constant $C$ can be seen as a smoothing factor, defining how fast the memory of older data points will be discarded. We express the decay constant as $C = \Delta T/T_{\mathrm{decay}}$, where $\Delta T$ is the data sampling interval and $T_{\mathrm{decay}}$ represents the time-length of the averaging scale. Strictly, the lower limit of time-averaging scale $T_{\mathrm{decay}}$ is defined by the sampling interval (i.e. $T_{\mathrm{decay}} \ge \Delta T$). In practice, it should be sufficiently large to allow a meaningful, stable representation of the signal statistical properties. The upper limit of $T_{\mathrm{decay}}$ depends on the type of analysed data, the time scale of the targeted transients, and the total length of the signal. In general, this parameter should be set separately for each particular case study, through a preliminary testing step.

Illustration of the adopted recursive scheme is provided in Fig. 2 on an example of two local earthquakes (Fig. 2a), and a 1 hr tectonic tremor record (Fig. 2c). Fig. 2(b) illustrates how the kurtosis function (eq. 7) reacts to P-wave arrival onsets, and shows that the decay tail of the local kurtosis maxima increases with $T_{\mathrm{decay}}$. The amplitude of the kurtosis peak does not depend on the amplitude of the signal itself—the small and the large event have comparable values of kurtosis—and is rather related to the information on how fast and significant the change in the statistics of the signal associated with the transient is (e.g. how impulsive is the arrival of the phase).

Fig. 2(d) shows a comparison of envelope functions (eq. 8) for different values of $T_{\mathrm{decay}}$ for the tectonic tremor seismogram. The envelope reacts to local signal energy and the time-averaging scale $T_{\mathrm{decay}}$ acts as a smoothing factor.
<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Seismogram of typical earthquakes recorded at a local station. The inset figure shows an enlarged view of preceding smaller event. (b) Kurtosis HOS CF of the record in panel (a) calculated using the recursive formulation (eq. 7), for different decay constants $T_{decay}$. (c) 1 hr raw seismogram of a typical energetic tectonic tremor recorded at a Hi-net station in Shikoku (southwestern Japan). (d) Envelope CF of the record in panel (e) calculated using the recursive RMS scheme (eq. 8) for different decay constants $T_{decay}$.</figcaption></figure>

### 2.1.2 Time-frequency representation and broad-band CF

Seismic signals are broad-band by their nature. These properties are better characterized through a time-frequency representation of the signal (i.e. a spectrogram), like the short-time Fourier transform, the continuous wavelet transform, or the S-transform (see Tary et al. 2014 for a review). A signal processing scheme combining the time-domain and frequency-domain analysis will potentially yield a more complete representation of the signal, by providing an estimate of both the time at which a particular transient occurs and its frequency characteristics, without requiring a-priori knowledge of its frequency content (e.g. Gibbons et al. 2008).

Figs 3(b) and (f) show an example of spectrograms (S-transform; Stockwell et al. 1996) calculated for the same signals shown in Fig. 2: a recording of several local earthquakes (Fig. 3a), and a 1 hr tectonic tremor seismogram (Fig. 3e). The spectrogram of the earthquake record (Fig. 3b) indicates the presence of an energy arrival associated with a smaller earthquake preceding the main event, observable on the original seismic record and detected by the kurtosis CF, as shown in Fig. 2(b). Some later energy can as well be identified on the spectrogram within the coda of the large earthquake. These signals, predominant at frequencies over 10.0 Hz, are small earthquakes hidden in the coda of the larger event (see inset in Fig. 3a), and are not detected by the kurtosis CF of the original, non-filtered record (Fig. 2b). Looking at the spectrogram of the tectonic tremor record (Fig. 3f), one can identify both the duration and predominant frequency range of 1.0–15.0 Hz corresponding to the main tremor event, as well as short duration broad-band energy spikes corresponding to local earthquakes (e.g. at ~3000 s).

In order to properly extract non-stationary and frequency-dependent statistical features of a signal $u(t)$, we seek to extend the HOS and envelope CFs (eqs 7 and 8) to the time-frequency transform of the signal ($U(t, f)$). Therefore, a time-frequency HOS CF ($CF_{HOS}(t, f)$) is constructed by calculating the kurtosis HOS of the time-frequency signal $U(t, f)$ along the time axis for each discrete frequency $f$. An example of such CF is shown in Fig. 3(c). The

time-frequency HOS CF clearly depicts late energy arrivals, corresponding to the smaller earthquakes in the coda of the large event, and appearing as local maxima at frequencies larger than 10.0 Hz (e.g. at ~130 s and ~155 s).

We also define a time-frequency envelope CF ($CF_{env}(t, f)$), which can be seen as a smoothing (through recursive RMS computation) along the time axis of the time-frequency signal $U(t, f)$. The effect of this transformation is to enhance significant energy peaks and reduce background noise (Fig. 3g).

As a final stage of the time-frequency analysis, a single broad-band CF ($CF^{TF}(t)$) can be composed from the time-frequency CF ($CF(t, f)$). The strategy for performing the composition can differ in terms of the selected frequency interval [$f_{min}$ : $f_{max}$], and operator. Here we use two alternative operators: the maximum and the RMS, defined respectively as:

$$CF^{TF}(t) = \max_{f} CF(t, f), \quad f \in [f_{min}, f_{max}] \qquad (9)$$

and

$$CF^{TF}(t) = \sqrt{\frac{\int_{f_{min}}^{f_{max}} (CF(t, f))^2 df}{f_{max} - f_{min}}}. \qquad (10)$$

Selecting one operator over another depends on the type of information to be included into the CF: the RMS operators will provide an averaged overview for the given frequency interval, while the maximum is generally a better choice for enhancing narrow-band energy arrivals hidden in the coda of larger events, or masked by a noisy background. The choice of frequency interval allows including or disregarding some parts of the information (e.g. low/high-frequency noise, impulsive/emergent transients, ...).

An example illustrating the influence of selected frequency range ($[f_{min}, f_{max}]$) is provided in Figs 3(d) and (h). Here, a RMS composition operator is applied. The example shows that, including the higher frequencies into $CF^{TF}(t)$ allows incorporating the energy arrivals from the later smaller events hidden in the coda of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a–d) Time–frequency representation applied to a seismogram of a M 5.2 local earthquake: (a) raw seismogram. The inset figures show the enlarged views of preceding smaller event, and an early aftershock in the coda of M 5.2 earthquake respectively; (b) spectrogram (S-transform) of the seismogram in panel (a); (c) time–frequency kurtosis $CF_{HOS}(t, f)$ of the spectrogram in panel (b); and (d) examples of broad-band CF ($CF^{TF}(t)$) composed from transformed spectrogram (c) in different frequency ranges. (e–h) Time–frequency representation applied to a seismogram of a tectonic tremor recorded at a Hi-net station in southwestern Japan: (e) 1 hr raw seismogram; (f) spectrogram (short-time Fourier transform) of the raw seismogram in panel (a); (g) time–frequency envelope $CF_{env}(t, f)$ of the spectrogram in (f); (h) examples of summary CF recomposed from transformed spectrogram (g) in different frequency ranges.</figcaption></figure>

large earthquake (Fig. 3a) that do not appear on the kurtosis CF of the original non-filtered record (Fig. 2b). Similarly, in the case of tectonic tremor (Fig. 3h), since most of its energy is contained in the lower frequency range, using high frequencies (≥20.0 Hz) for constructing the CF will result in eliminating from $CF^{TF}(t)$ the information on tremor’s energy characteristics, and enhancing the impact of short-duration transients (e.g. earthquakes).

### 2.1.3 Time-frequency representation using multiband filter (MBF) algorithm

The problem of estimating the spectrum of a time-limited signal is a non-trivial issue with multiple trade-offs, depending on the type of selected transform, and the characteristics of the signal. Moreover, this problem can become computationally challenging, especially when it is applied to large data sets.

An efficient alternative to the time-frequency decomposition is a time-domain analysis scheme making use of filter-banks. Here we follow the multiband filter (MBF) recursive scheme of Lomax *et al.* (2012), detailed in Appendix A1. The MBF algorithm provides

a time-frequency decomposition of the signal as a set of band-pass filtered time-series ($U(t, f_n)$). This is achieved by running the original record $u(t)$ through a predefined bank of narrow band-pass filters covering the interval of interest $[f_{min}, f_{max}]$, and having $n = 0, \dots, N_{band}$ central frequencies $f_n$ (Fig. A1). Filtering can be efficiently implemented through a cascade of simple low-pass and high-pass, one-pole recursive filters.

Examples of the MBF signals $U(t, f_n)$ for the two types of the seismic records, earthquake (Fig. 4a) and tectonic tremor (Fig. 4e), are given in Figs 4(b) and (f), respectively. Filtered traces in Fig. 4(b) are obtained by applying twelve logarithmically spaced filters with central frequencies covering the range of 0.02–50.00 Hz to the earthquake record shown in Fig. 4(a). For the case of tremor record (Fig. 4e), a filter-bank of twelve logarithmically spaced filters in the range of 2.00–40.00 Hz is used. Resulted time-frequency decomposed signals $U(t, f_n)$ are shown in Fig. 4(f).

Following the same procedure as that presented for the spectral time-frequency decomposition of the signal, we estimate the time-frequency CFs using the output of the MBF analysis. Once the set of filtered signal is generated, a set of CFs ($CF(t, f_n)$) is calculated by applying the recursive kurtosis (Fig. 4d), or alternatively recursive
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a–d) Illustration of MBF algorithm and broad-band CFs on an example of a local earthquake: (a) Raw broad-band seismogram. (b) Filtered traces. (c) Broad-band CF—CF<sup>TF</sup>(t) composed over all the frequencies according to eq. (9). (d) Kurtosis CFs of filtered signals. (e–h) Illustration of the MBF algorithm and resulted final CF on an example of a tectonic tremor recorded at a Hi-net station in southwestern Japan: (e) Raw 1 hr seismogram. (f) Filtered traces. (g) Broad-band CF—CF<sup>TF</sup>(t) composed over all the frequencies according to eq. (10). (h) RMS envelope CFs of filtered signals.</figcaption></figure>

envelope to each of the signals $U(t, f_n)$ (Fig. 4h). A single time dependent broad-band CF (CF<sup>TF</sup>(t)) can then be constructed using one of the operators discussed above (eqs 9 and 10). For the example in Fig. 4, a maximum operator (eq. 9) is used for the earthquake dataset, while a RMS summation operator (eq. 10) is used for the tremor signal.

In general, spectral and time-domain (MBF) time-frequency compositions should provide similar characterization of the signal, when the number of filters in the MBF filter-bank is equal to the number of scales in the spectrogram (see e.g. Mill & Brown 2005). A comparison between the CFs of raw record, the broad-band CFs obtained from spectral decomposition, and the broad-band CFs obtained from MBF analysis is shown in Figs 5(a) and (b) for the case of earthquake and tectonic tremor records, respectively. It can be observed that the spectral and MBF time-frequency representations provide similar results.

### 2.1.4 HOS derivative and generalized CF

HOS (including kurtosis) are asymmetric functions that ramp up quickly and decay with time, when computed over a sliding window (eq. 3) or with a recursive scheme (eq. 6; Fig. 6). Furthermore, the HOS maximum is delayed with respect to the onset of the transient, since a number of samples are required to provide stable statistics. To reduce this effect, we look at the positive derivative of CF<sub>HOS</sub>(t), defined as:

$$ \dot{\text{CF}}_{\text{HOS}}^{+}(t) = \begin{cases} \dot{\text{CF}}_{\text{HOS}}(t) & \text{if } \dot{\text{CF}}_{\text{HOS}}(t) \ge 0 \\ 0 & \text{elsewhere} \end{cases} \qquad (11) $$

where $\dot{\text{CF}}_{\text{HOS}}(t)$ is the time derivative of the (broad-band) HOS CF. This new function has the local maximum corresponding to the maximum gradient of kurtosis, which is closer to the onset of

the transient (Fig. 6). To account for timing uncertainty, we finally convolve the positive derivative of HOS CF with a Gaussian window, having the width equal to the width of the HOS CF (Fig. 6).

We now have all the elements to define a generic CF for our problem, incorporating both the HOS and envelope CFs:

$$ \text{CF}(t) = \begin{cases} \dot{\text{CF}}_{\text{HOS}}^{+}(t) * e^{-t^2/4\sigma^2} \\ \text{CF}_{\text{env}}(t) \end{cases} \qquad (12) $$

where '*' is the convolution operator and $\sigma = T_{\text{decay}}/2$. In the above equation, $\text{CF}_{\text{HOS}}(t)$ and $\text{CF}_{\text{env}}(t)$ are either simple or broad-band CFs.

Such a definition allows for a flexible signal processing scheme in which representation of the signal can be carried by selecting a single transformation that is best suited for the problem, or running the two transformations in parallel, if a more complex characterization of the signal is targeted. This can be especially useful for characterizing and classifying the seismic sources in the environments where accumulated energy is released through different mixed processes such as, for example, regular and slow earthquakes in subduction zones.

### 2.1.5 Example application of the signal processing scheme

The proposed processing scheme provides an efficient way of signal characterization that accounts for the broad-band and non-stationary nature of the recorded signals, as well as for the variability in size and energy release associated with the seismic sources. The output of this step is a generic broad-band CF (eq. 12) accounting for a priori unknown transient signals, and preparing data for a better performance in detection and location steps. In Fig. 1(b), we provide a summarized block-diagram of the signal processing steps presented throughout the discussions above. Application of these steps to the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Comparison of the kurtosis CF (black line) for a raw recording of an earthquake (grey line), and the broad-band kurtosis CFs composed from the transformed spectrogram (blue line), and the MBF analysis (red line). (b) Comparison of the RMS envelope CF (black line) of a raw recording of a tectonic tremor (grey line), and the broad-band envelope CFs composed from the transformed spectrogram (blue line), and the MBF analysis (red line).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Seismogram of a typical earthquake recorded at a local seismic station. (b) 50 s zoom on the M 5.2 earthquake record corresponding to the green window in panel (a). (c) Red line—broad-band MBF kurtosis HOS CF of record (b) estimated using eq. (7), and decay constant $T_{\text{decay}} = 1.0$ s; green line—derivative of the CF; blue line—generic broad-band MBF kurtosis CF, resulted from convolution of the CF positive derivative with the Gaussian window of $\sigma^2 = T_{\text{decay}}/2$ (eq. 11).</figcaption></figure>

entire seismic network (e.g. Figs 7a and d) will result in replacement of the original raw seismic records (Figs 7b and e) by their broad-band CFs (Figs 7c and f) extracting a particular feature of the seismic signal. This provides the desired CFs necessary for applying coherency-based detection and location methods to the local and regional wide aperture arrays. The estimated broad-band CFs serve as input time-series for the detection and location algorithm (Fig. 1a).

## 2.2 Space-time imaging for detection and location problem

### 2.2.1 Station-pair time-delay estimate (TDE) function

Let's consider a set of receivers {$r_i$}, $i = 1, \dots, N$ recording the signal generated by a source at location $\mathbf{q}$. The seismic field due to this seismic source associated with the traction $\mathbf{T}$ operating in the direction $j$ at $\mathbf{q}$ and time $t'$, recorded by the $i$th station ($i = 1, \dots, N$) can be expressed as (e.g. Aki & Richards 2002):

$$u_i(t) = \mathbf{G}_{ji}(t; \mathbf{q}, t') * T_j(\mathbf{q}, t') + \varepsilon_i(t) \quad (13)$$

where $\mathbf{G}_{ji}(t; \mathbf{q}, t')$ is the elastodynamic Green's function between the source and the $i$th station located at $r_i$, $\varepsilon_i(t)$ is additive noise, and '*' is the convolution operator.

A correct evaluation of eq. (13) requires detailed knowledge of the spatiotemporal parameters of the source, as well as of the propagation medium for accurate estimation of the Green's functions, which is, seldom the case. The time-frequency signal processing scheme presented in previous section has an effect of replacing the observed seismic signals $\mathbf{u}(t) = [u_1(t), \dots, u_N(t)]^T$ with a set of generic CFs $\mathbf{c}f(t) = [\mathbf{c}f_1(t), \dots, \mathbf{c}f_N(t)]^T$ (eq. 12; Fig. 5). Eq. (13) can be then modified as:

$$\mathbf{c}f_i(t) = \text{CF}[\mathbf{G}_{ji}(t; \mathbf{q}, t') * T_j(\mathbf{q}, t') + \varepsilon_i(t)] \quad (14)$$

where CF[$\cdot$] is the time-frequency CF operator corresponding to the processing scheme illustrated in Fig. 1(b).

The CF defined in eq. (12) is mainly sensitive to phase or energy arrival times, through its HOS and RMS components, respectively (see e.g. Fig. 5). We, therefore, simplify eq. (14), assuming that the propagation term of the CFs of the signal is a simple time shift $\tau_i(\mathbf{q})$, depending on the distance from the source to receiver and the effective velocity of the medium. We thus write eq. (14) as:

$$\mathbf{c}f_i(t + \tau_i(\mathbf{q})) = \tilde{\mathbf{s}}(\mathbf{q}, t) \quad (15)$$

where $\tilde{\mathbf{s}}(\mathbf{q}, t)$ is an effective source term, responsible of the phase or energy arrivals that we are characterizing (incorporating noise term—$\tilde{\varepsilon}_i(t)$), and $\tau_i(\mathbf{q}) = \tau_i(\mathbf{q}, r_i)$ is the traveltime from the source to the $i$th station, expressed as $\tau_i(\mathbf{q}) = \int_{\mathbf{q}}^{r_i} c^{-1}(\mathbf{r})d\mathbf{r}$, with $c(\mathbf{r})$ being an effective propagation velocity in the medium.

We next partition the receivers into pairs, and denote $\Omega$ the set of independent pairs. For a source $s$ at $\mathbf{q}$ and a pair $p$: $[r_i, r_j] \in \Omega$, the theoretical time difference of arrival (TDOA) is defined as

$$\Delta\tau_p(\mathbf{q}) = \tau_i(\mathbf{q}) - \tau_j(\mathbf{q}) \quad (16)$$

We now define the functions of TDE (TDE functions) at each time $t$, for each receiver pair $p$ and source $s$, as a measurement $\Psi_p(t, \Delta\tau)$ satisfying the following conditions:

$$\Psi_p(t, \Delta\tau) \in [0, 1] \quad (17)$$

$$\Delta\tau \in [-\Delta\tau_{p,\text{max}}, \Delta\tau_{p,\text{max}}] \quad (18)$$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a–c) Example of MBF signal characterization applied to IMAD post-seismic network in Chile: (a) Map view of central Chile; coloured triangles indicate the IMAD seismic stations. The colour of triangle corresponds to the institution to which the station belongs. White star shows the epicentre of the M 5.2 earthquake after Rietbrock *et al.* (2012). (b) Record section of UD component seismograms corresponding to the 3 min record of M 5.2 earthquake at the stations. Dark yellow line indicates the P-wave moveout. Grey triangle marks the record presented in Figs 2(a)–6(a), shown in grey. (c) Record section of broad-band kurtosis HOS CFs from MBF analysis in the frequency range 0.02–50.00 Hz. Dark yellow and grey lines indicate, respectively, the move-out of P-wave corresponding to M 5.2 event, and later small event in its coda. (d)–f) Example of MBF signal characterization scheme applied to Hi-net stations in southwestern Japan: (d) Map view of southwestern Japan; grey triangles indicate the Hi-net seismic stations. White stars show the epicentres of earthquakes during 2100–2200 JST, 2012 May 27, from JMA catalogue. White circle correspond to the position of the tectonic tremor activity during same period reported by NIED. (e) Record section of 1 hr EW component seismograms of the tremor activity recorded during 2100–2200 JST, 2012 May 27 at Hi-net stations. Red shadow lines mark the times corresponding to the two earthquakes from JMA catalogue. (f) Broad-band RMS envelope CFs, estimated by applying the MBF analysis in the frequency range 2.0–40.0 Hz to the records in panel (e).</figcaption></figure>
where $\Delta\tau_{p,\max}$ is the maximum time delay for a given pair of stations, depending on the distance between the stations and the velocity of the medium.

The measurements correspond therefore to a p-dimensional vector of TDE functions:

$$ \mathbf{y}(t, \Delta\tau) = f(\Delta\tau, \mathbf{c}f(t)) \quad (19) $$

with

$$ \mathbf{y}(t) = [\psi_1(t, \Delta\tau), \dots, \psi_{N_p}(t, \Delta\tau)] \quad (20) $$

where $N_p$ is the total number of station-pairs. Among possible TDE functions are properly normalized local (i.e. time-dependent) cross-correlation functions.

Each TDE $\mathcal{Y}_p(t, \Delta\tau)$ at station-pair $p$, defines a set of 3-D potential source locations that, in a case of constant velocity model, are given by hyperboloid objects (Font et al. 2004). Therefore, at each time $t$, the TDE function itself, which is a continuous function, can be the basis of a likelihood function. Selecting a spatial source location $\mathbf{q}$, and assigning the TDOA and TDE function to a pair of stations $p$, the likelihood position of the source is defined as:

$$ P(\mathbf{y}(t)|\tilde{\mathbf{s}}(\mathbf{q})) = \mathcal{Y}_p(t, \tau_p(\mathbf{q})) \in [0, 1] \quad (21) $$

where $P(\cdot|\cdot)$ represents a conditional likelihood, normalized in $[0, 1]$.

We refer to the likelihood position of the source expressed by eq. (21) as the pairwise Spatial Likelihood Function (SLF). Estimation of pairwise SLF relies on TDE functions, which, as mentioned previously, can be represented by any properly normalized local cross-correlation (LCC) function, as long as conditions defined in eqs (17) and (18) are satisfied. Thus, eq. (21) can be written as:

$$ P(\mathbf{y}(t)|\tilde{\mathbf{s}}(\mathbf{q})) = CC_p(t, \Delta\tau_p(\mathbf{q})) \quad (22) $$

### 2.2.2 Local cross-correlation (LCC)

LCC gives a time-dependent measure of similarity between the two signals by estimating their cross-correlation in the neighbourhood of each data point (e.g. Hale 2006a; Fomel 2007). This provides an important advantage when detection and location of multiple, possibly mixed seismic sources is targeted. Most commonly, estimation of LCC is performed in a sliding window (e.g. Birchfield & Gillmor 2002; Birchfield 2004), using the following formulation:

$$ \begin{aligned} CC_p(t, \Delta\tau_p(\mathbf{q})) &= \int_{t-W/2}^{t+W/2} c f_i(t' - \tau_i(\mathbf{q})) c f_i(t' + \tau_j(\mathbf{q})) dt' \\ &= \int_{t+\tau_j-W/2}^{t+\tau_j+W/2} c f_i(t' - \Delta\tau_p(\mathbf{q})) c f_i(t') dt'. \end{aligned} \quad (23) $$

Here, $W$ is the size of the selected time window.

Direct implementation of eq. (23), for each time sample, can lead to substantial demand of computational resources. We therefore use a different approach, based on the algorithm proposed by Hale (2006a), which achieves better efficiency by avoiding explicit windowing of the signals for each time sample, implied by eq. (23). Moreover, calculations can be implemented recursively.

The details on the LCC calculation scheme of Hale (2006a) are provided in Appendix A2. In general terms, the idea consists in computing LCC as convolution of a Gaussian filter and the product of the shifted signals. This results in an entirely recursive implementation, provided that Gaussian filtering is performed using recursive

scheme (e.g. Young & van Vliet 1995). For a given station-pair $p: [\mathbf{r}_i, \mathbf{r}_j]$, and corresponding CFs, this can be expressed as:

$$ \begin{aligned} \text{LCC}_p(t, \Delta\tau_p(\mathbf{q})) = & \int_{-\infty}^{+\infty} c f_i \left(t' - \frac{\Delta\tau_p(\mathbf{q})}{2}\right) c f_j \left(t' + \frac{\Delta\tau_p(\mathbf{q})}{2}\right) \\ & \times v(t - t', \Delta\tau_p(\mathbf{q})) dt' \end{aligned} \quad (24) $$

with $v(t - t', \Delta\tau_p(\mathbf{q}))$ being a Gaussian filter, defined as:

$$ v(t, l) = e^{-\Delta\tau_p(q)^2/4\sigma^2} e^{-t^2/\sigma^2} \quad (25) $$

where the variance $\sigma^2$ corresponds to the half-width of Gaussian function used for filtering. The choice of this parameter should reflect uncertainties incorporated into the velocity model, and in the method for computing the TDOA ($\Delta\tau_p(\mathbf{q})$), and thus should be set for each case separately. The LCC function estimated following this algorithm is a 2-D function of time $t$ and TDOA—$\Delta\tau_p(\mathbf{q}) \in [-\Delta\tau_{p,\max}, \Delta\tau_{p,\max}]$ (eq. 18); its values lie in the interval $[0, 1]$, since both kurtosis HOS and envelope CFs are always positive functions. Thus, LCC satisfies the condition in eq. (17).

An example of LCC calculated using this scheme for two synthetic traces composed of identical triangular signals, with different time-lags, is shown in Fig. 8. The 2-D local cross-correlation LCC($t, \Delta\tau$)—function of time $t$ and discrete time-lag $\Delta\tau$—has local maxima at the point with coordinates $[t_{\text{LCCmax}}, \Delta\tau_{\text{LCCmax}}]$ (see Fig. 8b). Once these coordinates are known, it is possible to identify the time samples of the cross-correlated signals that correspond to this $\text{LCCmax}$ as:

$$ t_1 = t_{\text{LCCmax}} - \frac{\Delta\tau_{\text{LCCmax}}}{2}, \quad \text{and} \quad t_2 = t_{\text{LCCmax}} + \frac{\Delta\tau_{\text{LCCmax}}}{2}. \quad (26) $$

We shall see that this property is useful for origin time estimation.

An example of LCC for a pair of actual signals is shown in Figs 9(a) and (b) using a MBF kurtosis CF, applied to the earthquake records from IMAD stations in Chile. The example illustrates how LCC is able to provide a continuous estimation in time of the time-lag between the two signals.

### 2.2.3 Imaging function for detection and location problem

Eqs (21) and (22) define the pairwise SLF consisting of weighted non-overlapping 3-D objects which, in the case of constant velocity model, are reduced to hyperboloids. The actual location of the source can be then obtained by combining the pairwise SLFs of all independent station-pairs in $p: [\mathbf{r}_i, \mathbf{r}_j] \in \Omega$. Among the combination operators of interest are those commutative, monotonic associative and bounded between $[0, 1]$. Here we use the summation operator, which satisfies all the above conditions, and express the combination of the pairwise SLFs at each time $t$ as:

$$ \begin{aligned} P(\mathbf{y}(t)|\tilde{\mathbf{s}}(\mathbf{q})) &= \frac{1}{N_p} \sum_{p \in \Omega} \mathcal{Y}_p(t, \Delta\tau_p(\mathbf{q})) \\ &= \frac{1}{N_p} \sum_{p \in \Omega} \text{LLC}_p(t, \Delta\tau_p(\mathbf{q})) \end{aligned} \quad (27) $$

where $N_p$ is the number of station-pairs used in summation.

The above equation defines the complete SLF in 4-D space-time coordinates. In order to simplify the implementation of the algorithm, we compute eq. (27) over a set of sliding time windows of size $W_p = \Delta\tau_{p,\max}$, equal to the maximum possible time delay for a given pair of stations ($p \in \Omega$). It is possible then to reduce the TDE function $\text{LCC}_p(t, \Delta\tau_p(\mathbf{q}))$ for each of the station-pairs to a single-variable function of time-lag ($\Delta\tau_p(\mathbf{q})$) over this time
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Example of LCC calculated for two normalized synthetic traces—signal 1 (green line) and signal 2 (yellow line). (b) Zoomed-in view at the portion of the LCC function and synthetic signals highlighted in panel (a); corresponding time-lag $\Delta\tau$, and picked times of the signals (arrival times) estimated from maximum LCC are marked by dotted lines. (c) $LCC_{MAX}$ function, corresponding to the TDE function used in the detection and location scheme, estimated by applying eq. (28).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) In black—transformed MBF kurtosis GFs corresponding to seismograms (plotted in grey) from two stations in central Chile. The MBF analysis is done in the frequency range 0.02–50.00 Hz, and kurtosis is calculated assuming $T_{decay} = 1.0$ s. The green shadowed window highlights the maximum of CFs corresponding to the *P*-wave arrival from the M 5.2 earthquake. (b) LCC of the CFs in panel (a). (b') Zoomed-in view showing the LCC function highlighted by window in panels (a) and (b), and the corresponding $LCC_{MAX}$ TDE function. (c) One horizontal and two vertical cross-sections through the centre of 3-D grid of theoretical TDOA for the stations 1 and 2, calculated assuming homogeneous model with the constant *P*-wave velocity of 6.0 km/s. (d) One horizontal and two vertical cross-sections passing through the hypocentre location of the M 5.2 (Rietbrock *et al.* 2012) earthquake, showing the projection of $LCC_{MAX}$ TDE function (b') on the theoretical differential traveltime grid (c).</figcaption></figure>
window $W_p$. Here we use the maximum operator to construct this single-variable LCC function:

$$LCC_{MAX,p}(\tilde{t}, \Delta\tau_p(\mathbf{q})) = \max\{LCC_p(t, \Delta\tau_p(\mathbf{q}))|t \\ \in [\tilde{t} - W_p/2, \tilde{t} + W_p/2]\} \quad (28)$$

where $[\tilde{t} - W_p/2, \tilde{t} + W_p/2]$ is the time window centred at $\tilde{t}$.

Going further with the simplification, we may use a single time window for all the station pairs, with size fixed to the largest possible time-lag:

$$W = \max\{\Delta\tau_{p,\max}|p \in \Omega\} \quad (29)$$

This reduces the complexity of TDE functions manipulation taking, however, into account the move-out corresponding to all possible source locations $\mathbf{q}$ in the 3-D volume of interest.

Fig. 8(c) provides an illustration of $LCC_{MAX,p}(\tilde{t}, \Delta\tau_p(\mathbf{q}))$ constructed over a given time window on an example case of two simple synthetic signals (Fig. 8c).

The maximum likelihood location of the source expressed by the normalized imaging function in eq. (27) at each time $t = \tilde{t}$ will then be:

$$P(\mathbf{y}(\tilde{t})|\tilde{\mathbf{s}}(\mathbf{q})) = \frac{1}{N_p} \sum_{p \in \Omega} LCC_{MAX,p}(\tilde{t}, \mathbf{q}) \quad (30)$$

where $LCC_{MAX,p}(\tilde{t}, \mathbf{q})$ corresponds to the LCC function defined in eq. (28) mapped (projected) onto TDOA estimated in a common spatial coordinate system, under the assumption of a propagating phase and a given velocity model.

When working with CFs calculated for real signals (Fig. 9b'), the $LCC_{MAX,p}(\tilde{t}, \mathbf{q})$ function will possibly contain multiple peaks, and the time delay for the true source may not correspond to the highest peak. The summation over the multiple pairs of stations (eq. 30) handles this problem, enhancing the SLFs corresponding to the sources that are coherent for most of the receivers, and attenuating those of incoherent signals.

The LCC-based estimation of the TDE functions proposed here shares the same foundation with the equal differential time (EDT) location method formulation used for arrival-time earthquake location (e.g. Font et al. 2004; Lomax 2005); or the TDE approach used for determining the location of sound source in acoustics (e.g. Birchfield & Gillmor 2002; Birchfield 2004). However, instead of picking the maximum of cross-correlation over a given time window, and thus estimating the value of EDT, we stack the entire maximum LCC function (eq. 30). This allows taking into account all the available information before deciding on a possible detection, resulting in more robust behaviour of the method. Such delaying of decision-making is the base of the principle of least commitment (e.g. see Birchfield & Gillmor 2002; Birchfield 2004).

In Appendix A3, the presented imaging function is compared with the traditional beamforming and amplitude stack (often used in migration), demonstrating the effectiveness of using a cross-correlation term in constraining the location likelihood.

## 2.3 Space-time detection and location of seismic sources

### 2.3.1 Spatial location of seismic sources

The formulation of the detection and location problem provided here (e.g. eq. 21) implies that theoretical TDOA for assumed phase(s) and a given velocity model, are estimated in a designed 3-D coordinate system for all station-pairs before the summation in eq. (30) is

performed. In fact, our SLFs are composed of the TDE functions projected (or mapped) onto the TDOA assigned to pairs of stations.

Calculation of TDOA is an independent step of our method that can be done by any available tool. Once defined, a 3-D grid covering the volume in which potential seismic sources are expected, one can calculate beforehand theoretical traveltimes for any possible source-station configuration, under an assumption of a velocity model (1-D, 2-D or 3-D), and propagating phase (e.g. P, S, or P and S). Here we use the Grid2Time routine of NonLinLoc program by Lomax (2005, 2008) based on the Eikonal finite-difference scheme of Podvin & Lecomte (1991).

A potential source location (SLF) is then estimated by mapping the TDE functions according to pre-calculated TDOAs, and summing them over the pairs of stations. This procedure of mapping implies the change of support for the TDE function from time-lag $\Delta\tau_p(\mathbf{q})$ to spatial coordinates $\mathbf{q} = (x, y, z)$. An illustration of this step for a single pair of stations is given in Fig. 9 (and Supporting Information Fig. S1). The example shows a TDE function (Fig. 9b') calculated for a time window containing the P-wave arrival from a local earthquake that is mapped on the TDOA calculated for a constant velocity model (Fig. 9c). Resulting SLF (Fig. 9d) corresponds to a hyperboloid passing through the hypocentre of the source, marked by the star (Fig. 9d). When performing the summation of these SLFs over the increasing number of station-pairs (Fig. 10), the intersections of hyperboloid objects will define a region in space representing the total likelihood of source location. This corresponds to the imaging function in eq. (30). After all the information has been taken into account, the location with the highest likelihood is selected as the estimate for the position of the seismic source. The procedure is repeated for each position of the sliding window, providing thus a framework for an automated technique of analysing the continuous data, yielding the information on the spatial and temporal evolution of the seismic energy sources.

### 2.3.2 Origin time estimation

The expression of LCC in eq. (24) is independent of the origin time of the earthquake. Thus, a traditional 4-D problem of hypocentre location is reduced to a search for a 3-D maximum likelihood spatial position of the source, provided by the imaging function in eq. (30). Once the position of the seismic source is determined, we estimate a posteriori the origin time $T_0$ as the average of the observed arrival times ($t_i$) minus the predicted traveltimes ($\tau_i(\mathbf{q})$) for the peaks of detection functions that contributed to the corresponding imaging function. Following Moser et al. (1992) this can be expressed as:

$$T_0(\mathbf{q}) = \frac{\sum_{i=1}^{N_p} w_i (t_i - \tau_i(\mathbf{q}))}{\sum_{i=1}^{N_p} w_i}, \quad (31)$$

where $N_p$ is the number of station pairs and $w_i$ are weights attributed to the observed arrival times $t_i$. The weighting coefficients for each station can be expressed as, for example, a mean of the pairwise SLFs to which corresponding station is contributing. However, in the current version we assume a homogeneous weighting equal to one for all the arrival times.

Here, the arrival times $t_i$ are estimated using eq. (26) for each pair of stations that contributes to the imaging function, corresponding to the maximum likelihood location of the source $\mathbf{q}$. Thus, we obtain $(N_p - 1)$ estimates of $t_i$ from the same hypocentre for each of the stations. Although not exploited here, this information has a potential of providing statistical characteristics and error estimations of the observed arrival times.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Examples of imaging functions (eq. 30) for the case of (a) two, (b) three, (c) four and (d) 22 stations.</figcaption></figure>

# 3 APPLICATION EXAMPLES

We further present several examples illustrating the application of the proposed detection and location method for regular earthquakes and tectonic tremor in subduction zones. An example illustrating the performance of the method on a synthetic dataset is provided in Appendix A4.

## 3.1 Detecting and locating earthquakes

As discussed previously (Section 2.1.2), a CF based on broad-band signal representation and kurtosis HOS, has an advantage of emphasizing short transients (e.g. small earthquakes hidden by stronger, narrow-band noise) that often would be missed by the traditional signal processing algorithms (e.g. using pre-defined band-pass filters, and/or STA/LTA as CFs). We already demonstrated the advantage of such broad-band CFs for enhancing signals associated with small earthquakes hidden in the coda of the larger event on an example record corresponding to a local M 5.2 earthquake (Fig. 3a; Rietbrock et al. 2012) recorded at the IMAD stations from Central Chile (Fig. 7).

To illustrate the performance of the method in detecting and locating earthquakes using continuous broad-band records from a regional network, we selected a subarray of 22 stations of the IMAD post-seismic network (Figs 7a and 11). The total length of analysed data is 3 min. First, broad-band kurtosis CFs are estimated applying the MBF filter algorithm assuming a bank of 12 filters with logarithmically spaced central frequencies covering range of 0.02–50.00 Hz. The decay constant for recursive kurtosis calculation, $T_{decay}$ (eq. 6), is set to 1.0 s and is assumed to be constant for all subbands of the filter. This assumption of the constant $T_{decay}$ is a reasonable simplification, since we are mostly interesting in picking precise onset time of the phase arrivals corresponding to the aftershock seismicity recorded at a regional scale, and dominant in high-frequency seismic energy. For a more general case, were both low frequency and high frequency energy arrivals should be accounted for, frequency-dependent $T_{decay}$ assumption might be more appropriate. The final generic broad-band HOS CFs (eq. 12) are shown in Fig. 11(d). For this example, only vertical components of the records are considered under the assumption of P-wave propagation. Theoretical P-wave traveltimes are calculated using homogeneous velocity model with the constant velocity of 6.0 km/s, on the grid of 200 × 200 × 50 km³. Next, 3-D maximum likelihood
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Detection and location example for regular earthquakes. (a) A horizontal and two vertical sections thorough the maximum of imaging function corresponding to a small local earthquake not identified in the catalogue of Rietbroek *et al.* (2012), and shown according to a selected colour scale. Green star indicates the result of location. White triangles correspond to the location of the IMAD station from central Chile. (b) Same as panels (a), but for the M 5.2 local earthquake. White star shows the hypocentre of the earthquake after Rietbroek *et al.* (2012). (c) Same as panels (a) and (b) for an early aftershock, not listed in the catalogue. (d) Transformed MBF summary kurtosis CFs (black lines) of the original UD component seismograms (grey lines) of the IMAD stations used in example. The green rectangle indicates the time windows corresponding to the time intervals of the imaging functions in panels (a)–(c). Vertical red bars show the theoretical arrival times for determined locations (green stars).</figcaption></figure>
location of the sources are obtained by stacking the TDE functions (corresponding to maximum LCCs of normalized CFs, projected on the theoretical TDOA grids) over sliding windows. The size of the sliding window, as well as the maximum time-lag for calculating LCC, is set to 27 s, according to the maximum possible time-lag for the given station configuration (eq. 29). Simultaneous detection and location are performed when the value of imaging function, normalized by the number of station-pairs, exceeds the threshold value of 0.7. Our method could identify three events within the 3 min of the analysed continuous data. Maximum likelihood locations, as well as the normalized imaging functions, are shown in Figs 11(a)–(c). In particular, the result corresponding to the M 5.2 earthquake is presented in Fig. 11(b). Obtained location is in good agreement with that of Rietbrock et al. (2012), shown by the white star: the difference between the two locations is of about 3 km in horizontal and 8 km in vertical direction. Larger hypocentral depth, found in our case, can be attributed to the simple homogenous velocity model that does not take into account shallow, low velocity structures. However, it should be noticed, the catalogue location has a relatively high likelihood of 0.75.

Figs 11(a) and (c) show the locations obtained for an event preceding the M 5.2 earthquake, and a small event following it, and hidden in the coda (see e.g. Fig. 2a). Red vertical lines in Fig. 11(d) mark the corresponding arrivals. These events could not be identified in the catalogue of Rietbrock et al. (2012).

## 3.2 Detecting and locating an energetic tectonic tremor

The detection and location approach is further tested on the Nankai subduction zone in southwestern Japan. This area is characterized by intense seismic activity, comprising a wide variety of phenomena observed by the dense Hi-net seismic network of NIED. The tectonic tremors in this area are currently routinely located by NIED using the hybrid technique of Maeda & Obara (2009), in combination with clustering method of Obara et al. (2010). Together with the catalogue of background seismic activity published by JMA, these offer a perfect case-study for evaluating the capability and potential of the proposed methodology to detect and locate the energy release associated with possibly overlapping seismic radiation from earthquakes and low-frequency tectonic tremors.

We first present an example of detecting and locating a high-energy part of the signal associated with a tectonic tremor using the MBF envelope CFs (Fig. 4g). The analysis follows the same steps as in previous example (see diagram in Fig. 1). We first estimate CFs corresponding to broad-band envelope CFs though the MBF algorithm assuming a bank of 12 filters with logarithmically spaced central frequencies covering range of 2.0–40.0 Hz. The decay constant for recursive RMS envelope calculation, $T_{decay}$ (eq. 8), is set to 0.5 s for the central frequency band of MBF filter-bank. For other frequencies, the decay constants are calculated using the following weighting scheme:

$$T_{decay}^{n} = T_{decay} \cdot f_{N_{band}/2} / f_{n} \quad (32)$$

This provides a frequency-dependent RMS envelope estimation (e.g. Fig. 4h). The generic broad-band envelope CF is further composed by applying the RMS operator (eq. 10) to the time-frequency MBF envelopes (Fig. 4g). The resulted detection functions for the horizontal EW components and a subset of 22 stations, selected from the overall dataset (Fig. 7f), are shown in Fig. 12. In this example we only analysed 20 min of recording corresponding to the most energetic part of the tremor. Next, 3-D maximum likelihood

locations are obtained by stacking the TDE functions projected on the theoretical TDOA grids following the steps in Fig. 1(a). Here, we work under the assumption of S-wave propagation, and calculate theoretical traveltime grids for a homogeneous velocity model with the constant velocity of 3.5 km/s (Maeda & Obara 2009). Only station-pairs with inter-station distance not exceeding 80 km were used. The size of the window and the maximum time-lag for LCC are set to 60 s, according to the maximum distance between the stations (eq. 29). The resulting 3-D maximum likelihood location of the source for a single time window, associated with the high-energy part of the tectonic tremor within the analysed time interval, is shown in Fig. 12(a). For comparison we also show the 1 hr averaged tectonic tremor location provided by the NIED catalogue (white star in Fig. 12a). Our location is compatible with that of NIED, taking into account that the 1 hr averaged location is associated with the tremor activity within an aperture of 7 km (Obara et al. 2010), and that the depth of the tremor is fixed to the subduction plate interface, at 30–40 km.

## 3.3 Detecting and locating complex events: tectonic tremors and local earthquakes

It may very often be the case when the tectonic tremor activity overlaps with the background regular seismicity. To provide an example of such application for the described method we selected a 1 hr record from Hi-net, NIED stations in Shikoku (southwestern Japan), covering the period of 2012 May 31 from 21:59:00 to 22:59:00 JST (Japan Standard Time). NIED catalogue indicates a rather low-energy tremor activity during this hour, located in the northeast part of Shikoku (blue circle in Fig. 13a). According to the JMA seismicity catalogue, during the same hour a M 1.2 local earthquake occurs in the central part of the region (blue star, Fig. 13a). Visual inspection of the corresponding record section (Fig. 13b) indicates that the tectonic tremor and earthquake are overlapping in time. Tremor activity starts before the earthquake, continuing after it for several minutes. On the record section, the local earthquake can be identified as a short transient with high amplitude at around 1000 s.

We applied the detection and location scheme to these 1 hr records. Envelope CFs were calculated using the same MBF filter-bank, and decay constant settings as presented above (eq. 32), applied to the two components of horizontal waveforms (combined NS and EW). The TDE functions, corresponding to the maximum LCC functions, were projected according to S-wave theoretical traveltime for a homogeneous velocity model with the constant velocity of 3.5 km/s. By running a sliding window of 60 s, with a shift of 3 s, across the records, a location corresponding to the maxima of the normalized imaging function is provided for every time window for which this maxima exceeds the value of the trigger, set to 0.875. To avoid including incoherent information into the detection and location problem by cross-correlating all the available station-pairs, we divided them in two subarrays: 22 southwest stations and 19 northeast stations. The final locations are obtained by combining the information from the two subarrays. The results are summarized in Figs 13(a) and (b). The distribution of colour-coded circles in Fig. 13(a), corresponding to the resulted trigger locations, indicates that the tremor activity, during which the local M 1.2 event occurs, corresponds to the cluster of locations about 10 km south from the earthquake’s hypocentre. That tremor event is not reported in the NIED catalogue. We could as well recover the tremor activity that was reported by NIED catalogue and located northeast.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Detection and location example for an energetic part of a tectonic tremor. (a) A horizontal and two vertical sections thorough the maximum of imaging function shown according to the selected colour scale. Green star indicates resulted location. White star shows the location of the tremor associated with the period of 2100–2200 JST, 2012 May 27, reported by NIED. Black triangles show the location of the Hi-net stations in southwestern Japan used in the example. (b) 1 hr MBF RMS envelope CFs corresponding to the EW component seismograms of the tremor activity (plotted in grey) recorded at selected Hi-net stations. The green shadow rectangles indicate the time window corresponding to the spatial stack in panel. (a) Vertical red bars show the theoretical arrival times for determined location (green star).</figcaption></figure>

Analysis of the temporal evolution of the results (Fig. 13b) indicates that the two tremor episodes are partially overlapping in time. Most probably the tremor located close to the local earthquake was missed by the routine procedure of NIED because this part of the record was disregarded due to the large amplitude arrivals corresponding to the earthquake. Since in our method records are analysed in continuous manner, it was possible to locate this cluster of tremor activity. Thus, detection and location schemes based on the analysis of the coherency of continuous seismic records across the stations, and making use of efficient signal processing algorithm, can potentially provide a more complex and revealing picture of seismic energy release processes and its space–time evolution. This example emphasizes the ability of the presented method to detect and locate overlapping seismic sources in complex environments.

# 4 DISCUSSIONS AND CONCLUSIONS

We presented a new method for detection and location of seismic sources that lays the ground towards the development of fully automated multiscale framework for extraction (i.e. detection, location and classification) of multiple sources with a wide range of dynamic and spatial scales.

The proposed method is a multistep procedure consisting of signal processing and detection and location schemes designed to

enhance coherent statistical features of the wave field recorded across the array of sensors. The employed signal processing scheme makes use of multiband filter (MBF) algorithm and provides a frequency-dependent representation of non-stationary statistical features, through higher-order statistics (HOS) and envelope-based characteristic functions (CFs). The two types of CF can be either used separately, or run in parallel in the detection and location scheme, depending on the type of information one wants to use, or the type of signals that are targeted. Kurtosis CFs can be successfully applied for analysing and extracting short transients, like earthquakes, possibly hidden within strong, narrow-band noise. For the analysis of seismic records corresponding to tectonic tremors, like those from Hi-net stations in southwestern Japan, the use of both kurtosis and envelope CF can provide interesting complementary information (e.g. detecting low-frequency earthquakes within tectonic tremor activity). Thus, the choice of transform used for estimating the CF will depend strongly on the desired features to be extracted from the signal, allowing to take into account the variability in size and energy release associated with the seismic sources. The signal processing scheme plays a double role of providing a representation of recorded seismic signal and extracting frequency-dependent information. The main goal of this step is to account for a-priori unknown transient signals, potentially associated with different sources, and to prepare data for a better performance in the detection and location part.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Detection and location example for mixed tectonic tremor and earthquake sources. (a) Map view of southwestern Japan, and two vertical cross-sections showing the spatial and temporal distribution of the two tremor episodes and a local earthquake detected and located by analysing the continuous recordings from Hi-net stations (triangles), covering period of 2200–2300 JST, 2012 May 31. The subarrays of 22 southwest and 19 northeast stations are shown by green and grey triangles, respectively. The circles correspond to the location of the maxima of 3-D stacked SLF function (value for the trigger was set to 0.875). The colour of the circles denotes the timing of the triggered time-windows relative to the beginning of the records. Blue circle shows the location of tremor activity reported by NIED for the analysed time period. Blue star corresponds to the epicentre of the earthquake indicated in JMA catalogue. (b) Horizontal, EW component, seismograms filtered in 2.0–40.0 Hz frequency band, for Hi-net stations exemplifying the analysed dataset. Coloured rectangles mark in time the corresponding triggered spatial locations. The colour of the window is related to the colour of the two subarrays of station in panel (a).</figcaption></figure>

The detection and location step of the method is divided in two parts. First, station-pair time delay estimate (TDE) functions are evaluated by cross-correlation of the transformed signals. Then, simultaneous detection and location is performed by estimating the array imaging function corresponding to sum of TDE functions mapped to the corresponding theoretical time difference of arrivals (TDOAs), calculated in a common 3-D coordinate system, under the assumption of propagating phases for a given velocity model.

Through the sequence of these procedures, the signal CFs are transformed into time-series of 3-D spatial images, representing the likelihood of each pixel to be the part of a seismic source. The final source location is then provided as corresponding to the maxima of the 3-D spatial time-series maps. Assumed imaging function allows taking into account all the available information before making the decision about event detection, resulting in a more robust behaviour of the method (e.g. Birchfield & Gillmor 2002; Birchfield 2004).
The capability of the method to detect and locate small events hidden in the coda of larger ones (e.g. during an aftershock sequence) is exemplified on the case of the IMAD post-seismic deployment in Chile. This example emphasizes the importance of the multiband signal characterization scheme and the efficiency of HOS in extracting low-energy transients embedded in narrow-band, stronger background noise. The capability and potential of the method to detect and locate the energy release associated with overlapping seismic radiation from earthquakes and low-frequency tectonic tremors is further tested on the Nankai subduction zone in southwestern Japan. Both examples confirm that the proposed methodology can provide significant improvement to existing automated methods.

In the form presented here, the method is based on the assumption of a single phase, either *P* (earthquake location) or *S* (tremor location), using vertical or horizontal components. It can be however easily extended to include both *P* and *S*-phase arrivals at the same time. Such developments represent an important near-future expansion that will be presented in the following studies.

Being based on TDE, thus having a cross-correlation term in the imaging function, the detection and location scheme provides accuracy similar to that of beamforming (Appendix A3). A main difference with standard beamforming is that the proposed method does not require a delay-and-sum process (i.e. performing an explicit grid-search over all the possible source locations), but rather a direct mapping of the TDE function onto the 3-D space (i.e. a change of support). This property, along with a fully recursive implementation, makes the approach computationally efficient. This, as well, implies the suitability of the method for future real-time applications and for a straightforward parallelization. The parallel version of the code (BackTrackBB) was used for the application example presented in Section 3.3. It required 22 min to produce the final locations shown in Fig. 13(a), on a machine with 16 cores (2 x Intel Xeon E5-2690 2.90 GHz) and 64 GB of RAM. The BackTrackBB code is distributed on open-source basis, and can be downloaded from http://backtrackbb.github.io.

The proposed methodology is not limited to large-scale regional seismic networks, but can be as well applied at much smaller phenomenological scales like, for example: fault zone microseismicity or artificially induced seismicity associated with mining or fluid injection activity.

## ACKNOWLEDGEMENTS

The authors express their thanks to Frank Krueger, Francesco Grigoli and an anonymous reviewer for their thorough review. We would also like to acknowledge Aline Gendrin and other researchers from Schlumberger Gould Research Centre in Cambridge, UK for their comments and suggestions for improving the manuscript. We thank R.-E. Plessix, and A. Lomax for scientific discussions during this work. We also thank students Pierre Romanet and Vegar Flatoy pointing out and correcting a number of typos and mistakes in the text and formulas. We are grateful to the IMAD team in Chile and M. El-Aissaoui for their continuous effort in providing high quality seismic data. We used the JMA earthquake, and the NIED tremor catalogues as well as the data of Hi-net from NIED and IMAD, Central Chile, seismic networks. This study has been partly funded by the Centre de Recherches sur le Stockage Géologique du CO2 (IPGP/TOTAL/SCHLUMBERGER/ADEME), and UniEarthS LabEx program of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). Part of the data analysis was im-

plemented using ObsPy (Megies et al. 2011). Figures were prepared using GMT (Wessel et al. 2013) and Matplotlib (Hunter 2007).

## REFERENCES

Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, 2nd edn, p. 700, University Science Books.

Aldridge, D.F., 1990. The Berlage wavelet, *Geophysics*, **55**(11), 1508-1511.

Allen, R.V., 1978. Automatic earthquake recognition and timing from single traces, *Bull. seism. Soc. Am.*, **68**, 1521-1532.

Allen, R.V., 1982. Automatic phase pickers: their present use and future prospects, *Bull. seism. Soc. Am.*, **72**, 225-242.

Baillard, C., Crawford, W.C., Ballu, V., Hibert, C. & Mangeney, A., 2014. An automatic kurtosis-based *P*- and *S*-phase picker designed for local seismic networks, *Bull. seism. Soc. Am.*, **104**, 394-409.

Baker, T., Granat, R. & Clayton, R.W., 2005. Real-time earthquake location using Kirchhoff reconstruction, *Bull. seism. Soc. Am.*, **95**(2), 699-707.

Beck, S. et al., 2014. Advancing subduction zone science after a big quake, *EOS, Trans. Am. geophys. Un.*, **95**, 193-194.

Beroza, G.C. & Ide, S., 2011. Slow earthquakes and nonvolcanic tremor, *Annu. Rev. Earth. Planet. Sci.*, **39**, 271-296.

Birchfield, S.T., 2004. A unifying framework for acoustic localization, in *Proceedings of the 12th European Signal Processing Conference (EUSIPCO)*, Vienna, Austria.

Birchfield, S.T. & Gillmor, D.K., 2002. Fast Bayesian acoustic localization, in *Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)*, Orlando, Florida.

Cesca, S. & Grigoli, F., 2015. Full waveforms seismological advances for microseismic monitoring, *Adv. Geophys.*, **56**, 169-228.

Das, I. & Zoback, M.D., 2013a. Long-period, long-duration seismic events during hydraulic stimulation of shale and tight-gas reservoirs—Part 1. Waveform characteristics, *Geophysics*, **78**, KS97-KS108.

Das, I. & Zoback, M.D., 2013b. Long-period long-duration seismic events during hydraulic stimulation of shale and tight gas reservoirs—Part 2. Location and mechanisms, *Geophysics*, **78**, KS109-KS120.

Drew, J., White, R.S., Tilmann, F. & Tarasewicz, J., 2013. Coalescence microseismic mapping, *Geophys. J. Int.*, **195**(3), 1773-1785.

Fomel, S., 2007. Local seismic attributes, *Geophysics*, **72**(3), A29-A33.

Font, Y., Kao, H., Lallemand, S., Liu, C.-S. & Chiao, L.-Y., 2004. Hypocentral determination offshore eastern Taiwan using the maximum intersection method, *Geophys. J. Int.*, **158**, 655-675.

Freiberger, W.F., 1962. An approximate method in signal detection, Technical Report No. 12, Department of the Army, Division of Applied Mathematics, Brown University.

Gajewski, D. & Tessmer, E., 2005. Reverse modeling for seismic event characterization, *Geophys. J. Int.*, **163**, 276-284.

Geiger, L., 1912. Probability method for the determination of earthquake epicenters from the arrival time only (translated from German), *Bull. St Louis Univ.*, **8**, 56-71.

Gentili, S. & Michelini, A., 2006. Automatic picking of *P* and *S* phases using a neural tree, *J. Seismol.*, **10**, 39-63.

Gibbons, S.J., Ringdal, F. & Kværna, T., 2008. Detection and characterization of seismic phases using continuous spectral estimation on incoherent and partially coherent arrays, *Geophys. J. Int.*, **172**, 405-421.

Grigoli, F., Cesca, S., Vassalo, M. & Dahm, T., 2013. Automated seismic event location by travel-time stacking: an application to mining induced seismicity, *Seismol. Res. Lett.*, **84**(4), 666-677.

Grigoli, F., Cesca, S., Amoroso, O., Emolo, A., Zollo, A. & Dahm, T., 2014. Automated seismic event location by waveform coherence analysis, *Geophys. J. Int.*, **196**, 1742-1753.

Hale, D., 2006a. An efficient method for computing local cross-correlations of multi-dimensional signals, *CWP Report*, 656.

Hale, D., 2006b. Recursive Gaussian filters, *CWP Report*, 546.

Hunter, J.D., 2007. Matplotlib: a 2D graphics environment, *Comput. Sci. Eng.*, **9**(3), 90-95.

Ide, S., 2010. Striations, duration, migration and tidal response in deep tremor, *Nature*, **466**, 356-359.
Ide, S., Beroza, G.C., Shelly, D.R. & Uchide, T., 2007. A scaling law for slow earthquakes, *Nature*, **447**, 76–79.

Ito, Y., Obara, K., Shiomi, K., Sekine, S. & Hirose, H., 2007. Slow earthquakes coincident with episodic tremors and slow slip events, *Science*, **315**(5811), 503–506.

Kao, H. & Shan, S.-J., 2004. The source-scanning algorithm: mapping the distribution of seismic sources in time and space, *Geophys. J. Int.*, **157**(2), 589–594.

Kao, H. & Shan, S.J., 2007. Rapid identification of earthquake rupture plane using Source-Scanning Algorithm, *Geophys. J. Int.*, **168**(3), 1011–1020.

Kuperkoch, L., Meier, T., Lee, J., Friederich, W. & EGELADOS Working Group. 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics, *Geophys. J. Int.*, **181**(2), 1159–1170.

Langet, N., Maggi, A., Michelini, A. & Brenguier, F., 2014. Continuous Kurtosis-based migration for seismic event detection and location, with application to Piton de la Fournaise Volcano, La Réunion, *Bull. seism. Soc. Am.*, **104**, 229–246.

Larmat, C., Montagner, J.-P., Fink, M., Capdeville, Y., Tourin, A. & Clévédé, E., 2006. Time-reversal imaging of seismic sources and application to the great Sumatra earthquake, *Geophys. Res. Lett.*, **33**, doi:10.1029/2006GL026336.

Liao, Y.C., Kao, H., Rosenberger, A., Hsu, S.K. & Huang, S.B., 2012. De- lineating complex spatiotemporal distribution of earthquake aftershocks: an improved Source-Scanning Algorithm, *Geophys. J. Int.*, **189**(3), 1753–1770.

Lokajíček, T. & Klíma, K., 2008. High-order statistics approach: automatic determination of sign and arrival time of acoustic emission signals, in *28th European Conference on Acoustic Emission*, pp. 82–87, European Working Group on Acoustic Emission (EWGAE).

Lomax, A., 2005. A reanalysis of the hypocentral location and related obser- vations for the Great 1906 California Earthquake, *Bull. seism. Soc. Am.*, **95**, 861–877.

Lomax, A., 2008. Location and Tectonics of the focal region of the California Earthquake of 18 April 1906, *Bull. seism. Soc. Am.*, **98**, 846–860.

Lomax, A., Michelini, A. & Curtis, A., 2009. Earthquake location, direct, global-search methods, in *Complexity in Encyclopedia of Complexity and System Science*, Part 5, pp. 2449–2473, Springer, doi:10.1007/978-0-387-30440-3.

Lomax, A., Satriano, C. & Vassallo, M., 2012. Automatic picker devel- opments and optimization: FilterPicker—a robust, broadband picker for real-time seismic monitoring and earthquake early-warning, *Seismol. Res. Lett.*, **83**, 531–540.

Maeda, T. & Obara, K., 2009. Spatiotemporal distribution of seismic en- ergy radiation from low-frequency tremor in western Shikoku, Japan, *J. Geophys. Res.*, **114**, B00A09, doi:10.1029/2008JB006043.

McCausland, W., Malone, S. & Johnson, D., 2005. Temporal and spatial occurrence of deep non-volcanic tremor: from Washington to northern California, *Geophys. Res. Lett.*, **32**, L24311, doi:10.1029/2005GL024349.

McMechan, G.A., Luetgert, J.H. & Mooney, H.M., 1985. Imaging of earth- quake sources in Long Valley Caldera, California, 1983, *Bull. seism. Soc. Am.*, **75**(4), 1005–1020.

Megies, T., Beyreuther, M., Barsch, R., Krischer, L. & Wassermann, J., 2011. ObsPy – what can it do for data centers and observatories? *Ann. Geophys.*, **54**, 47–58.

Mill, R. & Brown, G., 2005. *Auditory-Based Time-Frequency Representations and Feature Extraction Techniques for Sonar Processing*, Speech and Hearing Research Group.

Moser, T.J. van Eck, T. & Nolet, G., 1992. Hypocenter determination in strongly heterogeneous earth models using the shortest path method, *J. geophys. Res.*, **97**, 6563–6572.

Nadeau, R.M. & Dolenc, D., 2005. Nonvolcanic tremors deep beneath the San Andreas Fault, *Science*, **307**(5708), 389–389.

Nadeau, R.M. & McEvilly, T.V., 1999. Fault slip rates at depth from re- currence intervals of repeating microearthquakes, *Science*, **285**(5428), 718–721.

Obara, K., 2002. Nonvolcanic deep tremor associated with subduction in southwest Japan, *Science*, **296**(5573), 1679–1681.

Obara, K., Kasahara, K., Hori, S. & Okada, Y., 2005. A densely distributed high-sensitivity seismograph network in Japan: Hi-net by National Research Institute for Earth Science and Disaster Prevention, *Rev. Sci. Instrum.*, **76**, 021301, doi:10.1063/1.1854197.

Obara, K., Tanaka, S., Maeda, T. & Matsuzawa, T., 2010. Depth-dependent activity of non-volcanic tremor in southwest Japan, *Geophys. Res. Lett.*, **37**, L13306, doi:10.1029/2010GL043679.

Okada, Y., Kasahara, K., Hori, S., Obara, K., Sekiguchi, S., Fujiwara, H. & Yamamoto, A., 2004. Recent progress of seismic observation networks in Japan—Hi-net, F-net, K-NET and KiK-net, *Earth Planets Space*, **56**, xv–xviii.

Podvin, P. & Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, *Geophys. J. Int.*, **105**, 271–284.

Rietbrock, A., Ryder, I., Hayes, G., Haberland, C., Comte, D., Roecker, S. & Lyon-Caen, H., 2012. Aftershock seismicity of the 2010 Maule $M_w = 8.8$, Chile, earthquake: correlation between co-seismic slip models and aftershock distribution? *Geophys. Res. Lett.*, **39**, L08310, doi:10.1029/2012GL051308.

Rogers, G. & Dragert, H., 2003. Episodic tremor and slip on the Cascadia subduction zone: the chatter of silent slip, *Science*, **300**(5627), 1942–1943.

Saragiotis, C.D., Hadjileontiadis, L.J. & Panas, S.M., 2002. PAI-S/K: a robust automatic seismic P phase arrival identification scheme, *IEEE Trans. Geosci. Remote Sens.*, **40**, 1395–1404.

Stockwell, R.G., Mansinha, L. & Lowe, R.P., 1996. Localization of the complex spectrum: the S transform, *IEEE Trans. Signal Process.*, **44**(4), 998–1001.

Suda, N., Nakata, R. & Kusumi, T., 2009. An automatic monitoring system for non-volcanic tremors in southwest Japan, *J. Geophys. Res.*, **114**, B00A10, doi:10.1029/2008JB006060.

Tary, J.B., Herrera, R.H., Han, J. & van der Baan, M., 2014. Spectral estimation—What is new? What is next?, *Rev. Geophys.*, **52**, doi:10.1002/2014RG000461.

Wech, A.G. & Creager, K.C., 2008. Automatic detection and location of Cascadia tremor, *Geophys. Res. Lett.*, **34**, 1–5.

Wessel, P., Smith, W., Scharroo, R., Luis, J. & Wobbe, F., 2013. Generic map-ping tools: improved version released, *EOS, Trans. Am. geophys. Un.*, **94**, 409–410.

Withers, M., Aster, R., Young, C., Beiriger, J., Harris, M., Moore, S. & Trujillo, J., 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection, *Bull. seism. Soc. Am.*, **88**(1), 95–106.

Withers, M., Aster, R. & Young, C., 1999. An automated local and regional seismic event detection and location system using waveform correlation, *Bull. seism. Soc. Am.*, **89**(3), 657–669.

Xuan, R. & Sava, P.C., 2010. Probabilistic micro-earthquake location for reservoir monitoring and characterization, *Geophysics*, **75**, MA9–MA26.

Young, T.I. & van Vliet, L.J., 1995. Recursive implementation of the Gaussian filter, *Signal Process.*, **44**, 139–151.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.** (a) In black—multiband filtering (MBF) root mean square (RMS) envelope characteristic functions (CFs) corresponding to the seismograms of a tectonic tremor (plotted in grey) recorded at two Hi-net stations in southwestern Japan. The MBF analysis is done in the frequency range 2.0–40.0 Hz, and RMS envelopes are calculated using the variable $T_{decay}$ scaled with the frequency. The green shadowed window highlights the maximum of CFs corresponding to the most energetic part of the tremor during 1 hr. (b) LCC of the CFs in panel (a). (b') Zoomed-in view showing the LCC function highlighted by window in
panels (a) and (b), and the corresponding LCC<sub>MAX</sub> function. (c) One horizontal and two vertical cross-sections through the centre of 3-D grid of theoretical TDOAs for stations 1 and 2, calculated assuming homogeneous model with the constant S-wave velocity of 3.5 km/s. (d) One horizontal and two vertical cross-sections through the hypocentre location of the tremor reported by NIED for 1 hr period of 2100–2200 UT, 2012 May 27, and showing the projection of LCC<sub>MAX</sub> function (b') on the theoretical TDOA grid (c) (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggw071/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

# APPENDIXES

## A1 Recursive multiband filtering (MBF) algorithm

The MBF formulation implemented here is based on the multiband processing scheme of Lomax et al. (2012) used in the FilterPicker algorithm built upon the definition of recursive filters.

Filter recurrence relations provide a way to determine the output samples of a filter in terms of the input samples and the preceding output. Given a discrete sampled signal $u(t_i)$, simple, one-pole recursive high-pass and low-pass filters can be defined as follows:

$$ 
\begin{aligned}
U^{\text{HP}}(t_i) &= C^{\text{HP}}[U^{\text{HP}}(t_{i-1}) + u(t_i) - u(t_{i-1})] \\
U^{\text{LP}}(t_i) &= U^{\text{LP}}(t_i) + C^{\text{LP}}[u(t_i) - U^{\text{LP}}(t_{i-1})]
\end{aligned}
\quad \text{(A1.1)}
$$

where $C^{\text{HP}} = w/(w + \Delta T)$ and $C^{\text{LP}} = \Delta T/(w + \Delta T)$ are filter constants, $\Delta T$ is the sampling rate of the data, $w = 1/2\pi f$ is the corner period of the filter, and $f$, its corner frequency. The filter constants $C^{\text{HP}}$ and $C^{\text{LP}}$—defined by the corner frequency $f$— determine how fast will the output of the filter decay, and how strongly it will be influenced by small changes in the input. The main advantage of implementing the recursive scheme of (A1.1) is its computational efficiency since it involves only logical and arithmetical calculations.

A broad-band characterization of the signal $u(t_i)$ can be achieved by running the signal through a bank of recursive filters defined by eq. (A1.1). Following the scheme of Lomax et al. (2012), we define a filter-bank of narrow band-pass filters with central frequencies $f_n \in [f_{\text{min}}, f_{\text{max}}]$, $n \in [1, N_{\text{band}}]$ (Fig. A1). The number of frequency bands $N_{\text{band}}$ depends on the desired sampling of the frequency space: one can, for instance, use linearly or logarithmically spaced central frequencies $f_n$. The set of filtered signals $U_n(t_i)$ is obtained by applying a cascade of two high-pass and two low-pass, one-pole recursive filters (eq. A1.1). This is equivalent to a two-pole band-pass filter (Fig. A1), and can be expressed in following way:

$$ 
\begin{aligned}
U_n^{\text{HP1}}(t_i) &= C_n^{\text{HP}}[U_n^{\text{HP1}}(t_{i-1}) + u(t_i) - u(t_{i-1})] \\
U_n^{\text{HP2}}(t_i) &= C_n^{\text{HP}}[U_n^{\text{HP2}}(t_{i-1}) + U_n^{\text{HP1}}(t_i) - U_n^{\text{HP1}}(t_{i-1})] \\
U_n^{\text{LP1}}(t_i) &= U_n^{\text{LP1}}(t_i) + C_n^{\text{LP}}[U_n^{\text{HP2}}(t_i) - U_n^{\text{LP1}}(t_{i-1})] \\
U_n^{\text{LP2}}(t_i) &= U_n^{\text{LP2}}(t_i) + C_n^{\text{LP}}[U_n^{\text{LP1}}(t_i) - U_n^{\text{LP2}}(t_{i-1})].
\end{aligned}
\quad \text{(A1.2)}
$$

Each of the filters defined in the relation above is equivalent to a single-pole recursive infinite impulse response (IIR) filter. The broad-band representation of the signal consists in the set of signals $U_n(t_i) \equiv U_n^{\text{LP2}}(t_i)$ corresponding to the output of the last filter

<figure><img src="image_16.png" /><figcaption>Figure A1. An example of frequency response for a filter-bank composed of eleven 2 high-pass and 2 low-pass serial recursive filters.</figcaption></figure>

in eq. (A1.2). The overall computational time of MBF increases proportionally with $N_{\text{band}}$.

## A2 An efficient algorithm for estimating local cross-correlation (LCC)

The cross-correlation of two continuous signals $f(t)$ and $g(t)$ is defined by

$$ cc(l) = (f \star g)(l) \equiv \int_{-\infty}^{\infty} f(t)g(t+l)dt \quad \text{(A2.1)} $$

where $l$ denotes cross-correlation lag, and ‘★’ is the correlation operator. The centred definition of cross-correlation can be obtained by the change of variable $t \to t - l/2$:

$$ cc(l) = (f \star g)(l) \equiv \int_{-\infty}^{\infty} f(t - l/2)g(t + l/2)dt \quad \text{(A2.2)} $$

which is symmetric in lag variable $l$.

Eq. (A2.1) can be written for two discrete sampled signals, $f(t_i)$ and $g(t_i)$, as:

$$ cc(l) = (f \star g)(l) \equiv \sum_{i=-\infty}^{\infty} f(t_i)g(t_i+l) \quad \text{(A2.3)} $$

LCC has the advantage of allowing to account for variations in signal features that are not constant for the duration of the signal or the selected time window, allowing estimating the cross-correlation in the vicinity of a given point.

A straightforward implementation of LCC would imply that the signal should be first truncated/tapered to zero outside of the specified window, and then cross-correlated. The procedure will be then repeated for a number of times corresponding to the number of assumed (overlapping) windows. An estimation of the LCC of the two signals $f(t_i)$ and $g(t_i)$ will imply thus a computational complexity $O(N_s \times N_w \times N_l)$, where $N_w$ is the number of non-zero points in
the time widow, $N_s$ number of samples on which the windows are centred, and $N_l$ is the number of time lags. This computational complexity can be high, especially for large windows.

The formulation of Hale (2006a) provides an efficient way of calculating the LCC for every sample of the two signals. This means that the number of window for which the cross-correlation is calculated is equal to the number of samples in the signals. To obtain an LCC, a smooth window function $\omega(t)$ is applied to each of the signals $f(t)$ and $g(t)$ before the computation. Hale (2006a) uses the Gaussian windowing function, due to its specific properties allowing to simplify significantly the LCC estimation algorithm (for more details refer to the original paper).

Let's consider that two windowed signals $f(t, t')$ and $g(t, t')$, defined as follows:

$$ 
\begin{aligned}
f(t, t') &\equiv f(t)\omega(t' - t) \\
g(t, t') &\equiv g(t)\omega(t' - t)
\end{aligned}
\qquad (A2.4) 
$$

where $\omega(t' - t) = e^{-\frac{(t'-t)^2}{2\sigma^2}}$ is a Gaussian window with the half-width $\sigma$, centred at $t = t'$. Getting back to the definition (A2.1), the LCC (or the cross-correlation of the two windowed signals) is defined as:

$$ \text{LCC}(t', l) \equiv \int_{-\infty}^{\infty} f(t, t')g(t+l, t')dt. \qquad (A2.5) $$

The change of variables $t \to t - l/2$ leads to a centred form of cross-correlation:

$$ 
\begin{aligned}
\text{LCC}(t', l) &\equiv \int_{-\infty}^{\infty} f\left(t - \frac{l}{2}, t'\right) g\left(t + \frac{l}{2}, t'\right) dt \\
&= \int_{-\infty}^{\infty} f\left(t - \frac{l}{2}\right) g\left(t + \frac{l}{2}\right) \\
&\qquad \times \omega\left(t' - t + \frac{l}{2}\right) \omega\left(t' + t + \frac{l}{2}\right) dt \\
&= \int_{-\infty}^{\infty} f\left(t - \frac{l}{2}\right) g\left(t + \frac{l}{2}\right) v(t' - t, l) dt \quad (A2.6)
\end{aligned}
$$

where the Gaussian window $v(t, l)$ is defined as:

$$ v(t, l) \equiv \omega\left(t + \frac{l}{2}\right) \omega\left(t - \frac{l}{2}\right) = e^{-l^2/4\sigma^2} e^{-t^2/\sigma^2}. \qquad (A2.7) $$

Eq. (A2.6) shows that for any lag $l$ the local cross-correlation LCC($t', l$) can be calculated for all time samples $t = t'$ by convolving Gaussian filter $v(t, l)$ with the product of shifted signals $f(t - \frac{l}{2})$ and $g(t + \frac{l}{2})$.

Here we work with discrete sampled signals $f(t_i)$ and $g(t_i)$, and thus, we need to compute LCC($t', l$) for sampled integer values of $t'$ and $l$. Taking this into consideration, we re-write the discrete version of eq. (A2.6), avoiding the linear interpolation for the odd samples of integer time lags $l$ as:

$$ \text{LCC}(t, l) = \sum_{i=-\infty}^{\infty} \tilde{h}(t_i, l) v(t - t_i, l) \qquad (A2.8) $$

where $t_i$ is the current time sample, $v(t - t_i, l)$ is the Gaussian filter, and $\tilde{h}(t_i, l)$ is expressed as:

$$ \tilde{h}(t_i, l) \equiv (f(t_i - l_f)g(t_i + l_f) + f(t_i - l_f)g(t_i + l_f))/2 $$

with $l_f \equiv [l/2]$ and $l_g \equiv [l/2]$.

The above suggests that the LCC estimated following the algorithm of Hale (2006a) is a 2-D function of lag $l$ and time $t$ that can be computed in two steps:

(1) Estimate a sum of the sample by sample product of shifted signals $f$ and $g$ given by $\tilde{h}(t_i, l)$.

(2) Filter the sum with a Gaussian window.

With the recursive implementation of the Gaussian filter in eq. (A2.8) (Young & van Vliet 1995; Hale 2006b), the cost of computing the LCC is proportional to the number of samples to be filtered $N_s$. It is independent of the width of Gaussian window, meaning that the complexity of computing LCC($t, l$) for any single lag $l$ is $O(N_s)$. Therefore, the complexity of computing $N_l$ lags of $N_s$ LCC is $O(N_l \times N_s)$, the same as that for $N_l$ lags of a single non-LCC.

Without entering into much detail, we will mention that the recursive Gaussian filtering scheme (Young & van Vliet 1995) resumes to applying a cascade of two IIR filters (implemented recursively) to the input data. The corresponding filter constants are being defined depending of the value of $\sigma$ in the Gaussian kernel. The value of $\sigma$ is in general selected by trial-and-error, depending on the desired goal of the filtering, and the type of filtered signal. For detailed explanation please refer to the paper of Young & van Vliet (1995).

## A3 Comparison of TDE imaging function and beamforming

Traditional delay-and-sum beamforming provide the maximum likelihood location of the source by computing the energy of the reconstructed signal $u_i(t)$ as:

$$ \text{P}_{\text{beam}}(\mathbf{q}) = \int_{t_0-W}^{t_0+W} \left[ \sum_{i=1}^{N} u_i(t + \tau_i(\mathbf{q})) \right]^2 dt \qquad (A3.1) $$

After some algebraic manipulations, the above expression can be re-written in the following form:

$$ 
\begin{aligned}
\text{P}_{\text{beam}}(\mathbf{q}) = 2 \sum_{i=1}^{N} \sum_{j=i+1}^{N} \int_{t_0-W}^{t_0+W} u_i(t + \tau_i) u_j(t + \tau_j) dt \\
+ \sum_{i=1}^{N} \int_{t_0-W}^{t_0+W} u_i^2(t + \tau_i) dt
\end{aligned}
\qquad (A3.2) 
$$

The integral under the summation in the first term of the above equation represents the $CC_p(\Delta\tau_p(\mathbf{q}))$ cross-correlation of equation (e.g. eq. 24); while the second corresponds to the combined energy that can be attributed to a given source $\mathbf{q}$ in all the signals ($E_i(\tau_i(\mathbf{q}))$):

$$ \text{P}_{\text{beam}}(\mathbf{q}) = 2 \sum_{i=1}^{N} \sum_{j=i+1}^{N} CC_p(\Delta\tau_p(\mathbf{q})) + \sum_{i=1}^{N} E_i(\tau_i(\mathbf{q})) \qquad (A3.3) $$

Thus, the difference between the method presented in this paper, and the classical beamforming is the energy term. It can be demonstrated, however (Birchfield & Gillmor 2002; Birchfield 2004), that taking into account the total energy of the signal does not provide better constraint on location, acting more as a scaling factor. A synthetic example comparing beamforming with the TDE method based on LCC-functions proposed here (eq. 24) and with amplitude stacking is presented in Fig. A2. The figure shows how the cross-correlation term dominates in the beamforming image (Fig. A2a),
<figure><img src="image_17.png" /><figcaption>**Figure A2.** (a) One horizontal and two vertical cross-section illustrating the example of a beamforming energy stack calculated for the synthetic signals in panel (d). White star corresponds to the assumed location of the source. (b) Same as panel (a) but for the current method using the imaging function based on stacked TDE functions mapped on theoretical traveltime grids. (c) Same as panels (a) and (b) for the amplitude stack. (d) Synthetic signal corresponding to the epicentre location shown by white star in (a)–(c), calculated for the stations under consideration.</figcaption></figure>

which has a spatial focusing similar to the LCC image (Fig. A2b). In contrast, methods based on a simple amplitude stack provide less spatial constraint (Fig. A2c).

The main difference between the LCC-based methodology proposed here and standard beamforming is that our approach does not require a—computationally expensive—delay-and-sum process (i.e. performing an explicit grid-search over all the possible locations of the source), but rather relies on direct mapping of the LCC function into the 3-D space (through a change of function support). This, along with a fully recursive implementation, makes the approach computationally effective.

## A4 Synthetic tests

The performance of the method is tested on a synthetic dataset set-up for the station geometry of Hi-net in southwestern Japan and composed of 21 stations (Fig. 12a). We consider four test events located within the 3-D local XYZ grid with 1 km grid spacing defined as in Fig. 12a. An isotopic source is associated with each of

the four events. The test is designed to demonstrate the efficiency of both time-frequency signal processing and the space-time detection and location schemes.

The synthetic signal at each station corresponds to a Berlage wavelet (e.g. Aldridge 1990) given by:

$$w(t) = AH(t) t^n e^{-\alpha t} \cos(2\pi f_0 t + \phi_0), \quad (A4.1)$$

where $H(t)$ is the Heaviside unit step function ($H(t) = 0$ for $t \le t_0$ and $H(t) = 1$ for $t > t_0$); $A, n, \alpha$ and $\phi_0$ are non-negative adjustable parameters controlling the shape of the wavelet, $f_0$ is the predominant frequency, and $t_0$ corresponds to the onset time of the wavelet given by the propagation time from the location of the test source to the station. Berlage wavelet respects the condition of causality required by the physical phenomena of wave propagation, and provides a suitable approximation of the short transient signal corresponding to a phase arrival. The final test signal $S(t)$ is composed by superposition of a band-limited time-series of filtered white noise, having a flat amplitude spectrum in the frequency band of 0.5–30 Hz, and the wavelet $w(t)$. The SNR between the
<figure><img src="image_18.png" /><figcaption>**Figure A3.** (a) Synthetic signal used for the test composed of white noise filtered in the frequency range 0.5–30.0 Hz and Berlage wavelet (eq. A4.1). (b) Synthetic signal filtered in the frequency range 2.0–8.0 Hz. Red dashed line marks the onset of the Berlage wavelet. SNR for this frequency band is 3. (c) Berlage wavelet calculated using following constants: $A = 1$, $n = 2$, $\alpha = 20$, $\phi_0 = -\pi/2$ and $f_0 = 5$ Hz. The phase arrival time is $t_0 = 20$ s, representing the onset of the wavelet and marked by red dashed line. (d) Amplitude spectra of synthetic signal in panel (a). Grey dashed line marks the central frequency ($f_0 = 5$ Hz) of the Berlage wavelet.</figcaption></figure>

noise and the wavelet in the frequency range 2.0–8.0 Hz (around the predominant frequency of the wavelet) is set to a predefined value. We generate synthetic signals with SNR = [3.0; 2.0; 1.5; 1.0] for the test. This synthetic signals allow testing the efficiency of the time-frequency signal processing scheme using MBF approach, since there will be a narrow frequency band for which the signal of interest is significantly over the level of the noise. Fig. A3 illustrates the major features of the synthetic signal described above (assumed SNR = 3.0). The application of the MBF time-frequency signal processing scheme for extracting the information on timing ($t_0$) of the wavelet is shown in Fig. A4. The example compares the kurtosis CF calculated on the original time-series (Fig. A4e) and the broad-band kurtosis CF obtained through the MBF algorithm (Fig. A4d). Note the difference in kurtosis CFs at time $t_0 = 20$ s (onset time of the wavelet), and general improved SNR for the recomposed broad-band kurtosis CF. This illustrates the efficiency of MBF scheme in combination with HOS kurtosis CF for extracting the information on narrow-band transient signals even if the predominant frequency of the signal of interest is not known, and there exists at least a narrow band of frequencies for which signal of interest is over the level of the background noise.

We further assume a location of a test source (Fig. A5c) and generate simplified synthetic waveforms following the procedure described above for all stations. The Berlage wavelet will be recorded at station *i* at a time $t_{0i}$, corresponding to the traveltime from test source to the station. We consider a homogeneous constant velocity model (3.5 km/s) and S-wave propagation. Assumed value of the

threshold for imaging function used in declaring detection is 0.75. The methodology including MBF signal processing and space-time detection and location schemes is tested first on this single source location and synthetic signals with different SNR = [3.0; 2.0; 1.5; 1.0] (eq. A4.1). We could confirm that with the exception of the case with SNR = 1.0 the source was successfully located. This is mainly because broad-band kurtosis CFs have clear maxima around the arrival times of the signals generated by the test source. Retrieved maximum likelihood locations are consistent with the location of the test source within the limits not exceeding the grid spacing set to 1 km (Fig. A6). In the case of synthetic signals with SNR = 1.0, the maximum of the imaging function is below the detection level (Fig. A5d), thus no detection was triggered. Concerning the origin times, the absolute time shifts with the respect to the true origin time vary in the range 0.02–0.08 s, increasing with the decrease of SNR. An example of a retrieved location for synthetic time-series with SNR = 2.0 is shown in Figs A5(a)–(c). It should be noticed that the real location of the source corresponds to the likelihood (imaging function) exceeding value of 0.8. We compared as well these results with the locations obtained using kurtosis CF applied directly to the original synthetic records skipping the time-frequency MBF analysis step. In general, we could recover the location of the test source similar to those that used MBF kurtosis CFs, however with much smaller values of maximum stack and worse focusing (e.g. Max_stack = 0.61 for SNR = 2.0). Thus, having the maximum value of the imaging function below the trigger threshold, in the example shown in Figs A5(d) and (e), location of the source could not be determined.
<figure><img src="image_19.png" /><figcaption>**Figure A4.** Illustration of MBF algorithm on an example of synthetic signal in Fig. A3. (a) Traces filtered with filter-bank of narrowband filters covering the frequency range 0.5–45 Hz. (b) Original raw synthetic signal. Dashed line marks the onset of Berlage wavelet. (c) Kurtosis CFs of filtered signals. (d) Red line: broad-band kurtosis $CF-CF^{TF}(t)$ composed over all the frequencies according to eq. (9). Blue line: broad-band MBF CF, resulted from convolution of the kurtosis' (red line) positive derivative with the Gaussian window of $\sigma^2 = T_{decay}/2$ (eq. 11). (e) Same as panel (d), only directly on the original raw signal in panel (b).</figcaption></figure>

Finally, we tested the methodology on four test sources with different locations, assuming SNR = 2.0 for calculating synthetic time-series. The results are summarized in Fig. A6. The errors of recovered locations are around 570 m in horizontal directions, while in vertical they vary from 450 to 720 m. Errors in origin times are around 0.08 s. In general, the resolution of the retrieved locations is limited by the grid discretization (equal to 1 km in this case). The final refined locations (on the order of meters) are obtained interpolating the imaging function around the loca-

tion of the maxima by subdividing the 1 km sampled grid into smaller grid-cells. The most important point of the synthetic test is that in all of the above cases, the true location of the source lies in the area with the likelihood exceeding the value of the threshold.

The above tests confirm the robustness of our approach, and emphasize the advantage of the signal processing scheme combining the representation of the signal using HOS and time-frequency decomposition.
<figure><img src="image_20.png" /></figure>

**Figure A5.** (a)–(c) Test source detection and location example using our MBF signal representation and space–time detection and location schemes: (a) A horizontal and two vertical sections thorough the estimated location, corresponding to the maximum of imaging function shown according to a selected colour scale. White triangles indicate the location of Hi-net stations used in the test. (b) Transformed MBF broad-band kurtosis CFs (black lines) of the synthetic signals (grey lines) for the stations in panel (a). Vertical red bars show the theoretical arrival times for estimated location. Synthetic signal for each station is constructed following the scheme illustrated in Fig. A3, with SNR = 2.0, and the onset time of the Berlage wavelets corresponding to the traveltime from the assumed location of the source (at [−10.0, −30.0, 30.0]) and the stations. (c) Zoomed image of panel (a) comparing the imaging function (maxima provides the determined location of the source) and the actual assumed source location marked by white circle. (d,e) Same as panels (a) and (b) for the detection and location scheme without MBF analysis. The transformed CFs are calculated directly on synthetic waveforms.
<figure><img src="image_21.png" /><figcaption>**Figure A6.** (a) Map view comparing the determined locations (red points) recovered using the presented method (MBF signal processing and detection and location schemes), and the actual location of the test sources (blue dots). (b) Zoomed views of the map in panel (a) for each of the sources S1–S4. Distances between the recovered and actual locations of the sources in horizontal and vertical planes are shown in metre.</figcaption></figure>