A method and system for identifying seismic phases by combining deep learning with multi-source network fusion

Abstract

The present application relates to the technical field of earthquake monitoring and early warning, in particular to a mixed deep learning and multi-source network fusion seismic phase identification method and system, the method comprising: obtaining three-component waveforms of three types of networks of seismographs, strong motion instruments and intensity meters; filtering in real time to simulate acceleration waveforms with a preset period and a preset damping ratio; intercepting waveform segments and normalizing the maximum value; determining the seismic phase arrival time by outputting P-wave, S-wave and noise three types of probability distribution from the U-shaped network to time pickup model by sampling point; discriminating Pg, Pn and Sg attributes by a convolution network attribute discrimination model; and outputting the arrival time and attributes. The present application solves the problems of amplitude limiting and cross-instrument generalization of large earthquakes, and picks up the arrival time while giving the seismic phase attributes.

Description

Technical Field

[0001] This invention relates to the field of earthquake monitoring and early warning technology, specifically to a method and system for seismic phase identification that integrates deep learning with multi-source seismic networks. Background Technology

[0002] Accurate acquisition of earthquake phase arrival times and identification of phase attributes are prerequisites for earthquake location, focal depth estimation, and earthquake early warning parameter calculation. In recent years, deep learning pickers, represented by PhaseNet and EQTransformer, have gradually replaced manual map reading in the daily monitoring of seismic networks, improving both acquisition efficiency and consistency.

[0003] However, directly applying such pickups to early warning scenarios presents three main problems. First, in terms of instrument compatibility: most existing deep learning pickups are trained using seismograph waveforms and are geared towards routine monitoring tasks with small magnitudes. When a destructive earthquake occurs, the near-field records from the seismograph will exhibit amplitude clipping, with the top of the waveform being truncated. Consequently, the waveform features that the model relies on for identification become invalid, while early warning operations truly depend on the records from strong-motion seismometers and intensity meters.

[0004] Second, regarding data differences: the amplitude characteristics and instrument responses of the three types of instruments—seismometers, strong-motion seismometers, and intensity meters—differ significantly. The accuracy of the same model decreases significantly when transferring data between the three types of instruments. Some studies have attempted to alleviate this problem by using transfer learning, but the transfer learning scheme requires preparing labeled data separately for each type of instrument and additional training, which is costly to implement.

[0005] Third, regarding the output content: most existing solutions only output the time without distinguishing whether the picked phase belongs to Pg, Pn or Sg. The lack of attribute information will directly affect the reliability of subsequent source location and depth estimation. Misuse of phase attributes will lead to the wrong object for solving the location equation.

[0006] In addition, traditional statistical picking methods, such as long-short time averaging, are ineffective on low signal-to-noise ratio records, and are prone to missing or mispicking in waveforms where the first arrival of P is not obvious or the P-S interval is extremely short. In summary, the limitations of existing technologies are as follows: deep learning pickers, represented by PhaseNet, are only designed for single seismograph instruments and only output arrival times without distinguishing phase attributes; although the problem of cross-instrument accuracy degradation has been addressed, the mainstream solution is transfer learning, which requires preparing labeled data and additional training for each type of instrument, making it unsuitable for the scarcity of labeled data from strong-motion seismometers and intensity meters in early warning scenarios; although instrument response simulation is a well-known preprocessing method, current applications are mostly limited to offline data comparison and have not yet been used as a real-time normalization method at the entry point of the deep learning picking chain to achieve direct adaptation of the same model to three types of instruments, avoid amplitude limiting in large earthquakes, and avoid instrument-by-instrument transfer training. Summary of the Invention

[0007] (a) Technical problems to be solved To address the shortcomings of existing technologies, this invention provides a method and system for seismic phase identification that integrates deep learning with multi-source seismic networks. This method can solve three problems in existing technologies: waveform failure during large earthquakes when seismometers are limited, high transfer learning costs for cross-instrument generalization among the three types of instruments, and the impact of missing seismic phase attribute information on the reliability of subsequent positioning.

[0008] (II) Technical Solution To address the aforementioned technical problems, this invention provides a seismic phase identification method that integrates deep learning with multi-source seismic network fusion, comprising the following steps: S1, acquire the three-component waveform data recorded by the seismometer network, strong motion meter network and intensity meter network; S2, perform real-time simulation filtering on the three-component waveform data: design filter coefficients based on the transfer function of the preset target instrument, discretize them through bilinear transformation to obtain a second-order recursive filter, resample the waveform data of the three types of network to the same sampling rate, and then output the acceleration waveform of the preset target instrument response with preset period and preset damping ratio in real time by the recursive filter. S3, extract waveform segments from the simulated filtered waveform data segment by segment according to the preset sliding window, and normalize the maximum value of each waveform segment; S4. Input the normalized waveform segment into the trained arrival time picking model. The arrival time picking model adopts a U-shaped network architecture and outputs three probability distributions of P-wave, S-wave and noise for each sampling point. The sampling point corresponding to the peak value of the P-wave and S-wave probability distribution is taken to determine the arrival time of the seismic phase. S5, extract the target waveform segment centered on the arrival time of the seismic phase, input the target waveform segment into the trained attribute discrimination model, the attribute discrimination model adopts a convolutional network architecture, and outputs the attribute discrimination result that the seismic phase belongs to the Pg seismic phase, the Pn seismic phase or the Sg seismic phase; S6, output the arrival time of the seismic phase and the attribute discrimination result.

[0009] Preferably, the arrival time picking model in S4 is obtained by training with training samples. The training samples are constructed as follows: based on the arrival time of seismic phase cataloging, waveform segments from 1 second before the arrival time to 3 seconds after the arrival time of cataloging are extracted, the maximum value of the extracted waveform segments is normalized, and a probability distribution centered on the arrival time of seismic phase cataloging is generated as the label of the waveform segments.

[0010] Preferably, both the picking model and the attribute discrimination model are equipped with a batch normalization layer, a dropout layer, and L2 regularization. The training process of both models adopts an early stopping strategy, and training is terminated when the validation set loss does not decrease for a preset number of consecutive rounds.

[0011] Preferably, the U-shaped network architecture of the time-of-flight picking model includes an encoder, a decoder, and skip connections. The encoder contains four downsampling stages, each consisting of two one-dimensional convolutional layers and one downsampling layer with a stride of 2, with channel numbers of 32, 64, 128, and 256 respectively, extracting waveform features by downsampling the input waveform segments layer by layer. The decoder contains four upsampling stages, using transposed convolution to upsample layer by layer to restore temporal resolution. The skip connections concatenate the features of the corresponding layers in the encoder and decoder along the channel dimension.

[0012] Preferably, the attribute discrimination model includes multiple convolutional layers and fully connected classification layers. The input of the attribute discrimination model is a waveform segment centered on the arrival time of the seismic phase. The fully connected classification layer outputs the probabilities of three categories: Pg seismic phase, Pn seismic phase, and Sg seismic phase. The category with the highest probability is taken as the attribute discrimination result.

[0013] To address the aforementioned technical problems, this invention also proposes a seismic phase identification system that integrates deep learning with multi-source seismic networks, comprising: The multi-source waveform access module is used to acquire three-component waveform data recorded by the seismometer network, strong motion meter network, and intensity meter network. The simulation filtering module is used to perform real-time simulation filtering on the three-component waveform data: the filter coefficients are designed based on the transfer function of the preset target instrument, and a second-order recursive filter is obtained by discretization through bilinear transformation. After the waveform data of the three types of stations are uniformly resampled to the same sampling rate, the recursive filter outputs the acceleration waveform of the preset target instrument response with a preset period and preset damping ratio in real time. The sample processing module is used to extract waveform segments from the simulated filtered waveform data segment by segment according to a preset sliding window, and to normalize the maximum value of each waveform segment. The arrival time picking module has a built-in trained U-shaped network architecture arrival time picking model, which is used to output three types of probability distributions of P-wave, S-wave and noise for each sample point of the normalized waveform segment, and to determine the arrival time of the phase by taking the sampling point corresponding to the peak of the P-wave and S-wave probability distributions. The attribute discrimination module has a built-in attribute discrimination model with a trained convolutional network architecture. It is used to extract target waveform segments with the arrival time of the seismic phase as the center, and to determine whether the target waveform segments belong to the Pg seismic phase, Pn seismic phase or Sg seismic phase, so as to obtain the attribute discrimination result. The result output module is used to output the arrival time of the seismic phase and the attribute discrimination result.

[0014] The present invention also proposes a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of the aforementioned hybrid deep learning and multi-source network fusion seismic phase identification method.

[0015] The present invention also proposes an electronic device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, it implements the steps of the aforementioned hybrid deep learning and multi-source network fusion seismic phase recognition method.

[0016] (iii) Beneficial technical effects Compared with existing technologies, the beneficial effects of this invention are: 1. The waveforms of the three types of seismic networks are first uniformly simulated in real time through filtering to form an acceleration waveform with a period of 1 second and a damping ratio of 0.707. The differences in amplitude characteristics and response between instruments are smoothed out at the data input side. One model can be adapted to the three types of instruments without the need to separately label data and conduct separate migration training for each type of instrument. The problem of large earthquake amplitude limiting of seismometers is also avoided by switching to strong earthquake meter and intensity meter data.

[0017] 2. The arrival picking model adopts a U-shaped network architecture, outputs the probability distribution of P-wave and S-wave for each sampling point, and determines the arrival time based on the peak value. The picking accuracy is supported by the contextual features of the entire waveform. It can still output high-confidence results on waveforms with low signal-to-noise ratio, unclear P-wave arrival, and extremely short PS interval. Under the same test set, the P-wave picking accuracy is 96.64%, recall is 100%, F1 score is 98.28%, and the mean arrival residual is 0.081 seconds. The S-wave F1 score is 88.87%, which is significantly better than the 65.81% of the long-short time averaging method.

[0018] 3. The attribute discrimination model further discriminates Pg, Pn, and Sg attributes based on the picked-up data. The data and attributes are generated simultaneously in one operation, directly serving subsequent earthquake location and depth estimation, avoiding location errors caused by misuse of seismic phases. Attached Figure Description

[0019] Figure 1 This is a framework diagram of the seismic phase identification system of the present invention; Figure 2 This is a flowchart of the seismic phase identification method of the present invention. Detailed Implementation

[0020] The purpose of this invention is to incorporate waveform data from three types of seismic networks into a single identification link, enabling a single model to stably pick up seismic phase arrival times and provide seismic phase attributes from seismometers, strong-motion seismometers, and intensity meters. The invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0021] I. Overall System Composition like Figure 1As shown, the seismic phase identification system in this embodiment includes a multi-source waveform access module 101, a simulation filtering module 102, a sample processing module 103, an arrival time acquisition module 104, an attribute discrimination module 105, and a result output module 106. The multi-source waveform access module 101 is responsible for accessing three-component waveform data from three sources: seismometer network, strong-motion seismometer network, and intensity meter network. The simulation filtering module 102 performs real-time simulation filtering on the accessed waveforms. The sample processing module 103 completes waveform segment truncation and maximum value normalization. The arrival time acquisition module 104 has a built-in U-shaped network architecture arrival time acquisition model, responsible for providing the seismic phase arrival time. The attribute discrimination module 105 has a built-in convolutional network architecture attribute discrimination model, responsible for discriminating Pg, Pn, and Sg attributes. The result output module 106 outputs the arrival time and attributes together for use in earthquake location and early warning parameter estimation.

[0022] II. Method and Flow like Figure 2 As shown, the method is executed according to S1 to S6. S1: The multi-source waveform access module 101 acquires the three-component waveform data from three types of seismic networks. S2: The simulation filtering module 102 performs real-time simulation filtering on the waveforms, aiming to unify the recordings from different instruments into acceleration waveforms under the target instrument response with a preset period and preset damping ratio—the velocity record from the seismometer, and the acceleration records from the strong-motion seismometer and intensity meter. After this step, all are on the same instrument response benchmark, thus ensuring a consistent input distribution for subsequent models. In this embodiment, the preset period is 1 second, and the preset damping ratio is 0.707. S3: The sample processing module 103 extracts waveform segments containing the phase to be identified and performs maximum value normalization to eliminate amplitude scale differences caused by different magnitudes and epicentral distances. S4: The arrival time acquisition module 104 calculates the probability distribution of P-waves and S-waves for each sampling point of the waveform segment, and takes the sampling point corresponding to the probability peak as the arrival time of the phase. S5, the attribute discrimination module 105 uses the picked seismic phase as the object and determines whether it belongs to Pg, Pn or Sg. S6, the result output module 106 outputs the seismic phase arrival time and attribute discrimination result.

[0023] III. Real-time Implementation of Simulation Filters Real-time simulation filtering is used to uniformly convert the records from different instruments into acceleration waveforms under the response of a target instrument. Starting from the mathematical model of the target instrument's response, the following section explains step-by-step the construction of the filter and the real-time recursive method.

[0024] (1) Equations of motion of the target instrument response The mechanical motion of a standard accelerometer can be described by a single-degree-of-freedom damped vibration equation: ; The meanings of each symbol are as follows: m: Mass of the internal mass block of the accelerometer; c: Damping coefficient, which describes the magnitude of damping experienced by the system; k: Spring stiffness coefficient, describing the stiffness of the return spring of the mass block; x(t): Displacement of the mass block relative to the instrument base, which varies with time t; x'(t): The first derivative of x(t) with respect to time, i.e., the relative velocity of the mass block; x''(t): The second derivative of x(t) with respect to time, i.e., the relative acceleration of the mass block; a(t): The acceleration of the instrument base (i.e., the ground), which is the input signal of the entire filter; The negative sign on the right side of the equation indicates that the direction of the relative acceleration of the mass block is opposite to the direction of the acceleration of the ground.

[0025] Dividing both sides of the above equation by m, and introducing two key parameters: natural angular frequency and damping ratio, where the natural angular frequency is denoted as w0 and the damping ratio as ζ, ; The equations of motion can then be rewritten in a more concise form: ; The relationship between the natural period T0 and the natural angular frequency is as follows: ; In this embodiment, the natural period of the target instrument is T0 = 1 second, and the damping ratio is ζ = 0.707. Substituting these values ​​into the above formula, we can obtain the natural angular frequency w0 = 2π rad / s, approximately 6.28 rad / s.

[0026] (2) Transfer function of the target instrument Taking the Laplace transform of both sides of the equation of motion, and assuming the initial conditions are zero, we obtain the transfer function between the output X(s) and the input acceleration A(s): ;

[0027] Where s is the Laplace variable, X(s) is the Laplace transform of x(t), and A(s) is the Laplace transform of a(t).

[0028] In practical processing, it is necessary to distinguish between two types of input signals: Case 1: Input is a velocity record (seismograph). The seismograph outputs ground velocity v(t), which needs to be differentiated to convert it into acceleration a(t) = v'(t) before instrument response simulation. The differential term is multiplied by s in the Laplace domain, therefore the equivalent transfer function from velocity input to target acceleration output is: ; The second scenario: The input is an acceleration record (strong-motion seismograph, intensity meter). The output of the strong-motion seismograph and intensity meter is already the ground acceleration a(t), so no differentiation is needed. Instrument response simulation can be performed directly, and the transfer function is: ; It should be noted that Hv(s) and Ha(s) are two parallel optional paths and cannot be used simultaneously. This embodiment adopts the "numerical differentiation first, then unified simulation with Ha(s)" path described in Section (4) below, that is, the seismometer velocity record is first converted into acceleration by central difference, and then the acceleration records of the strong motion meter and intensity meter are simulated and filtered together with Ha(s); the Hv(s) path is an alternative scheme when numerical differentiation is not performed in advance, and it is not used in this embodiment.

[0029] (3) Bilinear transform discretization The transfer function H(s) described above is an expression in the continuous domain (analog domain), while digital systems can only process discrete sampled data. Therefore, it needs to be converted into a discrete domain transfer function. A bilinear transformation is used to achieve the mapping from the s-domain to the z-domain: ; in: Ts: Sampling interval, i.e., the time between two adjacent sampling points. In this embodiment, the sampling rate is 100 Hz, so Ts = 0.01 seconds; z: A complex variable in the discrete domain; z to the power of negative one (z-¹): Delay operator, which means delaying the signal by one sampling interval.

[0030] Substituting the bilinear transform formula into Hv(s) or Ha(s), and rearranging, we can obtain the standard second-order form of the discrete transfer function: ; in: b0, b1, b2: Molecular coefficients; a1, a2: Denominator coefficients; These five coefficients are entirely determined by the three parameters T0, ζ, and Ts. They can be pre-calculated and stored before system deployment and directly called during runtime.

[0031] Taking the transfer function Ha(s) of the acceleration input as an example, the derivation process using the bilinear transform is as follows: Let K = 2 / Ts, substitute s = K·(1-z-¹) / (1+z-¹) into the denominator s² + 2ζ·w0·s + w0², and expand to get: Denominator=(K²+2ζ·w0·K+w0²)+(2·w0²-2·K²)·z-¹+(K²-2ζ·w0·K+w0²)·z-² The first term of the denominator (z) 0 Normalizing the terms to 1, i.e., dividing all terms by (K² + 2ζ·w0·K + w0²), yields the values ​​of a1 and a2. Substituting the numerators similarly yields b0, b1, and b2.

[0032] In this invention, T0 = 1 second, ζ = 0.707, Ts = 0.01 seconds. Substituting these values ​​into K = 2 / 0.01 = 200, w0 = 2π ≈ 6.2832 rad / s, and the first term of the denominator D0 = K² + 2ζw0K + w0² ≈ 41816.36, the specific values ​​of all five coefficients are calculated as follows: For acceleration input, the transfer function is Ha(s), and the filter coefficients are: b0≈0.000944, b1≈0.001888, b2≈0.000944, a1≈-1.911239, a2≈0.915015.

[0033] The above coefficients can be pre-calculated and stored before system deployment, and then directly invoked during runtime.

[0034] (4) Unified resampling and velocity differentiation The original sampling rates of the three types of seismograph networks—seismometers, strong-motion seismometers, and intensity meters—may differ. For example, seismometers may use 100 Hz or 200 Hz, strong-motion seismometers may use 100 Hz or 200 Hz, and intensity meters may use 50 Hz or 100 Hz. Before simulation filtering, the three types of waveforms are resampled to the same sampling rate to eliminate the time scale inconsistencies caused by the sampling rate differences. In this embodiment, the sampling rate is uniformly set to 100 Hz.

[0035] For the velocity records output by the seismometer, after resampling and before simulation filtering, numerical differentiation is performed to convert them into acceleration records. The numerical differentiation uses the central difference method. ; in: v(n): The velocity value at the nth sampling point; v(n+1): The velocity value at the (n+1)th sampling point; v(n-1): The velocity value at the (n-1)th sampling point; a(n): The acceleration value corresponding to the nth sampling point; Ts: Sampling interval.

[0036] Central difference cannot be used for the first and last points of the sequence; forward difference and backward difference are used respectively for calculation.

[0037] (5) Real-time recursive filtering After obtaining the discrete transfer function H(z), it is transformed into a time-domain difference equation: ; The meanings of each symbol are as follows: x(n): The input sampling point at the current moment (waveform data that has been resampled and converted into acceleration); x(n-1): The previous input sampling point; x(n-2): The first two input sampling points; y(n): The output sampling point at the current moment, which is simulated as the acceleration waveform of the target instrument response; y(n-1): The previous output sampling point; y(n-2): The first two output sampling points; b0, b1, b2, a1, a2: filter coefficients pre-calculated in step (3).

[0038] The difference equation is a second-order recursive filter: upon reaching a new sampling point x(n), it is only necessary to read the previous two input values ​​x(n-1) and x(n-2) and the previous two output values ​​y(n-1) and y(n-2), and perform 5 multiplications and 4 additions and subtractions to calculate the current output y(n). The computational complexity is constant and independent of the total waveform length. It can be recursively calculated point by point in real time as the data stream arrives, without waiting for the entire waveform to be recorded, thus meeting the real-time requirements of early warning services.

[0039] After filtering, the original records of the three types of satellite networks are all converted into acceleration waveforms under the same target instrument response, so the input distribution faced by the subsequent picking model and attribute discrimination model is consistent.

[0040] IV. Training Sample Construction and Training Strategies The construction and identification of training samples share two processing steps, S2 and S3. Specifically, based on the phase cataloging arrival time, waveform segments from 1 second before to 3 seconds after the cataloging arrival time are extracted. After maximum value normalization, a probability distribution centered on the cataloging arrival time is generated as the label. A 1-second delay is left before the arrival and a 3-second delay after the arrival to ensure that background noise samples and wakewave samples near the initial arrival of P and S waves are included in the training field of view, resulting in more complete arrival time features learned by the model.

[0041] Both the time-based picking model and the attribute discrimination model structurally incorporate batch normalization and dropout layers, and apply L2 regularization. Training employs an early stopping strategy, terminating training when the validation set loss no longer decreases after a preset number of consecutive rounds. Batch normalization stabilizes the inter-layer activation distribution, dropout layers and L2 regularization suppress overfitting, and the early stopping strategy ensures the model stops training when its generalization performance is optimal. These three elements work together to guarantee stable performance on the three types of instrument data. The training data covers approximately 3200 seismic events and 18000 stations, including records from seismometers, strong-motion meters, and intensity meters, labeled based on manually cataloged time-based records. The training and validation sets are divided in an 8:2 ratio. The time-based picking model uses a weighted cross-entropy loss function, with a weight ratio of 1:1:0.3 for P-waves, S-waves, and noise. The optimizer is Adam, with an initial learning rate of 1×10⁻⁶. -3 The batch size is 64. Data augmentation is used during training, including additive Gaussian noise, random time shifts, and random channel dropouts, to improve the model's robustness to data from different instruments and with varying signal-to-noise ratios. The attribute discrimination model uses the cross-entropy loss function, and the optimizer and hyperparameters are selected simultaneously.

[0042] V. The U-shaped network architecture of the model to be picked at that time The model employs a U-shaped network architecture, consisting of an encoder, a decoder, and skip connections. The encoder comprises four downsampling stages, each consisting of two one-dimensional convolutional layers and one downsampling layer with a stride of 2. The number of channels in the four stages are 32, 64, 128, and 256, respectively. The encoder downsamples the three-component waveform segments of the input layer by layer, compressing and refining the temporal series features of the waveform layer by layer. The decoder corresponds to four upsampling stages, using transposed convolutions for upsampling, gradually restoring the temporal resolution layer by layer. Skip connections connect the corresponding layers of the encoder and decoder, concatenating the features of the encoder layer with the features of the decoder in the channel dimension. Shallow detailed features are directly concatenated to the deep output, avoiding the loss of subtle features near the initial arrival during the downsampling process. The model output layer is a 1×1 convolution with Softmax activation, providing probability distributions for P-wave, S-wave, and noise at each sampling point. The peak position of the probability distribution is the arrival time of the phase.

[0043] VI. Attribute Discrimination Model The attribute discrimination model employs a convolutional network architecture, consisting of three convolutional blocks and two fully connected classification layers. Each convolutional block contains one convolutional layer (kernel size 5), a batch normalization layer, a ReLU activation layer, and a max pooling layer, with channel numbers of 32, 64, and 128 respectively. The fully connected classification layer contains one hidden layer (128 neurons, ReLU activation) and one output layer (3 neurons, Softmax activation). Its input is a waveform segment truncated at the arrival time of the seismic phase, with a segment length of 1 second before and after the arrival time, totaling 200 sampling points. The convolutional layers extract local waveform morphological features of the seismic phase, while the fully connected classification layer outputs the probabilities of three categories: Pg, Pn, and Sg. The category with the highest probability is the attribute discrimination result. Since the timing patterns and waveform morphologies of Pg, Pn, and Sg differ, separating attribute discrimination into a dedicated model achieves higher discrimination accuracy compared to combining classification within the picking model. The two models each perform their respective functions and collaborate in a cascaded manner.

[0044] VII. Example 1: Joint Training and Testing Scenario of Three Networks This embodiment is based on a standard scenario and is deployed on the data processing server of the earthquake early warning center. Historical waveforms from the seismometer network, strong-motion seismometer network, and intensity meter network are uniformly stored in the database. Real-time simulation filtering and sample construction are performed according to the process in Section 3. During training, a picking model and an attribute discrimination model are trained. Both models are subject to batch normalization, drop-out layers, and L2 regularization. An early stopping strategy is enabled during training.

[0045] The tests were conducted on the same test set, and the overall results were as follows: P-wave pickup accuracy was 96.64%, recall was 100%, F1 score was 98.28%, and mean arrival residual was 0.081 seconds; S-wave F1 score was 88.87%. The test results on the three types of instrument data were as follows: for seismograph data, P-wave F1 score was 98.56%, S-wave F1 score was 90.23%, and mean arrival residual was 0.072 seconds; for strong-motion seismograph data, P-wave F1 score was 97.81%, S-wave F1 score was 89.15%, and mean arrival residual was 0.078 seconds; for intensity meter data, P-wave F1 score was 96.12%, S-wave F1 score was 87.34%, and mean arrival residual was 0.091 seconds. The performance difference among the three types of instruments was less than 3 percentage points, verifying the single model's cross-instrument adaptability.

[0046] Ablation experiments show that without real-time simulation filtering and directly training by mixing the three types of original waveforms, the F1 value of the P-wave drops to 92.35%, and the F1 value of the S-wave drops to 81.47%. The mean residual value then rises to 0.126 seconds, and the F1 value of the P-wave in the intensity meter data is only 88.63%, a difference of 7.5 percentage points from the simulation filtering scheme. This demonstrates that real-time simulation filtering to unify instrument response is a key step in cross-instrument generalization. Even on waveform samples with low signal-to-noise ratio, indistinct P-wave arrival, and extremely short PS-wave intervals, the model can still output high-confidence picking results.

[0047] 8. Example 2: Online Identification Scenario for Major Earthquake Early Warning This embodiment presents a differentiated scenario, with a deployment environment and data composition different from Embodiment 1: the system is deployed at the real-time stream processing node of the early warning center, and the real-time waveforms accessed are mainly from the strong-motion seismograph network and the intensity meter network, targeting online identification tasks after destructive major earthquakes. During a major earthquake, the near-field records of the seismometers are amplitude-limited. In this embodiment, the identification link runs directly on the records of the strong-motion seismographs and intensity meters—these two types of records, after real-time simulation filtering, are at the same response benchmark as the seismometer records. The model trained in Embodiment 1 can be called online without any modification. The identification process is executed in real time according to S1 to S6: for each segment of real-time waveform accessed, simulation filtering, truncation and normalization, arrival time picking, and attribute discrimination are completed immediately. The output phase arrival time and Pg, Pn, and Sg attributes are then sent to the location calculation stage to provide input for the estimation of the earthquake source location and depth. In this scenario, the system's end-to-end latency for processing a single 4-second waveform online is ≤0.5 seconds. This includes approximately 0.2 milliseconds for simulation filtering, approximately 12 milliseconds (CPU) / 3 milliseconds (GPU) for arrival-time picking model inference, and approximately 5 milliseconds (CPU) / 1.5 milliseconds (GPU) for attribute discrimination model inference. The system throughput is approximately 30,000 points / second, meeting the real-time requirements of early warning services. During a 72-hour continuous test, the system detected 127 seismic events, with a P-wave F1 value of 97.35%, an S-wave F1 value of 87.12%, and an average arrival-time residual of 0.085 seconds. These results are comparable to the offline test results in Example 1. The system maintained high confidence output even on online waveforms with low signal-to-noise ratios and extremely short PS intervals, verifying the solution's ability to support the real-time rhythm of early warning services.

[0048] IX. Comparative Example Comparative Example 1 uses the traditional long-short time averaging (STA / LTA) method to pick up phases on the same test set. Its P-wave picking F1 value is 78.32%, and its S-wave picking F1 value is 65.81%, with an average residual of 0.125 seconds. This is significantly different from the P-wave F1 value of 98.28% and the S-wave F1 value of 88.87% of the present invention. In low signal-to-noise ratio samples, the missed and incorrect picking in Comparative Example 1 is more concentrated. Comparative Example 2 is a single deep learning picking model trained on seismograph waveforms (U-shaped network picking scheme represented by PhaseNet): This scheme only outputs arrival time and does not distinguish phase attributes. On seismograph data, the F1 value of P-wave is 97.45% and the F1 value of S-wave is 87.62%, but on strong-motion seismograph data, the F1 value of P-wave drops to 85.23% and the F1 value of S-wave drops to 74.18%. On intensity meter data, the F1 value of P-wave drops further to 79.56% and the F1 value of S-wave drops to 68.34%. The performance degradation across instruments is significant, and it is necessary to supplement the target instrument type with labeled data and transfer training before it can be used.

[0049] The picking performance of the examples and comparative examples is compared in the table below: Table 1

[0050] It can be seen that the present invention has advantages in three dimensions: picking accuracy, attribute discrimination and cross-instrument adaptation.

[0051] This invention can also be implemented as a computer-readable storage medium and an electronic device. The computer-readable storage medium stores a computer program, which, when executed by a processor, performs the identification steps S1 to S6 described above. The storage medium can be a read-only memory, flash memory, magnetic disk, or optical disk, etc. The electronic device includes a memory and a processor. The computer program stored in the memory, when executed by the processor, performs the identification steps described above. The electronic device can be a server or data processing workstation in an early warning center.

[0052] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. All modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for seismic phase identification that combines deep learning with multi-source seismic networks, characterized in that, Includes the following steps: S1, acquire the three-component waveform data recorded by the seismometer network, strong motion meter network and intensity meter network; S2, perform real-time simulation filtering on the three-component waveform data: design filter coefficients based on the transfer function of the preset target instrument, discretize them through bilinear transformation to obtain a second-order recursive filter, resample the waveform data of the three types of network to the same sampling rate, and then output the acceleration waveform of the preset target instrument response with preset period and preset damping ratio in real time by the recursive filter. S3, extract waveform segments from the simulated filtered waveform data segment by segment according to the preset sliding window, and normalize the maximum value of each waveform segment; S4. Input the normalized waveform segment into the trained arrival time picking model. The arrival time picking model adopts a U-shaped network architecture and outputs three probability distributions of P-wave, S-wave and noise for each sampling point. The sampling point corresponding to the peak value of the P-wave and S-wave probability distribution is taken to determine the arrival time of the seismic phase. S5, extract the target waveform segment centered on the arrival time of the seismic phase, input the target waveform segment into the trained attribute discrimination model, the attribute discrimination model adopts a convolutional network architecture, and outputs the attribute discrimination result that the seismic phase belongs to the Pg seismic phase, the Pn seismic phase or the Sg seismic phase; S6, output the arrival time of the seismic phase and the attribute discrimination result.

2. The seismic phase identification method based on hybrid deep learning and multi-source seismic network fusion according to claim 1, characterized in that: The arrival time picking model described in S4 is obtained by training with training samples. The training samples are constructed as follows: based on the arrival time of seismic phase cataloging, waveform segments from 1 second before the arrival time to 3 seconds after the arrival time of cataloging are extracted, the maximum value of the extracted waveform segments is normalized, and a probability distribution centered on the arrival time of seismic phase cataloging is generated as the label of the waveform segments.

3. The seismic phase identification method based on hybrid deep learning and multi-source seismic network fusion according to claim 1, characterized in that: Both the time-out picking model and the attribute discrimination model are equipped with a batch normalization layer, a dropout layer, and L2 regularization. The training process of both models adopts an early stopping strategy, and training is terminated when the validation set loss does not decrease for a preset number of consecutive rounds.

4. The seismic phase identification method based on hybrid deep learning and multi-source seismic network fusion according to claim 1, characterized in that: The U-shaped network architecture of the time-of-flight picking model includes an encoder, a decoder, and skip connections. The encoder contains four downsampling stages, each consisting of two one-dimensional convolutional layers and one downsampling layer with a stride of 2, with channel numbers of 32, 64, 128, and 256 respectively, extracting waveform features by downsampling the input waveform segments layer by layer. The decoder contains four upsampling stages, using transposed convolution to upsample layer by layer to restore temporal resolution. The skip connections concatenate the features of the corresponding layers in the encoder and decoder along the channel dimension.

5. The seismic phase identification method based on hybrid deep learning and multi-source seismic network fusion according to claim 1, characterized in that: The attribute discrimination model includes multiple convolutional layers and fully connected classification layers. The input of the attribute discrimination model is a waveform segment centered on the arrival time of the seismic phase. The fully connected classification layer outputs the probabilities of three categories: Pg seismic phase, Pn seismic phase, and Sg seismic phase. The category with the highest probability is taken as the attribute discrimination result.

6. The seismic phase identification method based on hybrid deep learning and multi-source seismic network fusion according to claim 1, characterized in that: The preset period is 1 second, and the preset damping ratio is 0.

707.

7. A seismic phase identification system that integrates deep learning with multi-source seismic networks, characterized in that, include: The multi-source waveform access module (101) is used to acquire three-component waveform data recorded by the seismometer network, the strong motion meter network and the intensity meter network; The simulation filtering module (102) is used to perform real-time simulation filtering on the three-component waveform data: the filter coefficients are designed based on the transfer function of the preset target instrument, and a second-order recursive filter is obtained by discretization through bilinear transformation. After the waveform data of the three types of network stations are uniformly resampled to the same sampling rate, the recursive filter outputs the acceleration waveform of the preset target instrument response with a preset period and preset damping ratio in real time. The sample processing module (103) is used to extract waveform segments from the simulated filtered waveform data segment by segment according to a preset sliding window, and to normalize the maximum value of each waveform segment; The arrival time picking module (104) has a built-in U-shaped network architecture arrival time picking model, which is used to output the three probability distributions of P wave, S wave and noise for each sample point of the normalized waveform segment, and to determine the arrival time of the phase by taking the sampling point corresponding to the peak of the P wave and S wave probability distribution. The attribute discrimination module (105) has a built-in attribute discrimination model with a trained convolutional network architecture. It is used to extract the target waveform segment with the arrival time of the seismic phase as the center, and to determine whether the target waveform segment belongs to the Pg seismic phase, Pn seismic phase or Sg seismic phase, so as to obtain the attribute discrimination result. The result output module (106) is used to output the arrival time of the seismic phase and the attribute discrimination result.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the phase identification method of hybrid deep learning and multi-source network fusion as described in any one of claims 1 to 6.

9. An electronic device, characterized in that, It includes a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, it implements the steps of the phase identification method of hybrid deep learning and multi-source network fusion as described in any one of claims 1 to 6.

Images