Multi-task multi-station earthquake monitoring: An all-in-one seismic Phase picking, Location, and Association Network (PLAN)

The proposed multi-station multi-task PLAN (Fig. 1) employs a Graph neural network (GNN) (Kipf and Welling, 2016) as the backbone to integrate the four functional modules of waveform feature extraction, earthquake location, multi-station association, and a physics-informed multi-station phase picking (further details are provided in the Methods section). Compared with CNN, GNN is naturally suited for handling seismic data acquired from irregularly spaced stations whose quantity and location may vary (McBrearty and Beroza, 2023). In constructing our GNN, we define the graph nodes with seismic stations and define the feature vector for each node with the location and recoded seismic signal of the corresponding station. All nodes (corresponding to the seismic stations) are linked together and the linking weights are learnt during training to infer the relationships among the stations. We construct the GNN layers with TransformGConvs (Shi et al., 2020), which are designed based on an attention mechanism (Vaswani et al., 2017) to learn the dynamically linking weights among different stations. (details about TransformGConvs are provided in the Methods section). In addition, the graph nodes are not fixed so that the GNN could be adapted to variations in the station number and location.

Raw three-component seismic signals are input into the graph nodes. The front-end waveform feature extraction module, constructed as an encoder-decoder CNN and shared among the nodes, extracts their corresponding key features. The station feature extraction block, constructed as two MLPs and shared among the nodes, extract geographic features from the normalized input longitudes, latitudes, and elevations of the stations. The earthquake location module then concatenates the extracted waveform and geographic features and employs multiple TransformGConvs to aggregate these features from multiple nodes to predict the event depth and station-event offset. The predicted offsets and depth are further used to determine the event location by a triangulation algorithm (Yu and Segall, 1996). In the earthquake location module, we do not directly predict the event location but the station-event offsets and depth instead, because we will use them as input into a followed multi-station association module to estimate the time shifts needed to align the P and S arrivals.

This multi-station association module plays a key role in bridging the tasks of earthquake location and multi-station phase picking and introduces physical constraints between the two tasks. Prior to aggregating the waveform features from different stations for multi-station phase picking, the features corresponding to the same earthquake are required to be initially aligned or associated; otherwise, aggregation of unaligned features could mutually interfere and ultimately degrade the picking results. The multi-station phase-picking module includes a non-trainable physical layer, implemented with the Pytorch (Paszke et al., 2019) roll function, to shift and align the waveform features (from the decoder of the waveform feature extraction module) using the time shifts. Subsequently, multiple TransformGConvs in the phase-picking module aggregate the aligned waveform features to enhance the phase-picking features in the aligned space. Eventually, another physical layer unshifts the aggregated features back to the original space, followed by two convolutional layers to obtain the P/S-wave picks at all the stations.

Three regression loss functions are defined for the three modules of phase picking, association, and earthquake location and then combined to jointly train the entire network. Because all the modules are interconnected within the entire network, the training process finds an optimal network that could perform all the tasks both accurately and consistently. Moreover, after training, the multi-station association module could be detached from the network and utilized to calculate the S-P differential travel time with inputs of offsets and event depth. Further details on this module are provided in the section titled “Multi-station association module”.

We tested the proposed PLAN in two regions of Ridgecrest and Japan. For the Ridgecrest region (Fig. 2A), seismic recordings from 16 California Integrated Seismic Network stations within an epicentral distance of < 80 km were collected from July 4, 2019, to October 4, 2019, for a total of more than 71,000 M > 0 earthquakes. The data for Japan (Fig. 3A) included M > 2 earthquakes that occurred between January 1, 2011, and December 31, 2011, including the Mw 9.1 Tohoku sequence. We collected the 3-component High Sensitivity Seismograph Network (NIED Hi-net) (Obara et al., 2005; Aoi et al., 2020) data from over 35,000 events. Subsequently, the data were randomly divided into training, validation, and test sets (85%, 5%, and 10%, respectively) in both regions.

The number of stations corresponding to each event in the training samples varied, and the trained network can flexibly handle situations where the number of stations changes in actual data. Further, the distributions of the number of stations per event in the training and test sets were balanced. The results for the test sets in the two regions are presented in Figs. 5 and 6, respectively. To accommodate different range scales in the two study regions, we used different window lengths in two regions (30.72 s for Ridgecrest and 61.44 s for Japan) with the same sampling frequency (100 Hz).

To ensure a fair comparison with the existing phase-picking methods, we followed the same data preprocessing procedures used in previous studies (Zhu and Beroza, 2018; Mousavi et al., 2020). This involved normalizing the data by removing the mean and dividing by the standard deviation, and using a Gaussian-shaped target function as training labels for the P/S-phase arrival times. The use of Gaussian-shaped targets is effective in phase picking  (Zhu and Beroza, 2018; Mousavi et al., 2020), and in our study, the probabilities of P wave and S wave arrival times were set to 1 at the first arrival time and decreased to 0 within a 20 sample window before and after each phase arrival.

We compared the performance of PLAN with that of other established deep learning methods for earthquake phase picking (PhaseNet (Zhu and Beroza, 2018) and EQTransformer (Mousavi et al., 2020)) and location (Aggreated-GNN (van den Ende and Ampuero, 2020)). All of the methods were retrained on the same training set and evaluated on a common test set. As shown in the Ridgecrest application in Fig. 2B-2C, the performance of PLAN in phase picking was superior to that of the other two deep learning-based methods. Specifically, the residual distribution of the P-wave picks for PLAN was more concentrated than that of the other methods, indicating a higher overall accuracy. For S-wave picks, PLAN performed significantly better than EQTransformer because the distribution of PLAN was narrower whereas the difference in performance between PLAN and PhaseNet was relatively minor.

In terms of localization, our method (PLAN) outperformed the other approaches (Fig. 2D-2E; Table 1). The distribution of PLAN was notably more concentrated than that of Aggregated-GNN, particularly in terms of offset prediction. To further demonstrate the effectiveness of TransformGConv, we replaced the TransformGConv layers (Fig. S1) with GCN (Kipf and Welling, 2016), SAGE (Hamilton et al., 2017), and GATv2 (Brody et al., 2021) respectively. Among the various methods compared, our method yielded the best results in terms of offset, with an average error of 1.09 km and a standard deviation of 1.41 km. Furthermore, PLAN also outperformed Aggregated-GNN in terms of depth localization, regardless of whether it was based on GCN, GATv2, or TransformGConv. These results demonstrated the superiority of the proposed PLAN in location estimations.

Furthermore, we used three metrics of mPrecision, mRecall, and mF1 (described in the Methods section), to quantitatively evaluate the performance of the five methods (Table 2). In five of the six metric scores for the P-wave and S-wave picking results, our attention mechanism-based GNN method outperformed the other methods. The only exception was the mPrecision metric of P-wave picking, where the EQTransformer showed slightly higher scores than PLAN. Notably, even the simplified version of the multi-station phase-picking method, such as the SAGE-based PLAN, outperformed both the single station-based picking methods of EQTransformer and PhaseNet in mF1 scores for S-wave picking. This indicated that the phase-picking accuracy is significantly improved by multi-station picking, which effectively utilizes inter-station contextual information. We not only adjusted the time threshold while maintaining a constant picking probability for evaluation but also fixed the time threshold and changed the probability of picking to calculate and plot the precision-recall curves for four models (Fig. S2). Consistent with the aforementioned results, the PLAN model, based on TransformGConv, demonstrated superior performance in terms of F1, encompassing both P-wave and S-wave.

We retrained all the methods on the Japan training set for the evaluation. Compared to its performance in the Ridgecrest region, PLAN exhibited an even better performance in Japan (Fig. 3B-3C). Further, PLAN demonstrated a remarkably better performance than PhaseNet and EQTransformer for both P- and S-wave picks. The offset predicted by PLAN was notably more accurate than that predicted by Aggregated-GNN, with higher and narrower residuals (Fig. 3D-3E). In terms of depth estimation, although PLAN maintained a narrower residual distribution, the highest point of the distribution was shifted systematically, compared with the Aggregated-GNN method. Table 1 presents the comprehensive quantitative comparison of the results. Although the TransformGConv-based PLAN method did not demonstrate particular superiority in depth estimation, it excelled in offset estimation. Further, the GATv2-based PLAN showed the lowest depth error, indicating potential improvement of localization capabilities of the proposed PLAN.

Similar to the Ridgecrest example, we assessed the phase-picking performance of various models applied to the test data from Japan using mPrecision, mRecall, and mF1 metrics (Table 3) and precision-recall curves (Fig. S3). The TransformGConv-based PLAN model achieved superior results in terms of mRecall (95.14 for P-waves and 85.09 for S-waves) and mF1 (95.46 for P-waves and 86.72 for S-waves), whereas EQTransformer performed best in terms of mPrecision of P-waves and S-waves. TransformGConv-based PLAN demonstrated high mRecall scores, indicating that a large proportion of the samples containing P/S-waves were correctly detected. However, this was achieved at the expense of a slightly lower mPrecision compared to that of the EQTransformer, with some non-P/S-waves incorrectly classified as P/S-waves. The mF1 score provided a more comprehensive evaluation of the model performance, considering both the reduction in missed detections and the increase in correct detections. In this context, TransformGConv-based PLAN had the highest F1 score, indicating that it effectively reduced the missed detections of P/S-waves and increased the proportion of correct detections.

To simultaneously pick P/S-waves from multiple stations, PLAN utilizes GNNs, which typically aggregates raw signals received at different stations. Feature aggregation across multiple stations introduces inter-station constraints and enhances the features at each station. However, because of different travel times of the same source across different stations, directly aggregating the signals from multiple stations would deteriorate multi-station picking.

To address this issue, the proposed method employes a multi-station association module to estimate the time shifts as illustrated in Fig. 1. The input to this module is the offset of each station with respect to the event and its depth. The module output is the corrected time shifts of the P/S-wave for each station. To align the P/S-wave features, the correction criteria for the P/S-wave features were set at the 10 and 15 s in the Ridgecrest region and at the 20 and 32 s in Japan. Using these criteria, the multi-station association module was trained to estimate the corrections of P/S-wave for each station. These corrections are then used to align the waveform features, enabling the graph convolution to aggregate the features in a temporally aligned space. Consequently, the method could enhance or compensate for the features at each station by fusing the aligned features from other stations, allowing simultaneous and accurate multi-station picking.

The multi-station association module can be utilized independently after training. It converts the distance and depth information into arrival information and calculates the S-P differential travel time (Crotwell et al., 1999). Figure 4 illustrates the arrival time differences of the P/S-waves at various stations in Japan. Although the training process utilizes a maximum of 37 stations for a single event, the module can be adapted to cases with any number of stations (e.g., hundreds of stations shown in Fig. 4A - Fig. 4D) to estimate the P/S-wave arrival time differences for all stations. These results indicated that the module effectively enforced physical constraints based on time shifts within the overall network.

Several studies have shown that GNNs have the potential to deal with irregularly spaced stations for phase association and event localization (McBrearty et al., 2019b; van den Ende and Ampuero, 2020; Yano et al., 2021; McBrearty and Beroza, 2022b, 2023; Zhang et al., 2022; Bilal et al., 2022). Here, we build a graph-based network (Fig. 1) for mult-station earthquake monitoring. To utilize the GNN, we first need to change the data from the matrix format to the graph format and employ a graph-based representation of the stations, where each station is represented as a node in the graph and the three-channel data and the station location are used as the features of each node. In contrast to the current single-station processing methods (Ross et al., 2018b, a; Zhu and Beroza, 2018; Mousavi et al., 2019; Zhu et al., 2019; Pardo et al., 2019; Wang et al., 2019; Liu et al., 2020; Mousavi et al., 2020; Zhu et al., 2022a), which treat each three-channel data as an individual input sample, our approach inputs all the three-channel data received from multiple stations per event as a single sample. This allows for efficient aggregation of information from multiple stations during network training. As a result, the features of different stations could be effectively integrated using GNNs during the aggregation process.

In this study, we have evaluated various graph aggregation methods, including GCN (Kipf and Welling, 2016), GraphSAGE (Hamilton et al., 2017), GAT (Veličković et al., 2017), GATv2 (Brody et al., 2021), and TransformGCONV (Shi et al., 2020). Through this evaluation, we have determined that TransformGCONV, which is based on attention mechanism (Vaswani et al., 2017), is the most suitable module for the proposed PLAN. The message aggregation of TransformGCONV could be represented as:

$$\mathbf{x}_{i}^{\prime}=\mathbf{W}_{1}\mathbf{x}_{i}+\sum_{j\in\mathcal{N}(i)}\alpha_{i,j}\mathbf{W}_{2}\mathbf{x}_{j},$$

(1)

where $\mathbf{x}_{i}^{\prime}$ represents the aggregated features at the source node, and $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ represent the features of the source and distant nodes before aggregation, respectively. $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ are the trainable matrices. In addition, the attention coefficients $\alpha_{i,j}$ are computed via dot-product attention as follows:

$$\alpha_{i,j}={softmax}\left(\frac{\left(\mathbf{W}_{3}\mathbf{x}_{i}\right)^{\top}\left(\mathbf{W}_{4}\mathbf{x}_{j}\right)}{\sqrt{d}}\right),$$

(2)

where $\mathbf{W}_{3}$ and $\mathbf{W}_{4}$ are the trainable matrices. Similar to the attention mechanism(Vaswani et al., 2017), the source feature $\mathbf{x}_{i}$ and distant feature $\mathbf{x}_{j}$ are transformed into query vector and key vector, respectively, using $\mathbf{W}_{3}$ and $\mathbf{W}_{4}$. Compared to other graph aggregation methods, the use of the attention-based mechanism (equation 1) in TransformGCONV allows for a more fine-grained representation of the relationship between different stations, thereby improving the accuracy and efficiency of the proposed method.

Here, we design a multi-station multi-task network for simultaneous phase picking, association, and location. The network (Fig. S1) comprises four components: a waveform feature extraction module, an earthquake location module, a multi-station association module, and a physics-informed multi-station phase-picking module. Similar to previous deep-learning-based phase-picking approaches (Zhu and Beroza, 2018; Mousavi et al., 2020; Zhu et al., 2022a, c), we design an encoder to extract waveform features and a decoder to produce phase-picking results. However, to address the multi-station phase-picking problem, we introduce the GNN-based TransformGCONV for aggregating features from multiple stations.

Because aligned waveform features are easily used and aggregated in GNNs for multi-station phase picking, we do not employ it in the waveform feature extraction module (Fig. S1A), where the features are relatively shifted in time. Although we use a U-shape neural network for feature extraction to solve the phase-picking problem, it could be replaced with other single-station-based phase-picking networks, such as EQTransformer. No matter what type of network architecture is used, the features extracted from the middle of the network are input into the earthquake location module, and the structure of the final few layers of the network are modified for the purpose of multi-station phase picking. Additionally, the kernel size of all convolutional layers in the waveform feature extraction network is set to 7.

For the earthquake location module (Fig. S1B), we first extract features from the normalized coordinate information of the stations within the range of [0,1] through two fully connected layers (3-48-96). Simultaneously, the waveform features extracted from each station are further processed through several convolutional layers and then flattened. Subsequently, the position and waveform features are concatenated and passed through two fully connected layers (192-192-96). This fuses the position and waveform features at each station. The fused features are further aggregated among multiple stations by several GNN layers to predict the offsets of each station with respect to the event and its depth. Because there is only one depth parameter for each sample, we add a global average pooling before the output. In summary, this module allows the integration of both location and waveform information into the feature extraction process, which is crucial for accurate event localization.

Finally, in the physics-informed multi-station phase-picking module (Fig. S1D), we incorporate physics-motivated constraints of time alignment among waveforms corresponding to the same earthquake event. We first utilize a mulit-station association module (Fig. S1C) to calculate the relative alignment shifts between stations using the estimated offsets and the depth of the event. We then use the shifts to align the waveforms to a common time standard and aggregate the features across multiple stations in the phase-picking module. Subsequently, the aggregated features at each station are unaligned and fed to two layers of convolution to yield final P/S-wave picking results. This process leverages the physical information of the event location to improve the robustness and accuracy of the multi-station phase picking.

Our multi-task learning network model has three output results corresponding to phase picking, phase association, and earthquake localization. To train the model, we define three different loss functions for these three different tasks. For phase picking, instead of using the commonly used cross-entropy, we choose the mean square error (MSE) as the loss function, which is suitable for training in multi-task problems. To estimate offsets and depth, which is similar to event localization, we also use MSE as the loss function, as suggested by previous studies (van den Ende and Ampuero, 2020; Zhang et al., 2022). Finally, to calculate the P/S-wave shift, we define the loss function as follows:

$$\mathcal{L}_{\Delta p}=\sum_{i=1}^{n}\mid{CTime}_{p}-\left({label}_{p_{i}}+\Delta t_{p_{i}}\right)\mid,$$

(3)

$$\mathcal{L}_{\Delta s}=\sum_{i=1}^{n}\mid{CTime}_{s}-\left({label}_{s_{i}}+\Delta t_{s_{i}}\right)\mid,$$

(4)

where ${CTime}_{p}$ and ${CTime}_{s}$ represents the reference times where the P- and S-wave picks are aligned, respectively. Moreover, ${label}_{p_{i}}$ or ${label}_{s_{i}}$ represents the manually picked P/S-wave arrival time for each station, and $\Delta t$ represents the predicted P/S-wave shift. Finally, we combine the three types of loss functions to form the overall objective function:

$$L_{{total}}=\lambda_{1}\mathcal{L}_{{picking-p}}+\lambda_{1}\mathcal{L}_{{picking-s}}+\lambda_{2}\mathcal{L}_{\Delta p}+\lambda_{2}\mathcal{L}_{\Delta s}+\lambda_{3}\mathcal{L}_{{offset}}+\lambda_{3}\mathcal{L}_{{depth}}.$$

(5)

Here, we set the coefficients $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ to 1.

During the training process, the model was optimized using the ADAM (Kingma and Ba, 2014) method with an initial learning rate of 0.001, which is gradually decreased with a decay rate of 0.9 every 100 epochs. To enhance the training efficiency, we randomly selected 2048 events from the training set for each epoch, rather than using the entire data. The model was trained for a total of 2000 epochs with a batch size of 16, and the training process required approximately 24 h using 1 NVIDIA Tesla A100 GPU.

In previous studies (Mousavi et al., 2020; Zhu et al., 2022c), true positive phase picks were defined as those within 0.5 s of the predicted pick. The rest were counted as false positives. Nevertheless, owing to potential errors in the labels of the dataset, such statistical results based on a single threshold may not be reliable. Thus, to better evaluate the performance of algorithms, we introduce new metrics, mPrecision, mRecall, and mF1, which are calculated using multiple thresholds, following previous research(Zheng et al., 2022). The metrics are defined as:

$$\mathrm{mPrecision}=(\mathrm{Precision}@11+\mathrm{Precision}@12+\cdots+\mathrm{Precision}1@50)/40,$$

(6)

$$\mathrm{mRecall}=(\mathrm{Recall}@11+\mathrm{Recall}@12+\cdots+\mathrm{Recall}1@50)/40,$$

(7)

$$\mathrm{mF}1=(\mathrm{F}1@11+\mathrm{F}1@12+\cdots+\mathrm{F}1@50)/40.$$

(8)

where x@11, x@12, $\cdots$ , x@50 are Precision, Recall, or F1 metrics when the thresholds are 11, 12, $\cdots$, 50 samples (corresponding to 0.11 s, 0.12 s, $\cdots$, 0.5 s of time), respectively. These metrics, mPrecision, mRecall, and mF1 reward detectors with better picking results and, therefore, can more reasonably or fairly assess the performance of the different methods.