Seismic phase picking method based on time-frequency analysis and deep learning

Abstract

The present application relates to a seismic phase picking method based on time-frequency analysis and deep learning, belonging to the field of seismology, a deep neural network comprising an empirical mode decomposition module, a residual shrinkage module, an LSTM layer, a transformer layer, an up-sampling layer and a convolution layer is established, a linear rectifier function is used as an activation function, an Adam adaptive optimization function is used for back propagation, a binary cross-entropy loss function is used as a loss function, and the average absolute error function, error standard deviation, precision, recall and F1 score are used as evaluation indexes to obtain a deep neural network model. When there are differences in the frequencies of noise, P waves and S waves, empirical mode decomposition can effectively separate the three, thereby helping the network to better identify P waves and S waves. Even for wideband noise, the noise will be decomposed into different intrinsic mode functions, thereby reducing the noise content in the intrinsic mode function where the P wave and S wave are located. This time-frequency analysis method can effectively improve the noise immunity of the network.

Description

Technical Field

[0001] This invention relates to the field of seismology, and in particular to a seismic phase picking method based on time-frequency analysis and deep learning. Background Technology

[0002] Seismic phase picking is one of the most fundamental components of earthquake monitoring and early warning. Accurate picking of the arrival times of P-waves and S-waves not only improves the reliability of earthquake catalogs but also provides strong support for the establishment of velocity models for P-waves and S-waves. Furthermore, the arrival time significantly affects the accuracy of source location; even a 0.01-second error in P-wave picking can result in a source location error of tens of meters.

[0003] Currently, traditional P-wave and S-wave acquisition relies on manual picking by human analysts. However, with the increasing number of seismic stations, manual picking is insufficient to handle the massive amounts of seismic data. Researchers have developed automated algorithms, such as the Short-Time Average / Long-Time Average Ratio (STA / LTA) method and the Autoregressive-Akaike Information Criterion (AR-AIC) method. These automated algorithms identify the arrival times of P-waves and S-waves based on abrupt changes in amplitude, energy, and statistical characteristics in seismic records. However, these algorithms are susceptible to noise interference and require setting optimal initial parameters for different ground motions, making them inconvenient to apply.

[0004] In recent years, with the development of computer technology and artificial intelligence, deep learning has been widely used in phase picking. However, complex and ubiquitous noise still significantly reduces the picking accuracy of the model. Therefore, based on the advantage of time-frequency analysis in effectively separating noise, a deep learning model based on empirical mode decomposition is proposed to pick the arrival times of P-waves and S-waves. Summary of the Invention

[0005] The purpose of this application is to address the technical problems of low accuracy and insufficient noise resistance in phase picking methods by providing a deep learning method based on empirical mode decomposition.

[0006] A seismic phase picking method based on time-frequency analysis and deep learning includes the following steps:

[0007] Step 1: Collect seismic motion records and build a dataset:

[0008] Select known ground motion data, and use the P-wave and S-wave picked up from the ground motion data as the corresponding label values ​​of the earthquake record to obtain the phase picking dataset.

[0009] Step 2: Establish a deep neural network model incorporating the Empirical Mode Decomposition (EMD) module:

[0010] The deep neural network consists of an encoder and a decoder. The encoder's first part is an empirical mode decomposition module, which decomposes earthquake vibrations into multiple intrinsic mode functions (IMFs), each representing a different frequency component of the earthquake vibration. The second part is composed of multiple residual shrinking modules connected in series, which add a soft thresholding function to the traditional residual module to filter out background noise. The third part consists of two long short-term memory layers, used to extract temporal contextual information. The fourth part is a transformer module, which assigns weights to different positions in the sequence based on an attention mechanism to amplify key features.

[0011] The decoder consists of multiple upsampling layers and convolutional layers, which expand the feature map to the same length as the input.

[0012] Step 3: Training the deep neural network model:

[0013] The deep neural network is trained by using a loss function and a mean absolute error function to ensure training accuracy, thereby gradually reducing the loss function and obtaining the trained deep neural network model.

[0014] Step 4: Seismic Phase Picking

[0015] The three-component ground motion records are input into the deep neural network model trained in step three, and the output probability sequences corresponding to P-waves and S-waves are obtained. The location of the maximum value of the probability sequence is found, thereby realizing phase picking.

[0016] In step two, the empirical mode decomposition of the encoder includes the following steps:

[0017] 1) Identify all local maxima and local minima in the signal;

[0018] 2) Local maxima of interpolation form the upper envelope, and local minima of interpolation form the lower envelope;

[0019] 3) Calculate the mean of the upper and lower envelopes;

[0020] 4) Subtract the average value from the signal to generate a new signal; if this new signal meets the conditions to be an intrinsic mode function, it is considered an intrinsic mode function; an intrinsic mode function is a special type of signal that satisfies the same number of extrema as the number of zero crossings, and the average value of the upper and lower envelopes must be close to zero;

[0021] 5) Subtract this intrinsic mode function from the original signal, and repeat this process for the remaining signal to extract the next intrinsic mode function.

[0022] The encoder incorporates four max-pooling layers to reduce the input length. The encoder includes one empirical mode decomposition module, eight residual shrinking modules, four max-pooling layers, two LSTM layers, and one transformer layer. The empirical mode decomposition module, the first residual shrinking module, the first max-pooling layer, the second residual shrinking module, the second max-pooling layer, the third residual shrinking module, the third max-pooling layer, the fourth residual shrinking module, the fourth max-pooling layer, the fifth residual shrinking module, the sixth residual shrinking module, the seventh residual shrinking module, the eighth residual shrinking module, the first LSTM layer, the second LSTM layer, and the transformer layer are connected sequentially.

[0023] The decoder consists of four upsampling layers and five 1D convolutional layers, which are connected sequentially: the first upsampling layer, the first 1D convolutional layer, the second upsampling layer, the second 1D convolutional layer, the third upsampling layer, the third 1D convolutional layer, the fourth upsampling layer, the fourth 1D convolutional layer, and the fifth 1D convolutional layer. A batch normalization layer is added after each convolutional layer, and a linear rectified function is used as the activation function.

[0024] In step three, the loss function is the binary cross-entropy loss function, as shown in formula (2):

[0025]

[0026] It can measure the difference between the probability distribution of the model output and the label value; where BCELoss is the binary cross-entropy loss function; N is the number of sampling points in the probability sequence; y i p represents the probability value at the i-th position of the true probability sequence label. i This is the probability value at the i-th position in the probability sequence predicted by the neural network.

[0027] In step three, the deep neural network model is trained, and data augmentation is performed on the training data, including random translation of records and addition of Gaussian white noise. The data augmentation methods include random translation of records with a 99% probability and adding Gaussian white noise with zero mean and a standard deviation of 0.01 to 0.15 times the maximum amplitude of the record to the records with a 50% probability.

[0028] The batch size of the trained deep neural network model is 64, the number of training epochs is 25, the initial learning rate is 0.01, and the learning rate is reduced to 1 / 10 of the original when the loss on the validation set remains unchanged for 3 consecutive epochs.

[0029] In step one, each ground motion data point is standardized on each channel and used as the final network input data.

[0030] The standardized operation is shown in formula (1):

[0031]

[0032] Where: x * σ represents the standardized data; μ represents the mean of the data; σ represents the standard deviation of the data; and x represents the data before standardization.

[0033] The beneficial effects of this invention are:

[0034] This invention proposes a deep neural network based on time-frequency analysis for picking up the arrival times of P-waves and S-waves. The network consists of an encoder and a decoder, incorporating an empirical mode decomposition (EMD) module to convert seismic motions into intrinsic mode functions (EMFs). It also utilizes multiple residual contraction modules, LSTM, and transformer modules to improve the network's ability to learn spatial and temporal information from seismic motions. The main inventive aspect of this invention lies in adding an EMD module to the traditional neural network architecture. This module expresses the time-frequency characteristics of seismic motions by decomposing them into EMFs of different frequencies. When there are differences in the frequencies of noise, P-waves, and S-waves, EMD can effectively separate them, thus helping the network better identify P-waves and S-waves. Even for broadband noise, the noise is decomposed into different EMFs, thereby reducing the noise content in the EMFs containing P-waves and S-waves. This time-frequency analysis method can effectively improve the network's noise resistance.

[0035] Compared with current deep learning picking methods, the seismic phase picking method based on time-frequency analysis and deep learning in this invention has higher picking accuracy and better picking performance on noisy records with low signal-to-noise ratio. Attached Figure Description

[0036] Figure 1 The selected earthquake epicenters are shown in the image, with red circles representing the epicenters of the earthquakes.

[0037] Figure 2 This is a distribution map of magnitude and signal-to-noise ratio in the dataset.

[0038] Figure 3 This is a structural diagram of the residual shrinkage module of the network in the embodiment.

[0039] Figure 4 This is a schematic diagram of the deep neural network in the embodiment.

[0040] Figure 5 This example compares the P-wave and S-wave picking results of the deep neural network model and the current best picking method EQTransformer on the test set at different signal-to-noise ratios. Detailed Implementation

[0041] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0042] A seismic phase picking method based on time-frequency analysis and deep learning includes the following steps:

[0043] Step 1: Collect seismic motion records and build a dataset:

[0044] The P-waves and S-waves given in the database are picked up as the corresponding label values ​​of the seismic records to obtain the seismic phase picking dataset.

[0045] Step 2: Establish a deep neural network model incorporating the Empirical Mode Decomposition (EMD) module:

[0046] Deep neural networks are based on a fully convolutional network architecture, consisting of an encoder and a decoder. The encoder extracts features from the three-component seismic motion input to generate a high-dimensional feature map, and the decoder transforms the high-dimensional feature map into a one-dimensional probability sequence. The probability value of each point represents the probability that the point is reached by P-waves or S-waves.

[0047] The encoder's first part is an Empirical Mode Decomposition (EMD) module, which decomposes seismic motion into multiple intrinsic mode functions (EMFs), each representing a different frequency component of the seismic motion. The second part consists of multiple residual shrinking modules cascaded together, which add a soft thresholding function to the traditional residual module to filter out background noise. The third part comprises two Long Short-Term Memory (LSTM) layers for extracting temporal contextual information. The fourth part is a transformer module that assigns weights to different positions in the sequence based on an attention mechanism to amplify key features. Four max-pooling layers are added to the encoder to reduce the input length.

[0048] The decoder consists of multiple upsampling and convolutional layers, which expand the feature map to the same length as the input. The network uses Rectified Linearity (ReLU) as the activation function, and all convolutional layers are 1D convolutions. A batch normalization layer is added after each convolutional layer.

[0049] The loss function used is the binary cross-entropy loss function, which can measure the difference between the probability distribution of the model output and the label value.

[0050] The empirical mode decomposition of the encoder can be achieved through the following steps:

[0051] 1) Identify all local maxima and local minima in the signal.

[0052] 2) Local maxima of interpolation form the upper envelope, and local minima of interpolation form the lower envelope.

[0053] 3) Calculate the mean of the upper and lower envelopes.

[0054] 4) Subtract the average value from the signal to generate a new signal. If this new signal meets the conditions for becoming an eigenmode function, it is considered an eigenmode function. An eigenmode function is a special type of signal that satisfies a similar number of extrema as the number of zero crossings, and the average value of its upper and lower envelopes must be close to zero.

[0055] 5) Subtract this intrinsic mode function from the original signal, and repeat this process for the remaining signal to extract the next intrinsic mode function.

[0056] The final empirical mode decomposition decomposes the signal into five eigenmode functions and one remaining seismic signal component.

[0057] Step 3: Training the deep neural network model:

[0058] A deep neural network was trained using a seismic phase picking dataset. Data augmentation was performed on the training data, including recording random translations and adding Gaussian white noise. Training accuracy was ensured by using a loss function and a mean absolute error (MAE) function, allowing the loss function to gradually decrease, resulting in the trained deep neural network model.

[0059] Step 4: Seismic Phase Picking

[0060] The three-component ground motion records are input into the deep neural network model trained in step three, and the output probability sequences corresponding to P-waves and S-waves are obtained. The location of the maximum value of the probability sequence is found, thereby realizing phase picking.

[0061] Example 1

[0062] This embodiment of the seismic phase picking method based on time-frequency analysis and deep learning is implemented according to the following steps:

[0063] Step 1: Collect seismic motion records and build a dataset:

[0064] 1,029,421 three-component ground motion records were selected from the STEAD database, and the epicenter distribution is as follows: Figure 1 As shown, magnitude and signal-to-noise ratio are similar. Figure 2As shown. The P-wave and S-wave values ​​from the database are used as the label values ​​for the corresponding ground motion data. Based on the label values, a one-dimensional probability sequence with the same recording time length is generated, with the probability of the actual arrival point set to 1, linearly decreasing to 0 within a range of 10 surrounding sampling points. This yields the earthquake phase-picking dataset. For each ground motion data point, a standardization operation is performed on each channel, which serves as the final network input data.

[0065] The standardized operation is shown in formula (1):

[0066]

[0067] Where: x * σ represents the standardized data; μ represents the mean of the data; σ represents the standard deviation of the data; and x represents the data before standardization.

[0068] The seismic motion dataset was randomly divided into a training set, a validation set, and a test set, with a ratio of 17:1:2.

[0069] Step 2: Establish a deep neural network combining time-frequency analysis:

[0070] like Figure 4 The deep neural network shown includes an encoder and a decoder. The encoder comprises one Empirical Mode Decomposition (EMD) module, eight residual shrinking modules, four max-pooling layers, two LSTM layers, and one transformer layer. The EMD module, first residual shrinking module, first max-pooling layer, second residual shrinking module, second max-pooling layer, third residual shrinking module, third max-pooling layer, fourth residual shrinking module, fourth max-pooling layer, fifth residual shrinking module, sixth residual shrinking module, seventh residual shrinking module, eighth residual shrinking module, first LSTM layer, second LSTM layer, and transformer layer are sequentially connected. The decoder comprises four upsampling layers and five 1D convolutional layers. The first upsampling layer, first 1D convolutional layer, second upsampling layer, second 1D convolutional layer, third upsampling layer, third 1D convolutional layer, fourth upsampling layer, fourth 1D convolutional layer, and fifth 1D convolutional layer are sequentially connected. A batch normalization layer is added after each convolutional layer, and the ReLU function is used as the activation function, such as... Figure 3 As shown.

[0071] Step 3: Training the deep neural network model:

[0072] The deep neural network model was trained, and the training data was augmented, including randomly shifting the records with a 99% probability and adding Gaussian white noise with zero mean and a standard deviation of 0.01 to 0.15 times the maximum amplitude of the records with a 50% probability.

[0073] Training accuracy is ensured by using loss function and mean absolute error function. The batch size of the network model is 64, the number of training epochs is 25, the initial learning rate is 0.01, and when the loss of the validation set remains unchanged for 3 consecutive epochs, the learning rate is reduced to 1 / 10 of the original, thus obtaining a trained neural network model based on time-frequency analysis.

[0074] The loss function is shown in equation (2):

[0075]

[0076] Where: BCELoss is the binary cross-entropy loss function; N is the number of sampling points in the probability sequence; y i p represents the probability value at the i-th position of the true probability sequence label. i This is the probability value at the i-th position in the probability sequence predicted by the neural network.

[0077] Step 4: Seismic Phase Picking

[0078] The three-component ground motion records are input into the deep neural network model trained in step four, which outputs the arrival time probability sequences of P-waves and S-waves. The network uses the Adam adaptive optimization function for backpropagation, and the binary cross-entropy loss function is used as the loss function. The mean absolute error (MAE), standard deviation (SD), precision, recall, and F1 score are used as evaluation indicators to identify the location of the maximum value of the probability sequence, thereby realizing the pickup of the arrival time of P-waves and S-waves.

[0079] With the development of computer technology and deep learning, seismic data is gradually increasing, the number of layers in neural networks is also increasing, and the nonlinear fitting ability is becoming stronger. Various network models are used for seismic phase picking. However, these neural network methods often face insufficient noise resistance and struggle to handle low signal-to-noise ratio records. This invention compares these methods with current deep learning models; the comparison results are as follows... Figure 5 As shown in Table 1, the deep neural network model trained in this embodiment performs better. In the table, MAE and SD represent the mean absolute error and standard deviation of phase picking, respectively.

[0080] Table 1 Comparison of prediction results between the present invention and current deep learning models.

[0081]

[0082] In the table, EQTransformer represents the deep learning model mentioned in the paper published in 2020 entitled "Earthquake transformer—an attentive deep-learning model for simultaneous earthquake detection and phasepicking."

[0083] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and not restrictive.

Claims

1. A seismic phase picking method based on time-frequency analysis and deep learning, characterized in that, Includes the following steps: Step 1: Collect seismic motion records and build a dataset: Select known ground motion data, and use the P-wave and S-wave picked up from the ground motion data as the corresponding label values ​​of the earthquake record to obtain the phase picking dataset. Step 2: Establish a deep neural network model incorporating the Empirical Mode Decomposition (EMD) module: A deep neural network consists of an encoder and a decoder. The first part of the encoder is an empirical mode decomposition module, which decomposes earthquake motion into multiple intrinsic mode functions, each of which represents a different frequency component in the earthquake motion. The second part consists of multiple residual shrinking modules connected in series. The residual shrinking module adds a soft threshold function to the traditional residual module to filter out background noise. The third part consists of two long short-term memory layers, which are used to extract temporal contextual information; the fourth part is the transformer module, which assigns weights to different positions in the sequence based on the attention mechanism to amplify the features of key parts. Four max pooling layers are added to the encoder to reduce the input length; the encoder includes one empirical mode decomposition module, eight residual shrinking modules, four max pooling layers, two LSTM layers and one transformer layer; the empirical mode decomposition module, the first residual shrinking module, the first max pooling layer, the second residual shrinking module, the second max pooling layer, the third residual shrinking module, the third max pooling layer, the fourth residual shrinking module, the fourth max pooling layer, the fifth residual shrinking module, the sixth residual shrinking module, the seventh residual shrinking module, the eighth residual shrinking module, the first LSTM layer, the second LSTM layer and the transformer layer are connected in sequence; The decoder consists of multiple upsampling layers and convolutional layers, which expand the feature map to the same length as the input. Step 3: Training the deep neural network model: The deep neural network is trained using a loss function and a mean absolute error function to ensure training accuracy, gradually decreasing the loss function to obtain the trained deep neural network model. Data augmentation is performed on the training data, including random shifting of records and adding Gaussian white noise. The data augmentation methods include randomly shifting records with a 99% probability and adding Gaussian white noise with zero mean and a standard deviation of 0.01 to 0.15 times the maximum amplitude of the record with a 50% probability. Step 4: Seismic Phase Picking The three-component ground motion records are input into the deep neural network model trained in step three, and the output probability sequences corresponding to P-waves and S-waves are obtained. The location of the maximum value of the probability sequence is found, thereby realizing phase picking.

2. The seismic phase picking method based on time-frequency analysis and deep learning according to claim 1, characterized in that, In step two, the empirical mode decomposition of the encoder includes the following steps: 1) Identify all local maxima and local minima in the signal; 2) Local maxima of interpolation form the upper envelope, and local minima of interpolation form the lower envelope; 3) Calculate the mean of the upper and lower envelopes; 4) Subtract the average value from the signal to generate a new signal; if this new signal meets the conditions to be an intrinsic mode function, it is considered an intrinsic mode function; an intrinsic mode function is a special type of signal that satisfies the same number of extrema as the number of zero crossings, and the average value of the upper and lower envelopes must be close to zero; 5) Subtract this intrinsic mode function from the original signal, and repeat this process for the remaining signal to extract the next intrinsic mode function.

3. The seismic phase picking method based on time-frequency analysis and deep learning according to claim 1, characterized in that, In step two, the decoder includes four upsampling layers and five 1D convolutional layers, which are connected sequentially: the first upsampling layer, the first 1D convolutional layer, the second upsampling layer, the second 1D convolutional layer, the third upsampling layer, the third 1D convolutional layer, the fourth upsampling layer, the fourth 1D convolutional layer, and the fifth 1D convolutional layer. A batch normalization layer is added after each convolutional layer, and a linear rectified function is used as the activation function.

4. The seismic phase picking method based on time-frequency analysis and deep learning according to claim 1, characterized in that, In step three, the loss function is the binary cross-entropy loss function, as shown in formula (2): Official (2) It can measure the difference between the probability distribution of the model output and the label value; where BCELoss is the binary cross-entropy loss function; N is the number of sampling points in the probability sequence; y i p represents the probability value at the i-th position of the true probability sequence label. i This is the probability value at the i-th position in the probability sequence predicted by the neural network.

5. The seismic phase picking method based on time-frequency analysis and deep learning according to claim 1, characterized in that, In step three, the batch size of the trained deep neural network model is 64, the number of training epochs is 25, the initial learning rate is 0.01, and the learning rate is reduced to 1 / 10 of the original when the loss of the validation set remains unchanged for 3 consecutive epochs.

6. The seismic phase picking method based on time-frequency analysis and deep learning according to claim 1, characterized in that, In step one, each ground motion data point is standardized on each channel and used as the final network input data. The standardized operation is shown in formula (1): Official (1) Where: x * σ represents the standardized data; μ represents the mean of the data; σ represents the standard deviation of the data; and x represents the data before standardization.

Images