A focal mechanism joint inversion method based on polarity similarity and full waveform

Abstract

The application relates to the technical field of seismic information processing, and discloses a focal mechanism joint inversion method based on polarity similarity and full waveform, which comprises the following steps: obtaining a reference seismic data set, actual seismic observation data and a theoretical full waveform data set; preprocessing the data; constructing a loss function; taking a P-wave initial motion polarity window of the reference seismic data set and a theoretical P-wave initial motion polarity window as inputs of a twin neural network, taking P-wave initial motion polarity similarity as output, and training the twin neural network; performing dynamic time shift alignment according to observation waveforms and theoretical waveforms, and calculating total waveform residuals; introducing a scale matching factor, and constructing a joint objective function according to the total waveform residuals and a polarity constraint term; and performing global optimization by minimizing the joint objective function based on a neighborhood algorithm to obtain a focal mechanism parameter combination. Thus, fine inversion of focal mechanisms of small and medium earthquakes can be realized.

Description

Technical Field

[0001] This application relates to the field of earthquake information processing technology, and for example to a method for joint inversion of source mechanisms based on polarity similarity and full waveform. Background Technology

[0002] Earthquake disasters are characterized by their suddenness, wide destructive range, and long chains of secondary disasters. Earthquakes and the secondary disasters they trigger, such as landslides and collapses, can lead to casualties and economic losses. Focal mechanism solutions are crucial parameters for understanding regional tectonic stress fields, fault kinematics, and earthquake generation mechanisms. With the development of modern digital seismic observation technology and improved computing power, focal mechanism inversion methods have been continuously developed. Existing methods mainly utilize the initial motion polarity and waveform information of seismic waves, and achieve parameter inversion through linear or nonlinear methods, including: traditional initial motion polarity and amplitude ratio methods, full waveform fitting methods, and neighborhood algorithms.

[0003] However, in practical applications, especially in regions with high mountains and canyons and complex media structures, the relevant inversion work still faces challenges, and existing methods have uncertainties in the analysis of the fine structure of the regional stress field: (1) Inversion multiple solutions under complex media: For areas with huge topographic relief, the underground velocity structure changes drastically in the lateral direction. Relying solely on long-period waveform inversion (such as the CAP method) is easily affected by velocity model interference, resulting in multiple solutions or large errors in the source mechanism solution. It is difficult to lock the physical true solution by relying solely on waveform fitting. (2) Errors caused by polarity-constrained signal-to-noise ratio interference: Existing methods typically rely on manual visual interpretation, which introduces a degree of subjectivity. Furthermore, traditional inversion algorithms (such as HASH and FPFIT) often treat polarity as a binary classification pattern of positive or negative polarity. In cases of suboptimal signal-to-noise ratio, misjudgment of polarity at a single key station can lead to inversion results deviating from the true value. (3) Insufficient characterization of regional stress field: Existing studies are mostly based on the inversion of average stress field from relatively coarse source parameters, which cannot finely distinguish the subtle differences in stress state between the main rupture surface and the branch fracture.

[0004] It should be noted that the information disclosed in the background section above is only used to enhance the understanding of the background of this application, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] To provide a basic understanding of some aspects of the disclosed embodiments, a brief summary is given below. This summary is not intended as a general commentary, nor is it intended to identify key / important components or describe the scope of protection of these embodiments, but rather as a prelude to the detailed description that follows.

[0006] This disclosure provides a method for joint inversion of source mechanisms based on polarity similarity and full waveform, so as to achieve fine inversion of source mechanisms of small and medium-sized earthquakes.

[0007] In some embodiments, the joint inversion method for source mechanisms based on polarity similarity and full waveform includes: S10, acquiring a manually labeled reference seismic dataset, actual seismic observation data of the target seismic region, and a theoretical full waveform dataset; S20, performing a first preprocessing on the reference seismic dataset to obtain a P-wave initial motion polarity window; performing a second preprocessing on the theoretical full waveform dataset to obtain a theoretical P-wave initial motion polarity window and a theoretical S-wave window; pairing the P-wave initial motion polarity window of the reference seismic dataset and the theoretical P-wave initial motion polarity window to obtain a training dataset; performing a third preprocessing on the actual seismic observation data to obtain a P-wave initial motion polarity window of the actual observation dataset and a S-wave window of the actual observation dataset; S30, constructing a loss function based on the margin constraint of negative sample pairs and the similarity predicted by the Siamese neural network, wherein the negative sample pairs are waveforms with opposite polarities; using the P-wave initial motion polarity window of the reference seismic dataset and the theoretical P-wave initial motion polarity window as inputs to the Siamese neural network, and then... The initial motion polarity similarity is used as the output, and the loss function is used to train the Siamese neural network; S40, the trained Siamese neural network is evaluated; if the evaluation meets the preset requirements, the P-wave initial motion polarity window of the actual observation dataset and the theoretical P-wave initial motion polarity window are paired, and the similarity output by the Siamese neural network is used as the polarity constraint term in the joint inversion; S50, the P-wave initial motion polarity window of the actual observation dataset and the theoretical P-wave initial motion polarity window are dynamically time-shifted and the P-wave window fitting residual between the two waveforms is calculated, the S-wave window of the actual observation dataset and the theoretical S-wave window are dynamically time-shifted and the S-wave window fitting residual between the two waveforms is calculated, and the P-wave window fitting residual and the S-wave window fitting residual are added to obtain the total waveform residual; a scale matching factor is introduced, and a joint objective function is constructed according to the total waveform residual and the polarity constraint term; S60, the joint objective function is minimized based on the neighborhood algorithm to perform global optimization to obtain the source mechanism parameter combination.

[0008] The present disclosure provides a joint inversion method for source mechanisms based on polarity similarity and full waveform, which can achieve the following technical effects: 1. Introduce a Siamese neural network to extract the similarity between the P-wave initial motion polarity window in the actual observation dataset and the theoretical P-wave initial motion polarity window: Unlike directly using AI models to evaluate polarity, this method inputs the P-wave initial motion polarity windows from the actual observation dataset and the theoretical P-wave initial motion polarity windows into a Siamese neural network to obtain P-wave initial motion polarity similarity. It replaces the traditional binary polarity label with a continuous similarity value in the 0-1 range, which to some extent reduces the susceptibility of traditional polarity interpretation to noise and velocity model errors, as well as the inversion instability that may result from misjudgments at key stations. This provides a robust constraint method for focal mechanism inversion. Simultaneously, it reduces the subjectivity and efficiency bottlenecks of manual selection. Mathematically, it transforms the discrete polarity sign constraint into a continuous, differentiable, noise-resistant, and fault-tolerant soft constraint of similarity scores, which can improve the stability of focal mechanism inversion to a certain extent. 2. A joint inversion objective function is constructed by integrating the P-wave initial motion polarity similarity based on a Siamese neural network with the waveform matching residual: To address the challenge of complex velocity structures, a joint constraint method is proposed that utilizes both long-period waveform fitting and AI-identified high-frequency initial motion polarity. This complementary fusion of broadband information alleviates, to some extent, the period jump and multiple solutions problems in single-waveform inversion, improves the robustness and stability of focal mechanism inversion in tectonically complex regions, and enhances the accuracy of focal mechanism inversion for small and medium-sized earthquakes.

[0009] The above general description and the description below are exemplary and illustrative only and are not intended to limit this application. Attached Figure Description

[0010] One or more embodiments are illustrated by way of example with reference to the accompanying drawings. These illustrations and drawings do not constitute a limitation on the embodiments. Elements having the same reference numerals in the drawings are shown as similar elements. The drawings are not to be scaled. And wherein: Figure 1 This is a schematic diagram of a joint inversion method for source mechanisms based on polarity similarity and full waveform provided in an embodiment of this disclosure; Figure 2 This is a schematic diagram of the 1D-CNN twin neural network architecture provided in the embodiments of this disclosure; Figure 3 This is a schematic diagram of the input and output of the Siamese neural network provided in the embodiments of this disclosure; Figure 4 This is a schematic diagram showing the confusion matrix of the Siamese neural network model provided in this embodiment on the validation set, reflecting the recognition of Similar and Dissimilar samples. Figure 5 The output score distribution histogram of the Siamese neural network model provided in this embodiment shows the distribution characteristics of the sample prediction results in different score intervals; Figure 6This is a graph showing the changes in the training set and validation set loss values ​​and accuracy with the number of training rounds during the training process of the Siamese neural network model provided in this embodiment, used to characterize the convergence and generalization ability of the model; Figure 7 This is a schematic diagram of the changes in the joint objective function and the convergence of candidate points in the Luding main earthquake under the condition of using only waveform term constraints, provided in an embodiment of this disclosure; Figure 8 This is a schematic diagram of the polarity similarity change and candidate point convergence of the Luding main shock under the condition of using only waveform term constraints, provided in the embodiments of this disclosure; Figure 9 This is a schematic diagram of the waveform fitting residual change and candidate point convergence of the Luding main shock under the condition of using only waveform term constraints, provided in the embodiments of this disclosure; Figure 10 This is a schematic diagram of the change of the joint objective function and the convergence of candidate points in the Luding main earthquake under the condition of only using polarity similarity constraint, provided in the embodiments of this disclosure; Figure 11 This is a schematic diagram of the polarity similarity change and candidate point convergence of the Luding main earthquake under the condition of only using polarity similarity constraint, provided in the embodiments of this disclosure; Figure 12 This is a schematic diagram of the waveform fitting residual change and candidate point convergence of the Luding main shock under the condition of only using polarity similarity constraint provided in the embodiments of this disclosure; Figure 13 This is a schematic diagram comparing the focal mechanism solutions of the Luding mainshock under boundary weight conditions provided in the embodiments of this disclosure; Figure 14 This is a schematic diagram of the P / S waveform window fitting of the main shock in Luding provided in an embodiment of this disclosure; Figure 15 This is a schematic diagram comparing the focal mechanism solutions under the boundary weight conditions of the main aftershocks in Luding, provided in an embodiment of this disclosure; Figure 16 This is a schematic diagram of the P / S waveform fitting of the 4.5 magnitude aftershock in Luding provided in an embodiment of this disclosure. Detailed Implementation

[0011] To provide a more detailed understanding of the features and technical content of the embodiments of this disclosure, the implementation of the embodiments of this disclosure will be described in detail below with reference to the accompanying drawings. The accompanying drawings are for illustrative purposes only and are not intended to limit the embodiments of this disclosure. In the following technical description, for ease of explanation, several details are used to provide a full understanding of the disclosed embodiments. However, one or more embodiments may still be implemented without these details. In other cases, well-known structures and devices may be simplified in their depiction to simplify the drawings.

[0012] The terms "first," "second," etc., used in the specification and accompanying drawings of this disclosure are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate for the embodiments of this disclosure described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion.

[0013] Unless otherwise stated, the term "multiple" means two or more.

[0014] In this embodiment of the disclosure, the character " / " indicates that the objects before and after it are in an "or" relationship. For example, A / B means: A or B.

[0015] The term "and / or" describes an association between objects, indicating that three relationships can exist. For example, A and / or B means: A or B, or A and B.

[0016] The term "correspondence" can refer to an association or binding relationship. The correspondence between A and B means that there is an association or binding relationship between A and B.

[0017] Combination Figure 1 As shown, this disclosure provides a joint inversion method for source mechanisms based on polarity similarity and full waveform, including: S10: Obtain a reference seismic dataset with manual annotations, actual seismic observation data of the target seismic area, and a theoretical full waveform dataset.

[0018] The reference earthquake dataset is constructed based on observation data from the China Earthquake Science Array (ChinArray). It includes three-component waveform records and features manually picked Pg phases and P-wave first motion polarity labels for the Pg phases. It contains reliable P-wave arrival time and first motion polarity information, which can provide the polarity labels required for training the Siamese neural network.

[0019] The actual earthquake observation data comes from the records of relevant stations in the earthquake sequence of the target area and directly corresponds to specific earthquake events.

[0020] Various randomly generated source mechanism models are set up to calculate theoretical full waveform datasets under different station azimuth angles, epicentral distances, depths, magnitudes, and source time functions.

[0021] S20: Perform the first preprocessing on the reference earthquake dataset to obtain the P-wave initial motion polarity window of the reference earthquake dataset; perform the second preprocessing on the theoretical full waveform dataset to obtain the theoretical P-wave initial motion polarity window and the theoretical S-wave window; pair the P-wave initial motion polarity window of the reference earthquake dataset and the theoretical P-wave initial motion polarity window to obtain the training dataset; perform the third preprocessing on the actual earthquake observation data to obtain the P-wave initial motion polarity window of the actual observation dataset and the S-wave window of the actual observation dataset.

[0022] The focus of training data processing is on constructing a representative, balanced, and noise-resistant sample set, enabling the model to learn the polarity characteristics of P-wave initial motion under different signal-to-noise ratios and develop stable similarity discrimination capabilities. Therefore, the first preprocessing of the reference earthquake dataset specifically includes the following steps: Raw data preprocessing and noisy sample construction S21, based on the original waveform records in the seismic dataset, first performs mean averaging on each record to eliminate the influence of baseline offset.

[0023] Then, a fourth-order Butterworth bandpass filter is used to filter the original waveform from 1 Hz to 20 Hz to extract the main frequency components related to initial motion identification. Noise samples are uniformly processed using a fixed frequency band at this stage to ensure consistency in the noise database construction process.

[0024] After filtering, the waveform is downsampled from 100 Hz to 50 Hz. Then, standard deviation normalization is performed.

[0025] After the above processing, the noise samples are finally constructed, providing a foundation for the synthesis of subsequent low signal-to-noise ratio training samples.

[0026] Generate labeled polarity samples S22. Seismic waveforms with initial motion polarity labels are extracted from the reference seismic dataset to construct basic labeled samples. These labeled samples are then processed, with the filtering frequency range dynamically determined based on the event magnitude; the smaller the magnitude, the smaller the filtering frequency range. This allows the filtering frequency range to adjust smoothly with magnitude changes, thus more reasonably preserving frequency information related to the initial motion among samples of different magnitudes. Finally, the labeled samples are screened, with a signal-to-noise ratio threshold used to control sample quality.

[0027] S23. To enhance the model's adaptability to slight deviations during initial motion, the basic labeled samples are filtered and a small random time perturbation is added. After the time perturbation, the waveform is reduced from 100 Hz to 50 Hz, and a window of a preset number of sampling points is extracted as the input sample, using the downsampled center point as the reference. Optionally, the preset number is 64.

[0028] S24, to expand the number of samples and balance the positive and negative polarity samples, retains both the original and flipped forms of each input sample, forming expanded labeled samples. This processing effectively supplements samples of opposite polarity without changing the relative structure of the waveform, improving the balance of different categories in the training data.

[0029] After completing the above steps, all samples were uniformly normalized using standard deviation to ensure consistent network input scale. In addition to waveform and polarity labels, the data also simultaneously retained auxiliary information such as event latitude, longitude, depth, magnitude, station location, epicentral distance, and sample signal-to-noise ratio, providing a foundation for subsequent sample statistical analysis.

[0030] Constructing low signal-to-noise ratio synthetic samples S25. High signal-to-noise ratio (SNR) records are selected from the expanded labeled samples as relatively pure seismic signals. Then, a noise window is randomly extracted from the noise sample library. The relatively pure seismic signal and the noise window are superimposed according to a given target SNR to construct a low SNR synthetic sample. Optionally, the samples from the relatively pure seismic signal source must have an SNR of not less than 20 dB in decibel form.

[0031] Sample balancing and final training set formation S26. First, randomly mix high signal-to-noise ratio (SNR) samples, low SNR samples, relatively pure samples, and low SNR synthetic samples. Then, extract positive and negative polarity samples respectively. From each of the positive and negative polarity samples, extract a predetermined number of records to form a balanced training set. To meet the needs of subsequent classification tasks, negative polarity is encoded as 0, and positive polarity as 1. Optionally, the predetermined number of records is 2000. Let the set of negative polarity samples be "D" _"-", and the set of positive polarity samples be "D" _"+", then the training set can be written as: .

[0032] S27. After class balancing is completed, the samples in the balanced training set are randomly shuffled again to avoid positive and negative samples having a fixed order in the data arrangement. All samples are randomly shuffled to obtain the final training dataset.

[0033] Thus, in the training data preparation process, in addition to routine preprocessing, it is also necessary to combine steps such as sample screening, noise construction, low signal-to-noise ratio sample synthesis and class balancing to further organize and expand the original samples.

[0034] Actual seismic observation data no longer emphasizes sample augmentation, but focuses more on waveform quality control, consistency processing of data records from different stations, and effective information extraction suitable for focal mechanism inversion. In observation data processing, a unified waveform preprocessing workflow and necessary screening criteria are used to improve data comparability and inversion constraint capabilities. Therefore, the third preprocessing of actual seismic observation data specifically includes the following steps: Three-component waveform preprocessing and resampling S28. For each valid station, its three-component SAC waveform is read, in the order of vertical component, north-south component, and east-west component. Since different components may have mean shift, baseline drift, and discontinuous record boundaries in the original records, a unified basic preprocessing is performed on each of the three components before resampling. This includes mean removal, detrending, and applying a tapered window. Detrending reduces the impact of low-frequency drift and long-term baseline changes on subsequent analysis. Applying a tapered window helps reduce boundary effects caused by abrupt changes at both ends of the record, improving the stability of subsequent resampling and window truncation. Finally, the sampling interval and record length of the three components are checked for consistency. Only when the sampling parameters of the three components are consistent, and the length remains consistent after resampling, does the station record enter the subsequent data organization stage. This process ensures that the three-component data of the same station are strictly aligned on the time axis, providing a unified data foundation for polarity analysis and waveform fitting.

[0035] S29 involves resampling the processed three-component SAC waveform. Only when the sampling parameters of the three components are consistent, and the length remains consistent after resampling, does the station record proceed to the subsequent data organization stage. This process ensures that the three-component data from the same station are strictly aligned on the time axis, providing a unified data foundation for polarity analysis and waveform fitting. After three-component waveform preprocessing and resampling, waveform information is obtained.

[0036] At that time, parameter unification and seismic phase information organization will be implemented. S210 converts the absolute phase arrival times in the polarity file (the arrival time file extracted by PhaseNet) into arrival time parameters relative to the start time of the current waveform record, thus obtaining the phase arrival times for direct use in subsequent polarity window extraction and waveform fitting. During the specific window extraction process, if it is necessary to further determine the position of the phase in the discrete waveform, the relative arrival time can be converted into a sampling point index by combining the sampling interval dt. After the above processing, the phase arrival times retain their physical meaning in seconds and can be easily mapped to the sampling sequence, providing a unified parameter basis for polarity window truncation and waveform fitting under different sampling rate conditions.

[0037] Station Spatial Information Calculation S211 extracts the station's latitude and longitude information from the head segment of the three-component SAC waveform, and calculates the epicentral distance and azimuth from the event to the station based on the event's epicenter location, thus obtaining the station's spatial location information relative to the epicenter. Let the event latitude and longitude be respectively... The station's latitude and longitude are respectively The epicentral distance from the event to the station can then be obtained through spherical geodesy calculations. and azimuth To facilitate the use of a planar approximation to represent the relative positions of stations in subsequent inversion, the planar coordinates of stations relative to the epicenter are further decomposed into east-west and north-south components. Therefore, the planar coordinates of stations relative to the epicenter can be expressed as: , Where x represents the east-west offset and y represents the north-south offset. This processing allows the spatial location of each station to be used more intuitively in subsequent focal mechanism analysis and waveform fitting calculations.

[0038] Event filtering and final observation dataset determination S212: Construct observation data using waveform information obtained from S29, phase arrival time obtained from S210, and spatial location information obtained from S211.

[0039] The study focused on small to medium-sized earthquakes with a magnitude of at least 3.0. Based on this, the observational data were further screened by considering the completeness of the phase arrival time, the quality of the three-component waveforms, and the availability of station records for each event. Ultimately, multiple earthquake events were retained to form a dataset of actual observations for subsequent joint inversion of focal mechanisms.

[0040] A second preprocessing step is performed on the theoretical full waveform dataset. This second preprocessing includes: first, using the arrival times of the theoretical P-wave and theoretical S-wave to extract the theoretical full waveform dataset, obtaining the theoretical P-wave initial motion polarity window and the theoretical S-wave window. Then, for each record, the P-wave and S-wave are filtered using bandpass filters of different frequency bands to extract the main frequency components related to initial motion identification or waveform matching.

[0041] Combination Figure 2 As shown, this method constructs a deep one-dimensional convolutional Siamese neural network architecture to address the temporal characteristics of the P-wave waveform, comprising the following three core modules: (1) One-dimensional input layer Unlike the traditional three-component full waveform input method, to highlight the kinematic characteristics of the initial polarity, the network input is defined as a single-component time series within a fixed time window, i.e., the vertical Z-component waveform. The input sequence needs to be normalized to its standard deviation before entering the network. , in, Waveform window standard deviation A small constant introduced to prevent the denominator from being zero (e.g., taking...) ).

[0042] This operation mathematically separates the amplitude difference caused by the magnitude of the source energy and the geometric diffusion, allowing the network to focus on the undulation of the waveform rather than the absolute value of the energy, thereby concentrating computing power on the polarity characteristics of the waveform.

[0043] (2) Shared Feature Encoder For feature extraction from one-dimensional time series, the shared feature encoder is constructed using a deep one-dimensional convolutional Siamese network (1D-CNN Siamese Network). This network consists of two symmetrical branches sharing weights, each containing five convolutional modules. Each convolutional module comprises the following: One-dimensional convolutional layer (Conv1D): used to extract local morphological features and translation-invariant features of waveforms on the time axis. The number of channels is set to 64 and the convolutional kernel is set to 3.

[0044] Batch Normalization: This layer is used to normalize the distribution of features in the intermediate layers, reduce the influence of internal covariate shifts, and accelerate network convergence. The number of channels is set to 64, and the number of convolutional kernels is set to 3.

[0045] Activation function layer (ReLU): Used to introduce non-linear expressive power, with 64 channels and 3 convolution kernels.

[0046] MaxPooling1D: By downsampling and compressing the feature dimension, it reduces some high-frequency noise interference while retaining representative peak and trough features. The number of channels is set to 128 and the convolution kernel is set to 3.

[0047] Dropout: This is used to reduce the interdependence between nodes by randomly zeroing some weights or outputs of the hidden layers during the learning process, thereby achieving regularization of the neural network. The number of channels is 128, the convolution kernel is set to 3, and the dropout rate is set to 0.2.

[0048] After multi-layer feature extraction and flattening, the original waveform is mapped into a dense low-dimensional feature vector. .

[0049] (3) Metric Layer Obtaining the feature vectors of the left and right branches and Next, a suitable difference metric function needs to be selected. Considering that Euclidean distance is relatively sensitive to extreme outliers, a Manhattan distance-based function is adopted here. Style difference represents the difference between two feature vectors: .

[0050] Subsequently, the difference vector is input into a fully connected decision layer with a Sigmoid activation function, and the output is compressed to the (0,1) interval. The result is the global polarity similarity term Spol between the two waveforms.

[0051] S30: A loss function is constructed based on the margin constraints of negative sample pairs and the similarity predicted by the Siamese neural network, where negative sample pairs are waveforms with opposite polarities. The P-wave initial motion polarity windows from the reference earthquake dataset and the theoretical P-wave initial motion polarity windows are used as inputs to the Siamese neural network, and the P-wave initial motion polarity similarity is used as the output. The loss function is then used to train the Siamese neural network. Based on the above deep one-dimensional convolutional Siamese neural network architecture, S30 specifically includes the following steps: S31, Construct a loss function based on the margin constraints of negative sample pairs and the similarity predicted by the Siamese neural network: , in, For the total number of sample pairs, Let i be the true label of the i-th sample pair. The similarity predicted by the Siamese neural network. The negative sample margin parameter is set.

[0052] The loss function: For positive sample pairs, the loss term prompts the network to update its parameters, thus improving the predicted similarity. The target score is set as close to 1 as possible. For negative sample pairs, margin constraints are used to widen the distance between samples in the feature space, making the predicted score approach the preset boundary of 0. By combining data augmentation strategies such as sample mining and random perturbation of the waveform time axis, the network can learn polar feature representations that are robust to noise.

[0053] S32, the training dataset is input into a deep one-dimensional convolutional Siamese neural network architecture through a one-dimensional input layer.

[0054] S33, the training dataset includes a reference earthquake dataset P-wave initial motion polarity window and a theoretical P-wave initial motion polarity window, each with a corresponding theoretical polarity label assigned based on the theoretical travel time and polarity determination results. This theoretical polarity label is determined by calculating the ray exit angle using TauP and combining it with the radiation coefficient determination method.

[0055] Based on this, a theoretical polarity window is extracted centered on the arrival time of the theoretical P-wave, and it is matched with the observed polarity window for magnitude, thus forming the input sample pair for the Siamese network. Therefore, both the observed and theoretical windows have clear polarity labels before entering the Siamese neural network, enabling the construction of sample pairs with the same and opposite polarities.

[0056] S34, combined Figure 3 As shown, the P-wave initial motion polarity window of the reference earthquake dataset and the theoretical P-wave initial motion polarity window are used as inputs and fed into the first and second branches of the Siamese neural network, respectively.

[0057] S35, after the input convolutional features are extracted by two branches respectively, the feature vectors output by each branch are first flattened, and then mapped to embedded feature vectors by a 128-dimensional fully connected layer for subsequent dual-branch feature comparison.

[0058] During the feature comparison phase, the absolute difference between the feature vectors output by the two branches is calculated using the L1 distance layer.

[0059] S36 inputs the absolute difference obtained in S34 into a fully connected layer with a Sigmoid activation function, maps it to an embedded feature vector through the fully connected layer, and outputs the P-wave initial polarity similarity between 0 and 1.

[0060] The Siamese neural network was trained using a loss function with the Adam optimizer, a learning rate of 0.0001, and a batch size of 512. The model used mean squared error (MSE) as the loss function, ensuring that the output smoothly reflects the similarity between waveforms while simultaneously determining sample pair consistency. After 200 training epochs, the loss curves of the Siamese neural network on both the training and validation sets tended to converge stably.

[0061] Thus, embedding the trained Siamese neural network as a polarity similarity scorer into the inversion framework can demonstrate advantages over traditional methods to some extent. Traditional initial polarity constraints typically employ a rigid "0 / 1" discrete matching method, where consistent polarity is denoted as 1 and inconsistent polarity as 0. This type of constraint causes the joint objective function of the inversion to exhibit strong discontinuities in the parameter space, making the optimization algorithm more susceptible to misjudgments of noise from individual stations and prone to getting trapped near local minima.

[0062] In contrast, Siamese neural networks output continuous polarity similarity. This constraint not only quantifies the consistency of polarity characteristics but also provides smoother error guidance for the optimization algorithm. Even if there is a slight phase shift between the synthesized and observed waveforms due to the simplification of the one-dimensional velocity model, the network may exhibit a certain tolerance for such deviations due to its generalization ability learned from the data. The conversion from discrete constraints to continuous constraints improves the global search stability of source mechanism inversion in noisy environments to some extent.

[0063] S40, evaluate the trained Siamese neural network; if the evaluation meets the preset requirements, pair the P-wave initial motion polarity window of the actual observation dataset with the theoretical P-wave initial motion polarity window, and use the similarity of the Siamese neural network output as the polarity constraint term in the joint inversion. Specifically, this includes the following steps: S41. To evaluate the discriminative ability of the trained Siamese neural network on location samples, a confusion matrix is ​​used to perform statistical analysis on the test set results. The true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN) of the trained Siamese neural network on the independent test set are calculated.

[0064] S42 calculates multiple evaluation metrics based on true positive, false positive, true negative, and false negative results. These metrics include precision, recall, and F1 score. Details are as follows: , , .

[0065] S43, when all the above evaluation indicators are greater than the preset threshold, the trained Siamese neural network is determined to meet the preset requirements. Optionally, the preset threshold is 0.5.

[0066] S44. If the evaluation result of the trained Siamese neural network meets the preset requirements, the similarity output of the Siamese neural network is used as the polarity constraint term in the joint inversion.

[0067] Thus, the similarity output of the Siamese neural network, serving as a polarity constraint term in the joint inversion, helps to reduce the impact of noise interference and travel time bias on polarity discrimination, and provides a more stable polarity constraint for subsequent source mechanism solutions. Therefore, the polarity similarity output of the Siamese neural network can be used as a polarity term in the joint objective function to participate in the subsequent source mechanism parameter search.

[0068] S50: Dynamically time-shift align the P-wave initial motion polarity windows of the actual observation dataset and the theoretical P-wave initial motion polarity windows, and calculate the P-wave window fitting residual between the two waveforms. Similarly, dynamically time-shift align the S-wave windows of the actual observation dataset and the theoretical S-wave windows, and calculate the S-wave window fitting residual between the two waveforms. The total waveform residual is obtained by adding the P-wave window fitting residual and the S-wave window fitting residual. A scale matching factor is introduced, and a joint objective function is constructed based on the total waveform residual and the polarity constraint term. Specifically, the following steps are included: Waveform fitting residual term aligned with dynamic time shift S51, within a preset time range, calculate the first cross-correlation function of the waveforms of the actual observation dataset P-wave initial motion polarity window and the theoretical P-wave initial motion polarity window, and the second cross-correlation function of the waveforms of the actual observation dataset S-wave window and the theoretical S-wave window.

[0069] S52, the time shift corresponding to the point where the first cross-correlation function reaches its maximum value is used as the optimal alignment time shift for the waveforms of the actual observed P-wave initial motion polarity window and the theoretical P-wave initial motion polarity window. Similarly, the time shift corresponding to the point where the second cross-correlation function reaches its maximum value is used as the optimal alignment time shift for the waveforms of the actual observed S-wave window and the theoretical S-wave window. This step can be expressed by the following expression: .

[0070] S53: After time-shift correction is completed using S51 and S52, the difference between the waveforms of the actual observed P-wave initial motion polarity window and the theoretical P-wave initial motion polarity window is calculated to obtain the P-wave window fitting residual. The difference between the waveforms of the actual observed S-wave window and the theoretical S-wave window is then calculated to obtain the S-wave window fitting residual. The first waveform fitting residual is added to the second waveform fitting residual to obtain the total waveform residual.

[0071] The P-wave window fitting residual is expressed as: .

[0072] The S-wave window fitting residuals are expressed as follows: .

[0073] Therefore, by further using the L2 norm of the waveform in the actual observation dataset as a normalization factor to normalize the total waveform residual, a dimensionless waveform error can be obtained, which can reduce the dominant role of the absolute amplitude difference between different stations on the total error.

[0074] The total waveform residual is expressed as: .

[0075] in, δ and λ are the three parameters of the focal mechanism: This indicates the fault strike. The dip angle is the fault angle (dip). Rake angle; To maximize the joint objective function, the time shift corresponding to the peak value of the cross-correlation function is required here. This is the time offset for the trial calculation; This represents the upper limit of the time-shift search range; The P-wave observation waveform extracted from the i-th component of the k-th station; The S-wave observation waveform extracted from the i-th component of the k-th station; The composite theoretical seismic waveform of the i-th component at the k-th station; The residuals are the fitting values ​​for the P-wave window. The residuals are fitted to the S-wave window. This represents the total waveform residual; The number of valid stations participating in the inversion; Representing vectors Norm; The time shift of the maximum cross-correlation value of the i-component P-wave at k stations; The time shift of the maximum cross-correlation value of the i-component S-wave at k stations; The theoretical waveform of the i-component P-wave window from k stations after cross-correlation and time-shift alignment; The theoretical waveform of the i-component S-wave window of k stations after cross-correlation and time-shift alignment; For the first The first station The optimal time shift amount obtained by each component during the dynamic alignment process.

[0076] The denominator term, as a normalization scaling factor, can weaken the impact of amplitude magnitude differences between different events and stations on the residual value, thereby improving the comparability of waveform errors between different samples. The smaller the value, the better the theoretical waveform fits the actual observed waveform in terms of amplitude characteristics and time structure. This term mainly provides constraints on the radiation characteristics of candidate source mechanisms based on the morphology of P-waves and S-waves.

[0077] Polarity similarity constraint term based on Siamese neural network When the output value of the Siamese neural network is close to 1, it indicates that the theoretical polarity window and the observed polarity window have a high degree of consistency; when the output value is close to 0, it indicates that the two are quite different or have opposite polarity trends.

[0078] S54, calculate the arithmetic mean of the similarities of all valid stations, as the global polarity similarity term: , in, This is the global polarity similarity term; The number of valid stations participating in the inversion; k is the station number; This represents a completed twin neural network. and They represent the first The observation polarity window and theoretical polarity window for each station.

[0079] Global polarity similarity essentially transforms traditional discrete polarity labels into continuous similarity metrics, thus creating a constraint mechanism in joint inversion. Under this framework, even if individual stations exhibit some degree of discrimination ambiguity, it will not immediately cause excessive perturbation to the overall joint objective function, thereby helping to improve the joint inversion's tolerance to noise and local misjudgments.

[0080] Adaptive matching factor based on the spatial fluctuation scale of initial parameters S55, calculate the global standard deviation of waveform residuals and polarity similarity for all candidate models, where the candidate models are parametric models automatically generated during inversion. The scale matching factor is then: , in, Scale matching factor; This represents the global standard deviation of the waveform residual term; This represents the global standard deviation of the polarity term.

[0081] The standard deviation is used as the scaling metric primarily because it shares the same dimensions as the original objective term, making it easy to directly reflect its fluctuation range in the parameter space and to be used for subsequent scaling. The scaling factor serves to coordinate the numerical scales between different physical objective terms by introducing it into the polarity term. This makes the overall fluctuation range of the polarity term and the waveform residual term more comparable in scale, thereby improving the synergy between the two types of constraints in joint optimization.

[0082] S56, based on the total waveform residual, global polarity similarity term, and scale matching factor, construct a joint objective function: , in, This represents the total waveform residual; is the global polarity similarity term; W is a weighting parameter used to control the relative contribution of the total waveform residual and the global polarity similarity term to the joint objective function. For example, when A value of 0.5 indicates that, assuming scale matching is achieved, the total waveform residual and the global polarity similarity term participate in the construction of the joint objective function with the same weight. This design helps improve the adaptability and stability of the joint objective function under different event scales and station coverage conditions.

[0083] The weighting parameter is the ratio of the relative contribution of the total waveform residual to the global polarity similarity term in the joint objective function; the weighting parameter is determined in the following way: Based on the magnitude of the main shock M s The actual earthquake sequence is divided into three magnitudes: magnitude 1, magnitude 2, and magnitude 3, with magnitude 1 being greater than magnitude 2, and magnitude 2 being greater than magnitude 3. Optionally, when M... s A magnitude of ≥5.0 is considered the first earthquake; when 4.0 ≤M s When the magnitude is less than 5.0, it is considered the second magnitude; when the magnitude is less than 5.0, it is considered the second magnitude. s When the magnitude is less than 4.0, it is considered the third magnitude earthquake.

[0084] When the magnitude is first, the weighting parameter is the first ratio. When the magnitude is second, the weighting parameter is the second ratio. When the magnitude is third, the weighting parameter is the third ratio. The first ratio is greater than the second ratio, and the second ratio is greater than the third ratio. Optionally, the first ratio is... The second ratio is The third ratio is .

[0085] The main consideration for adopting the above-mentioned hierarchical strategy is that larger magnitude events usually have more complete waveform information and better signal-to-noise ratios, making them suitable for using waveform fitting and first motion polarity to jointly constrain the focal mechanism. Smaller magnitude events, on the other hand, are more significantly affected by station coverage, signal-to-noise ratio, and path effects, resulting in relatively weaker constraints from waveform fitting in the parameter space. In these cases, the contribution of first motion polarity to the stability of the mechanism solution is more prominent. Therefore, retaining a smaller proportion of waveform constraints for smaller earthquakes is primarily to utilize their auxiliary constraint effect on the mechanism parameter space without significantly amplifying waveform uncertainties. Through this hierarchical scheme, the neighborhood search and joint inversion ideas established earlier are uniformly applied to the entire Luding earthquake sequence, forming a catalog of focal mechanism solutions that can be used for subsequent mechanism analysis and regional stress studies.

[0086] In this way, the waveform fitting residuals and the similarity of the P-wave initial motion polarity are jointly embedded into the joint objective function using S50, comprehensively utilizing two types of constraint information in the same inversion process to improve the stability of the source mechanism parameter solution. Regarding polarity information processing, a Siamese neural network is introduced to measure the similarity between the initial motion polarity windows of the observed and theoretical waveforms, transforming the originally discrete polarity constraints into a continuous score form. This processing helps to reduce the impact of misjudgments from individual stations on the overall inversion results under low signal-to-noise ratio conditions and also provides smoother constraint information for subsequent parameter searches.

[0087] S60 uses a neighborhood algorithm to minimize the joint objective function for global optimization to obtain the combination of source mechanism parameters. Specifically, it includes the following steps: S61, Set the search range of the parameter space, where the parameters include: direction ,inclination and sliding angle In the 0th iteration (initialization phase), Latin Hypercube Sampling (LHS) is used to generate an initial set of candidate models.

[0088] S62 calculates the joint inversion joint objective function value based on the initial candidate model set, selects the parameter model with the smaller joint objective function value as the dominant seed, and the neighborhood algorithm gradually concentrates computational resources on the region with the lower joint objective function value through a cyclical iterative process of "Voronoi cell partitioning - dominant seed point selection - intra-cell resampling", eventually converging to the global optimum. It should be noted that the neighborhood algorithm is existing technology and will not be elaborated here.

[0089] S63: When the maximum number of iterations is reached or the joint objective function converges, the iteration stops and the results are output. The maximum number of iterations is set to 30, primarily to control the total computation time of batch inversion and prevent the iteration process of individual complex events from being extended indefinitely. When the relative change of the optimal joint objective function value over N consecutive generations is less than 10... -4 At this point, the joint objective function is considered to have converged sufficiently. Optionally, N can be 3.

[0090] The output results include two types: Optimal mechanism solution: The combination of source mechanism parameters that corresponds to the minimum value of the joint objective function during the entire iteration process, reflecting the optimal solution in the sense of the joint objective function.

[0091] Central Mechanism Solution and Uncertainty Interval: The posterior probability distribution is constructed from the set of models whose joint objective function values ​​are less than 1.2 times the optimal value in the last five iterations. The statistical mean of this set is calculated to obtain the central mechanism solution, and the standard deviation of each parameter is calculated to obtain the 95% confidence interval. Simultaneously, the stability of the inversion results is quantitatively assessed by calculating the Kagan angle (a standard indicator of the geometric difference between two source mechanism solutions) between the optimal mechanism solution and the central mechanism solution: a Kagan angle less than 5° indicates high reliability of the inversion results; a Kagan angle greater than 10° indicates significant ambiguity in the event, requiring further constraints based on other observational data.

[0092] Thus, by employing an adaptive global search strategy using the neighborhood algorithm in S60, the focal mechanism parameter space is gradually narrowed and focused on sampling. Compared to exhaustive search with a fixed step size, this method can concentrate more computational resources on low-mismatch regions while maintaining global search capabilities, thereby improving search efficiency and mitigating the interference of local extrema on the inversion results to some extent.

[0093] The present disclosure provides a joint inversion method for source mechanisms based on polarity similarity and full waveform, which can achieve the following technical effects: 1. Introduce a Siamese neural network to extract the similarity between the P-wave initial motion polarity window in the actual observation dataset and the theoretical P-wave initial motion polarity window: Unlike directly using AI models to evaluate polarity, this method inputs the polarity windows of the P-wave initial motion from the actual observation dataset and the theoretical forward-modeled waveform into a Siamese neural network to obtain the P-wave initial motion polarity similarity. It replaces the traditional binary polarity label with a continuous similarity value in the 0-1 range, which to some extent reduces the susceptibility of traditional polarity interpretation to noise and velocity model errors, as well as the inversion instability that may result from misjudgments at key stations. This provides a robust constraint method for focal mechanism inversion. Simultaneously, it reduces the subjectivity and efficiency bottlenecks of manual selection. Mathematically, it transforms the discrete polarity sign constraint into a continuous, differentiable, noise-resistant, and fault-tolerant soft constraint of similarity scores, which can improve the stability of focal mechanism inversion to a certain extent.

[0094] 2. A joint inversion objective function is constructed by integrating the P-wave initial motion polarity similarity based on a Siamese neural network with the waveform matching residual: To address the challenge of complex velocity structures, a joint constraint method is proposed that utilizes both long-period waveform fitting and AI-identified high-frequency initial motion polarity. This complementary fusion of broadband information alleviates, to some extent, the period jump and multiple solutions problems in single-waveform inversion, improves the robustness and stability of focal mechanism inversion in tectonically complex regions, and enhances the accuracy of focal mechanism inversion for small and medium-sized earthquakes.

[0095] Comprehensive performance analysis: To fully verify the performance of the inversion method proposed in this application, this embodiment verifies its stability and accuracy.

[0096] (1) Stability verification: Evaluate the Siamese neural network model on an independent test set. Figures 4 to 6 It is evident that the model demonstrates good recognition ability in classification tasks. Figure 4The confusion matrix on the test set is presented. Of the samples in the "Similar" class, 1887 were correctly identified and 161 were misclassified as "Dissimilar." Of the samples in the "Dissimilar" class, 1892 were correctly identified and 156 were misclassified as "Similar." Overall, the number of correctly classified samples significantly exceeded the number of misclassified samples, indicating that the model has good discriminative ability for both classes and does not show a significant bias towards either class. Based on the confusion matrix results, the overall classification accuracy of the model is approximately 92.3%, indicating that the model can effectively distinguish between similar and dissimilar samples.

[0097] Figure 5 The distribution characteristics of the model's output scores are shown in the figure. As can be seen, the samples are mainly concentrated at the extremes close to 0 and 1, while there are fewer samples in the intermediate range. This indicates that the model can provide relatively clear discrimination results for most samples, exhibiting strong polarization and reflecting good class separability.

[0098] Figure 6 This reflects the convergence of the model during training. As the number of training epochs increases, the loss values ​​on both the training and validation sets decrease rapidly and gradually stabilize, while the accuracy on both the test and validation sets continues to rise and eventually stabilizes at a high level. The overall trends of the training and validation set curves are consistent, with little difference between them, indicating that the model did not exhibit significant overfitting during training and possesses good convergence and generalization ability.

[0099] (2) Accuracy verification: Taking the Luding mainshock as an example, we examine the convergence characteristics of the candidate solutions in the parameter space when only waveform terms are used and when only polarity terms are used.

[0100] Figures 7 to 9 The candidate point convergence results are presented when only waveform constraints are used. It can be seen that the waveform term can form a relatively concentrated low-value region in the parameter space, indicating that it has a strong constraint ability on the overall radiation characteristics of the mainshock. Figures 10 to 12 The candidate point convergence results are presented when only polarity constraints are used. Compared with waveform constraints, the polarity term can more directly distinguish the geometric quadrants of fault nodal planes, but when used alone, it lacks constraints on the overall waveform morphology and focal depth, thus failing to form a stable dominant solution region in the parameter space comparable to that of the waveform term. Therefore, it is evident that using either the waveform term or the polarity term alone is insufficient to simultaneously satisfy both the requirements of overall waveform fitting and nodal plane geometric discrimination, which is precisely the main reason why joint constraints are used for the mainshock.

[0101] Figure 13The focal mechanism solutions for the Luding mainshock are presented under the following boundary conditions: using waveform fitting as the joint objective function (with a polarity similarity weight of 0), using polarity similarity as the joint objective function (with a waveform fitting weight of 0), and under the equilibrium weight (waveform fitting: polarity similarity = 0.5:0.5). The focal mechanism solutions obtained through joint inversion are compared with previous studies, and the Kagan angle is less than 10°.

[0102] Under joint constraints, the theoretical synthesized waveform corresponding to the final optimal mechanism solution of the mainshock maintains good consistency with the observed waveform. Figure 14 The waveform fitting results for each station during the mainshock event are presented. It can be seen that within the main P-wave and S-wave windows, the synthesized results can fit the main amplitude characteristics and phase structure of the observed waveforms well, indicating that the joint constraints do not significantly weaken the effective constraint of the waveform terms on the mainshock. Meanwhile, the introduction of the polarity term improves the ability to discriminate the geometric relationships of fault nodal planes, thus making the final mechanism solution more stable than the single constraint result. Overall, the mainshock results show that when the event has a high signal-to-noise ratio and sufficient station coverage, the waveform and polarity terms can complement each other, and the joint constraint inversion can improve the stability of the mechanism solution while maintaining the waveform fitting quality.

[0103] Inversion results: for For a medium-sized event, a final inversion is performed using joint constraints under fixed depth conditions. The final optimal mechanism solution is: .

[0104] Figure 15 The waveform fitting results corresponding to the optimal mechanism solution for this event are presented. It can be seen that there is good consistency between the theoretically synthesized waveform and the observed waveform, especially within the main seismic phase window, where the forward modeling results can effectively track the main morphological changes of the observed waveform. This indicates that in this event, the waveform term remains an effective source of constraint, and the introduction of polarity similarity did not come at the expense of waveform fitting quality. On the contrary, the joint constraints enhance the stability of the fault geometry parameter solution to some extent, making the inversion results more consistent with pure waveform constraints.

[0105] Figure 16 The beach ball model is presented as the focal mechanism solution for the main aftershocks in Luding under the boundary conditions: using waveform fitting as the joint objective function (with polarity similarity weight of 0), using polarity similarity as the joint objective function (with waveform fitting weight of 0), and under the equilibrium weight (waveform fitting: polarity similarity = 0.4:0.6). The focal mechanism solution obtained by joint inversion is compared with previous studies, and the Kagan angle is less than 10°.

[0106] In summary, moderate events in Ms4.5 represent a transitional type between mainshocks and minor earthquakes in terms of methodological applicability. On the one hand, the waveform recording quality and station coverage are still sufficient to support waveform fitting constraints; on the other hand, the addition of the polarity term can further improve the stability of the mechanism solution. Therefore, the inversion results for this event indicate that for moderate earthquake events with relatively good constraints, the inversion framework combining waveform and polarity constraints can be effective, and within a reasonable weight range, it can effectively balance waveform fitting quality and mechanism solution stability.

[0107] The above are merely preferred embodiments of this application and are not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

[0108] While the specific embodiments of the present invention have been described above, they are not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art without creative effort based on the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A joint inversion method for source mechanisms based on polarity similarity and full waveform, characterized in that, include: S10: Obtain a reference seismic dataset with manual annotations, actual seismic observation data of the target seismic area, and a theoretical full waveform dataset; S20, Perform a first preprocessing on the reference earthquake dataset to obtain the P-wave initial motion polarity window of the reference earthquake dataset; The theoretical full waveform dataset is subjected to a second preprocessing step to obtain the theoretical P-wave initial polarity window and the theoretical S-wave window; The P-wave initial motion polarity window of the reference earthquake dataset and the theoretical P-wave initial motion polarity window are paired to obtain the training dataset; The actual earthquake observation data is subjected to a third preprocessing step to obtain the P-wave initial motion polarity window and the S-wave window of the actual observation dataset. S30, a loss function is constructed based on the margin constraint of negative sample pairs and the similarity predicted by the Siamese neural network, wherein the negative sample pairs are waveforms with opposite polarities; the P-wave initial motion polarity window of the reference earthquake dataset and the theoretical P-wave initial motion polarity window are used as inputs to the Siamese neural network, the P-wave initial motion polarity similarity is used as output, and the Siamese neural network is trained using the loss function; S40, evaluate the trained Siamese neural network; if the evaluation meets the preset requirements, pair the P-wave initial motion polarity window of the actual observation dataset with the theoretical P-wave initial motion polarity window, and use the similarity output of the Siamese neural network as the polarity constraint term in the joint inversion. S50, dynamically time-shift align the P-wave initial motion polarity window of the actual observation dataset and the theoretical P-wave initial motion polarity window, and calculate the P-wave window fitting residual between the two waveforms; dynamically time-shift align the S-wave window of the actual observation dataset and the theoretical S-wave window, and calculate the S-wave window fitting residual between the two waveforms; add the P-wave window fitting residual and the S-wave window fitting residual to obtain the total waveform residual; introduce a scale matching factor, and construct a joint objective function based on the total waveform residual and the polarity constraint term; S60, based on the neighborhood algorithm, the joint objective function is minimized to perform global optimization in order to obtain the combination of source mechanism parameters.

2. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 1, characterized in that, In step S30, the loss function constructed based on the margin constraint of negative sample pairs and the similarity predicted by the Siamese neural network includes: ， in, The total number of sample pairs; Let i be the true label of the i-th sample pair; The similarity predicted by the twin neural network; The negative sample margin parameter is set.

3. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 1, characterized in that, In step S30, training the Siamese neural network includes: Input the P-wave initial motion polarity window of the reference earthquake dataset into the first branch, and input the theoretical P-wave initial motion polarity window into the second branch; The feature vectors of the first branch output and the second branch output are flattened, and then the absolute difference between the feature vectors of the two branches output is calculated using the L1 distance layer. The absolute difference is input into a fully connected layer with a Sigmoid activation function, mapped to an embedded feature vector by the fully connected layer, and the output is a polar similarity between 0 and 1.

4. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 3, characterized in that, In S50, The dynamic time-shift alignment of the P-wave initial motion polarity window in the actual observation dataset and the theoretical P-wave initial motion polarity window includes: Within a preset time range, calculate the first cross-correlation function of the waveforms of the actual observation dataset P-wave initial motion polarity window and the theoretical P-wave initial motion polarity window; The time shift corresponding to when the first cross-correlation function reaches its maximum value is taken as the corresponding optimal alignment time shift; The dynamic time-shift alignment of the S-wave window of the actual observation dataset and the theoretical S-wave window includes: Within a preset time range, calculate the second cross-correlation function of the waveforms of the actual observation dataset S-wave window and the theoretical S-wave window; The time shift corresponding to when the second cross-correlation function reaches its maximum value is taken as the corresponding optimal alignment time shift.

5. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 4, characterized in that, In step S50, the introduction of a scale matching factor and the construction of a joint objective function based on the total waveform residual and the polarity constraint term include: Calculate the arithmetic mean of the similarities of all valid stations, and use it as the global polarity similarity term: ， in, This refers to the global polarity similarity term; The number of valid stations participating in the inversion; k is the station number; This represents a completed twin neural network. and They represent the first The observation polarity window and theoretical polarity window for each station; Calculate the global standard deviation of the waveform residual term and polarity similarity term corresponding to all candidate models, then the scale matching factor is: ， in, The scale matching factor; This represents the global standard deviation of the waveform residual term; This represents the global standard deviation of the polarity term; Based on the total waveform residual, the global polarity similarity term, and the scale matching factor, a joint objective function is constructed: ， Where W is the weight parameter; The total waveform residual; This refers to the global polarity similarity term.

6. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 5, characterized in that, The time shift is expressed as: ； The P-wave window fitting residual is expressed as: ； The S-wave window fitting residual is expressed as: ； The total waveform residual is expressed as: ； in, The fault strike; The fault dip angle; The sliding angle; The time shift required to maximize the joint objective function; This is the time offset calculated on a trial basis; This represents the upper limit of the time-shift search range; The P-wave observation waveform extracted from the i-th component of the k-th station; The S-wave observation waveform extracted from the i-th component of the k-th station; The composite theoretical seismic waveform of the i-th component at the k-th station; The residuals are the fitting values ​​for the P-wave window. The residuals are the fitting values ​​for the S-wave window. The total waveform residual; The number of valid stations participating in the inversion; Representing vectors Norm; The time shift of the maximum cross-correlation value of the i-component P-wave at k stations; The time shift of the maximum cross-correlation value of the i-component S-wave at k stations; The theoretical waveform of the i-component P-wave window from k stations after cross-correlation and time-shift alignment; The theoretical waveform of the i-component S-wave window of k stations after cross-correlation and time-shift alignment; For the first The first station The optimal time shift amount obtained by each component during the dynamic alignment process.

7. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 6, characterized in that, The weighting parameter is the ratio of the relative contribution of the total waveform residual to the global polarity similarity term in the joint objective function; the weighting parameter is determined in the following way: Based on the magnitude of the main shock, the actual earthquake sequence is divided into magnitude 1, magnitude 2, and magnitude 3; When the magnitude is the first earthquake, the weighting parameter is the first ratio; When the magnitude is the second earthquake, the weighting parameter is the second ratio; When the magnitude is the third earthquake, the weighting parameter is the third ratio. Wherein, the first magnitude is greater than the second magnitude, the second magnitude is greater than the third magnitude; the first ratio is greater than the second ratio, and the second ratio is greater than the third ratio.

8. The method for joint inversion of source mechanisms based on polarity similarity and full waveform as described in claim 6, characterized in that, In step S40, the trained Siamese neural network is evaluated, including: Calculate the true positives, false positives, true negatives, and false negatives of the trained Siamese neural network on an independent test set; Based on the true positive, false positive, true negative, and false negative results, multiple evaluation indicators are calculated; these multiple evaluation indicators include: precision, recall, and F1 score. When all the evaluation metrics are greater than the preset threshold, the trained Siamese neural network is determined to meet the preset requirements.

9. A method for joint inversion of source mechanisms based on polarity similarity and full waveform according to any one of claims 1 to 8, characterized in that, In step S20, the first preprocessing of the reference seismic dataset includes: The waveforms in the reference seismic dataset are sequentially averaged, filtered, downsampled, and normalized to standard deviation to form noise samples; Seismic waveforms with initial motion polarity labels are extracted from the reference seismic dataset to construct a basic labeled sample. The base labeled samples are filtered and time perturbed. Then, based on the center point after downsampling, a window with a preset number of sampling points is extracted as the input sample. For each input sample, retain both the original and flipped forms to form an expanded labeled sample; High signal-to-noise ratio (SNR) records are selected from the expanded labeled samples as relatively pure seismic signals. Then, a noise window is randomly extracted from the noise sample, and the relatively pure seismic signal and the noise window are superimposed according to a given target SNR to construct a low SNR synthetic sample. High signal-to-noise ratio (SNR) samples, low SNR samples, relatively pure samples, and the low SNR synthesized samples are randomly mixed, and then positive and negative polarity samples are extracted separately; a preset number of records are extracted from each of the positive and negative polarity samples to form a balanced training set; The samples in the balanced training set are randomly shuffled.

10. A method for joint inversion of source mechanisms based on polarity similarity and full waveform according to any one of claims 1 to 8, characterized in that, The actual seismic observation data includes the three-component SAC waveforms of the stations; In step S20, the third preprocessing of the actual seismic observation data includes: The three-component SAC waveform of each valid station is processed by removing the mean, detrending, and applying a tapered window; The processed three-component SAC waveform is resampled to obtain waveform information; The absolute arrival time of the seismic phase is uniformly converted into an arrival time parameter relative to the start time of the current waveform record in order to obtain the seismic phase arrival time; The station's latitude and longitude information is extracted from the head segment of the three-component SAC waveform, and the epicentral distance and azimuth from the event to the station are calculated in combination with the location of the event's epicenter to obtain the spatial location of the station relative to the epicenter. Observational data are constructed using the waveform information, the arrival time of the seismic phase, and the spatial location information. The observation data are comprehensively filtered by combining the completeness of the phase arrival time, the quality of the three-component waveforms, and the availability of station records for each event, to obtain the actual observation dataset.

Images