A seismogram records how the ground moves and usually contains a rich mix of different seismic signals. They may be transient such as the ground motions caused by earthquakes, surface processes such as glacier sliding (Weaver & Malone 1979;Lipovsky et al. 2019) or landslides (Keefer 1984;Weichert et al. 1994) and human activities such as cars, trains, ships and machinery from factories (Schippkus et al. 2020). They may be more diffuse such as the microseism (Cessaro 1994), the seismic hum (Rhie & Romanowicz Many studies across sciences have focused on separating the diffuse ambient data from the transient signals. In event-based seismology, the diffuse field is often considered 'noise'. Improving the signal-to-noise ratio (SNR) of transient signals is important to characterize better the source of these signals (earthquakes, landslides, etc.). Typically, the impulsivity of the transient signals is parametrized in order to measure the origin time of the events and the traveltime to the receivers. The impulsivity is strongly affected by noise levels, so the main task of denoising transient signals aims to improve the signal impulsivity. In continuous-based seismology, the transient signals are often considered 'noise'. Proper use of the ambient-seismic field requires that the ambient field is diffused and equipartitioned (Sánchez-Sesma & Campillo 2006;Tsai 2010), without transient events. Several studies have focused on reducing the effects of transient signals in the ambient seismic noise to improve the accuracy of correlation seismology. Typical pre-processing techniques such as sign-bit or amplitude clipping, spectral whitening (Bensen et al. 2007), and network-based equitization (Bowden et al. 2015;Seydoux et al. 2017) are applied to suppress large-amplitude signals. These steps alter the linearity of the cross-correlation (Fichtner et al. 2020).
The most commonly used denoising method in seismology is Fourier-based spectral filtering. This approach assumes that the earthquake and ambient noise signals are well separated in the frequency domain. Although this technique has been proven to be effective in numerous cases, it can sometimes distort the signals and/or cause artefacts around the impulsive signals (Douglas 1997;Mousavi & Langston 2017). Moreover, earthquake and ambient noise signals often overlap in the frequency domain, and direct filtering may be challenging to separate them (Peterson et al. 1993;Mc-Namara et al. 2019). Denoising using time-frequency representations of the signal is another widely applied and effective technique to separate the earthquake and ambient noise signals when they overlap. Many innovative algorithms and methods have been developed, for example, using time-frequency transforms such as the Stockwell S-transform (Stockwell et al. 1996), the Radon transform, the wavepacket transform, the continuous wavelet transform, or others using f-x or f-k filtering, singular spectrum analysis, sparse transformbased denoising, which are extensively reviewed and discussed in Mousavi & Langston (2017). Most of these transform-based denoising methods achieve noise suppression through thresholding methods, that is, determining some hard- (Donoho & Johnstone 1994) or soft- (Chang et al. 2000) thresholds to separate seismic and noise signals. Although those transform-based methods are shown to be very practical and possible to be automated (Mousavi et al. 2016;Mousavi & Langston 2017); they still require manual intervention, and the parameter tuning is often performed using trial-and-error approaches.
With the recent leap of computational power, memory and data storage, machine learning (ML) has provided a diverse set of powerful tools in the geosciences (Bergen et al. 2019). Many ML algorithms are built to (1) automatically perform complex prediction task; (2) create a representation that approximates numerical simulations or captures relationships and (3) reveal new patterns, structures, or relationships from data (Bergen et al. 2019). ML algorithms are powerful in many different seismological tasks, including but not limited to waveform classification and earthquake detection (Li et al. 2018;Perol et al. 2018;Kong et al. 2019;Mousavi et al. 2019bMousavi et al. , 2020;;Beroza et al. 2021;Johnson et al. 2021), phase picking and association (Li et al. 2018;Meier et al. 2019;Zhu et al. 2019a;Liu et al. 2020;Mousavi et al. 2020;Walter et al. 2021), source location and characterization (Perol et al. 2018;Mousavi & Beroza 2019;Ren et al. 2020;van den Ende & Ampuero 2020;Kuang et al. 2021;Münchmeyer et al. 2021a;Zhou et al. 2021), earthquake early warning (Li et al. 2018;Münchmeyer et al. 2021b) and many others.
ML methods have also been developed to denoise and decompose the seismic data (Chen et al. 2016;Siahsar et al. 2017;Saad et al. 2018;Zhu et al. 2019b;Tibi et al. 2021). Saad et al. (2018) develop a stacked denoising auto-encoder to smooth the ambient noise and output a time-series for better arrival time detection, in a similar fashion to an impulsive filter (Allen 1978). Zhu et al. (2019b) develop a DeepDenoiser network, which applies deep learning to the short-time Fourier transform and focuses on the classification on pixels of the spectrogram of seismograms (Zhu et al. 2019b;Tibi et al. 2021). Recent work by Novoselov et al. (2022) shows that it is possible to separate signals in the time domain using a dual-path recurrent neural network.
This study aims to separate the waveforms of both earthquake and ambient noise signals through the ML network, thereby benefiting both earthquake and ambient-noise seismology. We develop a new multitask encoder-decoder ML network, which we name WaveDecompNet, to separate both types of signals simultaneously. We treat the problem as a time-series extrinsic regression (Tan et al. 2021) and directly work on the seismic data in the time domain. This setting can save human intervention by not tuning parameters for the time-frequency representation and provides potential for a near-real-time framework, given its lower computation cost. Our results show that the performance of the WaveDecompNet is notable even at a low SNR. Although both are working in the time domain, the architecture of WaveDecompNet is different from that proposed by Novoselov et al. (2022), and we expand beyond their work to apply WaveDecompNet to practical cases in seismology.
We take the Big Island of Hawai'i as a natural laboratory of the seismically complex environment. Oceanic islands record strong microseismic signals, which are the basis for ambient noise seismology (Longuet-Higgins 1950). The volcanic region also exhibits dynamic tectonics activities. We choose the broad-band seismic station IU.POHA at Pohakuloa, Hawai'i. The richness of this dynamical system presents particular challenges in monitoring the volcano-tectonic activities and the temporal evolution of structural changes from ambient-noise seismology. We apply WaveDecomp-Net to 1-month-long continuous 3-component seismic data recorded by IU.POHA, and show that WaveDecompNet can successfully separate the earthquake and ambient-noise waveforms. Many small local earthquakes, which were buried in the noisy raw waveforms, can be easily detected in the separated waveforms.
We further test the applicability of our technique by performing two standard single-station measurements. We first test the improvements in detecting transients by applying the standard short-time average (STA) through long-time average (LTA) method (Allen 1978) and find increased trigger rates and potentially improved picking accuracy. We also measure the coherence of single-station ambient-noise cross-correlation functions and find increased stability in the coda of these functions. Our results highlight the potential for WaveDecompNet to improve seismograms used in both eventand continuous-based seismology.
We design a multitask encoder-decoder network, WaveDecompNet, to handle the seismic data directly in time domain. The input to the network is 1-min-long, 3-component (east-west, north-south, and up-down) raw seismograms. The output of the network is two 1-min-long, 3-component seismograms as earthquake and noise (Fig. 1). In order to seek the generality of applications, we only apply minimum pre-processing, including demean and amplitude normalization of the waveforms.
WaveDecompNet is an encoder-decoder network, which is a popular network design in ML problems in generating dialogues (Serban et al. 2017), labelling semantic image (Badrinarayanan et al. 2015(Badrinarayanan et al. , 2017)), detecting image forgeries (Bappy et al. 2019) and predicting vehicle trajectory (Park et al. 2018). In seismology, many previous studies applied a similar network architecture, for instance, seismic phase picking (Mousavi et al. 2020), or denoise seismic signals in time-frequency domain (Zhu et al. 2019a). In this study, our WaveDecompNet consists of three major parts:
(i)The encoder branch takes the role of extracting useful, highlevel features from the seismic time-series. Through training with sufficient data and updating its parameters, the encoder aims to learn features of the input data that can help characterize the earthquake and ambient noise signals. We use 1-D convolutional layers with an increasing number of kernels to extract high-level features with a minimal number of parameters (Fig. 1). The stride of the convolution is adjusted to downsample the time-series along the time axis. We tested using of MaxPooling instead of convolutional strides but found poorer network performance. After each 1-D convolutional layer, batch normalization is applied to normalize the output to zero-mean and unit variance. Finally, a rectified linear unit function (ReLU) is used as the activation function for non-linearity. (ii)The decoder branch translates the learned features from the encoder branch and reconstructs the separated earthquake and noise time-series. A two-branch decoder block handles both the separated earthquake and noise waveforms individually. It performs better than a single branch network, which only outputs the earthquake signal and assumes that the noise is recovered from the difference between the original signal and the reconstructed earthquake signals. The branches are composed of 1-D transpose convolutional layers. In symmetry with the encoder block, the number of kernels gradually decreases, and we manually select the stride to incorporate the high-level features back into the time domain. We also apply batch-normalization and ReLU activation following each 1-D transpose convolutional layer like the encoder block. The parameters of the two branches are learned independently. (iii)The bottleneck blocks link the encoder and decoder branches. Their purpose is to learn the mapping relation between the encoderextracted features of the composite waveform (earthquakes plus noise) and the features of the separated earthquake and noise timeseries, respectively. The design of the bottleneck block greatly impacts the performance of the algorithm and is subject to investigation in this study.
To train the WaveDecompNet, we need two data sets: (1) earthquake and (2) ambient noise waveforms. Combining both, we can use the 'noisy' waveforms as the input of WaveDecompNet, and train the network towards the corresponding earthquake and noise waveforms (Fig. 1).
(i)The data set for earthquake waveforms comes from the STandford EArthquake Dataset (STEAD, Mousavi et al. 2019a) because of its broad coverage of global earthquakes. This data set was curated to provide many high signal-to-noise ratio waveforms of local (source-receiver distance less than 350 km) earthquakes and a set of 'noise' diffuse signals recorded globally. There are 234 526 samples of 3-component seismograms of ambient noise and 1 030 231 samples of 1-min 3-component seismograms associated with 450 000 earthquakes located in various regions of the world (Mousavi et al. 2019a). All the seismograms in the data set are raw waveforms where the instrumental response is preserved. We only keep a subset of STEAD earthquake data set with the highest quality (SNR >10 4 ) to approximate the noise-free ground-truth earthquake signal. The SNR in this study is defined as:
where ||S|| 2 and ||N|| 2 are the power of signal and noise, respectively. After this event filtering, we mostly have M2+ earthquakes in the STEAD data set. From this subset, we randomly picked 187 000 samples of STEAD earthquake waveforms with replacement, then lowpass-filter the waveform below 5 Hz and downsample them from 100 to 10 Hz. This high-frequency cut-off still provides sufficient bandwidth for earthquakes of magnitude 2 and above. Increasing the sampling rate can lower the resolvable magnitude range of earthquakes but also significantly increases the computational expenses.
In this study, we focus on comparing with a catalogue of earthquakes of magnitude 2 and greater in the later application to the continuous seismic data.
(ii)The data set for ambient noise waveforms combines a 'global noise' from STEAD noise waveforms with a 'local noise' from IU.POHA station. We randomly pick 100 000 3-component STEAD samples, low-pass filtering first, and then down-sample them to 10 Hz. For the 'local noise', we select noise waveforms from 2month-long continuous seismic data recorded by IU.POHA (from 1 June 2021 to 31 July 2021), and downsample the data to 10 Hz. The continuous data may contain known and unknown earthquakes. We remove the time windows containing large transient signals as follows. First, we calculate the amplitude median of the 3-component waveform. Then, we remove all the time windows with amplitudes larger than six times the amplitude median. By visual inspection, we confirm that all of the large transient signals can be removed while the features of ambient noise can be kept. Finally, we obtain 87 000 samples of 1-min ambient noise time-series.
To represent more time-series for better generalization, we apply two strategies of data augmentation (Zhu et al. 2020). First, we shift the earthquake signals with a random time from a uniform distribution of -30 to 60 s to handle the uncertain arrival time of earthquake signals in the time window during real application and allow for a redistribution of the weights in the encoder 1-D convolutional layers. We truncate the signals outside 0-60 s time window after shifting and fill with zeros to keep the constant window length of 1 min/600 samplings. Secondly, we randomly scale the amplitude of earthquake waveforms related to the noise signals to increase the range of SNR from 10 -2 to 10 4 . We linearly stack the shifted and scaled earthquake and ambient noise waveforms to generate the composite waveforms. We report that the resulting SNR of our training data is uniformly distributed in the log space.
We split the entire data set containing 187 000 samples of 3component waveforms into the training, validation and testing data sets with the proportion of 60, 20 and 20 per cent, respectively. Only The network consists of five main blocks: the encoder branch, two bottlenecks, and two decoder branches. The encoder and decoder branches contain seven 1-D convolutional and transpose convolutional layers. The layer parameters 'x kr y, stride=z' refer to x kernels with y features and stride of z. Each convolutional or transpose convolutional layer is followed by a batch normalization (BN) layer and a ReLU activation layer. The structure of the bottleneck block is tested and discussed in the main text for details. Six residual connections (skip-connection layer by summation) directly connect the encoder to the decoder to improve the convergence of training and prediction performance (Ronneberger et al. 2015;He et al. 2016;Zhu et al. 2019a). The red numbers in the brackets show how the tensor size varies between layers [number of batches × number of channels × length of signal sequence].
the training data set is used to train the network, update the model parameters, and minimize the loss function. The validation data set is used to validate the model and track over-fitting during training. The testing data set is used to evaluate the performance of models after training.
WaveDecompNet has a flexible structure and can utilize various types of bottleneck blocks. The bottleneck block aims at learning the mapping relation between features in the encoder and decoder, and those features are necessary to reconstruct the separated signals in the decoder blocks. There are multiple choices for the bottleneck structure in time-series analysis. In this study, we explore five of them and evaluate their impacts on model performance to choose the best structure for WaveDecompNet using the test data set:
(i)None: no specified bottleneck. The encoder and decoder are directly connected. The total number of trainable parameters in the network is 78 090. (ii)Linear: a linear regression layer between the encoder and the decoder. The total number of trainable parameters in the network is 86 410. (iii)LSTM: a bidirectional long-short-term memory (LSTM) layer between the encoder and the decoder (Hochreiter & Schmidhuber 1997). The total number of trainable parameters in the network is 178 442. (iv)Attention mechanism: a multihead attention layer between the encoder and the decoder (Vaswani et al. 2017). Position embedding is used to keep the sequence order information (Vaswani et al. 2017). We use a 4-head dot-product self-attention layer with a dropout probability of 0.2. Other numbers of heads were tested but did not significantly affect the results. The total number of trainable parameters in the network is 110 858. (v)Transformer: the standard transformer encoder layer made up of self-attention and feed-forward network (Vaswani et al. 2017). The transformer model has been shown to be a powerful tool in different seismological applications such as earthquake detection and phase picking (Mousavi et al. 2020), earthquake source characterization (Münchmeyer et al. 2021a) and early warning (Münchmeyer et al. 2021b). We only use one layer but find that adding more layers can significantly downgrade the model performance, which we attribute to insufficient training. The total number of trainable parameters in the network is 640 394.
We use the identical training and validation data set to train all these networks. We train and validate over the batch size of 128. We also test a batch size of 256 but found little change in the results. We require each network to be trained for 100 epochs. The learning rate quantifies how much the model parameters are updated in each epoch during training. We find that it impacts the stability of the model training. We explore different schemes to assign a learning rate and find that a variable learning rate with warm-up followed by an exponential decay with respect to epochs improves the stability of the training (Fig. S1). In the warm-up phase of the first 10 epochs, we linearly increase the learning rate from 10 -4 to 10 -3 . This learning rate warm-up has been shown to limit parameter divergence in training deep neural networks (Gotmare et al. 2018). After 10 epochs, the learning rate begins to decay with a power law of exponent of 0.95 (learning rate = 10 -3 × 0.95 epoch ) till the 50th epoch. Finally, we keep a constant learning rate of =10 -3 × 0.95 50 = 7.70 × 10 -5 until the terminating 100th epoch. The specific values of these steps are chosen by trial-and-error. We further inspect the training curves of each model to guarantee the stable convergence of each model's optimization (Fig. 2). We have also tested the training with an early stopper, which automatically terminates the model training if the validation loss has not been improved for a few training epochs (we choose 10 in this study). The final test loss values of models with early stopping are similar (Fig. S2) to those with fixed 100 training epochs shown in Fig. 2.
Normalization of the input data is necessary to stabilize the optimization. We use a standard scaler to normalize the time-series to zero-mean and unit variance. The loss error is the mean square error (MSE) function, and the training is performed using the Adam optimizer (Kingma & Ba 2014). Finally, we rescale the earthquake and noise signals using that same factor after the two decoder branches of WaveDecompNet. The overall training behaves properly for all models because the training loss curve keeps decreasing. The validation loss curve decreases and approaches the training loss curve as the training goes, and finally stays almost constant at some epoch, which indicates the convergence of optimization when the model is no longer improved.
LSTM and attention models achieve the minimum final loss value for both training and validation data sets (Figs 2c and d). The transformer model, although with more model parameters than any other model, shows a higher validation loss (Fig. 2e). The minimal None and Linear bottlenecks exhibit the highest final training loss values (Figs 2a and b). We also show the partial loss curves from individual branches of the earthquake and the noise waveforms, approximately half of the total loss.
After training, we test the trained models with the 37 400 samples of the test data set. The test data set is not included during the training process, so it can be used to evaluate the model performance. Figs 34show two examples with different types of noise signals: Fig. 3 contains the regional noise from the IU.POHA station, which is energetic mostly in low-frequency below 1 Hz and characteristic of microseism at a dominant period of 7 s, Fig. 4 contains the STEAD noise signal, which is broad-band and rich in high frequencies in a band that overlaps with the earthquake signal. Visually, the earthquake waveforms are relatively well recovered over a broad range of frequencies. In particular, it is able to decompose the signals with overlapping frequency content (Fig. 4b), which is often a challenge for filtering-based denoising methods (Douglas 1997;Mousavi & Langston 2017).
We also compare the test losses for different structures. The LSTM and Attention models have achieved the minimum mean test loss value of 0.0439 and 0.0452, respectively. The Transformer model has a mean test loss value of 0.0672. The None and Linear models have mean test loss values of 0.0717 and 0.0713, respectively. The distribution of test loss values is also shown in Fig. S3.
Furthermore, we inspect the waveform fitting for different models. For the same input X (composite waveform), we obtain the predicted output waveform ỹ and compare it with the ground truth waveform y. Other seismic denoising studies have reported improved SNR values as performance metrics (Zhu et al. 2019a;Tibi et al. 2021). Tibi et al. (2021) also use the signal-to-distortion ratio (SDR) metric (Nakajima et al. 2018), but the SDR metric is unbounded. Here, we calculate the explained variance (EV) score for both separated earthquake waveforms and ambient noise waveforms: EV = 1 -Var(y-ỹ) Var(y) , Var means variance of the time-series. The best possible EV score is 1.0, corresponding to perfect waveform reconstruction. An EV score of 0.0 means that no waveform has been reconstructed (ỹ = 0). A negative EV score means a false waveform reconstruction, for example, in the time window where there is no earthquake waveform (y ≈ 0), but the network reconstructs a spurious waveform (ỹ = 0).
The results are compared between bottleneck architectures in Fig. 5. All models can reconstruct both earthquake and noise, with over half of the tested samples that achieve a high EV score of around 1 (Fig. 5). The network with the LSTM bottleneck recovers most test samples with a high EV score of around 1 for the earthquake waveforms (Fig. 5a). There is a bimodal distribution in the EV scores for earthquake waveforms. All models show two peaks around EV scores of 1 and 0, especially when the earthquake-to-noise ratio is less than 1. The bimodal pattern is suggestive of the behaviours of this encoder-decoder network. In cases the network can recognize, the network can reconstruct the waveforms accurately. In cases the network can hardly recognize, the network tends to output a timeseries of zeroes rather than an artificial or a spurious waveform. The null time-series has an EV score of zero. All networks show similar performance in EV score for the noise waveforms (Fig. 5b). We find no obvious bimodal pattern in the ambient noise reconstruction, which indicates a lower likelihood of outputting strictly zero noise. There are, however, spurious reconstructions of the noise waveforms for lower noise-to-earthquake amplitude ratios or when transient signals dominate the time-series (Fig. 5b).
Moreover, we explore how the EV score varies with SNR for both earthquake (Fig. 5a) and ambient noise waveforms (Fig. 5b). First, all models present the same pattern: the EV score monotonically increases with the signal amplitudes, quantified by SNR for Figure 3. Example 1 of waveform separation using one of the bottleneck architecture (LSTM, see Section 2.3). (a) 3-component (E, east; N, north; Z, vertical from top to bottom) velocity seismograms normalized with the same scaling factor of maximum amplitude in each component. Left-hand panels: one STEAD earthquake waveform (in red) and IU.POHA local noise is stacked to get the noise input waveform (in black). Middle panels: comparison between the separated earthquake waveforms (blue) with the ground truth earthquake waveform (red). Right-hand panels: comparison of the separated noise waveform (blue) and ground truth noise waveform (red). (b) 3-component waveform Fourier amplitude spectra. Top panels: the spectrum of the input waveform is shown in black, with the ground truth earthquake spectrum (in red), and the separated earthquake spectrum (in blue). Bottom panels: the ground truth noise spectrum is shown in red, and the separated noise spectrum is shown in blue. earthquake waveforms or 1/SNR for ambient noise waveforms. This is expected because it is easier to learn the intrinsic features of the signal waveforms and reconstruct the earthquake and ambient noise signals if they have high amplitude in the original input waveform. All models perform similarly. Take the case of earthquake waveform as an example (Fig. 5a). For SNR >10 1 , all models can correctly separate the earthquake waveform almost perfectly (EV score ∼ 1). The performance of any model drops as the SNR decreases. For example, at SNR = 1, the median EV score of different models varies from 0.8 to 0.9, and LSTM and Attention have the best performance. The discrepancy between model performance is exacerbated at low SNR. For instance, it is visually difficult to extract the earthquake signal when the SNR is small, =10 -1 = 0.1. However, the LSTM and attention models can still achieve a median EV score of about 0.75, the None and Linear models have a median EV score of about 0.65, and the Transformer model has a median EV score of 0.60. For even smaller SNR values, =10 -1.8 = 0.02, all models tend to fail, with most EV scores close to 0, except the LSTM and Attention models still have a median EV score of about 0.35. The variance of the EV scores also increases with decreasing SNRs, indicating more uncertainties in the reconstructed waveforms.
Similar behaviours can be observed for the EV score of the ambient noise part, with larger amplitudes of noise yields to better model performance (Fig. 5b). However, we note some different features of ambient noise from the earthquake waveforms (Fig. 5a). For instance, LSTM works better than other models at a low noise-tosignal ratio. When the noise amplitude is low, LSTM returns mostly zero-sequence as noise (EV score ∼ 0), while other models produce artefacts with negative EV scores. Moreover, the Transformer model performs slightly better than the None and Linear models, especially at low SNRs.
To summarize, we evaluate the performance of models with different bottleneck models by applying them to the same testing data set. In general, the order of performance of our model is LSTM ≥ Attention > Transformer > None ≈ Linear based on the variation of EV score with SNR. The LSTM bottleneck outperforms other bottlenecks, slightly better than the Attention model, in reconstructing both earthquake and noise waveform, especially for a situation Downloaded from https://academic.oup.com/gji/article/231/3/1806/6651381 by University of Washington user on 03 November 2022 with low SNR (Fig. 5). In the following applications to the real data, we only focus on the WaveDecompNet with the LSTM bottleneck.
We now apply the trained WaveDecompNet to continuous timeseries. It is straightforward to apply the model to any continuous data, provided that it has the same sampling rate. We select continuous recordings at IU.POHA from 31 July 2021 to 1 September 2021. We first down-sample the three-component 1-month-long waveforms to 10 Hz. The typical data pre-processing steps to apply machine-learning models are (i) applying a sliding time window to obtain short time-series (1-min long, 600 samples here) and (ii) applying the data normalization by removing the mean and normalizing it with the standard deviation (STD) to have the zero-mean, unit variance time-series. The most intuitive order to apply these processing steps is (i), then (ii). We found that ordering (i) then (ii) leads to spurious effects when concatenating the 1-min waveforms into a 1-month-long waveform due to offset (means) and trends that rendered the application to continuous time-series unpractical. Instead, we experimented with the order of (ii) then (i) and found much better performance with few artefacts when stitching back the waveforms. Then we slide through the data with 1-min-long windows without overlap. We have also tested the sliding time window with overlap but found that multiple values at one time stamp often bring discontinuity and artefacts. Finally, we apply the WaveDecompNet to all 1-min-long windows, concatenate all ML-filtered windows, and scale back the 1-month-long time-series with the standard deviation and mean for both the earthquake/transient time-series and the noise time-series.
We show the results of decomposed transients and noise waveforms in Fig. 6. Most of the transient signals are well above the suppressed noise level (Fig. 6b). Furthermore, we validate that the new transient signals correspond approximately to real earthquakes from the regional M2.5+ earthquake catalogue within 5 • great arc distance to the station in the same month. The earthquakes are from the Hawaii Volcano Observatory (HOV) catalogue downloaded from the IRIS FDSN event server using Obspy (Beyreuther et al. 2010). We use the TauP module with IASP91 velocity model (Kennett & Engdahl 1991) to estimate the arrival time of these local events (green squares in Figs 6a and b). We find that many of these small local earthquakes, which have been hidden by the large ambient noise signals (Fig. 6a), can be separated and shown in the decomposed earthquake waveforms (Fig. 6b). A zoom-in example is also shown in Fig. 7. In the next section, we further test the validity and usefulness of these transient waveforms by applying an impulsivity filter that is standard for earthquake detection. The reconstructed noise waveforms exhibit more levelled, constant amplitudes throughout the month (Fig. 6c). This implies a persistent distribution of ambient noise sources, which are dominated by the 7 s period microseismic band (Fig. 3b). Some transient signals remain, especially in the coda part of large waveforms from teleseismic earthquakes. For additional evaluation of the usefulness of this network, we apply the single-station cross-correlation functions used in ambient-seismic noise monitoring in a later section.
Most of the previous denoising networks would construct the noise time-series from direct subtraction of the 'denoised' earthquake waveforms from the raw data (Zhu et al. 2019a;Novoselov et al. 2022). Our WaveDecompNet has two branches that learn features of the earthquake and noise waveforms, with only weak dependence since their only connection is through the skip-connection layers to the same encoder block (Fig. 1). Because the noise window is not the linear difference between the original and the transient/earthquake signal, we also investigate the waveform residuals and show them in Fig. 6(d). The residuals are obtained from subtracting the earthquake (Fig. 6b) and noise (Fig. 6c) waveforms from the input waveform 6a). In general, the amplitudes of the residual waveforms are quite small (about 20-30 times lower than the ambient noise in the standard deviation of waveform amplitude) except in the window between 10 August 2021 and 15 August 1. These are due to the teleseismic earthquakes (30-100 • great arc distance). No teleseismic waveform was used during the training, in part because the input data length is restricted to 1 min. Therefore, our model does not handle longer seismic periods at this stage, and the coda reconstruction of these long waveforms is imperfect. To illustrate this, we mark the P arrival, calculated from TauP with IASP91 earth model, of these large M5.5+ teleseismic earthquakes to illustrate that in Fig. 6. Including more teleseismic earthquake waveforms in the model training can potentially help to mitigate those large residuals, and we leave it as a future direction to explore.
In the following sections, we use two examples to show the improvement in the data quality after the separation of earthquake and noise waveforms. For the earthquake waveforms, we apply the STA to LTA trigger method (STA/LTA), which is routinely applied in many automatic earthquake detection (e.g. Grigoli et al. 2013Grigoli et al. , 2014;;Hansen & Schmandt 2015) and earthquake early warning (e.g. Chung et al. 2020). For the ambient noise waveforms, we apply the single-station cross-correlation, which is commonly used to monitor the near-surface processes (e.g. Wegler & Sens-Schönfelder 2007;De Plaen et al. 2016;Viens et al. 2018;Illien et al. 2021, etc.). We need to emphasize that in this study, we are not proposing any new method/algorithms for earthquake detection and ambient noise cross-correlation. Instead, we attempt to provide a tool to improve the data quality of both earthquake and noise waveforms. The decomposition of seismic signals can be easily combined with any state-of-the-art algorithms or techniques.
We apply a recursive STA/LTA to the continuous data (Vanderkulk et al. 1965;Allen 1978;Withers et al. 1998). This particular STA/LTA algorithm produces a decaying exponential impulse response, and that is a sharper impulse than the original STA/LTA algorithm (Withers et al. 1998). The settings of the STA/LTA parameters are chosen partly from Trnkoczy (2009) and also from trial-and-error tests. The short time window length is 2.0 s, the long time window length is 60.0 s, the on-threshold is 6.0, and the offthreshold is 2.0. Because we simultaneously run STA/LTA detection on 3-component waveforms, we perform a coincidence trigger with a threshold of 2, which means that a detection trigger occurs when the STA/LTA ratios of any two of three components exceed the on-threshold.
The STA/LTA time-series are a lot cleaner in the separated earthquake waveforms than in the original seismograms (Fig. 7), which manifests itself in two aspects. First, the increased signal-to-noise ratio of the separated earthquake waveforms improves the accuracy of the detection time automated by STA/LTA triggers. For the example shown in Fig. 7, the arrival time cannot be correctly picked in the raw data using automated STA/LTA thresholding detector (Fig. 7a). Nevertheless, with the noise separated by WaveDecomp-Net, the event arrival can be easily detected, and the accuracy of picking the first P arrivals can be improved compared with that estimated from the IASP91 velocity model (Fig. 7b). Secondly, we can detect many smaller signals (either smaller magnitude or more distant events) buried in the noise, which is suggested by the increased number of coincidence triggers in the separated waveform, from 38 in the original time-series to 667 in the separated waveforms. The wavefield separation increases the detection by a factor of 17.6 in the current STA/LTA settings. This ratio varies with the choice of the threshold. It is 2.6 (478 events in the raw waveform versus 1223 events in the separated waveforms) for a coincidence threshold of 1. It is 5.0 (19 events in the raw waveforms versus 95 events in the separated waveforms) for a coincidence threshold of 3. Furthermore, the STA/LTA triggers from the separated earthquake waveforms exhibit a greater coincidence with the catalogued local M2+ earthquakes (Fig. S4a).
The routine application of STA/LTA often requires a specific filtering of the data before calculating the STA/LTA ratio. The choice of frequency band can significantly alter the STA/LTA ratio. For example, we test the STA/LTA method on the continuous waveforms filtered at 0.5 Hz and above. The STA/LTA ratio is in general lower than that of the WaveDecompNet-separated waveforms (Fig. S5) so as the number of detected events (Fig. S4b). Tuning the parameters of STA/LTA detector is not the scope of this study, but our results suggest that using WaveDecompNet before STA/LTA would increase the quality of the picks.
Single-station cross-correlations of ambient noise signals are related to the zero-offset Green's function (Claerbout 1968;Draganov et al. 2009;Saygin et al. 2017;Clayton 2020). Monitoring phase changes in the single-station measurements have enabled the monitoring of changes in the near-surface environment that occur during earthquakes (Wegler & Sens-Schönfelder 2007;Viens et al. 2018), volcanic unrest (De Plaen et al. 2016), and to monitor shallow hydrology (Illien et al. 2021). Here, we do not attempt to verify that the single-station correlation is proportional to the Green's function. Instead, we evaluate the temporal stability of the singlestation cross-correlations of ambient noise signals decomposed by WaveDecompNet.
We calculate all nine components of the cross-correlation tensor. We select 1-min-long windows of waveforms and pad them with zeroes from 600 samples to 2048 samples (204.8 s). We then follow the spectral method from Viens et al. (2020) to calculate the ambient noise single-station cross-correlation function (CCF):
where i, j corresponds to components (E, N and Z), âi is the Fourier transform of the i-component waveform, * represents the complex conjugate. The operator F -1 ( • ) is the inverse Fourier transform, and | • | is a spectral whitening operator. The amplitude spectrum whitening process uses a running mean of eight samples in the frequency domain as in conventional processing (Bensen et al. 2007).
We substack the CCFs every 6 hr to evaluate their coherence over time. We show the causal part (positive lag) of CCFs in different frequency bands: high frequency (HF) 2.0-4.0 Hz (Fig. 8); medium frequency (MF) 1.0-2.0 Hz (Fig. S6) and low frequency (LF) 0.1-1.0 Hz (Fig. S7). Each figure shows the single-station cross-correlations from the original raw waveforms and ones obtained from the decomposed ambient noise waveforms. We also stack all of the CCFs to form a reference month-stack, from 31 July 2021 to 1 September 21. We calculate the correlation coefficient between each 6-hr substack and the month-stack reference waveform and show them in Fig. 9. We use the stability of cross-correlation as a success metric of ambient seismic noise recovery.
Because of transient earthquake signals in the raw waveforms, we can see large fluctuations (mostly reduced amplitudes) in the crosscorrelation functions at all frequency bands, especially for the E-Z, N-Z, Z-E and Z-N functions (Figs 8,S6 and S7). These fluctuations in the correlation functions and drops in their coherence arise from transient signals in the original time-series, mostly the local earthquake signals (green dots in Figs 8, S6 and S7). However, many of these fluctuations disappear when using the decomposed noise signals to calculate the CCFs. Furthermore, some of the coda phases that are weak in the original CCFs appear clearly in the correlation functions built from the separated noise. These coda phases potentially correspond to discontinuities in seismic wave speed beneath the seismic station. With the transient earthquake signals removed, WaveDecompNet can potentially help better constrain the velocity structure underneath the seismic station and at least improve the time-dependent monitoring of Earth structure. Additional work is necessary to verify the nature of these coda phases and whether they can be related to Earth structure.
In general, the improvements are greatest in the HF band (Fig. 9c) and weakest in the LF band (Fig. S7). This could be due to the fact that removing the small earthquakes that are rich in high frequency content has the most significant effect. While the improvements are visible at most components and frequencies, some components are deteriorated, for example, the N-Z component and Z-E component in MF (Figs 9b and S6). Further analysis is needed to understand and improve on these specific cases.
We develop a machine-learning-based model, WaveDecompNet, to decompose earthquake and ambient noise signals from raw seismic data. We combine the STEAD and local ambient noise to form the data set to train and test WaveDecompNet. Our network consists of one encoder branch, two decoder branches, and two bottlenecks. We systematically explore the performance of models using different types of bottlenecks, and we find that the network using the LSTM bottleneck has the best performance.
One possible explanation for the different performances between bottleneck architectures is the difference in model complexity. On the one hand, the bottlenecks might be too simple. None and Linear models have fewer parameters and higher training loss (about 0.07) than other models, implying that they may not be enough to understand the internal features of the seismograms properly. On the other hand, the Transformer model is more complex than the other models with the most parameters, but its mean test loss (0.067) is higher than LSTM (0.044) and attention (0.045) models. The Transformer model may be too complex for the current problem and may be improved as more diverse data are included in the training. The LSTM and Attention models share a similar overall complexity and achieve the lowest test loss and the most stable training. The similarity and systematically low values of the training and validation losses for these models also indicate that those two models have already correctly 'learned' the features in the training data sets (Figs 2c and d). It is interesting to note that the LSTM model performs slightly better than the Attention and Transformer models, especially on the ambient noise part with a higher EV score (Fig. 5b). This may imply that the sequential information is essential for reconstructing the waveforms and that LSTM performs better for short sequential series in bottlenecks. We speculate that the feature extraction of the encoder block suffices in representing the temporal sequencing in the bottleneck layer, thus the limitations of LSTMs that long memory is not long enough become no longer important.
We test the performance of WaveDecompNet on observed continuous data. We apply the trained model to a 1-month-long continuous seismic data at IU.POHA, and successfully separate the corresponding local earthquake and noise signals, except the long-duration teleseismic signals that are not included in the training data set. Additionally, we apply an automated transient detector (STA/LTA) and an established ambient-noise seismology monitoring method to the separated earthquake and noise signals, respectively. Our results show that the quality of both decomposed earthquake and noise signals has been improved significantly: we can obtain more STA/LTA triggers and a highly coherent ambient-noise correlation function.
There are some limitations to the model. First, it is only trained on local earthquake waveforms (<350 km). The lack of teleseismic waveforms in the training data set, especially those from large earthquakes, may decrease the performance of WaveDecompNet when handling the time windows with teleseismic earthquake waveforms. While we extract the general patterns of the teleseismic earthquake waveforms correctly, there remain large residuals in the ambient noise waveforms and residual waveforms. Secondly, we only include the local noise from a single island station, IU.POHA. We also test with other stations from the Hawaiian Volcano Observatory and find the network trained from IU.POHA can still successfully separate the earthquake, but the coherence of ambient noise worsens, and there are higher waveform residuals. The expansion of WaveDecompNet to other seismic stations may highlight some of the limitations of supervised deep learning: these methods tend to interpolate within the data space used for training. Therefore, continued training with other stations and other data may be required for optimal performance. A good direction to improve the network performance is to include additional and different types of seismic data. For example, we can include more regional data sets of earthquake waveforms such as these of Italy (Michelini et al. 2021), Chile (Münchmeyer et al. 2020) or southern California (Southern California Earthquake Data Center, SCEDC 2013), and some available teleseismic data for better separation of earthquake waveforms. We should also include the ambient noise waveforms from other stations and regions for a specific regional or global ambient noise study.
Future developments may involve the integration of multiple stations. The combination of multiple stations to combine the automated triggered events help reduce the false (non-tectonic) detections. It also helps locate the event and build a complete earthquake catalogue. Furthermore, modifying the network to add more stations may help improve the stability of the inter-station cross-correlations, which in turn can be used for better Earth imaging. An additional effort could compare more systematically against other methods (Zhu et al. 2019a;Novoselov et al. 2022) in a similar way to the phase picking comparison exercise of Münchmeyer et al. (2022). This task will require modifying the input layers for different window lengths and performing a systematic continued or retraining on the same data set. Finally, further investigations of the residual waveforms would be interesting because they may include different types of seismic signals beyond earthquakes or ambient noise, such as volcanic or non-volcanic tremors or instrumental glitches.
The STEAD data set is available at https://github.com/smousavi0 5/STEAD. The continuous seismic data from IU.POHA (IU: doi: 10.7914/SN/IU) are downloaded using Obspy (available at https://gi thub.com/obspy/obspy/wiki). PyTorch machine learning framework (https://pytorch.org) is to build and train the network. The early stopping module is from https://github.com/Bjarten/early-stoppin g-pytorch. The local earthquake mentioned in Fig. 7 is from http: //service.iris.edu/f dsnws/event/1/query?eventid=11465634. All the codes to reproduce this work are hosted on Github at https://github.c om/yinjiuxun/WaveDecompNet-paper, WaveDecompNet is hosted on https://github.com/yinjiuxun/WaveDecompNet.
Downloaded from https://academic.oup.com/gji/article/231/3/1806/6651381 by University of Washington user on 03 November 2022