Maximum covariance direction method for unconventional seismic sources

SUMMARY

We propose a new array covariance matrix analysis method, named ‘maximum covariance direction method’ (MCD method), to detect and locate unconventional seismic sources of weak signals without clear onsets. The MCD method builds a normalized-covariance matrix of time-shifted seismic waveforms recorded in a seismic array and determines the existence of source based on the consistency of the maximum covariance direction with the theoretical prediction. Synthetic tests demonstrate effectiveness of the MCD method in detecting one and multiple isolated sources with low signal-to-noise ratios. As a data application, we study 1-hr long-period tremors (LPTs) around Aso Volcano of Japan in 2014 November 24. A total of 26 LPTs are detected near the Naka-dake first crater of Aso Volcano, with the uncertainties of source location of about 7 km. Using the recorded background noise at the seismic stations, we show that the MCD method can detect LPTs even when the LPT signals are buried in the background noise and become indiscernible in the seismic data. Unlike traditional methods that employ the coherent features of seismic signals for source detection, the MCD method places emphases on both the coherence of seismic signals and consistency of the direction of the coherent signals from a potential source location. The synthetic tests and data application indicate that the MCD method provides a good alternative to other traditional methods for detecting and locating unconventional seismic sources, with a major improvement of avoiding source misidentification in the presence of strong incoherent signals.

Time-series analysis, Statistical seismology, Volcano seismology, Volcano monitoring

1 INTRODUCTION

Many types of natural and anthropogenic seismic sources do not register strong seismic signals with clear onsets. Such seismic signals include microseisms generated by the atmosphere–ocean system (e.g. storms over ocean) interacting with the solid Earth (e.g. Longuet-Higgins 1950; Bromirski et al. 2013; Gerstoft & Bromirski 2016), tectonic tremors originated from slow fault slipping at plate interfaces (e.g. Obara 2002; Rogers & Dragert 2003; Shelly et al. 2007), volcanic tremors related to volcano activities (e.g. Ogiso & Yomogida 2012; Wassermann 2012; Kawakatsu & Yamamoto 2015) and some anthropogenic seismic signals (e.g. Díaz et al. 2017; Li et al. 2018). The characteristic of those signals makes it difficult for the traditional methods to detect and locate the seismic sources.

Over the years, many array processing methods have been developed to detect and locate those sources based on the coherency of the signal across a seismic array (e.g. Schmidt 1986; Goldstein & Archuleta 1987; Goldstein & Archuleta 1991; Johnson & Dudgeon 1993; Rost & Thomas 2002), with the most notable ones being the beamforming method and the multiple signal classification (MUSIC) technique. Beamforming method includes the time-domain beamforming method (e.g. Rost & Thomas 2002) and the frequency-domain beamforming method (e.g. Rost & Thomas 2002; Gerstoft et al. 2006; Gerstoft & Tanimoto 2007). Specifically, the beamforming method defines a beamformer output and determines the source positions based on the maximal peaks of the beamformer output, with the time-domain method using the stack of the time-shifted waveforms for a hypothetic source location as the beamforming output and the frequency-domain method using the average of the elements in the cross-spectral density matrix of the waveforms as the beamformer output. The beamforming method has been used to locate the sources related to microseisms (e.g. Gerstoft et al. 2006; Gerstoft & Tanimoto 2007; Gerstoft et al. 2008; Zhang et al. 2009; Zhang et al. 2010; Liu et al. 2016; Nishida & Takagi 2016), tectonic tremors (e.g. Ghosh et al. 2009) and train-generated seismic signals (e.g. Li et al. 2018). The MUSIC method determines source positions based on the largest peaks of the spatial spectrum, defined as the inverse of the projection of the theoretical steering vector from a potential source location onto the noise subspace of the eigenvectors of the covariance matrix constructed from seismic recordings of an array (e.g. Schmidt 1986; Goldstein & Archuleta 1987). The MUSIC technique has been adopted to study volcanic tremors (e.g. Goldstein & Chouet 1994; Chouet et al. 1997; Saccorotti et al. 2004) and tectonic tremors (e.g. Ueno et al. 2010). In addition to the above two methods, several array covariance matrix analysis methods were proposed utilizing various properties of the matrix. For example, Seydoux et al. (2016) detected volcanic tremors in the Piton de la Fournaise volcano based on the spectral width of the covariance matrix eigenvalue distribution, Soubestre et al. (2018) identified clustered tremor sources based on the similarity of the first eigenvectors of daily array covariance matrices and the theory that the first eigenvector represents the most coherent part of the wavefield and may characterize the dominant seismic source (Wagner & Owens 1996) and Soubestre et al. (2019) determined the daily tremor location by the maximal value at zero lag-time of the stacked time-shifted cross-correlation envelopes between the first-eigenvector components over potential source locations.

In this study, we develop a new array covariance matrix analysis method, named the maximum covariance direction method (MCD method), to detect and locate unconventional isolated seismic sources of weak signals without clear onsets. The MCD method searches all possible sources in space and determines the existence of isolated seismic sources based on the consistency of the first eigenvector of the covariance matrix (i.e. the maximum covariance direction) with the theoretical prediction from the potential source locations. Unlike traditional methods that employ the coherent features of seismic signals for source detection, the MCD method places emphases on both the coherence of seismic signals and consistency of the direction of the coherent signals from a potential source location, and avoids misidentification of seismic sources. We present the method, theoretical justification and synthetic tests in Section 2, application of the method to the detection of the long-period tremors (LPTs) in Aso Volcano region of Japan in Section 3, and comparisons with other methods and difference of the covariance matrix construction in the time and frequency domains in Section 4.

2 MCD METHOD

The MCD method builds an |$N \times N\,\,$|normalized-covariance matrix of N time-shifted seismic waveforms recorded in a seismic array of N stations, and determines the existence of seismic sources based on the consistency of the maximum covariance direction with the theoretical prediction.

Seismic waveform recorded by station n is denoted as |${X_n}( t ).$| For each potential location |$( {x,y} )$| as a possible isolated source, |${T_n}( {x,y} )$| is the theoretical travel time from the potential source location |$( {x,y} )$| to station|${\rm{\,\,}}n$|⁠. The time-shifted waveform at station n can be written as

$$\begin{eqnarray*} {\bar X_n}\left( {x,y,t} \right) = {X_n}{\rm{}}\left[ {t + {T_n}\left( {x,y} \right)} \right]. \end{eqnarray*}$$

(1)

An |$N \times N$| normalized-covariance matrix |${R_{nm}}( {x,y,f} )$| at frequency f is constructed based on the seismic data in an array of N stations within an averaging time window L:

$$\begin{eqnarray*} {R_{nm}}\left( {x,y,f} \right) = \frac{{{{\langle {U_n}\left( {x,y,f} \right)U_m^\dagger \left( {x,y,f} \right)\rangle }_L}}}{{\left| {{{\langle {U_n}\left( {x,y,f} \right)U_m^\dagger \left( {x,y,f} \right)\rangle }_L}} \right|}}{}, \end{eqnarray*}$$

(2)

where |${U_n}( {x,y,f} )\,\,$| is the complex Fourier spectrum of |${\bar X_n}( {x,y,{\rm{t}}} ),$|m and n are the indexes of station pair, |${\langle \,\, \rangle _L}$| denotes the averaged Fourier cross-spectra matrices over M consecutive subwindows within the averaging time window L, |$\dagger $| denotes the conjugate operator and |$| {\,\,} |$| denotes the absolute value of a complex number.

Let |${\lambda _{\max}}( {x,y,f} )$| be the maximum eigenvalue of |${\rm{R}}( {{\rm{x}},{\rm{y}},f} ){\rm{\,\,}}$|and |${v_{\max}}( {x,y,f} )$| be the first eigenvector corresponding to the maximum eigenvalue |${\lambda _{\max}}( {x,y,f} ).$| We define a ‘source likelihood function’ |$C( {x,y,f} )$| for the potential source location |$( {x,y} )$|⁠:

$$\begin{eqnarray*} C\left( {x,y,f} \right) = \left| {{v_{\mathrm{ ref}}} \cdot {v_{\max}}\left( {x,y,f} \right)\,\,} \right|, \end{eqnarray*}$$

(3)

where vector |${v_{\mathrm{ ref}}} = [ {\frac{1}{{\sqrt N }}\,\,\,\,\,\,\frac{1}{{\sqrt N }}\,\, \ldots \,\,\frac{1}{{\sqrt N }}} ]\,\,\,\,$| is a ‘reference vector’, and |$\cdot $| represents a dot product of two vectors. We will show in Section 2.1, if the array signals are generated by only one isolated seismic source located at |$( {x,y} )$|⁠, all elements of |$R( {x,y,f} )$| are 1, and |${v_{\max}} = [ {\frac{1}{{\sqrt N }}\,\,\,\,\,\,\frac{1}{{\sqrt N }}\,\, \ldots \,\,\frac{1}{{\sqrt N }}} ]\,\,$| is the only eigenvector whose eigenvalue is non-zero. The source likelihood function |$C( {x,y,f} )$| thus provides a likelihood that the seismic data contain coherent signals from the potential source location |$( {x,y} )$|⁠. |$C( {x,y,f} )$| exhibits a peak at the true source location with a value of 1, and low values at other locations. In real seismic data because of possible presence of multiple isolated sources, background noise and inaccurate velocity model, |$C( {x,y,f} )$| is less than 1 even at the true source location.

We determine the source existence based on |$C( {x,y,f} )$|⁠. We search all the potential locations of isolated source in the averaging time window L. A source is detected when |$C( {x,y,f} ) > {C_{\mathrm{ cri}}}$|⁠, with the location of the source determined to be at the location |${\rm{\,\,}}( {x,y} )$| where |$C( {x,y,f} )$| is the maximum |${C_{\max}}$| and the uncertainty of the source region is determined to be the region where |${C}( {x,y,f} ) > 0.95 \cdot {C_{\max}}$|⁠. We name this method ‘maximum covariance direction method’ (MCD method) (Fig. 1).

Figure 1.

Workflow of the MCD method.

|${C_{\mathrm{ cri}}}\,\,$|is determined based on a criterion that the probability of noise being detected as a coherent signal is lower than a certain percentage (⁠|${\rm{0}}{\rm{.01\,\,per\,\,cent}}$| in this paper). In practice, we determine |${C_{\mathrm{ cri}}}$| based on the results of a certain number of repeated tests on noise. For example, for a target percentage of |${\rm{0}}{\rm{.01\,\,per\,\,cent\,\,}}$|of misclassifying noise as a coherent signal at frequency |${f_0}$|⁠, we perform 10 000 tests on noise and obtain 10 000 source likelihood functions|$\,\,{C_{\mathrm{ noise}}}( {{f_0}} )$|⁠. |${C_{\mathrm{ cri}}}{\rm{\,\,}}$| is set to be 1.2 times of the maximal value of those 10 000 |${C_{\mathrm{ noise}}}( {{f_0}} )$| functions, a value that would ensure the probability of detecting noise as a coherent signal is less than 1 out of 10 000 (⁠|${\rm{0}}{\rm{.01\,\,per\,\,cent}}$|⁠).

Note that only the first eigenvector is used to detect and locate the source in the MCD method. By adopting this approach, the MCD method ensures that the first eigenvector in the data, which represents the direction of most coherent energy in the data for the source location being tested, is consistent with the theoretical prediction from the test source location, and thus minimizes the possibility of source misidentification. At the same time, the method also preserves the capability of detecting multiple sources by only using the first eigenvector as the detection criterion. When the MCD method searches a true source location, the energy of other sources may appear in other secondary eigenvectors and is ignored in the first eigenvector search. However, when the source searching is at the true position of the other sources, the energy of the other sources would appear in the first eigenvector and be identified by the MCD method.

In the following subsections, we present theoretical background and synthetic validation of the MCD method for a single and multiple source detections. We also discuss the effectiveness of the MCD method in dealing with noise, inaccurate velocity model and sources with different durations.

2.1 MCD detection for one single isolated source

Let us start with a single isolated harmonic source of a frequency |${f_0}$| without noise, with the source time function represented as

$$\begin{eqnarray*} S\left( t \right) = A\left( t \right){\rm{exp}}\left[ {2\pi i{f_0}t + i\emptyset \left( t \right)} \right], \end{eqnarray*}$$

(4)

where |$A( t )$| is the source strength and |$\emptyset ( t )$| is the phase of signal, both randomly varying with time.

When the potential location |$( {x,y} )$| is at the true location of source |$( {{x_s},{y_s}} )$|⁠, the theoretical travel time |${T_n}( {x,y} )$| from source |$( {{x_s},{y_s}} ){\rm{\,\,}}$| to station n is equal to the observed travel time |${\tau _n}( {x,y} )$| from source |$( {{x_s},{y_s}} ){\rm{\,\,}}$| to station n. Thus, the time-shifted seismogram |${\bar X_n}( {{x_s},{y_s},t} )$| at station n is presented as

$$\begin{eqnarray*} {\bar X_n}\left( {{x_s},{y_s},t} \right) &=& {B_n}\left( {t - {\tau _n}\left( {{x_s},{y_s}} \right) + {T_n}\left( {{x_s},{y_s}} \right)} \right) \nonumber \\ && \times \, {\rm{exp}}[2\pi i{f_0} \cdot \left( {t - {\tau _n}\left( {{x_s},{y_s}} \right) + {T_n}\left( {{x_s},{y_s}} \right)} \right) \nonumber \\ && +\, i\emptyset \left( {t - {\tau _n}\left( {{x_s},{y_s}} \right) + {T_n}\left( {{x_s},{y_s}} \right)} \right)] \nonumber \\ &=& {B_n}\left( t \right){\rm{exp}}\left[ {2\pi i{f_0}t + i\emptyset \left( t \right)} \right] \end{eqnarray*}$$

(5)

where |${B_n}( t )$| is the product of source strength |$A( t )$| and the geometric spreading factor of the wave propagation at station n.

The complex Fourier spectrum of |${\bar X_n}( {{x_s},{y_s},t} )$| in the time window from t to |$t + \delta t$| at frequency |${f_0}$| is presented as |$\,\,{U_n}( {{x_s},{y_s},{f_0}} )$|

$$\begin{eqnarray*} {U_n}\left( {{x_s},{y_s},{f_0}} \right) &=& \mathop \int \limits_t^{t + \delta t} {B_n}\left( \tau \right)\exp \left[ {2\pi i{f_0}t + i\emptyset \left( \tau \right) - 2\pi i{f_0}t} \right]\mathrm{ d}\tau \nonumber \\ &=& \mathop \int \limits_t^{t + \delta t} {B_n}\left( \tau \right)\exp \left[ {i\emptyset \left( \tau \right)} \right]\mathrm{ d}\tau \nonumber \\ &=& {W_{n,j}} + i{X_{n,j}} \end{eqnarray*}$$

(6)

where |$\delta t$| is the length of the subwindow and is set to be twice of the time period at frequency f, j is the truncated integer value of |$t/\delta t$|⁠, representing the |${j{\mathrm{ th}}}$| subwindow ranging from 1 to M =|$\,\,L/\delta t$|⁠, |${W_{n,j}}$| and |${X_{n,j}}$| are the real and imaginary parts of the complex Fourier spectrum.

The element in the normalized-covariance matrix |${R_{nm}}( {{x_s},{y_s},f} )$| for stations m and n at frequency |${f_0}$| is given by

$$\begin{eqnarray*} {R_{nm}}\left( {{x_s},{y_s},{f_0}} \right) &=& \frac{{{{\langle {U_n}\left( {{x_s},{y_s},{f_0}} \right)U_m^\dagger \left( {{x_s},{y_s},{f_0}} \right)\rangle }_L}}}{{\left| {{{\langle {U_n}\left( {{x_s},{y_s},{f_0}} \right)U_m^\dagger \left( {{x_s},{y_s},{f_0}} \right)\rangle }_L}} \right|}} \nonumber \\ &=& \frac{{{S_{mn}}/M\mathop \sum \nolimits_{j = 1}^{j = M} \left( {{W_{n,j}}^2 + {X_{n,j}}^2} \right)}}{{\left| {{S_{mn}}/M\mathop \sum \nolimits_{j = 1}^{j = M} \left( {{W_{n,j}}^2 + {X_{n,j}}^2} \right)} \right|}} = 1, \end{eqnarray*}$$

(7)

where |${S_{mn}}$| is the ratio of the geometric spreading factor of the wave propagation at station m to that at station n.

At frequency |${f_0}$|⁠, eq. (7) is independent of station and |${R_{nm}}( {{x_s},{y_s},{f_0}} )$| are all equal to 1. In this case, there is only one positive non-zero eigenvalue |${\lambda _{\max}}\,\,( {{x_s},{y_s},{f_0}} ) = \,\,N$| with the associated eigenvector |$\,\,{v_{\max}}\,\,( {{x_s},{y_s},{f_0}} ) = \,\,{v_{ref}} = [ {\frac{1}{{\sqrt N }}\,\,\,\,\,\,\frac{1}{{\sqrt N }}\,\, \ldots \,\,\frac{1}{{\sqrt N }}} ]\,\,.$| We obtain the source likelihood function |${\rm{C\,\,}}( {x,y,{f_0}} ) = {\rm{\,\,}}1$| only at the true source location |$( {{x_s},{y_s}} )$|⁠. |${\rm{C\,\,}}( {x,y,{f_0}} ) = {\rm{\,\,}}1$| means that the coherent seismic signals are detected in the seismic array, and the seismic signals come exactly from the direction of source location |$( {{x_s},{y_s}} )$|⁠.

Our synthetic test set-up includes a synthetic single harmonic source located at |$( { - 70^\circ {\rm{E}},33^\circ {\rm{N}}} )$| and a seismic array of 72 seismic stations evenly distributed along a circle of |$5^\circ {\rm{\,\,}}$|centred at the source (Fig. 2a). The source has a single frequency of 0.17 Hz and lasts 3600 s (Figs 2b and c). Each element in the normalized-covariance matrix |$R( {x,y,{f_0}} )$| is calculated by stacking the cross-spectra of time-shifted seismic waveforms recorded by two stations in the seismic array with a subwindow |${\rm{\,\,}}\delta t\,\, = {\rm{\,\,}}12{\rm{\,\,s}}$|⁠. Using the same station set-up and processing procedure, we perform 10 000 tests on synthetic noise, and determine the source criterion |${C_{\mathrm{ cri}}}$| to be 0.46 based on 1.2 times of the maximal source likelihood function of the noise tests (Fig. 3).

Figure 2.

A synthetic single isolated harmonic source setting. (a) Location of a synthetic seismic source (red cross) and an array of 72 seismic stations (triangles) evenly distributed along a circle of |$5^\circ {\rm{\,\,}}$| centred at the seismic source with the seismograms at some example stations (red triangles) shown in panel (b). (b) 1-hr seismograms at 0.17 Hz in some example stations (red triangles in panel a), with the enlarged seismograms in the time window between red dashed lines shown in the bottom figure.

Figure 3.

Source criterion determination. Source criterion |${C_{\mathrm{ cri}}}$| (red line mark) is set to be 1.2 times of the maximal value of 10 000 |${C_{\mathrm{ noise}}}( {f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| functions (black dots) through 10 000 tests on noise.

We apply the MCD method to identify and locate the synthetic source. The source likelihood function |$C( {{{x}},{{y}},{f_0}} )$| exhibits one peak at the true source location with a value of 1 (Fig. 4a). The uncertainty of the source region determination is about 3 km, based on the region with the source likelihood functions larger than 95 per cent of the maximal value (Fig. 4b).

Figure 4.

Determination of a single isolated source. (a) Source likelihood functions |$C( {x,y,f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| as a function of potential source location |$( {x,y} )$| for a single isolated harmonic source with a signal duration of 3600 s (seismogram examples shown in Fig. 2b). Grey surface marks the value of |${C_{\mathrm{ cri}}} = \,\,0.46,$| and black triangles denote the 72 seismic stations whose synthetic seismograms are used to locate the source. (b) The source likelihood functions in the enlarged area of the blue box in (a). Black star indicates the position of detected source (i.e. location with the maximal source likelihood function) and black contour marks the uncertainty of the source region determination (i.e. 95 per cent of the maximal source likelihood function). Red crosses in (a) and (b) denote the position of true seismic source.

We now add noise |$N( t )$| to the seismic signals with different signal-to-noise ratios (⁠|$\mathrm{ SNRs}$|⁠). The |$\mathrm{ SNR}$| is defined as

$$\begin{eqnarray*} \mathrm{ SNR}\,\, = \frac{{\sqrt {\mathop \smallint \nolimits_0^L {{\left[ {{X^f}\left( t \right)} \right]}^2}\mathrm{ d}t} }}{{\sqrt {\mathop \smallint \nolimits_0^L {{\left[ {{N^f}\left( t \right)} \right]}^2}\mathrm{ d}t} }}, \end{eqnarray*}$$

(8)

where L is the averaging time window length of 3600 s, |${X^f}( t )$| is the signal |$X( t )$| at frequency f, and |${N^f}( t )$| is the random noise at frequency |$f.$| The source becomes undetectable by the method when SNR decreases to 1/13 and |${\rm{C}}( {{x_s},{y_s},{f_0}} )$| becomes lower than the source criterion |${\rm{\,\,}}{C_{\mathrm{ cri}}}$| (Fig. 5a). An example of detection at |$\mathrm{ SNR}\,\, = {\rm{\,\,}}1/12$| is shown in Fig. 5(b) and an example of non-detection at |$\mathrm{ SNR}\,\, = {\rm{\,\,}}1/13$| is shown in Fig. 5(c).

Figure 5.

Event detection of a single isolated source at various noise levels. (a) Source likelihood functions |$C( {x,y,f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| at the source location as a function of SNR for a single isolated harmonic source with a signal duration of 3600 s in Fig. 2(a). Red line marks the value of |${C_{\mathrm{ cri}}} = \,\,0.46$|⁠, red dots represent the cases of two noise levels with the maximal source likelihood functions shown in panels (b) and (c). (b and c) Left panels: same as Fig. 4(a), except for the source likelihood functions with |$\mathrm{ SNR}\,\, = \,\,1/12$| in (b) and |$\mathrm{ SNR}\,\, = \,\,1/13\,\,$|in (c). Right panels: same as Fig. 4(b), except for the enlarged area of the blue box in the left panels.

2.2 MCD detection for multiple isolated sources

Let us set q isolated harmonic sources of a frequency |${f_0}$| without noise located at|${\rm{\,\,}}( {{x_{s1}},{y_{s1}}} )$|⁠, …, |$( {{x_{sq}},{y_{sq}}} )$|⁠, with the source time function of source i at location |$( {{x_{si}},{y_{si}}} )$| represented as:

$$\begin{eqnarray*} {S_i}\left( t \right) = {A_i}\left( t \right){\rm{exp}}\left[ {2\pi i{f_0}t + i{\emptyset _i}\left( t \right)} \right], \end{eqnarray*}$$

(9)

where |${A_i}( t )$| and |${\emptyset _i}( t )$| are the source strength and phase of signal for the source, both randomly varying with time.

When the potential location |$( {x,y} )$| is at the true source location at |$( {{x_{si}},{y_{si}}} )$|⁠, the time-shifted seismogram at station n can be written as

$$\begin{eqnarray*} \bar{X}_n(x_{si}, y_{si}, t) &=& B_{ni}(t) {\rm exp}\left[2{\pi}if_{0}t + i {\emptyset_i} (t)\right] \nonumber \\ && +\, \mathop \sum \limits_{k\,\, = \,\,1}^q {B_{nk}}\left( {{\rm{t}} - {\tau _{ni}} + {T_{nk}}} \right) \nonumber \\ && \times \, \exp \left[ 2\pi i{f_0} \cdot \left( t - {\tau _{ni}} + {T_{nk}} \right) \right. \nonumber \\ && \left. + \, i{\emptyset _k}\left( {t - {\tau _{ni}} + {T_{nk}}} \right) \right] \,\,\left( {k \ne i} \right), \end{eqnarray*}$$

(10)

where |${B_{ni}}( {\rm{t}} )$| is the product of source strength |${A_i}( t )$| and the geometric spreading factor of the wave propagation from source i to station n, |${\emptyset _i}( t )$| is the phase of signal for source i, |${\tau _{ni}}$| is the observed travel time from source i to station n, and |${T_{nk}}$| is the theoretical travel time from source k to station n. |${B_{nk}}( {\rm{t}} )$| and |${\emptyset _k}( t )$| are the same definitions as |${B_{ni}}( {\rm{t}} )$| and |${\emptyset _i}( t )$| except for source k.

The first term in the right-hand side of eq. (10) is the time-shifted seismogram associated with source i, and the second term consists of the time-shifted waveforms associated with other |$q - 1$| sources. If source i is much stronger than other sources|${\rm{\,\,}}( {{A_i} \gg {A_k},k \ne i} )$|⁠, the normalized-covariance matrix |${\rm{R}}( {{x_{si}},{y_{si}},{f_0}} )$| is a nearly all-one matrix and the source likelihood function |${\rm{C}}( {{x_{si}},{y_{si}},{f_0}} )$| will be peaked with a value close to 1. If the relative source strength ratio of |${A_i}$| to |${A_k}$| decreases, |${\rm{R}}( {{x_{si}},{y_{si}},{f_0}} )$| would deviate from an all-one matrix and |${\rm{C}}( {{x_{si}},{y_{si}},{f_0}} )$| may still be peaked with a value smaller than 1. In that case, if the relative source strength ratio of |${A_i}$| to |${A_k}$| is not too small |${\rm{\,\,}}( {k \ne i} )$| and |${\rm{C}}( {{x_{si}},{y_{si}},{f_0}} ) > {C_{cri}}$| can be satisfied, source i can still be detected and located.

We perform a synthetic test using four isolated sources without noise, located at|${\rm{\,\,}}( { - 73^\circ {\rm{\,\,E}},33^\circ {\rm{\,\,N}}} ),{\rm{\,\,\,\,}}( { - 67^\circ {\rm{\,\,E}},33^\circ {\rm{\,\,N}}} ),$||$( { - 70^\circ {\rm{\,\,E}},36^\circ {\rm{\,\,N}}} )$| and |$( { - 70^\circ {\rm{\,\,E}},30^\circ {\rm{\,\,N}}} )$| with equal source strength (Fig. 6a). Using the station set-up and processing procedure in Section 2.1, we perform source detections in an averaging time window of 3600 s (Figs 6b and c). Four peaks of source likelihood function |${{C}}( {{{x}},{{y}},{f_0}} )$| satisfy the source criterion |${{C}}( {{{x}},{{y}},{f_0}} ) > {C_{\mathbf{ cri}}}$| and appear at the true locations of all four sources (Fig. 7a), with an uncertainty of the source region determination of about 3 km for each source (Fig. 7b).

Figure 6.

Four synthetic isolated harmonic sources setting. Same as Fig. 2, except for four synthetic isolated harmonic sources setting.

Figure 7.

Determination of four isolated sources. Same as Fig. 4 except for the four sources setting in Fig. 6.

We add random noise with different SNRs to the synthetic signals of four sources. As the SNR decreases, |${{C}}( {x,y,{f_0}} )$| decreases (Fig. 8a). |${{C}}( {{{x}},{{y}},{f_0}} )\,\,$| exhibits four peaks at the four true source locations when |$\mathrm{ SNR}{\rm{\,\,}} \ge 1/5$| (Fig. 8b for an example with |$\mathrm{ SNR}\,\, = \,\,1/5$|⁠), and becomes lower than |${\rm{\,\,}}{C_{\mathrm{ cri}}}$| at one of the true source positions with |$\mathrm{ SNR}\,\, = \,\,1/6$| and at all the source positions when |$\mathrm{ SNR}\,\, = \,\,1/9$| (Fig. 8c).

Figure 8.

Event detection of four isolated sources at various noise levels. Same as Fig. 5 except for the four sources setting in Fig. 6.

2.3 Effect of using inaccurate velocity model in MCD source detection

We consider the effect of using an inaccurate velocity model in source detection, taking a single isolated source as an example. We adopt the same source and station configuration in Section 2.1. When the potential location |$( {x,y} )$| is at the true location of source|${\rm{\,\,}}( {{x_s},{y_s}} )$| but with an inaccurate velocity model in source detection, we have the time-shifted seismogram at station n:

$$\begin{eqnarray*} {\bar X_n}\left( t \right) = {B_n}\left( {t + {e_n}} \right){\rm{exp}}\left[ {2\pi i{f_0} \cdot \left( {t + {e_n}} \right) + i\emptyset \left( {t + {e_n}} \right)} \right], \end{eqnarray*}$$

(11)

where |${e_n}$| is a random travel time error at station n due to the usage of an inaccurate velocity model in the detection, and is assumed to be bounded within |$\pm D \cdot T$| where T is the period of signal and D is a positive constant number representing the level of theoretical travel time error caused by the inaccuracy of the velocity model.

Synthetic tests indicate that the source likelihood function |$C( {{x_s},{y_s},{f_0}} )$| decreases rapidly with increasing D (Fig. 9a). The MCD method is capable of detecting the source when |$D \le 0.3{\rm{\,\,}}$| with |${{C}}( {x,y,{f_0}} ){\rm{\,\,}}$| greater than |${C_{\mathrm{ cri}}}$| (Fig. 9a) and the maximum of |${{C}}( {x,y,{f_0}} )$| appearing at the synthetic source position |$( {{x_s},{y_s}} )$| (Figs 9b and c).

Figure 9.

Event detection of a single isolated source with an inaccurate velocity model. (a) Source likelihood functions |$C( {x,y,f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| at the source location as a function of travel time error D (see the text for explanation) for a single isolated harmonic source with a signal duration of 3600 s in Fig. 2(a). Red line marks the value of |${C_{\mathrm{ cri}}} = \,\,0.46$|⁠, red dots represent the cases of using two inaccurate models with the maximal source likelihood functions shown in (b) and (c). (b and c) Left-panels: same as Fig. 4(a), except for the source likelihood functions with|$\,\,D\,\, = {\rm{\,\,}}0.3$| in (b) and |$D\,\, = {\rm{\,\,}}0.35$| in (c). Right panels: same as Fig. 4(b), except for the enlarged area of the blue box in the left panels.

2.4 Detection of sources with different durations

We apply the MCD method to detect synthetic single isolated harmonic sources (Fig. 2a) with different durations. We set three types of source with durations of 1800 s (Fig. 10a), 900 s (Fig. 10b) and 600 s (Fig. 10c), respectively. These sources repeat after a time window that is equal to their duration (top panels, Fig. 10). Using the same station set-up and processing procedure in Section 2, we obtain |${C_{\mathrm{ cri}}} = {\rm{\,\,}}0.46$|⁠, 0.45 and 0.47 for the sources with durations of 1800, 900 and 600 s, respectively.

Figure 10.

Detection of sources with durations of (a) 1800 s, (b) 900 s and (c) 600 s. In each subfigure, the panels from top to bottom are 1-hr seismograms at 0.17 Hz in some example stations (red triangles in Fig. 2a), the corresponding source likelihood functions at the source location in each average time window, the source likelihood functions as a function of potential source location |$( {x,y} )$|⁠, and the source likelihood functions at the source location as a function of SNR. Red line in the top panel marks the averaging time window used in the third and bottom panels. Red line in the second and bottom panels mark the source criterion under each signal duration. Black star, red cross and black contour in the third panels indicate the position of detected source, the true synthetic seismic source position and the detected source uncertainty (i.e. 95 per cent of the maximal source likelihood function), respectively.

In the absence of background noise, the method accurately detects sources and their repeated patterns, regardless of the source duration (second to third panels, Fig. 10). However, in the presence of background noise, the method can detect the source of a longer duration with a lower SNR. The detection limits of SNRs are |$1/11,\,\,1/9$| and |$1/5$| for the sources with durations of 1800, 900 and 600 s (bottom panels, Fig. 10), respectively. The seismic data with a source of longer duration contains a longer time span of coherent signals, making the source detectable with a lower SNR.

3 MCD METHOD DETECTION OF LPTS AROUND ASO VOLCANO, JAPAN

We apply the MCD method to detect the sources of LPTs around the Naka-dake first crater of Aso Volcano, Japan (Figs 11a and b). The volcanic tremors around Aso Volcano have been continuously observed since 1930s (e.g. Sassa 1935; Kaneshima et al. 1996; Kawakatsu et al. 2000; Kawakatsu & Yamamoto 2015). The sources of LPTs are interpreted to be composed of an inclined tensile crack with a strike nearly parallel to the chain of craters, located a few hundred meters southwest of the active crater at depths between 1 and 1.8 km (Yamamoto et al. 1999; Legrand et al. 2000).

Figure 11.

Long-period tremor (LPT) detection around Aso Volcano in Japan. (a) Aso Volcano (red dot) and 18 F-net seismic stations (black triangles) used in the study. (b) Elevation in the area marked by the white box in (a) based on the ASTER Global Digital Elevation Model v2. (c) Source likelihood functions |$C( {x,y,f\,\, = {\rm{\,\,}}0.067{\rm{\,\,Hz}}} )$| in the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24 (denoted by black line in Fig. 12b) in the area marked by the white box in (a). Green star and black contour indicate the position of detected LPT (i.e. location with the maximal source likelihood function) and the uncertainty of the source region determination (i.e. 95 per cent of the maximal source likelihood function). Brown crosses in (b) and (c) represent the position of the Naka-dake first crater of Aso Volcano.

We use the vertical components of ground velocity recorded by 18 broad-band seismic stations of the Full Range Seismograph Network of Japan (F-net) (Okada et al. 2004). The selected stations are distributed within distances smaller than 3° from the Naka-dake first crater (Fig. 11a). We adopt this example to illustrate the method application, as the event detection can be readily verified from the prominent seismic signals in the data and subsequent noise effects can be studied. Note that the LPTs around Aso Volcano can be clearly identified from the waveforms bandpass filtered between 10 and 30 s (Fig. 12a) and determined to be composed mainly of Rayleigh waves from their propagation velocity in the vertical component seismic recordings (Fig. 12a). Seismic data and spectra exhibit prominent energy around a period of 15 s (Fig. 12b, and an example station TKDF in Fig. 12c).

Figure 12.

1-hr vertical components of ground velocity and corresponding spectral characteristics. (a and b) 1-hr vertical components of ground velocity at 18 F-net seismic stations (black triangles in Fig. 11a), recorded from 15:00:00 to 16:00:00 UTC, 2014 November 24, filtered in a period band of 10–30 s in (a), and 13–17 s in (b). (c) Spectrogram of seismogram at station TKDF. Each seismogram in (a) is aligned based on the distance from the Naka-dake first crater of Aso Volcano to the seismic station, and the red arrows indicate LPT signals. In panel (b), red lines denote the 100 s averaging time windows of detected LPTs at a period of 15 s and black line marks the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24 in which the seismic recordings are used to locate the LPT sources in Fig. 11 and left panels of Fig. 16, blue line marks the 34th averaging time window from 15:55:00 to 15:56:40 UTC, 2014 November 24 in which the seismic recordings are used to locate the seismic sources in right panels of Fig. 16. Red line in (c) marks the frequency 0.067 Hz (corresponding to a period of 15 s).

We analyse 1-hr continuous data around a narrow period of 15 s from 15:00:00 to 16:00:00 UTC, 2014 November 24, as an example of data application. A typical LPT signal lasts less than 100 s (Fig. 12a). We choose a 100 s averaging time window. For each 100 s averaging time window, the seismic data are analysed by the MCD method with a subwindow|${\rm{\,\,}}\delta t\,\, = {\rm{\,\,}}30{\rm{\,\,s}}$|⁠. The source likelihood function |$C( {x,\,\,y,\,\,f} )$| is calculated at |$f\,\, = {\rm{\,\,}}0.067$| Hz (corresponding to a period of 15 s) over an area of longitudes between |$130.935^\circ {\rm{\,\,E}}$| and |$131.235^\circ {\rm{\,\,E}}$|⁠, and latitudes between |$32.735^\circ {\rm{\,\,N}}$| and |$33.035^\circ {\rm{\,\,N}}$|⁠, with a horizontal source searching grid size of |$0.003^\circ $|⁠.

We choose 1-hr seismic recordings started from 09:00:00 to 10:00:00 UTC, 1 January 2014 as the background noise. The seismic recordings in that time window are noise-like with no visible Rayleigh waves (recording at an example station YTYF in Fig. 13b). We use the selected background noise to determine the source criterion |${C_{\mathrm{ cri}}}$|⁠. Using the same station set-up and processing procedure in detecting the 1-hr LPT sources, we obtain the average of 36 maximal source likelihood functions at |$f\,\, = {\rm{\,\,}}0.067$| Hz as 0.61, extracted from the maximal source likelihood function at the 10 000 potential source searching grids for each 100 s averaging time window in the 1-hr background noise recordings. We thus determine the source criterion |${C_{cri}}$| to be 0.73 at 0.067 Hz based on 1.2 times of the average of maximal source likelihood functions of 0.61.

Figure 13.

Detection of scaled-down LPTs. (a) 1-hr vertical component ground velocity at seismic station YTYF (black triangle in Fig. 11a), recorded from 15:00:00 to 16:00:00 UTC, 2014 November 24, bandpass filtered in a period band of 10–30 s. (b) 1-hr background noise from 09:00:00 to 10:00:00 UTC, 2014 January 1 at YTYF, bandpass filtered in a period band of 10–30 s. (c–h) Left panels: seismograms obtained by multiplying the signal seismogram in (a) by a scaling factor S and superimposing it with the background noise in (b). S value and the maximal amplitude of the scaled-down seismogram are labelled at the left top of each trace. Right panels: the corresponding source likelihood functions |$C( {x,y,f\,\, = \,\,0.067\,\,{\rm{Hz}}} )$| and signal-to-noise ratios (⁠|$\mathrm{ SNRs}$|⁠) in each 100 s averaging time window with the red lines marking the source criterion of |${C_{\mathrm{ cri}}} = \,\,0.73$|⁠. The black lines mark the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24.

We detect 26 LPT sources around the Naka-dake first crater of Aso Volcano (Fig. 12b). A Rayleigh-wave velocity of 3.25 km s⁻¹ is used in the source detection. This propagation velocity is chosen based on the maximal source likelihood function and seismogram alignment of the time-shifted LPT arrivals in vertical components of ground velocity. We try propagation velocities ranging from 2.00 to 4.00 km s⁻¹ with an increment of 0.05 km s⁻¹, with the largest peak value of the source likelihood function obtained at a velocity of 3.25 km s⁻¹ (Supporting Information Fig. S1). The estimated source positions (green stars in Fig. 11c) are close to the crater (brown cross in Fig. 11c), and the uncertainty of the source region determination is about 7 km (area marked by black contour in Fig. 11c). Our method also avoids any detections in the seismic data of a teleseismic event that occurred elsewhere (Sandanbata et al. 2015), despite its strong signals around minute 55 (Fig. 12).

To test its effectiveness of the method in the presence of noise, we apply the MCD method to scaled-down signals with the background noise. We generate a new ‘data’ by multiplying the signals by a scaling factor and superimposing them with the background noise; we then apply the method to the new ‘data’ to perform event detection (Fig. 13). When the scaling factor is 1.0, the LPT signals in the new ‘data’ are clear and 20 LPTs are detected with |$SNR$| larger than 1.5 (Fig. 13c). With decreasing of the scaling factor (Figs 13c–g) and thus the SNR, the value of the source likelihood function |${{C}}( {x,\,\,y,\,\,f\,\, = {\rm{\,\,}}0.067{\rm{\,\,Hz}}} )$| decreases and lesser events are detected in the hour (Figs 13c–g) with similar uncertainties (Fig. 14). When the scaling factor decreases to 0.2, all LPT signals become indiscernible in the ‘data’, but the method still detects six LPTs with a minimum SNR of 0.47 (Fig. 13g). No events are detected when the scaling factor reaches 0.1 (Fig. 13h). These examples demonstrate the capability of the MCD method to detect seismic sources with signals lower than the background noise.

Figure 14.

An example of scaled-down LPT detection. Source likelihood functions |$C( {x,y,f\,\, = {\rm{\,\,}}0.067{\rm{\,\,Hz}}} )$| in the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24 (denoted by black lines in Figs 13c–h) with the scaling factors from (a) 1.0 to (f) 0.1. In each subfigure, the brown cross represents the position of the Naka-dake first crater of Aso Volcano, the green star indicates the detected LPT position and the black contour marks the uncertainty of the source region determination (i.e. 95 per cent of the maximal source likelihood function).

To further test the robustness of LPT source detection, we perform source detection on 1-hr synthetics with multiple 0.067 Hz LPT sources located at the Naka-dake first crater of Aso Volcano, using the same station configuration and procedure in the real data processing (Supporting Information Fig. S2). The synthetic sources can be resolved completely in the MCD method, with a similar location uncertainty of about 7 km (Supporting Information Fig. S3). We also scale down the synthetic signals with synthetic background noise to test the detectability of the source embedded in noise. With decreasing of the scaling factor (Supporting Information Figs S4c–g) and thus the SNR, the value of the source likelihood function |${{C}}( {x,\,\,y,\,\,f\,\, = {\rm{\,\,}}0.067{\rm{\,\,Hz}}} )$| decreases and lesser events are detected in the hour (Supporting Information Figs S4c–g) with similar uncertainties of 7 km (Supporting Information Fig. S5), resembling that in the real data processing. The above synthetic tests support the robustness of LPT source detections around Aso Volcano by the MCD method.

4 DISCUSSION

In this section, we compare the MCD method with the beamforming method (Rost & Thomas 2002; Gerstoft et al. 2006; Gerstoft & Tanimoto 2007), MUSIC (Schmidt 1986; Goldstein & Archuleta 1987) and maximal first-eigenvector method (for discussion purpose, we use this name for the method of Soubestre et al. (2018, 2019). We also discuss the practical difference between constructions of the covariance matrix in the time and frequency domains.

4.1 Comparisons of the MCD method with other methods

We discuss the comparisons of the MCD method with the beamforming method, MUSIC, and maximal first-eigenvector methods using the same source parameters and station configuration in Sections 2.1 and 2.2.

The time-domain beamforming method could detect one isolated source (left-hand panel, Fig. 15a), but fails to detect four isolated sources (right-hand panel, Fig. 15a). The frequency-domain beamforming method is able to detect all the true sources for either one isolated source (left-hand panel, Fig. 15b) or four isolated sources (right-hand panel, Fig. 15b), but with several small side lobes locating on the non-source regions as the method cannot suppress unreal sources with relatively high signal similarities in some of seismic stations (Fig. 15b). The MUSIC method could detect one isolated source (left-hand panel, Fig. 15c). However, for the four-source setting, while high peaked spatial spectrum emerges at all the four true source positions, relatively high spatial spectrum also appears in the non-source positions (right-hand panel, Fig. 15c). The maximal first-eigenvector method fails to detect the sources with a subwindow |$\delta t\,\, = \,\,12\,\,s$| used in the other methods (Fig. 15d), but could detect all the sources for either one isolated source or four isolated sources with a subwindow |$\delta t\,\, = \,\,300\,\,s$| (Fig. 15e).

Figure 15.

Source detection of the beamforming method, MUSIC and maximal first-eigenvector method for a single isolated harmonic source and four isolated harmonic sources. (a) The normalized beamformer output |${b_0}( {x,y,f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| in a time-domain beamforming method. (b) The normalized beamformer output |$b( {x,y,f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| in a frequency-domain beamforming method. (c) The normalized spatial spectrum |$P( {x,y,f\,\, = \,\,0.17\,\,{\rm{Hz}}} )$| in MUSIC, derived using the largest eigenvalue for the single source and four largest eigenvalues for the four isolated sources. (d) The network response at 0.17 Hz in the maximal first-eigenvector method. A subwindow |$\delta t\,\, = \,\,12\,\,s$| is used in (a)–(d). (e) Same as (d), except a subwindow |$\delta t\,\, = \,\,300\,\,s$| is used. In each subfigure, the red cross denotes the true seismic source position, the black star indicates the position of detected source (maximal peaks of the output) and the black triangles denote the 72 seismic stations whose synthetic seismograms are used to locate the sources.

A major improvement of the MCD method over the other methods lies on its capability to avoid misidentification of incoherent strong signals as a potential source existence. This is clearly demonstrated in the detection results of various methods in the time window of the teleseismic signals of the F-net stations in Section 3. All the other methods detect the LPT source as the MCD method does, with the detected LPT source positions distributed in less than 5 km apart between the methods (Fig. 16a). However, in the time window of teleseismic signals starting from minute 55, the MCD method detects no seismic event (top panel in Fig. 16b), while the beamforming and MUSIC methods misattribute the teleseismic signals to a source in a region less than 10 km away from the crater (second to fourth panels in Fig. 16b) and the maximal first-eigenvector method shows large values at zero lag-time of the stacked time-shifted cross-correlation envelopes in the source searching region (bottom panel in Fig. 16b).

Figure 16.

Comparison of the MCD, beamforming, MUSIC and maximal first-eigenvector methods for (a) an LPT event detection and (b) a teleseismic event misdetection around Aso Volcano. In each subfigure, the panels from top to bottom are the source likelihood function |$C( {x,y,f\,\, = \,\,0.067\,\,{\rm{Hz}}} )$| in the MCD method, the normalized beamformer output |${b_0}( {x,y,f\,\, = \,\,0.067\,\,{\rm{Hz}}} )$| in a time-domain beamforming method, the normalized beamformer output |$b( {x,y,f\,\, = \,\,0.067\,\,{\rm{Hz}}} )$| in a frequency-domain beamforming method, the normalized spatial spectrum |$P( {x,y,f\,\, = \,\,0.067\,\,{\rm{Hz}}} )$| in MUSIC derived using the largest eigenvalue and the network response at |$0.067\,\,{\rm{Hz}}$| in the maximal first-eigenvector method. The LPT signals and teleseismic signals are in 21.67–23.33 min (denoted by black line in Fig. 12b) and 55–56.67 min (denoted by blue line in Fig. 12b) during the 1-hr of 15:00:00–16:00:00 UTC, 2014 November 24, respectively. A subwindow |$\delta t\,\, = \,\,200\,\,s$| is used in the maximal first-eigenvector method, while a subwindow |$\delta t\,\, = \,\,30\,\,s$| is used for all the other methods. In each subfigure, the brown cross represents the position of the Naka-dake first crater of Aso Volcano, the green star indicates the detected seismic event position (maximal peaks of the output) and the black contour in the top panel of (a) marks the uncertainty of LPT source region determination in the MCD method (i.e. 95 per cent of the maximal source likelihood).

The MCD method determines the source existence based on the consistency of the maximum covariance direction with the theoretical prediction, while the other methods determine the source existence based on the maximal values of a particular coherent data feature. Thus, the MCD method places emphases on both the coherence of seismic signals and consistency of the direction of the coherent signals from a potential source location, while the other methods only place emphases on the coherence of seismic signals. When a coherent signal is observed for a real source, all methods are capable of identifying the source because all the data features from those methods would exhibit a maximal value at the source location. However, when some strong seismic signals from other sources appear in the seismic recordings (e.g. the teleseismic signals discussed above), the beamforming and MUSIC methods would misidentify the signals as a source at the non-source region based on the maximal values of coherent data features, leading to some misidentification of source. This is also true for the maximal first-eigenvector method that determines the source existence based on a data feature (characteristic of the first eigenvector) similar to the MCD method. It could misattribute the teleseismic signals to a local source or sources close to the crater based on the large values at zero lag-time of the stacked time-shifted cross-correlation envelopes (bottom panel in Fig. 16b). Besides, the maximal first-eigenvector method requires a relative longer subwindow than the other methods, as the subwindow needs to contain coherent signals across the seismic array (see an example of failed detection of the four isolated sources with a subwindow |$\delta t\,\, = \,\,12\,\,s$| in the right-hand panel of Fig. 15d). As a result, the maximal first-eigenvector method usually has poorer temporal and spatial resolutions than the other methods, as a long subwindow would not be able to separate sources that are close in time and the interference of signals from other sources in a long subwindow would deteriorate the coherence of the current source signals across the seismic array (see a LPT detection with a subwindow |$\delta t\,\, = \,\,200\,\,s$| in the bottom-panel of Fig. 16a when some signals of another LPT source are also present). Moreover, using cross-correlation envelopes could also be a contributor to the reduced spatial resolution of the maximal first-eigenvector method (e.g. the bottom-panel of Fig. 16a).

4.2 A note on constructing the covariance matrix in the time and frequency domains in the MCD method

Some studies suggested that seismogram alignment in the time domain before the construction of covariance matrix is an important pre-processing procedure in detecting and locating the seismic source (e.g. Goldstein 1988; Goldstein & Archuleta 1991). In the MCD method, the covariance matrix is also constructed by the time-shifted seismic recordings in the time domain. Mathematically, the covariance matrix could be constructed equivalently in the frequency domain, i.e. performing a phase-shifting to the data spectrum in the frequency domain. In practice, there are some shortcomings in the approach of the frequency domain. For a complete construction of the covariance matrix in the frequency domain, the averaging time window used for the construction is required to be long enough to encompass the signals of the source in all seismic stations. The frequency-domain construction will encounter problems when multiple sources occur within a time span that is shorter than the signal time span of the studied source across the seismic array. In that case, using a long averaging time window for the construction of the covariance matrix would include the signals from other sources that contaminate the current source signals, whereas using a short averaging time window would miss the signals of the studied source in some stations and yield artifacts or failure in source detection.

We show an example of location difference between constructing the covariance matrix in the time and frequency domains for a LPT signal in Aso Volcano region that spans a time window of 200 s across the seismic array (black line in Fig. 12b). The MCD method using the time-domain construction of covariance matrix can detect the LPT source with an averaging time window of 100 s (Fig. 11c), while the MCD method using the frequency-domain construction fails to detect the source as no sufficient signals are included in the averaging time window of 100 s (Supporting Information Fig. S6a). Increasing the averaging time window does not resolve the issue. With arrivals of another LPT signal starting to contaminate the current LPT signal in an increased averaging time window (Fig. 12), the MCD method using the frequency-domain construction fails to detect the current LPT source with a longer averaging time window (Supporting Information Fig. S6b).

5 CONCLUSION

We propose a new array covariance matrix analysis method, the MCD method, to detect and locate unconventional seismic sources of weak signals without clear onsets. The MCD method builds a normalized-covariance matrix of time-shifted seismic waveforms recorded in a seismic array, and determines the source existence based on a source likelihood function defined to measure the consistency of the maximum covariance direction with the theoretical prediction. A source is detected when the source likelihood function exceeds a source criterion that the probability of noise being detected as a coherent signal is lower than a certain percentage. Synthetic tests demonstrate effectiveness of the method in detecting one and multiple sources with low signal-to-noise ratios. As an example of data application, we employ the MCD method to detect the LPTs around Aso Volcano of Japan during one hour in 2014 November 24. A total of 26 LPT sources are detected near the Naka-dake first crater of Aso Volcano, with the uncertainties of source location of about 7 km. With the recorded background noise at the stations, we show that the MCD method is capable of detecting LPT sources even when the LPT signals are buried in the background noise and become indiscernible in the seismic data. Unlike traditional methods that employ the coherent features of seismic signals for source detection, the MCD method places emphases on both the coherence of seismic signals and consistency of the direction of the coherent signals from a potential source location. Synthetic tests and data application show that the MCD method provides a good alternative to other traditional methods for detecting and locating unconventional seismic sources, with a major improvement of avoiding source misidentification in the presence of strong incoherent signals.

SUPPORTING INFORMATION

Figure S1. Source likelihood functions |$C( {x,y,{{f\,\,}} = {{\,\,}}0.067{\rm{\,\,Hz}}} )$| at the source location as a function of assumed Rayleigh-wave velocity. Red line marks the value of source criterion |$C_{\mathrm{ cri}} = {\rm{\,\,}}0.73$|⁠. Detection is for the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24 (donated by black line in Fig. 12b). Rayleigh-wave velocity that produces the maximum source likelihood function is determined to be the best-fitting velocity (grey vertical dashed line, 3.25 km s⁻¹).

Figure S2. 1-hr vertical components of synthetic 0.067 Hz ground velocity at 18 F-net seismic stations (black triangles in Fig. 11a) from 15:00:00 to 16:00:00 UTC, 2014 November 24. Black line marks the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24 in which the synthetics are used to locate the synthetic LPT source in Fig. S3.

Figure S3. Synthetic LPT source detection around Aso Volcano in Japan. Same as Fig. 11(c), except that detection uses the synthetics in the 14th averaging time window from 15:21:40 to 15:23:20 UTC, 2014 November 24 (denoted by black line in Fig. S2) with a source criterion of |${C_{\mathrm{ cri}}} = \,\,0.84$|⁠.

Figure S4. Detection of scaled-down synthetic LPTs. Same as Fig. 13, except that the synthetic seismograms are used.

Figure S5. An example of scaled-down synthetic LPT source detection. Same as Fig. 14, except that the synthetic seismograms are used (denoted by black lines in Figs S4c–h).

Figure S6. LPT source detection around Aso Volcano by the MCD method using frequency-domain construction in (a) 21.67–23.33 min and (b) 21.67–25 min from 15:00:00 to 16:00:00 UTC, 2014 November 24. In each subfigure, the brown cross represents the position of the Naka-dake first crater of Aso Volcano, and the green star indicates the detected LPT source position.

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.

ACKNOWLEDGEMENTS

We thank the editor Lapo Boschi, Gilberto Saccorotti and an anonymous reviewer for their constructive reviews that have significantly improved the paper. This research is supported by the National Natural Science Foundation of China under grants NSFC41674045 and NSFC41874002, and by the Key Research Program of Frontier Sciences of Chinese Academy of Sciences (QYZDY-SSW-DQC020).

DATA AVAILABILITY

The National Research Institute for Earth Science and Disaster Resilience, Yamagata, Japan provides the F-net broad-band seismic data (accessible at http://www.fnet.bosai.go.jp/). The F-net broad-band seismic data are freely obtained using the Python package HinetPy (https://seisman.github.io/HinetPy/), and processed using the ObsPy (Beyreuther et al. 2010, Megies et al. 2011, Lion et al. 2015). The ASTER Global Digital Elevation Model v2 is provided by Japan’s Ministry of Economy, Trade, and Industry (METI) and U.S. National Aeronautics and Space Administration (NASA) at https://doi.org/10.5067/ASTER/ASTGTM.002. Figures are plotted using the Generic Mapping Tools (Wessel et al. 2013).