Unsupervised seismic phase pickup method based on auto-encoder

Abstract

The invention relates to the technical field of seismology and artificial intelligence, in particular to an unsupervised seismic phase pickup method based on an auto-encoder, which comprises the following steps: preprocessing an original three-component seismic waveform, and removing noise through band-pass filtering, mean removal and baseline correction; determining an optimal time window length through reconstruction precision evaluation, dividing the optimal time window length into time windows through a fixed non-overlapping segmentation method, extracting features by using an auto-encoder in combination with 2D convolution and position coding, obtaining intra-window and inter-window attention weights through a double-branch attention mechanism, and constructing a physical constraint function containing intra-window and inter-window loss; according to the technical scheme, an unsupervised objective function is established through loss minimization, a feature sensitivity matrix is calculated, an optimization equation is constructed in combination with gradient and a second derivative, parameters are solved through back propagation, abnormal scores are calculated to recognize seismic phases, and iteration is carried out to convergence, manual annotation is not needed, the precision is improved through double-branch attention and physical constraints, the real-time performance is good, and the method is suitable for large-scale popularization and application. The method is suitable for earthquake early warning.

Description

Technical Field

[0001] The present invention relates to the field of seismology and artificial intelligence technology, and in particular to an unsupervised seismic phase picking method based on an autoencoder. Background Art

[0002] Assessing the accuracy of earthquake phase picking is fundamental to ensuring the quality of earthquake catalogs. The main challenge lies in accurately identifying the arrival times of P and S waves at seismic stations. Traditionally, seismologists manually determine arrival times using established principles, such as prioritizing P waves over S waves and body waves over surface waves due to their faster propagation speeds. However, with the advancement of disaster reduction efforts and the deployment of high-density seismic networks, manual picking is no longer sufficient to process the growing amount of waveform data, especially in high-risk earthquake areas and critical infrastructure areas. Compared with P waves, the precise picking of S-wave arrival times is more challenging and places a more stringent test on phase picking algorithms. This difficulty stems from the emergent nature of S waves, which originate from the tail of P waves and are often masked by preceding energy. Nevertheless, S-wave arrival times are crucial for earthquake early warning systems, as these systems rely on the time difference between the arrival of P and S waves to trigger rapid response mechanisms.

[0003] Recent advances in machine learning have spawned powerful automated data extraction tools. Deep neural networks and Transformer models are now adept at extracting complex features from waveform data, performing exceptionally well in tasks such as denoising and phase recognition. These models can automatically extract high-dimensional features from seismic waveforms and learn complex nonlinear relationships between input and output. However, these methods rely on supervised learning, which requires large amounts of high-quality labeled seismic data. Generating such labels is labor-intensive and time-consuming. This reliance limits their applicability in many real-world engineering scenarios, particularly where labeled data is scarce or unavailable. Summary of the Invention

[0004] (1) Technical problems solved

[0005] In view of the shortcomings of the existing technology, the present invention provides an unsupervised seismic phase picking method based on an autoencoder.

[0006] (2) Technical solution

[0007] To achieve the above object, the present invention provides the following technical solution: an unsupervised seismic phase picking method based on an autoencoder of the present invention comprises the following steps:

[0008] 1) Data preprocessing of the original three-component seismic waveform using bandpass filtering, de-averaging, and baseline correction to remove instrument response and environmental noise;

[0009] 2) Segment the preprocessed waveform using various window lengths, and obtain the optimal time window length W through reconstruction accuracy evaluation, which serves as the standard window size for subsequent segmentation;

[0010] 3) A fixed-length non-overlapping segmentation method is used to divide the continuous seismic records into L time window panes to avoid redundant information between adjacent windows;

[0011] 4) The time window is embedded and encoded through the autoencoder model, and the spatiotemporal features are extracted by convolution operation and position encoding is added to obtain the feature representation of the time window;

[0012] 5) Convert the time window features to the multi-scale attention domain through the dual-branch attention mechanism to obtain the intra-window attention weight and inter-window attention weight;

[0013] 6) Construct a physical constraint loss function based on the original waveform, the reconstructed waveform, and the dual-branch attention weights, including intra-window loss and inter-window loss, and calculate the root mean square of the reconstruction error;

[0014] 7) Minimize the reconstruction loss and physical constraint loss, establish an unsupervised learning objective function, and use the dual-branch attention mechanism to calculate the feature sensitivity matrix to represent the contribution weight of the time window to phase detection;

[0015] 8) Combine the sensitivity matrix and attention weight to calculate the gradient and second-order derivative of the loss function, introduce the length normalization factor as a preconditioning operator, and construct a multi-scale parameter optimization equation;

[0016] 9) The network parameter update amount is solved by the back-propagation algorithm, the reconstructed waveform is obtained using the decoder, and the anomaly score is calculated to identify the arrival time of the seismic phase;

[0017] 10) Repeat steps 3) to 9) until the maximum number of training rounds is reached or the reconstruction loss converges, and obtain the unsupervised phase picking model.

[0018] Preferably, in step 1), the seismic data preprocessing first performs filtering on the original three-component seismic waveform, assuming that the original seismic waveform is x(time)=[x EW (time), x NS (time), x UD (time)] T , where time = 1, 2, ..., TIME, TIME represents the time length of the record. The preprocessing operation can be expressed as:

[0019] x processed (time)=Filter(x(time)-μ);

[0020] Where μ is the DC component and Filter(·) is the bandpass filtering operation.

[0021] Further preferably, in step 2), the determination of the optimal window length is achieved by evaluating the reconstruction accuracy, and the candidate window length set Windows = {8, 16, 32, 64} is set, and for each window length Window index ∈Windows calculates the reconstruction signal-to-noise ratio, which is calculated as follows:

[0022]

[0023] Where sNR(·) represents the ratio of useful signal power to noise power, and argmax(·) represents the input point at which the function output value reaches its extreme value.

[0024] Again preferably, in step 3), the time window segmentation operation segments the continuous earthquake record of length TIME into non-overlapping time windows, the cnt-th time window is defined as:

[0025]

[0026] Among them, ch∈{EW,NS,UD} represents three components, and the three-component time windows are stacked to form a Pane cnt , W optimal The calculation method comes from step 3):

[0027]

[0028] Preferably, in step 4), the feature embedding process first extracts spatial-temporal features through 2D convolution, and then adds position encoding:

[0029]

[0030] Among them, Pos cnt is the learnable position embedding vector, d model is the model dimension.

[0031] Further preferably, in step 5), the implementation of the dual-branch attention mechanism first obtains the query Query, key Key and value Value and the alternative value scale parameter σ through linear transformation, assuming a dual-branch network with Num layers of attention mechanism and the output of the num-1 layer is Z (num-1) , for the numth layer, there are:

[0032]

[0033] Length normalization factor L f In (·), for the cnt-th time window:

[0034]

[0035] The calculation formula for attention within the window is:

[0036]

[0037] Among them, Query num ,Key num ,Value num Represents the query, key, and value results from the early part of step 6), d model is the dimension of the input representation of the multi-head attention mechanism, L f (·) is the length normalization factor, and soft max is the normalized exponential function.

[0038] For the inter-window attention mechanism, first calculate the distance bias matrix. For any two distances postion1 and postion2, the calculation method of the bias matrix is ​​Bias num (postion1, postion2) are as follows:

[0039]

[0040] Among them, σ num From the formula, assume a two-branch network with Num layer attention mechanism. For the calculation result of the num layer, L f (cnt) comes from the formula length normalization factor L f In (·), for the calculation result of the cnt-th time window, after flattening the time dimension, the inter-window attention is calculated as:

[0041]

[0042] in, For the flattened Query num ,Key num ,Value num Matrix, d model is the dimension of the input representation of the multi-head attention mechanism, and soft max is the normalized exponential function.

[0043] During the feature fusion and reconstruction process, the primary fusion feature representation Hybrid is obtained after the dual-branch attention mechanism at the bum layer num :

[0044]

[0045] in, They are the intra-window attention mechanism and the inter-window attention mechanism, and their calculation method follows the previous part of step 6).

[0046] The final latent representation of the cnt-th time window It is expressed as follows:

[0047]

[0048] Among them, LayerNorm is application layer normalization, ReLU is rectified linear unit, is the weight of this layer, is the bias of this layer.

[0049] Final potential representation Can be converted into a reconstructed single-pane waveform through nonlinear transformation

[0050]

[0051] where Θ(·) represents a learnable nonlinear projection module, is the projection parameter shared by all time windows, and the final reconstructed waveform By connecting all get.

[0052] Again, preferably, in step 6), the seismic signal is inherently multi-scale: the arrival of P-waves and S-waves manifests as sharp amplitude mutations, while the subsequent propagation exhibits long-term energy decay and dispersion. To address the limitations of a single loss function, a dual loss mechanism is designed:

[0053] A. In-window loss is used to capture local mutation sensitivity. For each reconstructed time window, With the original time window Pane cnt , calculate the local reconstruction error:

[0054]

[0055] Where Length is the number of non-overlapping time windows, and its calculation method follows step 4).

[0056] To further enhance the sensitivity to local changes, the attention-weighted KL divergence loss is combined. The loss within the window is defined as:

[0057]

[0058] in, It's Loss reconstruction The normalized result of The calculation result comes from the calculation formula of attention within the window, kintra Scaling factor to balance the divergence contribution;

[0059] B. Inter-window loss is used to capture the consistency of long-range propagation, and inter-window anomaly captures the time evolution of earthquake energy during propagation. For the original time window Pane cnt and reconstruction time window Calculate the global reconstruction error:

[0060]

[0061] Where Length is the number of non-overlapping time windows, and its calculation method follows step 4).

[0062] This formula emphasizes the continuity of seismic wave propagation. The significant changes in the reconstruction error pattern between adjacent time windows indicate coherent phase arrival rather than random noise fluctuations. The inter-window attention map captures cross-temporal window dependencies by directly reconstructing the error alignment:

[0063]

[0064] in, and Bias num They are the inter-window attention force output and distance bias matrix of the th layer respectively, and their generation methods follow the calculation formula of the intra-window attention and the calculation formula of the bias matrix respectively.

[0065] Preferably, in step 7), the proposed dual-branch attention mechanism calculates the feature sensitivity matrix, which is characterized by the structured deviation of the background noise according to the earthquake phase arrival. When the model attempts to reconstruct the time window containing the phase mutation, it will cause the local reconstruction error peak. For the cnt-th reconstructed window Its abnormal score is recorded as anomaly cnt ;

[0066]

[0067] Among them, Loss int ra，cnt Reflecting the local inconsistency captured by the attention mechanism, Loss int er,cnt Considering the mismatch of long-range propagation patterns, since the earthquake phase not only causes sudden local changes but also induces coherent energy propagation in time, it is necessary to reconstruct the energy E[Error cnt ] Correction score M cnt Scaling to emphasize significant structural deviations, anomalies cnt The time windows where the correction errors are large and the baseline energy is non-negligible are highlighted, which is consistent with the physical characteristics of earthquake phase mutations.

[0068] Further preferably, in step 8), in order to detect anomalies without artificial thresholds, a data-driven method is adopted, and ANOMALY={anomaly1, anomaly2, ....} is anomaly. cnt The set of output results, define the target abnormality rate γ, define the target abnormality rate Threshold = Q 100-γ / 2 (ANOMALY), where Q 100-γ / 2 (·) represents the 100-γ / 2th percentile, and the binary detection mask is defined as:

[0069]

[0070] Where cnt represents the cnt-th time window. cnt The calculation method follows step 8), Threshold comes from the formula in step 9), and if...otherwise is a conditional judgment.

[0071] The above formula generates a sparse sequence of markers that mark candidate P-wave or S-wave intervals, and the continuously activated time windows can be post-processed to determine the precise arrival time.

[0072] Again preferably, step 9) and step 10) are the general process of neural network training, and the parameter update of the complete training cycle adopts the Adam optimizer, and the update equation is:

[0073]

[0074] Among them, par inter and par intra Represents Loss inter and Loss intra The weight parameter, Loss inter and Loss intra The calculation method follows The loss definition formula of the window is: The inter-window attention map captures the cross-time window dependency formula by directly reconstructing the error alignment. θ represents the set of all parameters in the neural network. Specifically, θ train_epoch and θ train_epoch+1 Respectively represent the parameter sets of the model during the train_epoch and train_epoch+1 rounds of training, α represents the learning rate, and the training process repeats steps 3) to 9) until one of the following termination conditions is met: Or train_epoch ≥ EPOCH, where ζ represents the minimum allowed error and EPOCH represents the maximum number of rounds of model training.

[0075] (3) Beneficial effects

[0076] Compared with the prior art, the present invention provides an unsupervised seismic phase picking method based on an autoencoder, which has the following beneficial effects:

[0077] This technical solution combines a dual-branch Transformer attention mechanism with a physical constraint loss function for seismic phase detection, significantly improving phase picking accuracy and robustness. In feature extraction, it implements local mutation detection and long-range propagation modeling based on intra-window / inter-window attention theory, improving computational efficiency and real-time performance while ensuring detection accuracy.

[0078] This technical solution redefines seismic phase picking as an anomaly detection problem. It uses reconstruction error and attention weights to derive a scoring mechanism for describing time series anomalies. Phase recognition is performed using data-driven quantile thresholds. This allows for accurate detection of P- and S-wave arrival times without manual labeling, overcoming the limitations of traditional supervised learning methods, which rely heavily on large amounts of high-quality labeled data.

[0079] In addition, by introducing dual physical constraints that consider local features within the window and continuity of propagation between windows, the convergence speed of unsupervised learning is accelerated and the detection stability is improved. To this end, we propose an unsupervised algorithm to solve this problem. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 It is a flowchart of the present invention and a schematic diagram of the system structure block diagram;

[0081] Figure 2 This is a schematic flow chart of step 2) to step 6) of the present invention; DETAILED DESCRIPTION

[0082] 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.

[0083] See also Figure 1-2 The present invention provides an unsupervised seismic phase picking method based on an autoencoder, comprising the following steps:

[0084] 1) Data preprocessing of the original three-component seismic waveform using bandpass filtering, de-averaging, and baseline correction to remove instrument response and environmental noise;

[0085] 2) Segment the preprocessed waveform using various window lengths, and obtain the optimal time window length W through reconstruction accuracy evaluation, which serves as the standard window size for subsequent segmentation;

[0086] 3) A fixed-length non-overlapping segmentation method is used to divide the continuous seismic records into L time window panes to avoid redundant information between adjacent windows;

[0087] 4) The time window is embedded and encoded through the autoencoder model, and the spatiotemporal features are extracted by convolution operation and position encoding is added to obtain the feature representation of the time window;

[0088] 5) Convert the time window features to the multi-scale attention domain through the dual-branch attention mechanism to obtain the intra-window attention weight and inter-window attention weight;

[0089] 6) Construct a physical constraint loss function based on the original waveform, the reconstructed waveform, and the dual-branch attention weights, including intra-window loss and inter-window loss, and calculate the root mean square of the reconstruction error;

[0090] 7) Minimize the reconstruction loss and physical constraint loss, establish an unsupervised learning objective function, and use the dual-branch attention mechanism to calculate the feature sensitivity matrix to represent the contribution weight of the time window to phase detection;

[0091] 8) Combine the sensitivity matrix and attention weight to calculate the gradient and second-order derivative of the loss function, introduce the length normalization factor as a preconditioning operator, and construct a multi-scale parameter optimization equation;

[0092] 9) The network parameter update amount is solved by the back-propagation algorithm, the reconstructed waveform is obtained using the decoder, and the anomaly score is calculated to identify the arrival time of the seismic phase;

[0093] 10) Repeat steps 3) to 9) until the maximum number of training rounds is reached or the reconstruction loss converges, and obtain the unsupervised phase picking model.

[0094] How Data Preprocessing Works

[0095] Noise filtering mechanism: bandpass filtering, de-averaging, baseline correction: perform basic data preprocessing on the original three-component seismic waveform data to remove instrument response and environmental noise interference. The bandpass filtering range is 0.1-10Hz, and polynomial baseline correction is performed to ensure data quality. Assume that the original seismic waveform is x(time) = [x EW (time),x NS (time), x UD (time)] T , where time = 1, 2, ..., TIME, TIME represents the time length of the record. The preprocessing operation can be expressed as:

[0096] x processed (time)=Filter(x(time)-μ);

[0097] Where μ is the DC component and Filter(·) is the bandpass filtering operation.

[0098] How optimal time window length determination works

[0099] Reconstruction signal-to-noise ratio (SNR) optimization: It is achieved through reconstruction accuracy evaluation. Let the candidate window length set Windows = {8, 16, 32, 64}, and for each window length Window index ∈Windows calculates the reconstruction signal-to-noise ratio, which is calculated as follows:

[0100]

[0101] Where SNR(.) represents the ratio of useful signal power to noise power, and argmax(·) represents the input point at which the function output value reaches its extreme value.

[0102] How Non-Overlapping Time Window Split works

[0103] Redundancy elimination strategy:

[0104] The time window splitting operation splits the continuous earthquake record of length TIME into non-overlapping time windows, the cnt-th time window is defined as:

[0105]

[0106] Among them, ch∈{EW,NS,UD} represents three components, and the three-component time windows are stacked to form a Pane cnt :

[0107]

[0108] How feature embedding and encoding work

[0109] Spatiotemporal feature extraction: Spatiotemporal features of time windows are extracted by 2D convolution, combined with learnable position encoding (POS cnt ) marks the time sequence, the formula is:

[0110]

[0111] Among them, Z cnt For feature representation, position encoding enhances sequence temporal information;

[0112] Among them, Pos cnt is the learnable position embedding vector, d model is the model dimension.

[0113] How the dual-branch attention mechanism works

[0114] In-window attention (local mutation capture): obtain the query, key, value, and alternative value scale parameter σ through linear transformation. Assuming a two-branch network with Num-layer attention mechanism, for the num-th layer, there is:

[0115]

[0116] Length normalization factor L f In (·), for the cnt-th time window:

[0117]

[0118] The calculation formula for attention within the window is:

[0119]

[0120] Among them, Query num ,Key num ,Value num Represents the query, key, and value results from the early part of step 6), d model is the dimension of the input representation of the multi-head attention mechanism, L f (·) is the length normalization factor, and soft max is the normalized exponential function.

[0121] Inter-window attention (long-range propagation modeling): Introducing the distance bias matrix (Bias) to capture the temporal continuity of waveform propagation: For any two distances position1 and position2, the calculation method of the bias matrix Bias num (postion1, postion2) are as follows:

[0122]

[0123] Among them, σ num From the formula, assume a two-branch network with Num layer attention mechanism. For the calculation result of the num layer, L f (cnt) comes from the formula length normalization factor L f In (·), for the calculation result of the cnt-th time window, after flattening the time dimension, the inter-window attention is calculated as:

[0124]

[0125] in, For the flattened Query num ,Key num ,Value num Matrix, d modelis the dimension of the input representation of the multi-head attention mechanism, and soft max is the normalized exponential function.

[0126] During the feature fusion and reconstruction process, the primary fusion feature representation Hybrid is obtained after the dual-branch attention mechanism at the numth layer num :

[0127]

[0128] in, They are the intra-window attention mechanism and the inter-window attention mechanism, and their calculation method follows the previous part of step 6).

[0129] The final latent representation of the cnt-th time window It is expressed as follows:

[0130]

[0131] Among them, LayerNorm is application layer normalization, ReLU is rectified linear unit, is the weight of this layer, is the bias of this layer.

[0132] Final potential representation Can be converted into a reconstructed single-pane waveform through nonlinear transformation

[0133]

[0134] where Θ(·) represents a learnable nonlinear projection module, is the projection parameter shared by all time windows, and the final reconstructed waveform By connecting all get;

[0135] Dual-branch attention integrates local and global features to improve the ability to identify the starting point and propagation path of the earthquake phase.

[0136] How the physical constraint loss function works

[0137] Seismic signals are inherently multi-scale: the arrival of P- and S-waves manifests as sharp amplitude changes, while subsequent propagation exhibits long-term energy decay and dispersion. To address the limitations of a single loss function, a dual loss mechanism is designed:

[0138] In-window loss (local sensitivity): For each reconstructed time window Calculate the local reconstruction error:

[0139]

[0140] Where Length is the number of non-overlapping time windows, and its calculation method follows step 4).

[0141] To further enhance the sensitivity to local changes, the attention-weighted KL divergence loss is combined. The loss within the window is defined as:

[0142]

[0143] in, It's Loss reconstruction The normalized result of The calculation result comes from the calculation formula of attention within the window, k intra Scaling factor to balance the divergence contribution;

[0144] Inter-window loss (propagation consistency): The inter-window anomaly captures the time evolution of earthquake energy during propagation. cnt and Calculate the global reconstruction error:

[0145]

[0146] Where Length is the number of non-overlapping time windows, and its calculation method follows step 4).

[0147] This formula emphasizes the continuity of seismic wave propagation. The significant changes in the reconstruction error pattern between adjacent time windows indicate coherent phase arrival rather than random noise fluctuations. The inter-window attention map captures cross-temporal window dependencies by directly reconstructing the error alignment:

[0148]

[0149] in, and Bias num They are the inter-window attention force output and distance bias matrix of the th layer respectively, and their generation methods follow the calculation formula of the intra-window attention and the calculation formula of the bias matrix respectively.

[0150] How the anomaly score calculation works

[0151] Structured deviation identification: The dual-branch attention mechanism calculates the feature sensitivity matrix, which is represented by the structured deviation of the background noise according to the earthquake phase arrival. When the model tries to reconstruct the time window containing the phase mutation, it will cause the local reconstruction error peak. For the cnt-th reconstructed window Its abnormal score is recorded as anomaly cnt ;

[0152]

[0153] Among them, Loss int ra，cnt Reflecting the local inconsistency captured by the attention mechanism, Loss int er，cnt Considering the mismatch of long-range propagation patterns, since the earthquake phase not only causes sudden local changes but also induces coherent energy propagation in time, it is necessary to reconstruct the energy E[Error cnt ] Correction score M cnt Scaling to emphasize significant structural deviations, anomalies cnt The time windows with large correction errors and non-negligible baseline energy are highlighted, which are consistent with the physical characteristics of earthquake phase abrupt changes;

[0154] When the earthquake phase arrives, the local reconstruction error surges and the long-range propagation consistency is destroyed, resulting in anomalies. cnt Significantly increased.

[0155] How data-driven thresholding works

[0156] Quantile adaptive threshold: To detect anomalies without artificial thresholds, a data-driven approach is used, where ANOMALY = {anomaly1, anomaly2, ....} is the anomaly. cnt The set of output results, define the target abnormality rate γ, define the target abnormality rate Threshold = Q 100-γ / 2 (ANOMALY), where Q 100-γ / 2 (.) represents the 100-γ / 2th percentile, and the binary detection mask is defined as:

[0157]

[0158] Where cnt represents the cnt-th time window. cnt The calculation method follows step 8), Threshold comes from the formula in step 9), and if...otherwise is a conditional judgment.

[0159] The above formula generates a sparse sequence of markers that mark candidate P-wave or S-wave intervals, and the continuously activated time windows can be post-processed to determine the precise arrival time;

[0160] Avoid manual parameter adjustment and adaptively adapt to the noise levels of different earthquake scenarios.

[0161] How parameter optimization and model training work

[0162] Preconditioned multi-scale optimization: Steps 9) and 10) are the general process of neural network training. The parameter update of the complete training cycle uses the Adam optimizer, and the update equation is:

[0163]

[0164] Among them, par inter and par intintra Represents Loss inter and Loss intra The weight parameter, Loss inter and Loss intra The calculation method follows The loss definition formula of the window is: The inter-window attention map captures the cross-time window dependency formula by directly reconstructing the error alignment. θ represents the set of all parameters in the neural network. Specifically, θ train_epoch and θ train_epoch+1 Respectively represent the parameter sets of the model during the train_epoch and train_epoch+1 rounds of training, α represents the learning rate, and the training process repeats steps 3) to 9) until one of the following termination conditions is met: Or train_epoch ≥ EPOCH, where ζ represents the minimum allowed error and EPOCH represents the maximum number of rounds of model training.

[0165] Detailed workflow

[0166] Data preprocessing

[0167] The original three-component seismic waveform data are band-pass filtered, de-averaged and baseline corrected to remove instrument response and environmental noise interference.

[0168] Window Splitting

[0169] The preprocessed seismic waveform data is segmented into various window lengths, and the optimal time window length W is obtained through reconstruction accuracy evaluation.

[0170] Fixed-length non-overlapping segmentation

[0171] The continuous seismic records are divided into L time window panes using a fixed-length non-overlapping segmentation method to avoid the problems of redundant information in adjacent windows and high computational complexity.

[0172] Feature Embedding and Encoding

[0173] The segmented time windows are feature embedded and encoded under the initial autoencoder model. Spatial-temporal features are extracted through operations such as convolution and position encoding is added to obtain feature representations of all time windows.

[0174] Dual-branch attention mechanism

[0175] The time window feature representation is converted to a multi-scale attention field through a dual-branch attention mechanism to obtain the intra-window attention weight and inter-window attention weight.

[0176] Physical Constraint Loss Function

[0177] The physical constraint loss function is established using the original waveform data, reconstructed waveform data and dual-branch attention weights, including intra-window loss and inter-window loss, and the root mean square of the reconstruction error is calculated.

[0178] Unsupervised learning objective function

[0179] The unsupervised learning objective function is established by minimizing the reconstruction loss and physical constraint loss, and the feature sensitivity matrix is ​​calculated using the dual-branch attention mechanism to reflect the contribution weights of different time windows to phase detection.

[0180] Parameter optimization equation

[0181] The sensitivity matrix and attention weight are used to calculate the gradient and second-order derivative information of the loss function, and the length normalization factor is added as a preconditioning operator to construct a preconditioned multi-scale parameter optimization equation.

[0182] Backpropagation algorithm

[0183] The network parameter update amount is solved through the back-propagation algorithm, the reconstructed waveform is obtained using the autoencoder decoder, and the anomaly score is calculated to identify the arrival time of the seismic phase.

[0184] Training process

[0185] Repeat steps 3) to 9) until the maximum number of training rounds is reached or the reconstruction loss converges, and the final unsupervised phase picking model is obtained.

[0186] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. An unsupervised seismic phase picking method based on an autoencoder, characterized in that: The following steps are involved: 1) Data preprocessing of the original three-component seismic waveform using bandpass filtering, de-averaging, and baseline correction to remove instrument response and environmental noise; 2) Segment the preprocessed waveform using various window lengths, and obtain the optimal time window length W through reconstruction accuracy evaluation, which serves as the standard window size for subsequent segmentation; 3) A fixed-length non-overlapping segmentation method is used to divide the continuous seismic records into L time window panes to avoid redundant information between adjacent windows; 4) The time window is embedded and encoded through the autoencoder model, and the spatiotemporal features are extracted by convolution operation and position encoding is added to obtain the feature representation of the time window; 5) Convert the time window features to the multi-scale attention domain through the dual-branch attention mechanism to obtain the intra-window attention weight and inter-window attention weight; 6) Construct a physical constraint loss function based on the original waveform, the reconstructed waveform, and the dual-branch attention weights, including intra-window loss and inter-window loss, and calculate the root mean square of the reconstruction error; 7) Minimize the reconstruction loss and physical constraint loss, establish an unsupervised learning objective function, and use the dual-branch attention mechanism to calculate the feature sensitivity matrix to represent the contribution weight of the time window to phase detection; 8) Combine the sensitivity matrix and attention weight to calculate the gradient and second-order derivative of the loss function, introduce the length normalization factor as a preconditioning operator, and construct a multi-scale parameter optimization equation; 9) The network parameter update amount is solved by the back-propagation algorithm, the reconstructed waveform is obtained using the decoder, and the anomaly score is calculated to identify the arrival time of the seismic phase; 10) Repeat steps 3) to 9) until the maximum number of training rounds is reached or the reconstruction loss converges, and obtain the unsupervised phase picking model.

2. The unsupervised seismic phase picking method based on an autoencoder according to claim 1, characterized in that: In step 1), the seismic data preprocessing first performs filtering on the original three-component seismic waveform, assuming that the original seismic waveform is x(time)=[x EW (time), x NS (time), x UD (time)] T , where time = 1, 2, ..., TIME, TIME represents the time length of the record. The preprocessing operation can be expressed as: x processed (time)=Filter(x(time)-μ); Where μ is the DC component and Filter(·) is the bandpass filtering operation.

3. The unsupervised seismic phase picking method based on an autoencoder according to claim 2, characterized in that: In step 2), the optimal window length is determined by evaluating the reconstruction accuracy. Let the candidate window length set Windows = {8, 16, 32, 64}, and for each window length Window index ∈Windows calculates the reconstruction signal-to-noise ratio, which is calculated as follows: Where SNR(·) represents the ratio of useful signal power to noise power, and argmax(·) represents the input point at which the function output value reaches its extreme value.

4. The unsupervised seismic phase picking method based on an autoencoder according to claim 3, characterized in that: In step 3), the time window segmentation operation segments the continuous earthquake record of length TIME into non-overlapping time windows, the cnt-th time window is defined as: Among them, ch∈{EW,NS,UD} represents three components, and the three-component time windows are stacked to form a Pane cnt , W optimal The calculation method comes from step 3):

5. The unsupervised seismic phase picking method based on an autoencoder according to claim 4, characterized in that: In step 4), the feature embedding process first extracts spatial-temporal features through 2D convolution and then adds positional encoding: Among them, Pos cnt is the learnable position embedding vector, d model is the model dimension.

6. The unsupervised seismic phase picking method based on an autoencoder according to claim 5, characterized in that: In step 5), the implementation of the dual-branch attention mechanism first obtains the query Query, key Key and value Value and the alternative value scale parameter σ through linear transformation. Assuming a dual-branch network with Num layers of attention mechanism and the output of the num-1 layer is Z (num-1) , for the numth layer, there are: Length normalization factor L f In (·), for the cnt-th time window: The calculation formula for attention within the window is: Among them, Query num ,Key num ,Value num Represents the query, key, and value results from the early part of step 6), d model is the dimension of the input representation of the multi-head attention mechanism, L f (·) is the length normalization factor, and softmax is the normalized exponential function. For the inter-window attention mechanism, first calculate the distance bias matrix. For any two distances postion1 and postion2, the calculation method of the bias matrix is ​​Bias num (postion1, postion2) are as follows: Among them, σ num From the formula, assume a two-branch network with Num layer attention mechanism. For the calculation result of the num layer, L f (cnt) comes from the formula length normalization factor L f In (·), for the calculation result of the cnt-th time window, after flattening the time dimension, the inter-window attention is calculated as: in, For the flattened Query num , Kery num ,Value num Matrix, d model is the dimension of the input representation of the multi-head attention mechanism, and softmax is the normalized exponential function. During the feature fusion and reconstruction process, the primary fusion feature representation Hybrid is obtained after the dual-branch attention mechanism at the bum layer num : in, They are the intra-window attention mechanism and the inter-window attention mechanism, and their calculation method follows the previous part of step 6). The final latent representation of the cnt-th time window It is expressed as follows: Among them, LayerNorm is application layer normalization, ReLU is rectified linear unit, is the weight of this layer, is the bias of this layer. Final potential representation Can be converted into a reconstructed single-pane waveform through nonlinear transformation where Θ(·) represents a learnable nonlinear projection module, is the projection parameter shared by all time windows, and the final reconstructed waveform By connecting all get.

7. The unsupervised seismic phase picking method based on an autoencoder according to claim 6, characterized in that: In step 6), seismic signals are inherently multi-scale: the arrival of P and S waves manifests as sharp amplitude changes, while subsequent propagation exhibits long-term energy decay and dispersion. To address the limitations of a single loss function, a dual loss mechanism is designed: A. In-window loss is used to capture local mutation sensitivity. For each reconstructed time window, With the original time window Pane cnt , calculate the local reconstruction error: Where Length is the number of non-overlapping time windows, and its calculation method follows step 4). To further enhance the sensitivity to local changes, the attention-weighted KL divergence loss is combined. The loss within the window is defined as: in, It's Loss reconstruction The normalized result of The calculation result comes from the calculation formula of attention within the window, k intra Scaling factor to balance the divergence contribution; B. Inter-window loss is used to capture the consistency of long-range propagation, and inter-window anomaly captures the time evolution of earthquake energy during propagation. For the original time window Pane cnt and reconstruction time window Calculate the global reconstruction error: Where Length is the number of non-overlapping time windows, and its calculation method follows step 4). This formula emphasizes the continuity of seismic wave propagation. The significant changes in the reconstruction error pattern between adjacent time windows indicate coherent phase arrival rather than random noise fluctuations. The inter-window attention map captures cross-temporal window dependencies by directly reconstructing the error alignment: in, and Bias num They are the inter-window attention force output and distance bias matrix of the th layer respectively, and their generation methods follow the calculation formula of the intra-window attention and the calculation formula of the bias matrix respectively.

8. The unsupervised seismic phase picking method based on an autoencoder according to claim 7, characterized in that: In step 7), the proposed dual-branch attention mechanism calculates the feature sensitivity matrix. According to the structured deviation of the earthquake phase arrival as background noise, when the model tries to reconstruct the time window containing the phase mutation, it will lead to the local reconstruction error peak. For the cnt-th reconstructed window Its abnormal score is recorded as anomaly cnt ; Among them, Loss intra，cnt Reflecting the local inconsistency captured by the attention mechanism, Loss inter,cnt Considering the mismatch of long-range propagation patterns, since the earthquake phase not only causes sudden local changes but also induces coherent energy propagation in time, it is necessary to reconstruct the energy E[Error cnt ] Correction score M cnt Scaling to emphasize significant structural deviations, anomalies cnt The time windows where the correction errors are large and the baseline energy is non-negligible are highlighted, which is consistent with the physical characteristics of earthquake phase mutations.

9. The unsupervised seismic phase picking method based on an autoencoder according to claim 8, characterized in that: In step 8), in order to detect anomalies without artificial thresholds, a data-driven approach is adopted, and ANOMALY = {anomaly1, anomaly2, ....} is anomaly. cnt The set of output results, define the target abnormality rate γ, define the target abnormality rate Threshold = Q 100-γ / 2 (ANOMALY), where Q 100-γ / 2 (·) represents the 100-γ / 2th percentile, and the binary detection mask is defined as: Where cnt represents the cnt-th time window. cnt The calculation method follows step 8), Threshold comes from the formula in step 9), and if...otherwise is a conditional judgment. The above formula generates a sparse sequence of markers that mark candidate P-wave or S-wave intervals, and the continuously activated time windows can be post-processed to determine the precise arrival time.

10. The unsupervised seismic phase picking method based on autoencoder according to claim 9, characterized in that: Steps 9) and 10) are the general process of neural network training. The parameter update of the complete training cycle uses the Adam optimizer, and the update equation is: Among them, par inter and par intra Represents Loss inter and Loss intra The weight parameter, Loss inter and Loss intra The calculation method follows The loss definition formula of the window is: The inter-window attention map captures the cross-time window dependency formula by directly reconstructing the error alignment. θ represents the set of all parameters in the neural network. Specifically, θ train_epoch and θ train_epoch+1 Respectively represent the parameter sets of the model during the train_epoch and train_epoch+1 rounds of training, α represents the learning rate, and the training process repeats steps 3) to 9) until one of the following termination conditions is met: Or train_epoch ≥ EPOCH, where ζ represents the minimum allowed error and EPOCH represents the maximum number of rounds of model training.

Images