Seismic waveform data multi-task conjoint analysis method

Abstract

The invention discloses a seismic waveform data multi-task conjoint analysis method, and belongs to the technical field of seismic data analysis. The method comprises the following steps: acquiring a three-component seismic waveform and inputting into an encoder to extract a waveform feature vector; performing position coding on the waveform feature vector; obtaining a task vector according to the task identifier, and carrying out splicing and attention fusion on the task vector and the position-coded waveform feature vector; and inputting the fusion features into a parallel multi-task decoder, and obtaining the seismic type, the seismic phase and the initial movement direction at the same time. According to the method, the convolutional neural network, the recurrent neural network and the transformer are combined, so that the time-frequency characteristics and the sequence characteristics of the seismic data can be better utilized, and Pg, Sg, Pn and Sn seismic phase detection, P-wave initial motion direction judgment and event type judgment work can be completed at the same time on the basis of ensuring the detection precision.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data analysis, and more particularly to a multi-task joint analysis method for seismic waveform data. Background Art

[0002] In the field of earth science, seismic waveform data is the core basis for studying the internal structure of the earth, earthquake location, source mechanism inversion and earthquake early warning.

[0003] Currently, existing deep learning models for seismology have achieved quite good detection results. However, they focus on a single task and lack the accuracy to identify complex seismic phases (such as Pn, Pg, Sn, and Sg). For example, models such as PhaseNet (a seismic phase picking algorithm based on the U-Net architecture proposed by Zhu et al., 2018) and LPPN (a P-wave and S-wave phase picking model proposed by Yu et al., 2022) are only designed for P / S wave detection.

[0004] In addition, Ross et al. (2018) constructed a convolutional neural network model for determining the arrival time and initial motion polarity of the P wave and DiTingMotion (a convolutional neural network model constructed by Zhao et al., 2023) can only determine the initial motion polarity and cannot output it synchronously with other seismic phase detection.

[0005] Therefore, there is still a need for a general model that can contain enough information and fully explore the time-frequency characteristics of single-station waveform data for earthquake analysis. Summary of the Invention

[0006] In view of this, in order to at least partially solve the above technical problems, the present invention provides a multi-task joint analysis method for seismic waveform data.

[0007] In order to achieve the above object, the present invention adopts the following technical solutions:

[0008] A multi-task joint analysis method for seismic waveform data, comprising:

[0009] Obtain three-component seismic waveform and input it into the encoder to extract waveform feature vector;

[0010] Position encoding of waveform feature vectors;

[0011] Obtain a task vector based on the task identifier, concatenate the task vector with the position-encoded waveform feature vector, and perform attention fusion;

[0012] The fused features are input into a parallel multi-task decoder to obtain the earthquake type, earthquake phase and initial motion direction respectively.

[0013] In an optional embodiment, the encoder includes a plurality of convolutional residual modules connected in sequence, and the convolutional residual module includes:

[0014] The first one-dimensional convolution layer and the second one-dimensional convolution layer,

[0015] The output of the second one-dimensional convolutional layer is skip-connected to the input of the first one-dimensional convolutional layer.

[0016] Preferably, the convolution kernel size of the first one-dimensional convolution layer is 5, the step size is 2, and the number of filters is the preset base number of the current level; the convolution kernel size of the second one-dimensional convolution layer is 5, the step size is 1, and the number of filters is the same as that of the first one-dimensional convolution layer.

[0017] In an optional embodiment, a position encoder is used to perform position encoding on the waveform feature vector, and the position encoder is a single-layer bidirectional GRU network model.

[0018] In an optional embodiment, a task encoder is used to obtain a task vector according to a task identifier, and the task encoder includes a first one-hot encoding module, a second one-hot encoding module, and a task vector generation module;

[0019] A first one-hot encoding module is used to map the task identifier into a first binary vector, where the vector dimension is consistent with the total number of task categories;

[0020] The second one-hot encoding module is used to perform high-dimensional expansion on the first binary vector to generate a second binary vector;

[0021] The task vector generation module is used to map the second binary vector into a continuous task vector through a fully connected layer.

[0022] In an optional embodiment, the concatenated features are input into a Transformer model to obtain fused features.

[0023] In an optional embodiment, the multi-task decoder includes a first task decoder and a second task decoder, wherein:

[0024] The first task decoder is used to output the earthquake type,

[0025] The second task decoder is used to output the earthquake phase and initial motion direction.

[0026] In an optional embodiment, the first task decoder is a multilayer perceptron constructed by a multi-layer fully connected network.

[0027] The second task decoder includes a plurality of cascaded decoding units, each of which includes a one-dimensional upsampling module and a one-dimensional convolutional layer connected thereto.

[0028] In an optional embodiment, the second task decoder is further used to extract additional waveform features, reconstruct and output the seismic waveform.

[0029] The present invention discloses a multi-task joint analysis method for seismic waveform data. Compared with the existing technology, it can simultaneously complete Pg, Sg, Pn and Sn seismic phase detection, P wave initial motion direction determination and event type judgment. As a pre-trained model, it can be applied to other tasks by modifying the encoder or feature vector.

[0030] Specifically, the advantages of this application include:

[0031] 1. It integrates three modules: seismic phase detection, onset detection and event classification, which can comprehensively cover and effectively detect various seismic phases and onsets;

[0032] 2. The combination of convolutional neural networks, recurrent neural networks and transformers can better utilize the time-frequency and sequence characteristics of seismic data, thereby ensuring the accuracy of detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0034] Figure 1 This is a flow chart of the multi-task joint analysis method for seismic waveform data of the present invention;

[0035] Figure 2 Schematic diagram of the encoder structure of the present invention;

[0036] Figure 3 Schematic diagram of the structure of the position encoder of the present invention;

[0037] Figure 4 This is a schematic diagram of the structure of the task encoder of the present invention;

[0038] Figure 5 This is a flowchart of the application of the Transformer model of the present invention;

[0039] Figure 6 This is a structural diagram of the first task decoder of the present invention;

[0040] Figure 7 This is a structural diagram of the second task decoder of the present invention;

[0041] Figure 8Graphs showing the initial movement determination accuracy of the model at different resolutions of the present invention, wherein (a) shows the initial movement direction determination accuracy of the model at the highest resolution, and (b) shows the initial movement determination accuracy of the model at all resolutions;

[0042] Figure 9 This is a comparison chart of the seismic phase picking accuracy of the model under different signal-to-noise ratios of the present invention;

[0043] Figure 10 This is the earthquake classification test result diagram of the present invention. The left diagram is a model retrained based on data from Inner Mongolia, and the right diagram is a transfer learning model. DETAILED DESCRIPTION

[0044] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0045] The following description sets forth many specific details to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those described herein, and those skilled in the art may make similar generalizations without violating the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0046] The present invention discloses a multi-task joint analysis method for seismic waveform data. This method not only addresses the existing problem of phase-picking models being unable to distinguish more complex seismic phases and simultaneously perform onset detection and event classification, but also ensures detection accuracy. The present invention can be widely applied to fields such as earthquake early warning, earthquake detection, and earthquake cataloging.

[0047] In this embodiment, the multi-task joint analysis method of earthquake waveform data is as follows: Figure 1 , including the following steps:

[0048] Obtain three-component seismic waveform and input it into the encoder to extract waveform feature vector;

[0049] Position encoding of waveform feature vectors;

[0050] Obtain a task vector based on the task identifier, concatenate the task vector with the position-encoded waveform feature vector, and perform attention fusion;

[0051] The fused features are input into a parallel multi-task decoder to obtain the earthquake type, earthquake phase and initial motion direction respectively.

[0052] In one embodiment, in order to better extract seismic waveform features, three components of the seismic waveform are used as input, wherein each component data is normalized using the maximum value, and no other preprocessing is included.

[0053] In one embodiment, since seismic waveform data contains more high-frequency information, a multi-layer convolutional neural network with translation invariance and stretch invariance is used to construct an encoder to process the original waveform; in this embodiment, the encoder includes a plurality of convolutional residual modules connected in sequence, wherein the convolutional residual module includes: a first one-dimensional convolutional layer and a second one-dimensional convolutional layer, wherein the output of the second one-dimensional convolutional layer is jump-connected to the input of the first one-dimensional convolutional layer.

[0054] The first one-dimensional convolutional layer is used to align the dimension of the input data to be consistent with the output dimension of the second one-dimensional convolutional layer.

[0055] As a preferred embodiment, the encoder structure refers to Figure 2 The convolution kernel size of the first convolution layer is 5, the step size is 2, and the number of filters is the preset cardinality of the current level; the convolution kernel size of the second one-dimensional convolution layer is 5, the step size is 1, and the number of filters is the same as the first one-dimensional convolution layer.

[0056] In one embodiment, the waveform features lack position-related information, so a position encoder constructed by a single-layer bidirectional GRU network model is used (structured as Figure 3 ) is used to encode the feature position information, forming a waveform feature vector with position information.

[0057] In one embodiment, since the obtained features need to output information such as earthquake type, a task encoder is used to obtain a task vector according to the task identifier to provide additional features for subsequent analysis.

[0058] In this embodiment, the task encoder includes a first one-hot encoding module, a second one-hot encoding module, and a task vector generation module; Figure 4 ,

[0059] A first one-hot encoding module is used to map the task identifier into a first binary vector, where the vector dimension is consistent with the total number of task categories;

[0060] The second one-hot encoding module is used to perform high-dimensional expansion on the first binary vector to generate a second binary vector;

[0061] The task vector generation module is used to map the second binary vector into a continuous task vector through a fully connected layer.

[0062] In one embodiment, the waveform feature vector with position information and the task vector are vector-concatenated, and the concatenated features are input into the Transformer model (multi-head attention model) to obtain fusion features, so as to comprehensively consider the features of all waveforms and output the feature vector. The process is referred to as Figure 5 .

[0063] In one embodiment, the decoding part constructs a parallel multi-task decoder, including a first task decoder and a second task decoder, wherein:

[0064] The first task decoder is a multi-layer perceptron constructed from a multi-layer fully connected network, which is used to output the earthquake type.

[0065] The second task decoder includes multiple cascaded decoding units, each comprising a one-dimensional upsampling module and a connected one-dimensional convolutional layer for outputting earthquake phases and first-motion directions. Preferably, the one-dimensional upsampling module uses linear interpolation or zero-padding interpolation to upsample the input features by a factor of 2. The one-dimensional convolutional layer has a convolution kernel size of 5 and a stride of 1. The number of filters decreases stepwise according to a preset ratio. The first-stage decoding unit has 1024 filters, while the final-stage decoding unit has K filters, where K is the number of task-related parameters.

[0066] Different decoders can complete different tasks. In this embodiment, four main decoders are added during pre-training.

[0067] The first task decoder structure is as follows Figure 6 ,The final number of features K is specified by different tasks. The task decoder corresponds to the task encoding vector and outputs the earthquake type. When there are a total of 7 earthquake types, K = 7 is specified.

[0068] The second task decoder structure is as follows Figure 7 , it is built based on upsampling layers (using interpolation for upsampling) and convolutional neural network layers (the two are combined and called transposed convolution). Finally, the number of filters K is specified by different tasks, which corresponds to the waveform feature vector with position, so the output length is the same as the waveform length.

[0069] Specifically, when used to output seismic phase types, K = 5, including Pg, Sg, Pn, Sn, and noise;

[0070] When used to output the initial movement direction, K=2.

[0071] In a preferred embodiment, the second task decoder is further used to extract additional waveform features, reconstruct and output the seismic waveform. In this case, K=3, indicating three channels of the waveform.

[0072] In the above embodiment, the first three decoders use cross entropy as the loss function to encode manually labeled data, and the fourth decoder is an autoencoder, which is used to extract more features of the waveform. The purpose of adding this decoder is to allow the pre-trained feature vector to contain more waveform information in addition to the first three encoders.

[0073] The four decoders are trained simultaneously to ensure that the trained model can contain the information of all encoders.

[0074] Although the second, third, and fourth decoders in this application have the same convolutional layer structure, the output layer designs of the different decoders differ significantly. The purpose of the phase decoder is to identify multiple different seismic phases. The loss function of the phase decoder is optimized for phase classification, and ultimately outputs the probability distribution of Pg, Sg, Pn, Sn, and noise. The initial motion decoder is optimized for initial motion direction prediction and has only two output channels. The output of the fourth decoder has the same structure as the input three-component seismic waveform and therefore has three output channels.

[0075] To verify the multi-task joint analysis method of the present invention, the performance of phase type and first motion detection was evaluated:

[0076] Phase detection evaluation uses standard evaluation metrics, including precision (P), recall (R), F1-score (F1-Score), and the mean and standard deviation of the time residuals between the arrival times picked up by the model and the manually annotated arrival times. Three types of samples are defined: true positive (TP), false positive (FP), and false negative (FN). TP refers to a sample in which the error between the manually annotated arrival time and the model-picked arrival time is less than 1 second. FP refers to a sample in which the error is greater than 0.5 seconds. FN refers to a sample in which the model does not pick up the model but is within 1 second of the manually annotated arrival time.

[0077] The definitions of Precision, Recall, and F1 score are as follows:

[0078]

[0079]

[0080]

[0081] Furthermore, the statistical results of the picking accuracy of different seismic phases are shown in Table 1;

[0082] Table 1

[0083]

[0084] Table 1 shows that the model has high detection accuracy for the Pg and Sg phases. The Pg phase has the highest precision, recall, and F1 score, reaching 0.860, 0.866, and 0.863, respectively. The Sg phase has slightly lower precision, recall, and F1 score, at 0.838, 0.846, and 0.842, respectively. The Pn and Sn phases have lower detection accuracy, but similar to the Pg and Sg phases, the Pn phase has higher accuracy than the Sn phase. The Pn phase has a precision, recall, and F1 score of 0.757, 0.743, and 0.750, respectively. These indicators demonstrate that the model has a certain degree of accuracy and stability in detecting the Pn phase. The Sn phase has the lowest accuracy of all phase detections and exhibits the widest error distribution, with a mean error of -55.77 ms and a standard deviation of 1386.36 ms. This is primarily due to the Sn phase being the smallest in the dataset, resulting in a lack of sufficient data samples during training to capture its waveform characteristics. Furthermore, the Sn phase signal is typically weaker than the Pg, Sg, or Pn phases and is more susceptible to noise, making detection more difficult. This can affect the model's performance in real-world data applications. However, the model achieved a precision of 0.383, a recall of 0.377, and an F1 score of 0.380 for Sn phase detection, demonstrating that the model still possesses basic performance for Sn phase detection.

[0085] Initial motion detection evaluation calculates the proportion of the total number of samples correctly identified as "upward" by the model when the initial motion direction of the test sample is "upward" to the total number of samples with an initial motion direction of "upward". Similarly, the proportion of samples with an initial motion direction of "downward" is calculated.

[0086] The above results are visualized using a confusion matrix. Since the dataset used has the clarity of the initial movement annotated, two independent tests were conducted.

[0087] First, samples marked with "I" clarity in the dataset have the highest clarity, and this data was independently tested (16,590 samples) to evaluate the model's performance under optimal conditions. Second, to better reflect the data diversity in actual application scenarios and comprehensively evaluate the model's applicability under different signal-to-noise ratios, data containing all clarity levels (including unmarked clarity "M" and the more moderate clarity "E") was further tested (48,400 samples). The test results are as follows: Figure 8 As shown in the figure, (a) is the initial motion direction determination accuracy of the model at the highest definition, and (b) is the initial motion determination accuracy of the model at all definitions;

[0088] like Figure 8As can be seen, when testing only on the highest-resolution "I" dataset, the model demonstrated high performance, achieving 88.54% and 88.36% accuracy for identifying initial motion with upward (U) and downward (D) polarity, respectively. This result demonstrates that the model is capable of accurately determining the direction of initial motion when the data is high-resolution. When the evaluation scope was expanded to include samples with "M" and "E" resolutions, the model's performance declined slightly, but its accuracy remained high, reaching 83.38% for upward (U) polarity and 80.28% for downward (D) polarity. This demonstrates that the model can effectively identify initial motion even when subjected to varying resolutions and in the presence of more noise, demonstrating its generalizability.

[0089] Furthermore, in order to comprehensively and accurately evaluate the performance of the model under different signal-to-noise ratios, the signal-to-noise ratio is divided into four ranges: less than 0dB, 0dB-1dB, 1dB-5dB, and greater than 5dB. In each signal-to-noise ratio range, the corresponding precision, recall, and F1 score are calculated. Figure 9 ,

[0090] Depend on Figure 9 As can be seen, the model's performance improves with increasing signal-to-noise ratio. When the signal-to-noise ratio is greater than 5dB, the model's F1 scores for detecting Pg and Sg both exceed 0.9, the F1 score for detecting Pn reaches above 0.85, and the F1 score for detecting Sn also reaches above 0.75. The signal-to-noise ratio primarily affects the model's recall, with a smaller impact on the precision. At low signal-to-noise ratios, the model's precision for both Pg and Sg is above 0.8. Furthermore, the model is more susceptible to the signal-to-noise ratio when detecting Pn and Sn. This is primarily because the number of Pn and Sn samples in the dataset is smaller than that of the Pg and Sg phases, making the model more susceptible to the signal-to-noise ratio when detecting Pn and Sn.

[0091] In order to evaluate the event classification performance of the analysis model of this application, a classification test was conducted based on the data from Inner Mongolia. The total number of data was 419, of which 251 earthquake events were used for training and 168 earthquake events were used for testing. There are three main types of earthquake events in this region: natural earthquakes (EQ, 103 events), blasting (EP, 244 events) and collapse (SS, 72 events). The test is divided into two models. The first one is directly trained without transfer learning; the second one is transferred through a general pre-training model. The two models were uniformly iterated 501 times, and the learning rate was fixed at 1e-5. Strictly speaking, earthquake classification is not a new task, and there are also labels for classification problems in the pre-training process. However, due to the extreme imbalance of samples in direct training (natural earthquakes account for 99.5%), this makes the original model predict 100% of natural earthquakes, so it is necessary to use local data for transfer learning. The test results are shown in Figure 10 The left picture is a model retrained based on data from Inner Mongolia, and the right picture is a transfer learning model.

[0092] The original training model achieved a classification accuracy of 67.5%, and after transfer learning, the classification accuracy reached 92.2%. Before transfer learning, all collapse events were classified as blasting events. This was due to the high number of blasting events and the uneven distribution of samples, resulting in all collapses being identified as blasting events. Meanwhile, many blasting and collapse events were identified as natural earthquakes. After transfer learning, the accuracy was significantly improved, with all natural earthquakes successfully identified and the recognition accuracy for both collapses and blasting also improved. However, some collapse events were still identified as blasting events, indicating that some collapse events may share similarities with blasting events based on waveforms.

[0093] Currently, many deep neural networks for seismic phase detection, such as PhaseNet, EQTransformer, and LPPN, are retrained on manually annotated datasets. However, due to the limitations of the training data (all three used the STEAD dataset), these models suffer from accuracy degradation in practical applications. The maximum length of the STEAD dataset is 60 seconds, which means limited noise data before the Pg phase. However, the CSNCD dataset has at least 100 seconds of waveforms. This application selected 275 continuous waveforms from May 21, 2021 and May 22, 2022, including the Ms6.4 Yangbi mainshock and its aftershocks, to test the model and PhaseNet. The statistics are shown in Table 2.

[0094] Table 2

[0095]

[0096] As can be seen from Table 2, for the Pg and Sg phases, the recall rate of the model of the present application is 5% higher than that of PhaseNet, which means that the model of the present application can detect more phases. More training samples may enable the model of the present application to recognize more complex waveforms, thereby enabling it to recognize more phases. In addition, at a higher recall rate, the model of the present application detects fewer phases in continuous data. This means that the model will greatly reduce false detections. The amount of Pg and Sg detected is more balanced, which is beneficial for phase correlation. Since the length of the model input waveform is 10240, which is 102.39s for 100Hz sampling rate data, the model can use more noise data and phases to distinguish the phase type.

[0097] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.

[0098] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A multi-task joint analysis method for seismic waveform data, characterized in that: Obtain three-component seismic waveform and input it into the encoder to extract waveform feature vector; Position encoding of waveform feature vectors; Obtain a task vector based on the task identifier, concatenate the task vector with the position-encoded waveform feature vector, and perform attention fusion; The fused features are input into a parallel multi-task decoder to obtain the earthquake type, earthquake phase and first motion direction at the same time.

2. The analysis method according to claim 1, characterized in that The encoder includes a plurality of convolutional residual modules connected in sequence, and the convolutional residual module includes: The first one-dimensional convolution layer and the second one-dimensional convolution layer, The output of the second one-dimensional convolutional layer is skip-connected to the input of the first one-dimensional convolutional layer.

3. The analysis method according to claim 1, characterized in that The waveform feature vector is positionally encoded using a position encoder, which is a single-layer bidirectional GRU network model.

4. The analysis method according to claim 1, characterized in that Obtaining a task vector according to the task identifier using a task encoder, the task encoder comprising a first one-hot encoding module, a second one-hot encoding module, and a task vector generation module; A first one-hot encoding module is used to map the task identifier into a first binary vector, where the vector dimension is consistent with the total number of task categories; The second one-hot encoding module is used to perform high-dimensional expansion on the first binary vector to generate a second binary vector; The task vector generation module is used to map the second binary vector into a continuous task vector through a fully connected layer.

5. The analysis method according to claim 1, characterized in that The concatenated features are input into the Transformer model to obtain fused features.

6. The analysis method according to claim 1, characterized in that The multi-task decoder includes a first task decoder and a second task decoder, wherein, The first task decoder is used to output the earthquake type, The second task decoder is used to output the earthquake phase and initial motion direction.

7. The analysis method according to claim 6, characterized in that The first task decoder is a multi-layer perceptron constructed with a multi-layer fully connected network.

8. The analysis method according to claim 6, characterized in that The second task decoder includes a plurality of cascaded decoding units, each of which includes a one-dimensional upsampling module and a one-dimensional convolutional layer connected thereto.

9. The analysis method according to claim 6, characterized in that The second task decoder is also used to extract additional waveform features, reconstruct and output the seismic waveform.

Images