A method, system and electronic device for extracting seismic magnetic field anomalies

Abstract

The present disclosure belongs to the field of geomagnetic station earthquake anomaly extraction, and is an earthquake magnetic field anomaly extraction method, system and electronic equipment, comprising: constructing a three-dimensional non-negative tensor data body; adopting a spatial weighted non-negative tensor decomposition method to decompose the three-dimensional non-negative tensor data body, extracting R characteristic components, each characteristic component containing a frequency factor matrix, a time factor matrix and a station contribution factor matrix; calculating the proportion of the total energy of a target frequency band in the frequency factor matrix of each characteristic component in the total energy in the entire frequency range, and selecting the characteristic component with the largest proportion as the earthquake-related characteristic component; and based on the time factor matrix of the earthquake-related characteristic component, adopting an over-limit threshold method to extract an earthquake anomaly point. The present disclosure can retain and utilize all measured data to study earthquakes, and effectively detect earthquake anomalies by obtaining more relevant components of earthquake activity.

Description

Technical Field

[0001] This disclosure pertains to the field of seismic anomaly extraction from geomagnetic stations, specifically a method, system, and electronic device for extracting seismic magnetic field anomalies. Background Technology

[0002] Earthquake disasters are characterized by their suddenness, destructive power, and unpredictability, bringing immense calamities to human society. Earthquake prediction is recognized worldwide as one of the scientific challenges. However, numerous studies have revealed that precursory anomalies related to earthquake gestation can appear in the lithosphere, atmosphere, and ionosphere before an earthquake. Monitoring these precursory anomalies and studying their pre-earthquake behavior patterns is of great significance for the research and development of earthquake prediction. Among the many anomalous precursors, electromagnetic precursors are considered one of the research directions most likely to achieve breakthroughs in earthquake prediction. Numerous studies have reported pre-earthquake electromagnetic anomalies in the DC-HF frequency band, including observations from ground stations and satellites.

[0003] Undoubtedly, ground-based stations are closer to the earthquake source than satellites, resulting in less attenuation of pre-earthquake electromagnetic signals as they propagate to the stations, and making it easier to capture some localized, small-scale anomalous signals. However, given the complexity of the Earth's electromagnetic environment, ground-based stations are subject to not only local interference (such as lightning, high-voltage direct current transmission lines, trams, or other anthropogenic interference) but also global interference (including geomagnetic storms, solar activity, etc.). Studies based on a single station cannot accurately reflect pre-earthquake electromagnetic anomaly changes. Some studies have attempted to extract pre-earthquake electromagnetic anomalies by fusing data from multiple stations, such as principal component analysis and multi-channel singular spectrum analysis. Using such methods, common features from data from multiple observation platforms can be extracted, effectively reducing the randomness introduced by a single observation platform and increasing the reliability of the results. Hattori et al. used principal component analysis to process observation data from multiple geomagnetic stations and extracted geomagnetic fusion anomalies before the 2000 Izu Islands earthquake. Yu Zining et al. conducted multi-channel singular spectrum analysis on observation data from multiple borehole strain stations, removing periodic signals and noise from the dataset to obtain seismic signals. They then extracted strain anomalies before multiple earthquakes using a topological network. In fact, geomagnetic anomalies are essentially space-time coupled events; the above method only considers the temporal variation of the signal and lacks utilization of spatial information. Existing research shows that when anomalous signals from the source region propagate to stations at different epicentral distances, the anomalous energy decreases with increasing epicentral distance. Treating all stations at different epicentral distances equally may be detrimental to the extraction of pre-earthquake electromagnetic anomalies.

[0004] Accordingly, considering the spatial propagation characteristics of anomalous signals from the seismic source area, a method for extracting pre-seismic geomagnetic field anomalies based on spatially weighted non-negative tensor decomposition is proposed. This method aims to separate signal components with clear spatial attenuation characteristics originating from the seismic source area from the background field, thereby potentially more directly relating them to the physical processes of the seismic source. Summary of the Invention

[0005] This disclosure provides a method, system, and electronic device for extracting seismic magnetic field anomalies, which solves the problem that traditional data fusion methods lack prior information when processing geomagnetic stations, resulting in poor decomposition effects.

[0006] A method for extracting seismic magnetic field anomalies according to a first aspect of this disclosure includes:

[0007] Preprocessing of the observation data from K stations yields continuous raw data;

[0008] Discrete wavelet transform is performed on the raw data of each station every day to remove the low-frequency background field and obtain low-frequency de-filtered data;

[0009] Multi-channel empirical mode decomposition was performed on the low-frequency data from all stations. High-frequency noise was removed based on the kurtosis and energy entropy of each component after decomposition, and the denoised reconstructed signal was obtained.

[0010] Time-frequency analysis was performed on the denoised reconstructed signal of each station to obtain its two-dimensional time-frequency amplitude spectrum. The spectrum was then stacked along the third dimension according to the epicentral distance of the stations to construct a three-dimensional non-negative tensor data volume.

[0011] The three-dimensional non-negative tensor data volume is decomposed using the spatial weighted non-negative tensor decomposition method, and R feature components are extracted. Each feature component contains a frequency factor matrix, a time factor matrix, and a station contribution factor matrix.

[0012] Calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each characteristic component to the total energy of the entire frequency range, and select the characteristic component with the largest proportion as the earthquake-related characteristic component;

[0013] Based on the time factor matrix of the earthquake-related feature components, earthquake anomalies are extracted using the over-threshold method.

[0014] In some embodiments of the first aspect of this disclosure, the three-dimensional non-negative tensor data volume is decomposed using a spatially weighted non-negative tensor decomposition method to extract R feature components, including:

[0015] Obtain geomagnetic station parameters and construct a geomagnetic station epicentral distance weight vector based on the distance between each station and the potential earthquake source;

[0016] Based on the epicentral distance weight vector, an objective function for spatial weighted nonnegative tensor decomposition is constructed; wherein, the objective function includes the sum of a data fitting term and a spatial regularization term, the data fitting term is used to calculate the error between the three-dimensional nonnegative tensor data volume and the reconstructed tensor after decomposition, and the spatial regularization term is used to constrain the difference between the station contribution factor matrix to be solved and the epicentral distance weight vector through the regularization coefficient.

[0017] The objective function is solved iteratively. In each iteration, the frequency factor matrix, time factor matrix, and station contribution factor matrix are updated alternately. When updating the station contribution factor matrix, the epicentral distance weight vector is introduced for regularization constraint.

[0018] Determine whether the current iteration meets the preset convergence condition or reaches the maximum number of iterations; if so, stop the iteration and output the frequency factor matrix, time factor matrix, and station contribution factor matrix after spatial regularization constraint.

[0019] In some embodiments of the first aspect of this disclosure, the construction of the epicentral distance weight vector of the geomagnetic station specifically includes: calculating the spatial distance of each geomagnetic station relative to the potential earthquake source based on the spatial location of the station's latitude and longitude and the potential earthquake source region, and constructing a weight function that varies with distance based on the spatial distance to generate the corresponding epicentral distance weight vector.

[0020] In some embodiments of the first aspect of this disclosure, the following formula is used as the update rule for the station contribution factor matrix:

[0021] ,

[0022] in, This represents the pattern matrix representation of the third tensor. For time factor matrix, For frequency factor matrix, This represents element-wise multiplication. Contribution factor matrix for stations;

[0023] The convergence condition is:

[0024] ,

[0025] in, It is the convergence tolerance. Describe the objective function For factor matrix The gradient.

[0026] In some embodiments of the first aspect of this disclosure, the proportion of the total energy of the target frequency band in the frequency factor matrix of each feature component to the total energy of the entire frequency range is calculated, and the feature component with the largest proportion is selected as the earthquake-related feature component, including:

[0027] Construct a frequency proportion coefficient, which is used to characterize the ratio of the total energy in the target frequency band to the total energy in the entire frequency range in the frequency factor matrix;

[0028] Based on the numerical value of the frequency proportion coefficient, the feature vector groups representing the feature components in the decomposition result are sorted in descending order to generate a sorted feature vector sequence.

[0029] The feature component ranked first is identified as the target component most relevant to seismic activity and is used as the seismic-related feature component.

[0030] In some embodiments of the first aspect of this disclosure, seismic anomalies are extracted using an out-of-limit threshold method based on the time factor matrix of the seismic correlation feature components, including:

[0031] Calculate the abnormal amplitude threshold Statistical analysis is performed on the time eigenvectors in the time factor matrix, and their mean is calculated. and standard deviation And according to the preset coefficients Constructing abnormal amplitude thresholds ;

[0032] Calculate the duration of the anomaly: count the number of points in the time feature vector that consecutively meet the anomaly conditions. And set a lower threshold for the abnormal duration window. and upper limit threshold ;

[0033] Jointly identify outliers: when the amplitude of a data point in the time feature vector exceeds the outlier amplitude threshold. And the number of consecutive outliers to which the data point belongs satisfy If the time is right, it is determined to be an earthquake anomaly point.

[0034] In some embodiments of the first aspect of this disclosure, the target frequency band ranges from 0.01 Hz to 0.05 Hz.

[0035] A seismic magnetic field anomaly extraction system according to a second aspect embodiment of this disclosure includes:

[0036] The preprocessing module is used to preprocess the observation data from K stations to obtain continuous raw data;

[0037] The low-frequency removal module is used to perform discrete wavelet transform on the raw data of each station every day to remove the low-frequency background field and obtain low-frequency removed data.

[0038] The denoising module is used to perform multi-channel empirical mode decomposition on the low-frequency data of all stations, remove high-frequency noise based on the kurtosis and energy entropy of each component after decomposition, and reconstruct the denoised reconstructed signal.

[0039] The tensor construction module is used to perform time-frequency analysis on the denoised reconstructed signal of each station to obtain its two-dimensional time-frequency amplitude spectrum, and stack it along the third dimension according to the epicentral distance of the station to construct a three-dimensional non-negative tensor data volume.

[0040] The decomposition module is used to decompose the three-dimensional non-negative tensor data volume using the spatial weighted non-negative tensor decomposition method, and extract R feature components. Each feature component contains a frequency factor matrix, a time factor matrix, and a station contribution factor matrix.

[0041] The feature selection module is used to calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each feature component to the total energy of the entire frequency range, and select the feature component with the largest proportion as the earthquake-related feature component.

[0042] The anomaly extraction module is used to extract earthquake anomalies based on the time factor matrix of the earthquake-related feature components using an out-of-limit threshold method.

[0043] An electronic device according to a third aspect of this application includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements a method for extracting seismic magnetic field anomalies.

[0044] Compared with the prior art, the beneficial effects of this disclosure are as follows: by introducing the constraint of the spatial location of the station into the optimization process, this disclosure is good at separating the signal components with clear spatial attenuation characteristics originating from the source area from the background field, thereby potentially more directly related to the physical process of the source and improving the robustness of geomagnetic anomaly extraction.

[0045] This disclosure allows for the retention and utilization of all measured data for earthquake research, while simultaneously obtaining components more relevant to seismic activity for effective seismic anomaly detection. Attached Figure Description

[0046] Figure 1 A flowchart of a method provided for an embodiment of this disclosure;

[0047] Figure 2 A schematic diagram of the raw data from the seven stations provided in this embodiment of the disclosure;

[0048] Figure 3This is a schematic diagram of the original data after removing low frequencies, provided in an embodiment of the present disclosure;

[0049] Figure 4 A schematic diagram of the reconstructed data after removing high-frequency noise, provided in an embodiment of this disclosure;

[0050] Figure 5 A schematic diagram of a third-order nonnegative tensor data volume constructed from a two-dimensional time-frequency amplitude spectrum provided in an embodiment of this disclosure;

[0051] Figure 6 A schematic diagram of the frequency factor matrix, time factor matrix, and station contribution factor matrix obtained by spatial weighted nonnegative tensor decomposition according to embodiments of this disclosure;

[0052] Figure 7 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this disclosure. Detailed Implementation

[0053] To make the objectives, technical solutions, and advantages of this disclosure clearer, the following detailed description is provided in conjunction with embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this disclosure.

[0054] To address the problem that existing methods for extracting pre-earthquake electromagnetic anomalies only consider the temporal changes of signals and lack the utilization of spatial information, resulting in poor decomposition effects due to a lack of prior information when processing geomagnetic stations, this disclosure proposes a method, system, and electronic device for extracting seismic magnetic field anomalies. This method constructs a third-order non-negative tensor data volume from the two-dimensional time-frequency amplitude spectrum matrices of multiple stations. It then decomposes the third-order non-negative tensor data volume using spatially weighted non-negative decomposition to obtain a frequency factor matrix, a time factor matrix, and a station contribution factor matrix. Based on the frequency energy proportion in the frequency factor matrix, seismic-related components are selected, and seismic anomalies are extracted from the time factor matrix of these components using an out-of-limit threshold method. Daily seismic magnetic field anomalies are accumulated, and their deviation from the background fitted line is used to detect seismic anomalies. This disclosure can retain and utilize all measured data for earthquake research, while simultaneously obtaining components more relevant to seismic activity for effective seismic anomaly detection.

[0055] See Figure 1 The flowchart shown illustrates a method for extracting seismic magnetic field anomalies. This embodiment of the disclosure provides a method for extracting seismic magnetic field anomalies, comprising:

[0056] S101, preprocess the observation data from K stations to obtain continuous raw data;

[0057] S102, Perform discrete wavelet transform on the raw data of each station every day to remove the low-frequency background field and obtain low-frequency removed data;

[0058] S103. Perform multi-channel empirical mode decomposition on the low-frequency data of all stations, remove high-frequency noise based on the kurtosis and energy entropy of each component after decomposition, and reconstruct the denoised reconstructed signal.

[0059] S104. Perform time-frequency analysis on the denoised reconstructed signal of each station to obtain its two-dimensional time-frequency amplitude spectrum, and stack it along the third dimension according to the epicentral distance of the station to construct a three-dimensional non-negative tensor data volume.

[0060] S105, the three-dimensional non-negative tensor data volume is decomposed using the spatial weighted non-negative tensor decomposition method, and R feature components are extracted. Each feature component contains a frequency factor matrix, a time factor matrix, and a station contribution factor matrix.

[0061] S106, calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each feature component to the total energy of the entire frequency range, and select the feature component with the largest proportion as the earthquake-related feature component.

[0062] S107, Based on the time factor matrix of the earthquake-related feature components, the earthquake anomaly points are extracted using the over-threshold method.

[0063] In step S101, for the observation data of K stations, the data of the Z component is selected for processing. In order to reduce the influence of artificial noise, the data at midnight (LT00:00-04:00) is selected. The data marked 99999 in the data is removed, and then interpolation is performed to obtain the preprocessed original data.

[0064] In one example, taking the 1Hz vertical component data from seven stations near the epicenter of the Madoi earthquake on May 21, 2021 (UTC time) as an example, the data from midnight (local time: 00:00-04:00) is selected for processing. After removing invalid data from the original data, the missing data is interpolated using the cubic spline interpolation method. The cubic spline interpolation formula is as follows:

[0065] In the interval [ ]superior:

[0066] ,

[0067] in, The y-coordinate is the known data. This represents the curve starting from point ( , ) Tangent direction at the start It is half the curvature. Reflects from point arrive The average rate of change of curvature c This represents the interpolation result.

[0068] After data interpolation, the original waveforms of the seven stations each day are obtained. See [link / reference]. Figure 2 As shown in the figure, the original waveform data of stations 1-7 after preprocessing are shown in the figure.

[0069] To remove the background trend of the magnetic field, the original data was subjected to Discrete Wavelet Transform (DWT) to remove low-frequency background. See [link to relevant documentation]. Figure 3 As shown. From Figure 3 As can be seen, after the original waveforms of each station are processed by discrete wavelet transform, the low-frequency trend signals in the original waveforms are removed.

[0070] Discrete wavelet transform (DWT) continuously filters and downsamples the original signal, decomposing it into low-frequency and high-frequency components. By repeatedly decomposing the low-frequency components, different frequency components of the signal can be obtained. The formula for DWT is shown below:

[0071] ,

[0072] ,

[0073] In the formula, This is magnetic field data for a single station. It is a low-pass filter. For high-pass filters, It is the first The approximation coefficients of the layer decomposition represent the low-frequency information of that layer. It is the first The detail coefficients of the layer decomposition represent the high-frequency information of that layer. , indicating the first The low-frequency subband coefficient of the layer, This represents the discrete-time index of the coefficients in the current layer. This is the index of the filter coefficients, ranging from 0 to K, where K is the length of the filter.

[0074] The high-pass and low-pass filters are determined by the wavelet basis used.

[0075] In one example, the db4 wavelet basis is used, specifically represented as:

[0076] Low-pass filter (scaling function coefficients): , , , ,

[0077] The high-pass filter (wavelet function coefficients) are obtained from the low-pass coefficients through the orthogonal mirror filter relationship: ,

[0078] Based on the shape and frequency characteristics of the original data, the data is decomposed into 6 levels using discrete wavelet transform, and the low-frequency background is removed by removing the approximation coefficients of the last level.

[0079] In step 103, since high-frequency noise still exists at the stations, multi-channel empirical mode decomposition (MEMD) is used to decompose the low-frequency data from all stations to remove it, thus obtaining the intrinsic mode functions (IMF components) of each station. Multi-channel empirical mode decomposition can be expressed as:

[0080] ,

[0081] This represents the original multivariable signal. It is the first multivariable eigenmode functions of order, This is the first In the eigenmode function of order 1, the 1st order 2 is... The components of each channel. It is a multivariate residual, and the third variable in the residual signal is... The components of each channel.

[0082] The reconstructed signal is obtained after removing high-frequency components. In a specific example, seven stations are processed together to obtain the IMF components of each station. Then, the kurtosis and energy entropy of each IMF component of each station are calculated. The formula for calculating the energy entropy is as follows:

[0083] ,

[0084] ,

[0085] in, The energy for each IMF component, where b is a constant, typically taken as 2. This represents the entropy value. No. The proportion of the energy of each IMF component to the total energy of the entire signal.

[0086] The formula for calculating kurtosis is as follows: ,

[0087] N is the total number of data points. Overall mean Meaning indicates the first The values ​​of each data sample point This represents the overall mean of the data.

[0088] In one example, the IMF components with kurtosis > 3 and energy entropy < 0.3 were reconstructed to obtain the denoised reconstructed signal. The signals from each station after removing high-frequency noise are shown below. Figure 4 As shown.

[0089] The denoised and reconstructed signal is then subjected to synchronous wavelet-squeezed transform (WSST) to obtain the two-dimensional time-frequency amplitude spectrum of each station. Subsequently, The two-dimensional time-frequency amplitude spectra are stacked according to the epicentral distance of each station along the direction perpendicular to their third dimension to obtain a third-order non-negative tensor, which is called a three-dimensional non-negative tensor data volume, such as... Figure 5 (As shown). The three dimensions of this three-dimensional non-negative tensor data volume are time, frequency, and station observation channel, respectively. This construction not only simultaneously characterizes the spatial and frequency characteristics of the station magnetic field data, but also reflects the structural information of the two-dimensional characteristics of different observation data, that is, the three-dimensional characteristics of the data.

[0090] Synchronous wavelet compression transform has high time-frequency resolution, which can more clearly display the time-frequency distribution characteristics of magnetic field data, and is beneficial for the fusion of anomalous information from different observation data. This method first obtains the instantaneous frequency based on the coefficients of the wavelet transform, and then compresses the values ​​around any frequency to the current frequency to improve the time-frequency resolution of the result.

[0091] frequency nearby Values ​​compressed to The formula above is:

[0092] ;

[0093] ;

[0094] ;

[0095] In the formula, For magnetic field data Wavelet coefficients obtained after discrete wavelet transform As a scale factor, The translation factor is... For the mother wavelet used in discrete wavelet transform, the embodiments of this disclosure select the Morlet wavelet. The instantaneous frequency is obtained through wavelet coefficients. This represents the current compression scale factor.

[0096] A tensor is a multidimensional array. First-order tensors are vectors, second-order tensors are matrices, and third-order and higher-order tensors are higher-order tensors. The tensors described in this embodiment are higher-order tensors. The order and dimension of a tensor are the same, also called its modulus. order tensor , No. The size of the module is Tensors Modulus vector, i.e. Fiber bundles refer to its Multiple vectors in the modulus direction. A subarray of a tensor is formed by fixing some of its index subscripts while changing the rest, where a tensor slice is a matrix formed when only two subscripts change.

[0097] The rank of a tensor is the minimum number of rank-tensors required to represent a tensor as the sum of rank-tensors. (Any) order tensor The rank of is denoted as .matrix sum matrix The Kronecker product is defined as:

[0098]

[0099] The size of the new matrix is , Belongs to matrix In OK The elements of the column.

[0100] matrix sum matrix The Khatri-Rao product is defined as: ,

[0101] The size of the new matrix is And there are .

[0102] A rank-tensor is a tensor that can be represented by the outer product of one set of vectors. Rank Tensor It can be by vectors The outer product representation of , and this representation is unique:

[0103] , Indicates the first Eigenvectors at order (pattern).

[0104] Among them, symbols The outer product of vectors is represented by the following element-wise expression: , . Indicates the first vectors The first in Each component.

[0105] In step S104, the third-order nonnegative tensors for each day are decomposed using spatially weighted nonnegative tensor decomposition. This includes:

[0106] Obtain geomagnetic station parameters and construct a geomagnetic station epicentral distance weight vector based on the distance between each station and the potential earthquake source;

[0107] Based on the epicentral distance weight vector, an objective function for spatial weighted nonnegative tensor decomposition is constructed; wherein, the objective function includes the sum of a data fitting term and a spatial regularization term, the data fitting term is used to calculate the error between the three-dimensional nonnegative tensor data volume and the reconstructed tensor after decomposition, and the spatial regularization term is used to constrain the difference between the station contribution factor matrix to be solved and the epicentral distance weight vector through the regularization coefficient.

[0108] The objective function is solved iteratively. In each iteration, the frequency factor matrix, time factor matrix, and station contribution factor matrix are updated alternately. When updating the station contribution factor matrix, the epicentral distance weight vector is introduced for regularization constraint.

[0109] Determine whether the current iteration meets the preset convergence condition or reaches the maximum number of iterations; if so, stop the iteration and output the frequency factor matrix, time factor matrix, and station contribution factor matrix after spatial regularization constraint.

[0110] The objective function comprises two terms, with the data fitting term belonging to Non-Negative Tensor Factorization (NTF). This decomposition allows for the reflection of the main local features of the original data using a small amount of data. Using NTF to decompose the constructed third-order time-spectrum tensor accurately separates the seismic signals within the dataset. Furthermore, the introduction of coefficient components reflects the contribution of observation data from different stations to the anomalies.

[0111] Nonnegative tensor decomposition represents the original tensor as the sum of a finite number of rank-one tensors, i.e., the sum of the outer products of a finite number of vectors. For tensors of order three and above, this representation is usually unique. This depends on the number of eigencomponents in the decomposition. When the rank is equal to that of the original tensor, the decomposition is an exact decomposition, called "rank decomposition". However, the rank of a tensor cannot actually be calculated, so an approximate decomposition is usually performed. After decomposition, different multidimensional local features of the original higher-order tensor can be obtained.

[0112] A third-order tensor was constructed from geomagnetic amplitude spectrum data recorded by multiple stations. Nonnegative tensor decomposition represents the original tensor as the sum of a finite number of tensors of rank 1, expressed as:

[0113] ,

[0114] In the formula, It is denoted as rank, representing the number of eigencomponents in the decomposition. It is the factor matrix of the first pattern, i.e., the frequency factor matrix. This is the factor matrix of the second mode, namely the time factor matrix. This is the factor matrix for the third mode, namely the station contribution factor matrix. Represents the outer product of vectors. , , , which represents the size of each dimension of the tensor. , , For each pattern, Represents the rank of the decomposition.

[0115] Nonnegative tensor decomposition first randomly initializes the nonnegative matrix, calculates the difference between the original tensor and the decomposition result, and then obtains the decomposition matrix by iteratively optimizing to minimize the objective function. The objective function using only nonnegative tensor decomposition can be expressed as: ,

[0116] in , , ,

[0117] The iterative optimization update rule for nonnegative tensor decomposition follows alternating least squares, is derived through gradient descent, and guarantees nonnegativity. Specifically, it is expressed as:

[0118] ;

[0119] ;

[0120] ;

[0121] in The first character represents element-wise multiplication, and the second character represents element-wise division. The representation of the various mode matrixings of a tensor is shown. To distinguish the various mode matrixings of a tensor at different steps, the following will be used: Represented as the pattern matrix of the first tensor, This is represented as a pattern matrix of the second tensor. This is represented as a pattern matrix of the third tensor.

[0122] After each iteration, the Karush-Kuhn-Tucker criterion (KKT criterion) is used to determine whether the result has converged. When the result satisfies the KKT conditions or the number of iterations reaches the upper limit, the iteration stops and the final decomposition matrix is ​​obtained.

[0123] The convergence conditions are as follows:

[0124] ;

[0125] ;

[0126] ;

[0127] In the formula, It is the convergence tolerance. Describe the objective function For the factor matrix respectively The gradient. This is determined when the result satisfies the KKT conditions or the number of iterations reaches its upper limit. When the iteration stops, the result is obtained. , and Three decomposition matrices.

[0128] However, when decomposing the data fitting term, it relies solely on the statistical properties of the data, treating all station data equally. This means that interference signals such as global geomagnetic storms and trans-ionospheric disturbances (TIDs), due to their spatial breadth and consistency, may exhibit similar behavior across all stations, easily being decomposed into the same or a few components and confused with seismic anomalies. Furthermore, if some stations are severely contaminated by noise, their outliers may also significantly affect the global loss function.

[0129] The spatial regularization term introduces spatial weighting into geomagnetic data. It uses prior knowledge of physical space (i.e., the spatial relationship between stations and potential seismic source areas) as constraints, injecting it into the decomposition process that relies solely on statistical data characteristics. This makes the decomposed components not only discriminative in the time-frequency domain but also more physically meaningful in the spatial domain, thus more accurately separating out the real anomalous signals related to earthquake gestation.

[0130] The method, which includes decomposing data fitting terms and spatially regularized objective functions, is named Spatially Weighted Non-Negative Tensor Factorization (SW-NTF) in this embodiment. SW-NTF focuses on stations near the epicenter, guiding the decomposition process towards a physically more reasonable solution space. Specifically, it constructs a weighted vector of geomagnetic station epicentral distances: ,in Indicates the first The distance between each monitoring station and the epicenter is weighted, with smaller distances carrying greater weight. The objective function is then:

[0131] ,

[0132] In the formula, These are regularization coefficients. It's important to note that the spatial regularization term only applies to the station contribution factor matrix. .

[0133] The update rules for frequency factor matrix A and time factor matrix B are the same as above, due to the station contribution factor matrix... Since there is a spatial regularization term, the denominator of the multiplicative update rule needs to be increased by the derivative (positive value) of the spatial regularization term, and the numerator needs to be increased by the negative value of the derivative of the spatial regularization term. Final station contribution factor matrix. The update rules are as follows:

[0134] ,

[0135] Finally, the convergence condition under the KKT criterion is replaced with:

[0136] ,

[0137] Finally, we obtain three non-negative matrices after decomposition. (Frequency factor matrix) (Time Factor Matrix) (Station contribution factor matrix). Indicates the number of characteristic components in the decomposition. In Each eigenvector represents the frequency characteristics of each feature component of the geomagnetic data. In Each eigenvector represents the change of the feature component over time. In Each eigenvector represents the contribution of each station in the eigencomponents.

[0138] In one example, the spatially weighted nonnegative tensor decomposition method is used to decompose the tensors constructed from the seven stations each day. The parameters are set as follows: number of feature components R=3, regularization coefficient... =0.1, iteration count of 200, and convergence tolerance under the KKT criterion. The value is 0.00001, resulting in three non-negative matrices after decomposition. (Frequency factor matrix) (Time Factor Matrix) (Station contribution factor matrix). The three eigenvectors in the data represent the frequency characteristics of each feature component of the geomagnetic data. The three eigenvectors in the diagram represent the changes of each feature component over time. The three eigenvectors in the diagram represent the contribution of each station to these eigencomponents.

[0139] Calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each characteristic component to the total energy across the entire frequency range. Select the characteristic component with the largest proportion as the seismic-related characteristic component, including:

[0140] Construct a frequency proportion coefficient, which is used to characterize the energy proportion of each eigenvector in the frequency factor matrix within the target frequency band;

[0141] Based on the numerical value of the frequency proportion coefficient, the feature vector groups representing the feature components in the decomposition result are sorted in descending order to generate a sorted feature vector sequence.

[0142] The feature component ranked first is identified as the target component most relevant to seismic activity and is used as the seismic-related feature component.

[0143] In one example, based on the frequency proportion coefficient Sort each feature component in descending order to obtain , and and will As an earthquake-related component These represent the time vector and station contribution vector for the seismic components. The final decomposition result is as follows: Figure 6 As shown.

[0144] Spatial weighting strategies naturally suppress the contribution of noise signals to key components. This is because spatially weighted nonnegative tensor decomposition seeks patterns that are prominent in high-weight regions and weak in low-weight regions, which is precisely the characteristic of local seismic anomalies. Figure 6 It can be seen that the earthquake-related component is present in the frequency vector. The energy is mainly concentrated at 0.01 Hz, indicating that anomalies possibly related to earthquakes have been extracted, and the station contribution vector... In the middle, the contributions of several stations close to the epicenter were similar. Additionally, the earthquake-related components... The energy is mainly concentrated at 0.05Hz and 0.15Hz, while the frequency component of the third characteristic component is mainly concentrated at 0.45Hz. This demonstrates better signal separation performance.

[0145] In one embodiment, in order to extract anomalies more accurately, the criteria for judging seismic electromagnetic anomalies usually follow two principles: anomaly amplitude and duration.

[0146] Based on the time factor matrix of the earthquake-related feature components, earthquake anomalies are extracted using an out-of-limit threshold method, including:

[0147] (1) Calculate the abnormal amplitude threshold Statistical analysis is performed on the time eigenvectors in the time factor matrix, and their mean is calculated. and standard deviation And according to the preset coefficients Constructing abnormal amplitude thresholds Standard deviation The calculation formula is as follows:

[0148] ,

[0149] in, It is the total number of data points. It is each data point, It is the overall mean.

[0150] Calculate the duration of the anomaly: count the number of points in the time feature vector that consecutively meet the anomaly conditions. And set a lower threshold for the abnormal duration window. and upper limit threshold ;

[0151] Jointly identify outliers: when the amplitude of a data point in the time feature vector exceeds the outlier amplitude threshold. And the number of consecutive outliers to which the data point belongs satisfy If the time is right, it is determined to be an earthquake anomaly point.

[0152] This disclosure also provides a seismic magnetic field anomaly extraction system for implementing the above-described method. The steps in the above method can be used to explain the seismic magnetic field anomaly extraction system, including:

[0153] The preprocessing module is used to preprocess the observation data from K stations to obtain continuous raw data;

[0154] The low-frequency removal module is used to perform discrete wavelet transform on the raw data of each station every day to remove the low-frequency background field and obtain low-frequency removed data.

[0155] The denoising module is used to perform multi-channel empirical mode decomposition on the low-frequency data of all stations, remove high-frequency noise based on the kurtosis and energy entropy of each component after decomposition, and reconstruct the denoised reconstructed signal.

[0156] The tensor construction module is used to perform time-frequency analysis on the denoised reconstructed signal of each station to obtain its two-dimensional time-frequency amplitude spectrum, and stack it along the third dimension according to the epicentral distance of the station to construct a three-dimensional non-negative tensor data volume.

[0157] The decomposition module is used to decompose the three-dimensional non-negative tensor data volume using the spatial weighted non-negative tensor decomposition method, and extract R feature components. Each feature component contains a frequency factor matrix, a time factor matrix, and a station contribution factor matrix.

[0158] The feature selection module is used to calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each feature component to the total energy of the entire frequency range, and select the feature component with the largest proportion as the earthquake-related feature component.

[0159] The anomaly extraction module is used to extract earthquake anomalies based on the time factor matrix of the earthquake-related feature components using an out-of-limit threshold method.

[0160] This disclosure also provides an electronic device. Figure 7 An exemplary structural diagram of the electronic device is shown. See also Figure 7 The electronic device includes a processor 701 and a memory 702. The memory 702 stores a computer program, which, when run by the processor 701, causes the processor to execute a method for extracting seismic magnetic field anomalies.

[0161] The processor 701 may be a central processing unit (CPU), a graphics processing unit (GPU), or other form of processing unit with data processing capabilities and / or instruction execution capabilities.

[0162] The memory 702 may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. The volatile memory may include, for example, random access memory (RAM) and / or cache memory. The non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc. One or more computer program instructions may be stored on the computer-readable storage medium, and the processor 701 may execute the program to implement a method for extracting seismic magnetic field anomalies.

[0163] Depending on the specific application, the electronic device may also include any other suitable components. For example, the electronic device may also include a communication component 703.

[0164] In addition to the methods and systems, embodiments of this disclosure also provide a computer program product comprising a computer program that, when run by a processor, causes the processor to perform the steps of a seismic magnetic field anomaly extraction method according to embodiments of this disclosure. The computer program product can be written in any combination of one or more programming languages ​​to perform the operations of embodiments of this disclosure. The programming languages ​​include object-oriented programming languages ​​such as Java and C++, as well as conventional procedural programming languages ​​such as C or similar languages. The program code can be executed entirely on a user's computing device, partially on a user's computing device, as a standalone software package, partially on a user's computing device and partially on a remote computing device, or entirely on a remote computing device or server.

[0165] Furthermore, embodiments of this disclosure also provide a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, causes the processor to perform steps in a method for extracting seismic magnetic field anomalies. The computer-readable storage medium may be any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may, for example, include, but is not limited to, electrical, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatuses, or devices, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: electrical connections having one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0166] The above description is merely a preferred embodiment of this disclosure and is not intended to limit this disclosure. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this disclosure should be included within the scope of protection of this disclosure.

Claims

1. A method for extracting seismic magnetic field anomalies, characterized in that, The method includes: Preprocessing of the observation data from K stations yields continuous raw data; Discrete wavelet transform is performed on the raw data of each station every day to remove the low-frequency background field and obtain low-frequency-removed data; Multi-channel empirical mode decomposition was performed on the low-frequency data from all stations. High-frequency noise was removed based on the kurtosis and energy entropy of each component after decomposition, and the denoised reconstructed signal was obtained. Time-frequency analysis was performed on the denoised reconstructed signal of each station to obtain its two-dimensional time-frequency amplitude spectrum. The signals were then stacked along the third dimension according to the epicentral distance of the stations to construct a three-dimensional non-negative tensor data volume. The three-dimensional non-negative tensor data volume is decomposed using the spatial weighted non-negative tensor decomposition method, and R feature components are extracted. Each feature component contains a frequency factor matrix, a time factor matrix, and a station contribution factor matrix. Calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each feature component to the total energy of the entire frequency range, and select the feature component with the largest proportion as the earthquake-related feature component; Based on the time factor matrix of the earthquake-related feature components, earthquake anomalies are extracted using the over-threshold method.

2. The method for extracting seismic magnetic field anomalies according to claim 1, characterized in that, The three-dimensional nonnegative tensor data volume is decomposed using the spatially weighted nonnegative tensor decomposition method to extract R feature components, including: Obtain geomagnetic station parameters and construct a geomagnetic station epicentral distance weight vector based on the distance between each station and the potential earthquake source; Based on the epicentral distance weight vector, an objective function for spatial weighted nonnegative tensor decomposition is constructed; wherein, the objective function includes the sum of a data fitting term and a spatial regularization term, the data fitting term is used to calculate the error between the three-dimensional nonnegative tensor data volume and the reconstructed tensor after decomposition, and the spatial regularization term is used to constrain the difference between the station contribution factor matrix to be solved and the epicentral distance weight vector through the regularization coefficient. The objective function is solved iteratively. In each iteration, the frequency factor matrix, time factor matrix, and station contribution factor matrix are updated alternately. When updating the station contribution factor matrix, the epicentral distance weight vector is introduced for regularization constraint. Determine whether the current iteration meets the preset convergence condition or reaches the maximum number of iterations; if so, stop the iteration and output the frequency factor matrix, time factor matrix, and station contribution factor matrix after spatial regularization constraint.

3. The method for extracting seismic magnetic field anomalies according to claim 2, characterized in that, The construction of the epicentral distance weight vector for geomagnetic stations specifically includes: calculating the spatial distance between each geomagnetic station and the potential seismic source region based on the station's latitude and longitude and the spatial location of the potential seismic source region, and constructing a weight function that varies with distance based on the spatial distance to generate the corresponding epicentral distance weight vector.

4. The method for extracting seismic magnetic field anomalies according to claim 2, characterized in that, The following formula is used as the update rule for the station contribution factor matrix: ， in, This represents the pattern matrix representation of the third tensor. For time factor matrix, For frequency factor matrix, This represents element-wise multiplication. Contribution factor matrix for stations; The convergence condition is: ， in, It is the convergence tolerance. Describe the objective function For factor matrix The gradient.

5. The method for extracting seismic magnetic field anomalies according to claim 2, characterized in that, Calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each characteristic component to the total energy across the entire frequency range. Select the characteristic component with the largest proportion as the seismic-related characteristic component, including: Construct a frequency proportion coefficient, which is used to characterize the ratio of the total energy in the target frequency band to the total energy in the entire frequency range in the frequency factor matrix; Based on the numerical value of the frequency proportion coefficient, the feature vector groups representing the feature components in the decomposition result are sorted in descending order to generate a sorted feature vector sequence. The feature component ranked first is identified as the target component most relevant to seismic activity and is used as the seismic-related feature component.

6. The method for extracting seismic magnetic field anomalies according to claim 2, characterized in that, Based on the time factor matrix of the earthquake-related feature components, earthquake anomalies are extracted using an out-of-limit threshold method, including: Calculate the abnormal amplitude threshold Statistical analysis is performed on the time eigenvectors in the time factor matrix, and their mean is calculated. and standard deviation And according to the preset coefficients Constructing abnormal amplitude thresholds ; Calculate the duration of the anomaly: count the number of points in the time feature vector that consecutively meet the anomaly conditions. And set a lower threshold for the abnormal duration window. and upper limit threshold ; Jointly identify outliers: when the amplitude of a data point in the time feature vector exceeds the outlier amplitude threshold. And the number of consecutive outliers to which the data point belongs satisfy If the time is right, it is determined to be an earthquake anomaly point.

7. The method for extracting seismic magnetic field anomalies according to claim 5, characterized in that, The target frequency band ranges from 0.01 Hz to 0.05 Hz.

8. A seismic magnetic field anomaly extraction system, characterized in that, include: The preprocessing module is used to preprocess the observation data from K stations to obtain continuous raw data; The low-frequency removal module is used to perform discrete wavelet transform on the raw data of each station every day to remove the low-frequency background field and obtain low-frequency removed data. The denoising module is used to perform multi-channel empirical mode decomposition on the low-frequency data of all stations, remove high-frequency noise based on the kurtosis and energy entropy of each component after decomposition, and reconstruct the denoised reconstructed signal. The tensor construction module is used to perform time-frequency analysis on the denoised reconstructed signal of each station to obtain its two-dimensional time-frequency amplitude spectrum, and stack it along the third dimension according to the epicentral distance of the station to construct a three-dimensional non-negative tensor data volume. The decomposition module is used to decompose the three-dimensional non-negative tensor data volume using the spatial weighted non-negative tensor decomposition method, and extract R feature components. Each feature component contains a frequency factor matrix, a time factor matrix, and a station contribution factor matrix. The feature selection module is used to calculate the proportion of the total energy of the target frequency band in the frequency factor matrix of each feature component to the total energy of the entire frequency range, and select the feature component with the largest proportion as the earthquake-related feature component. The anomaly extraction module is used to extract earthquake anomalies based on the time factor matrix of the earthquake-related feature components using an out-of-limit threshold method.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method as described in any one of claims 1 to 7.

Images