Satellite magnetic field data earthquake anomaly detection method and system

Abstract

The application belongs to the technical field of satellite ionosphere magnetic field earthquake anomaly extraction, and particularly relates to a satellite magnetic field data earthquake anomaly detection method and system, which comprises the following steps: preprocessing satellite magnetic field observation data and converting the same into time-frequency representation to obtain a non-negative time-frequency amplitude matrix; taking the non-negative time-frequency amplitude matrix as a decomposition object, decomposing the non-negative time-frequency amplitude matrix into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template by using a deep expandable convolution non-negative matrix decomposition; extracting a local influence component related to an earthquake based on the decomposition result, performing anomaly orbit determination on the local influence component related to the earthquake; counting the number of anomaly orbits per day and accumulating the number, and determining an earthquake anomaly according to the deviation degree of the accumulated result from a background fitting straight line. The method can suppress global disturbance interference without relying on removing data of strong solar activity or strong geomagnetic activity dates, and improve data utilization rate and detection stability.

Description

Technical Field

[0001] This application belongs to the field of space ionospheric satellite magnetic field seismic anomaly extraction technology, specifically a method and system for detecting seismic anomalies in satellite magnetic field data. Background Technology

[0002] Earthquake anomalies typically refer to unusual phenomena occurring before and after an earthquake in its epicenter and surrounding area. Typical manifestations include abnormal changes in geophysical quantities such as geomagnetism, geoelectricity, and gravity. Existing research usually involves acquiring observational data of these physical quantities and analyzing their anomalous characteristics to identify earthquake anomalies and conduct correlation studies.

[0003] Existing observation methods mainly include ground-based station observations and satellite telemetry. Satellite monitoring has advantages such as wide coverage and the ability to acquire large-scale continuous observations, and existing research has shown that satellite observation can capture anomalies related to seismic activity, thus attracting widespread attention. Currently, research on seismic anomalies based on satellite magnetic field data is relatively concentrated. A common processing procedure involves selecting nighttime orbital data to reduce the influence of solar radiation, and further removing data from periods with high solar and geomagnetic activity indices to reduce interference from solar activity and geomagnetic disturbances on the magnetic field sequence; then, anomaly detection is directly performed on the remaining data after screening. However, this traditional processing method often relies on removing data from periods of "strong solar / strong geomagnetic activity," which can easily reduce the amount of usable data. Furthermore, data gaps or sample biases may occur during periods of frequent solar or geomagnetic activity, affecting the continuity and stability of anomaly detection. In addition, existing methods often perform anomaly detection without effectively separating mixed influence components, allowing local earthquake-related disturbances to be masked by global disturbances or background changes, reducing the specificity of detection. These methods suffer from problems such as the need for data removal, difficulty in highlighting local earthquake components, and insufficient detection stability. Summary of the Invention

[0004] This application provides a method and system for detecting seismic anomalies in satellite magnetic field data, which solves the problems of difficulty in highlighting local components of earthquakes and insufficient detection stability.

[0005] The first aspect of this application provides a method for detecting seismic anomalies in satellite magnetic field data, including: The satellite magnetic field observation data is preprocessed and converted into a time-frequency representation to obtain a non-negative time-frequency amplitude matrix; The non-negative time-frequency amplitude matrix is ​​used as the decomposition object. Using depthwise unfoldable convolutional non-negative matrix decomposition, the non-negative time-frequency amplitude matrix is ​​decomposed into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template. The two-dimensional time-frequency convolution template is used to characterize several basis features of the non-negative time-frequency amplitude matrix, and the non-negative weight matrix is ​​used to reflect the amplitude of each basis feature as time changes. Based on the decomposition results, earthquake-related local impact components are extracted, and abnormal trajectories are determined for these components. The number of abnormal orbits is counted daily and accumulated. The degree of deviation of the accumulated results from the fitted straight line relative to the background is used to determine the earthquake anomalies.

[0006] Furthermore, the satellite magnetic field observation data is preprocessed and represented in a time-frequency format to obtain a non-negative time-frequency amplitude matrix, including: Subtracting the CHAOS-6 magnetic field model data from the satellite magnetic field observation data yields the residual data, and the first difference of the residual data is then calculated to obtain the difference data. A short-time Fourier transform is performed on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the satellite magnetic field observation data.

[0007] Furthermore, a short-time Fourier transform is performed on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the satellite magnetic field observation data, including: The time-frequency transformation of the differential data is performed by short-time Fourier transform to obtain a time-frequency matrix that includes amplitude and phase information; The amplitude information of the time-frequency matrix is ​​extracted to construct a non-negative time-frequency amplitude matrix.

[0008] Furthermore, earthquake-related local impact components are extracted based on the decomposition results, and abnormal orbit determination is performed on the earthquake-related local impact components. This includes using the energy ratio criterion to screen the earthquake-related local impact components from the decomposition results, and determining abnormal orbits of the earthquake-related local impact components based on the over-limit rate.

[0009] Furthermore, using depthwise expandable convolutional nonnegative matrix factorization, the nonnegative time-frequency amplitude matrix is ​​decomposed into a combination of a nonnegative weight matrix and a two-dimensional time-frequency convolutional template, including: Construct an objective function that minimizes the reconstruction error and the weighted sum of the L1 norm of the multilayer activation coefficients based on the β-divergence. Based on the objective function, the iterative update process is unfolded into a deep network containing multiple cascaded network layers. Each layer corresponds to an update operation along the negative gradient direction of the objective function, transforming the iterative process into an end-to-end trainable network structure. The forward propagation of each layer performs reconstruction calculations, updates the non-negative weight matrix, and updates the two-dimensional time-frequency convolution template. After all network layers, the final two-dimensional time-frequency convolution template and non-negative weight matrix are output.

[0010] Furthermore, the forward propagation of each layer performs reconstruction computation, updates the non-negative weight matrix, and updates the 2D time-frequency convolution template, including: The time-frequency matrix of the current layer is reconstructed by performing a one-dimensional convolution with the two-dimensional time-frequency convolution template of the previous layer and the non-negative weight matrix. Based on the gradient of the reconstruction error with respect to the non-negative weight matrix, combined with the learnable step size parameter and the non-negative activation function, the non-negative weight matrix is ​​updated. Based on the gradient of the reconstruction error with respect to the 2D time-frequency convolution template, and combined with the learnable stride parameter and non-negative activation function, the 2D time-frequency convolution template is updated.

[0011] Furthermore, the energy ratio criterion is used to screen the earthquake-related local influence components from the decomposition results, and the abnormal trajectories of the earthquake-related local influence components are determined based on the exceedance rate, including: By calculating the ratio of energy in the study area to the energy of the entire orbit in each row vector of the non-negative weight matrix obtained by non-negative time-frequency amplitude matrix decomposition, the component with the largest energy ratio is selected as the earthquake-related local influence component. The root mean square of the average energy level of the reaction data is used as the threshold parameter. Anomalies are judged based on whether the value of the largest component is greater than k times the root mean square of the entire orbit. Anomalies are defined as orbits that only appear in the study area and whose corresponding flags indicate that the anomalies are not caused by known ionospheric activity.

[0012] Furthermore, the number of anomalous orbits is counted daily and accumulated. Seismic anomalies are determined based on the degree of deviation of the accumulated results from the background fitted straight line, including: The daily number of abnormal orbits and the daily total number of orbits were standardized by deviation to obtain the cumulative normalized result of the daily number of abnormal orbits and the cumulative normalized result of the daily total number of orbits. The cumulative normalized results of the total number of orbits for each day are subjected to least squares linear fitting to obtain the background fitting line; The Sigmoid function is used to fit the cumulative normalized result of the number of daily abnormal orbits and the deviation of the fitted line from the background.

[0013] A second aspect of this application provides a satellite magnetic field data seismic anomaly detection system for performing a satellite magnetic field data seismic anomaly detection method, including: The data preprocessing and time-frequency conversion module is used to preprocess satellite magnetic field observation data and convert it into a time-frequency representation to obtain a non-negative time-frequency amplitude matrix. The feature decomposition module is used to decompose the non-negative time-frequency amplitude matrix as the decomposition object, and to decompose the non-negative time-frequency amplitude matrix into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template using depthwise unfoldable convolution non-negative matrix decomposition; wherein, the two-dimensional time-frequency convolution template is used to characterize several basis features of the non-negative time-frequency amplitude matrix, and the non-negative weight matrix is ​​used to reflect the amplitude of each basis feature changing over time. The abnormal trajectory determination module is used to extract earthquake-related local influence components based on the decomposition results and determine the abnormal trajectories of the earthquake-related local influence components. The anomaly identification module is used to count the number of abnormal orbits on a daily basis and accumulate them. Based on the degree of deviation of the accumulated results from the background fitted straight line, the seismic anomaly is determined.

[0014] Compared with the prior art, the advantages of this application are as follows: This application decomposes the non-negative time-frequency amplitude matrix of data using deep, expandable convolutional non-negative matrix factorization. The decomposition process maintains non-negativity constraints, and the convolutional operation explicitly represents the duration, waveform shape, and local propagation characteristics of seismic anomalies along the latitudinal (time) axis. In the scenario of anomaly signal identification in seismogenic zones, it more naturally characterizes the event template and the pattern structure of the occurrence time and location, while the decomposition results have good physical interpretability. Its network structure enables end-to-end training, retaining the algorithmic structure and physical interpretability of traditional matrix factorization while expanding the data-driven parameter tuning capabilities of deep learning. Simultaneously, it can suppress global disturbance interference without relying on removing data from dates with strong solar or geomagnetic activity, improving data utilization and detection stability. It can retain and utilize all measured data for earthquake research, effectively using components more relevant to seismic activity for seismic anomaly detection. This overcomes the problems of reduced data volume, insufficient continuity, and easy masking of seismic-related components caused by existing technologies that filter nighttime data, remove high solar / geomagnetic activity data, and directly perform anomaly detection on the remaining data. Attached Figure Description

[0015] Figure 1 A flowchart of a satellite magnetic field data seismic anomaly detection method provided in this application embodiment; Figure 2 This is a schematic diagram of an earthquake research area provided in an embodiment of this application; Figure 3 The original data of the Y component of the satellite magnetic field provided in this application, the residual data curve after subtracting the CHAOS-6 model, and the curve after first-order difference are shown in (a); the solid line in (a) is the original data curve, (b) is the residual data curve after subtracting the CHAOS-6 model from the original data, and (c) is the first-order difference data curve of the residual data. Figure 4 A graph showing the short-time Fourier transform time-frequency amplitude matrix results of satellite Y-component magnetic field differential data provided in this application embodiment; Figure 5The diagram shows the two-dimensional time-frequency convolution template and non-negative weight matrix obtained by depth-expandable convolution non-negative matrix decomposition provided in the embodiments of this application. (a) Non-negative weight matrix, (b) Two-dimensional time-frequency convolution template corresponding to the first component hs1, (c) Two-dimensional time-frequency convolution template corresponding to the second component hs2, and (d) Two-dimensional time-frequency convolution template corresponding to the third component hs3. Figure 6 The temporal reconstruction results of the sorted depth-expandable convolutional nonnegative matrix factorization provided in the embodiments of this application are shown in the following diagrams: (a) is the temporal component corresponding to the first component hs1, (b) is the temporal component corresponding to the second component hs2, and (c) is the temporal component corresponding to the third component hs3. Figure 7 The graph shows the result of fitting a straight line between the cumulative daily abnormal trajectory and the background, as provided in the embodiments of this application. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of this application 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 application.

[0017] This application aims to provide a method and system for detecting seismic anomalies using satellite magnetic field data, so as to effectively separate the local seismic-related effects from the global solar / geomagnetic effects and background components while preserving as much observation data as possible, thereby improving the pertinence and reliability of seismic anomaly detection.

[0018] To achieve the above objectives, the core idea of ​​this application is as follows: preprocess and time-frequency represent the satellite magnetic field observation data, take the obtained non-negative time-frequency amplitude matrix as the decomposition object, use depth-expandable convolutional non-negative matrix decomposition to achieve multi-component decoupling, and on this basis extract earthquake-related local influence components to carry out abnormal orbit determination and daily-scale statistical detection.

[0019] The depthwise expandable convolutional nonnegative matrix factorization is used to decompose an arbitrarily high-dimensional nonnegative time-frequency amplitude matrix into a combination of a nonnegative weight matrix and a two-dimensional time-frequency convolutional template, making the convolutional reconstruction matrix approximate the original nonnegative time-frequency amplitude matrix. The two-dimensional time-frequency convolutional template is used to characterize several basis features of the original nonnegative time-frequency amplitude matrix, and the nonnegative weight matrix is ​​used to reflect the amplitude of each basis feature changing over time. Since the matrix elements before and after decomposition remain nonnegative, the decomposition result has clear physical interpretability, facilitating the separation of local anomalies from global perturbations and background changes.

[0020] In this embodiment, a short-time Fourier transform is performed on the sequence of satellite magnetic field (e.g., Y component) after model subtraction and difference to construct a non-negative time-frequency amplitude matrix. This non-negative time-frequency amplitude matrix is ​​then subjected to depthwise expandable convolutional non-negative matrix decomposition to obtain local and global influence components. Earthquake-related local influence components are then screened based on criteria such as energy ratio, and anomalous orbits are determined by the exceedance rate. Finally, daily-scale cumulative statistics are performed on the anomalous orbits, and an earthquake anomaly indicator is constructed based on the degree of deviation from the relative background fitting trend, thus achieving earthquake anomaly detection and visualization output.

[0021] See Figure 1 As shown, a method for detecting seismic anomalies in satellite magnetic field data is described, the method comprising: S101, preprocess the satellite magnetic field observation data and convert it into a time-frequency representation to obtain a non-negative time-frequency amplitude matrix; S102, taking the non-negative time-frequency amplitude matrix as the decomposition object, using depthwise unfoldable convolutional non-negative matrix decomposition, the non-negative time-frequency amplitude matrix is ​​decomposed into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template; wherein, the two-dimensional time-frequency convolution template is used to characterize several basis features of the non-negative time-frequency amplitude matrix, and the non-negative weight matrix is ​​used to reflect the amplitude of each basis feature changing with time. S103, Based on the decomposition results, extract the earthquake-related local influence components and determine the abnormal trajectory of the earthquake-related local influence components; S104: The number of abnormal orbits is counted daily and accumulated. The earthquake anomalies are determined based on the degree of deviation of the accumulated results from the background fitted straight line.

[0022] In step S101, the satellite magnetic field observation data is preprocessed and represented in a time-frequency format to obtain a non-negative time-frequency amplitude matrix, including: Subtracting the CHAOS-6 magnetic field model data from the satellite magnetic field observation data yields the residual data, and the first difference of the residual data is then calculated to obtain the difference data. A short-time Fourier transform is performed on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the satellite magnetic field observation data.

[0023] First, the background field calculated using the CHAOS-6 magnetic field model data is obtained by subtracting the original satellite magnetic field observation data from the background field data. This removes the larger background values ​​in the geomagnetic field amplitude, yielding residual data. Then, the first-order difference of the residual data is taken to remove the remaining smaller background values, resulting in the change in the magnetic field data. This result can be used for subsequent operations such as short-time Fourier transform. The formula is expressed as: , , in This refers to the Y-component magnetic field data from the original satellite magnetic field observation data. The background field calculated for the CHAOS-6 magnetic field model. For the first The time corresponding to each data point For the first The time corresponding to each data point The length of each track data, This refers to the preprocessed result of the first-order difference of the residual data, i.e., the difference data.

[0024] In one embodiment, constructing the non-negative time-frequency amplitude matrix corresponding to the magnetic field data by performing a short-time Fourier transform on the differential data includes: The time-frequency transformation of the differential data is performed by short-time Fourier transform to obtain a time-frequency matrix that includes amplitude and phase information; The amplitude information of the time-frequency matrix is ​​extracted to construct a non-negative time-frequency amplitude matrix, which serves as the input matrix for depthwise expandable convolutional non-negative matrix decomposition (DFD). DFD requires that every element in the input matrix be non-negative. Since the Y-component magnetic field data is single-channel data and cannot form a matrix, and the preprocessing results contain negative values, a short-time Fourier transform is used to perform a time-frequency transformation on the data to construct a non-negative time-frequency amplitude matrix, which serves as the input matrix for DFD. Furthermore, the short-time Fourier transform has relevant physical meaning; its result reflects the time-frequency characteristics of the data, i.e., the changes in the amplitude and phase distribution of the data frequency over time. Here, the changes in the data frequency amplitude distribution over time are used to construct a non-negative matrix for decomposition, which can then be further decomposed based on the time-frequency amplitude characteristics of the data.

[0025] Suppose the time series of each track data after preprocessing is as follows: , in The length of each track data, For position is +1 orbital data.

[0026] The discrete form of the short-time Fourier transform is: , , , in For time series data, For window functions, Indicates conjugate. This represents the conjugate of the window function, and the length of the window function is... Step size is , The total number of windows Indicates the index of the window. Indicates frequency index, Representing the imaginary unit, we finally obtain The complex matrix is ​​its time-frequency matrix: , Each element in the matrix is ​​a complex number , For the real part, The imaginary part contains amplitude and phase information, and its amplitude is... Phase is .

[0027] Obtain its amplitude information The non-negative time-frequency amplitude matrix is: .

[0028] In step S102, the non-negative time-frequency amplitude matrix is ​​decomposed into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template using depthwise expandable convolutional non-negative matrix decomposition, including: Build with -The objective function is to minimize the divergence, the reconstruction error, and the weighted sum of the L1 norm of the multilayer activation coefficients; Based on the objective function, the iterative update process is unfolded into a deep network containing multiple cascaded network layers. Each layer corresponds to an update operation along the negative gradient direction of the objective function, transforming the iterative process into an end-to-end trainable network structure. The forward propagation of each layer performs reconstruction calculations, updates the non-negative weight matrix, and updates the two-dimensional time-frequency convolution template. After all network layers, the final two-dimensional time-frequency convolution template and non-negative weight matrix are output.

[0029] Specifically, the time-frequency amplitude matrix of the magnetic field data is decomposed using depth-expandable convolutional non-negative matrix decomposition to obtain the corresponding two-dimensional time-frequency convolutional template and the corresponding non-negative weight matrix of the frequency distribution characteristics, thus separating the local influence component caused by earthquakes from the global influence component caused by solar activity and geomagnetic activity.

[0030] Depth-expandable convolutional nonnegative matrix factorization (DEF) can decompose any high-dimensional nonnegative matrix into two low-dimensional nonnegative matrices, such that the product of the two low-dimensional nonnegative matrices approximately approximates the original high-dimensional nonnegative matrix. One of the low-dimensional nonnegative matrices is called the basis matrix, and the other is called the weight matrix. The basis matrix can be viewed as a set of foundations for the original high-dimensional nonnegative matrix, with its column vectors representing the characteristics of the data. The weight matrix is ​​the projection of the high-dimensional nonnegative matrix onto the basis vector matrix space, and the row vectors of the weight matrix represent the amplitude of the corresponding column vectors of the basis matrix over time. Since all elements in the original high-dimensional nonnegative matrix are nonnegative, the elements in the two decomposed low-dimensional nonnegative matrices are also nonnegative, making the decomposition results more physically meaningful. Based on this property, the nonnegative time-frequency matrix corresponding to magnetic field data can be used as the original matrix for decomposition, yielding several time-frequency distribution features of the data and their amplitudes in the time domain. This allows for the separation of local influence components caused by earthquakes from global influence components caused by solar and geomagnetic activity, without removing data from days with strong solar and geomagnetic activity, while simultaneously obtaining components more relevant to earthquake activity.

[0031] The mathematical model for depthwise expandable convolutional nonnegative matrix factorization is as follows: ; in Let the non-negative time-frequency amplitude matrix be denoted as . A nonnegative matrix, It is a two-dimensional time-frequency convolution template. for The non-negative weight matrix, Represents the frequency dimension size and the number of decomposed features. The choice satisfies ; "Indicates a one-dimensional convolution operation along the time dimension, The total number of bases for decomposition.

[0032] use - Minimize the divergence between the observation matrix and the reconstructed matrix, and introduce a weighted sum of the L1 norm of the multilayer activation coefficients as the optimization objective. The objective function is: , in, This is a non-negative time-frequency amplitude matrix, which here represents the observation matrix. For the reconstructed time-frequency matrix, For the sparse regularization weights of the activation coefficients of all expanded layers, Including convolution kernel parameters and all layers sparse terms Used to constrain the activation range of each mode, for - Divergence, For the objective function, when The KL divergence is expressed as: , in, For frequency index, For time frame index.

[0033] The basic idea of ​​Deeply Expandable Convolutional Nonnegative Matrix Factorization (DNF) is to expand the finite-step iterative process of traditional nonnegative matrix factorization into a deep network containing multiple cascaded network layers. Each network layer corresponds to one update of the convolutional nonnegative matrix factorization, thereby introducing learnable parameters and enabling end-to-end training driven by data while preserving the algorithm's structure and interpretability. Assuming the number of expanded layers is L, the parameters of layer 0 are initialized. For the first layer Network forward propagation consists of the following three parts: The first part is the reconstruction calculation, which uses the two-dimensional time-frequency convolution template of the previous layer and the non-negative weight matrix to perform one-dimensional convolution to reconstruct the time-frequency matrix of the current layer. The calculation formula is shown in the figure.

[0034] , in, for The reconstructed time-frequency matrix output by the layer. for The two-dimensional time-frequency convolution template output by the layer, for The non-negative weight matrix output by the layer.

[0035] The second part involves updating the non-negative weight matrix. Based on the gradient of the reconstruction error with respect to the non-negative weight matrix, and combined with the learnable stride parameter and non-negative activation function, the non-negative weight matrix is ​​updated. The update operator of the depthwise expandable convolution non-negative matrix factorization is parameterized and generalized to a learnable operator after depthwise expansion, resulting in the update formula shown below. ReLU is a non-negative activation function. For the first The learnable step size or scaling parameter of a layer. Represents all learnable hyperparameters in this layer: ; in, for The non-negative weight matrix output by the layer, for The non-negative weight matrix output from layer -1 This represents the gradient of the reconstruction error with respect to the non-negative weight matrix.

[0036] The third part concerns the update of the 2D time-frequency convolution template. Based on the gradient of the reconstruction error with respect to the 2D time-frequency convolution template, combined with the learnable stride parameter and non-negative activation function, the 2D time-frequency convolution template is updated. Learnable iterations are performed on the 2D time-frequency convolution template, and the update formula for the 2D time-frequency convolution template is shown in the following equation. For the first The learnable step size parameter of the layer, For the corresponding set of learnable hyperparameters, It is a non-negative activation function. Through inter-layer recursion, it is finally obtained at the th... The final output of the layer is shown in the following formula.

[0037] , , in, for The two-dimensional time-frequency convolution template output by the layer, To reconstruct the gradient of the error with respect to the two-dimensional time-frequency convolution template, for The two-dimensional time-frequency convolution template output by the layer, for The non-negative weight matrix output by the layer.

[0038] The iterative update process of depthwise expandable convolutional nonnegative matrix factorization is a convolutional version of multiplicative update. The multiplicative update process involves convolution and deconvolution. For example, the KL divergence objective function is shown in the following equation: , Here For transposed convolution, where , This represents the Hadamard element-wise product. For the first The learnable exponential parameter of the layer, It is a small positive number with numerical stability. When The time-frequency degenerates into a standard multiplication update. Furthermore, the 2D time-frequency convolution template is not directly optimized, but rather the original 2D time-frequency convolution template is optimized. It can be shared across layers or independent within each layer. To ensure non-negativity, the two-dimensional time-frequency convolution template is parameterized using Softplus transform. For unconstrained real-valued parameters, the parameters For numerical stability, the Softplus function is shown below: , , in The slope parameter is used to explicitly guarantee nonnegativity while maintaining trainability. These updates are essentially a series of deterministic nonnegative iterative operations on the nonnegative weight matrix and the two-dimensional time-frequency convolution template until convergence.

[0039] In one example, the rank r (number of extracted patterns) of the decomposition result is set to 3. Number of unfolded layers. The number of layers corresponding to unfolding the iterative process into a deep network is set to 200, and its size reflects the number of iterations to approach the optimal decomposition in each forward propagation. Convolutional time kernel length. The length is set to 9, which controls the support range of each convolutional template in the time (dimensional) direction, directly determining the event duration scale that the model can represent. The epoch for a complete forward and backward propagation process through the entire dataset is set to 500. The optimizer used for model training employs the Adam optimizer (Adaptive Moment Estimation), which dynamically adjusts the learning rate based on gradient-based first-moment and second-moment estimations. The initial learning rate is set to 0.001. The learning rate determines the magnitude of weight updates in each step of the Adam optimizer, significantly impacting the convergence speed and stability of the training process.

[0040] The above process yields the two-dimensional time-frequency convolution template and non-negative weight matrix for each track.

[0041] In step S103, earthquake-related local influence components are extracted based on the decomposition results, and abnormal orbit determination is performed on the earthquake-related local influence components. This includes using the energy ratio criterion to screen the earthquake-related local influence components from the decomposition results, and determining abnormal orbits of the earthquake-related local influence components based on the over-limit rate.

[0042] Based on the energy ratio, the local impact component caused by the earthquake is selected. Anomaly trajectory judgment is then performed on this local impact component using the over-limit rate method. This involves calculating the ratio of energy within the study area to the energy of the entire orbit in each row vector of the non-negative weight matrix obtained from non-negative matrix decomposition. The component with the largest energy ratio is selected as the local impact component caused by the earthquake. The root mean square (RMS) parameter, which reflects the average energy level of the data, is then used as a threshold parameter. Anomaly point judgment is based on whether the value of the largest component is greater than k times the RMS of the entire orbit. Here, k is empirically selected, and anomalies are considered to occur only within the study area, and the corresponding marker indicates that the anomaly is not caused by known ionospheric activity. Since the impact of earthquake events on the data is local, it only affects the data within the earthquake study area of ​​the orbit. In this case, the data energy is mainly concentrated in the earthquake-affected area. Solar activity and geomagnetic activity, on the other hand, are global events and therefore affect the data of the entire orbit. In this case, the data should be more evenly distributed across the entire orbit. Therefore, the component with the largest energy ratio between the study area and the entire orbit is selected, as this component's data energy is mainly concentrated within the earthquake study area and is most likely to reflect the local impact of the earthquake.

[0043] The formula for calculating the energy ratio is: , , , in Indicates the starting position of the data within the study area. This indicates the endpoint of the data within the study area. Multiple energy ratios can be obtained through calculation; the component with the largest energy ratio is selected as the earthquake-related local influence component. For the first The energy ratio of each component For the first The component in the first Weight values ​​on each time frame.

[0044] Anomaly detection is performed on selected components using the over-limit rate method. The root mean square (RMS) of the data, reflecting its average energy level, is chosen as the threshold parameter. Anomalies should significantly exceed this value and persist for a period of time. Since the original data has already undergone smoothing during the short-time Fourier transform, outlier detection can be directly performed based on whether the value of the largest component is greater than k times the RMS of the entire trajectory. A time series of length N is calculated. The root mean square formula is: , The threshold is Points exceeding a threshold are considered outliers. Furthermore, outliers are considered to occur only within the study area, and orbits whose corresponding flags indicate that the anomaly is not caused by known ionospheric activity are considered anomalous orbits.

[0045] In one embodiment, the number of anomalous orbits is counted daily and accumulated. Seismic anomalies are determined based on the degree of deviation of the accumulated results from the background fitted straight line, including: The daily number of abnormal orbits and the daily total number of orbits were standardized by deviation to obtain the cumulative normalized result of the daily number of abnormal orbits and the cumulative normalized result of the daily total number of orbits. The cumulative normalized results of the total number of orbits for each day are subjected to least squares linear fitting to obtain the background fitting line; The Sigmoid function is used to fit the cumulative normalized result of the number of daily abnormal orbits and the deviation of the fitted line from the background.

[0046] First, the number of anomalous orbits each day is accumulated to obtain a cumulative result. Then, a least-squares linear fit of all daily orbit counts is used as a background fitting line. The degree to which the cumulative daily anomalous orbit count deviates from the background fitting line indicates the change process of seismic anomalies. Because the results obtained after judging all orbits that meet the research conditions as anomalous orbits cannot intuitively show their relationship with earthquakes, the number of anomalous orbits each day is accumulated, and a least-squares linear fit of all daily orbit counts is used as a background fitting line. The degree to which the cumulative daily anomalous orbit count deviates from the background line reflects the impact of earthquakes on the data; the greater the deviation, the greater the impact.

[0047] The cumulative number of abnormal orbits each day is as follows: , in This is the cumulative total of all tracks per day. For the first The cumulative number of all orbits for the study days.

[0048] Then, the total number of all tracks for each day is accumulated, and the result is: , in, This is the cumulative total of all tracks per day. For the first The cumulative number of all orbits for the study days.

[0049] Since the number of all orbits per day differs significantly from the number of abnormal orbits per day, deviation standardization is used to normalize both.

[0050] The cumulative normalized result of the number of daily abnormal orbits was calculated. for: , in, For the first The cumulative normalized result of the number of anomalous orbits over the study days.

[0051] Cumulative normalized result of the total number of orbits per day for: , in, For the first The cumulative normalized result of the total number of orbits for the study days.

[0052] Cumulative normalized result of the total number of orbits for each day Perform a least-squares linear fit, and let the fitted line be: ; Where parameters and The estimated solution formula is: ; ; in for , For the number of days, for The first in One element; for , for The first in One element; for The first in Each element and The first in Multiplying by each element for The average value, for The average value; the parameters of the linear equation were obtained through calculation. and ; The degree to which the cumulative curve of daily anomalous orbits deviates from the background fitted straight line is used to display the earthquake change process. The Sigmoid function formula for calculating the degree of deviation from the linear fit is as follows: , in, This is the maximum value of the Sigmoid function (usually 1). To control the steepness of the curve, the slope of the curve is controlled. The center or threshold of the Sigmoid function, i.e., where the curve is located. The value at that location is Half of it.

[0053] On the other hand, embodiments of this application provide a satellite magnetic field data seismic anomaly detection system for performing the satellite magnetic field data seismic anomaly detection method of the above embodiments, including: The data preprocessing and time-frequency conversion module is used to preprocess satellite magnetic field observation data and convert it into a time-frequency representation to obtain a non-negative time-frequency amplitude matrix. The feature decomposition module is used to decompose the non-negative time-frequency amplitude matrix as the decomposition object, and to decompose the non-negative time-frequency amplitude matrix into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template using depthwise unfoldable convolution non-negative matrix decomposition; wherein, the two-dimensional time-frequency convolution template is used to characterize several basis features of the non-negative time-frequency amplitude matrix, and the non-negative weight matrix is ​​used to reflect the amplitude of each basis feature changing over time. The abnormal trajectory determination module is used to extract earthquake-related local influence components based on the decomposition results and determine the abnormal trajectories of the earthquake-related local influence components. The anomaly identification module is used to count the number of abnormal orbits on a daily basis and accumulate them. Based on the degree of deviation of the accumulated results from the background fitted straight line, the seismic anomaly is determined.

[0054] In one example, considering a magnitude 7.8 earthquake that occurred on April 16, 2016, the Y-component magnetic field data (1Hz) of the Swarm Alpha satellite constellation is used. The study area of ​​the earthquake is selected, centered on the epicenter. The rectangular area defined by longitude and latitude is the study area. The study period is from 90 days before the earthquake to 30 days after the earthquake, i.e., from January 17 to May 16, 2016. The epicenter, the affected area, and the study area are shown below. Figure 2 As shown, the pentagram represents the epicenter.

[0055] Read and input Y-component magnetic field data of Swarm Alpha satellite from January 17 to May 16, 2016. (1Hz), and invalid data caused by measurement errors of the instrument itself is removed based on the flag bits in the dataset. Due to the special characteristics of satellite data measurement, namely that the measured values ​​vary with both time and space, and the need to subsequently separate the local effects caused by earthquakes from the global components caused by solar and geomagnetic activity, geomagnetic latitude is used. The orbits within the range and passing through the earthquake research area are used as the research objects for data processing.

[0056] Taking data from April 8, 2016 as an example, the residual data curves of the original magnetic field data minus the CHAOS-6 model data, and the first-order difference data curves are as follows: Figure 3 As shown, Figure 3 In (a), the curve of the original magnetic field data is shown. It can be seen that the amplitude of the original magnetic field data is very large, and it is impossible to observe its changes. It is difficult to directly decompose or extract anomalies. Figure 3 (b) in the figure represents the curve of the residual data obtained by subtracting the CHAOS-6 model data from the original magnetic field data. Figure 3 (c) in the figure represents the curve of the first-order difference data of the residual data. The results show that the original magnetic field data has a very large amplitude. Subtracting the magnetic field model can remove the background field with a large amplitude, but the residual data still has a certain background amplitude, which needs further removal. Therefore, the first-order difference is calculated on the residual data to remove the remaining background amplitude and obtain the change in the magnetic field data. The final result shows a very small background amplitude and prominent singular values.

[0057] The time-frequency amplitude matrix result of the short-time Fourier transform is shown in the figure below. Figure 4 As shown, the horizontal axis represents geographical latitude, the vertical axis represents frequency, and the color indicates the amplitude. Therefore, it can be seen that single-channel magnetic field data can be transformed into... The matrix, which consists of non-negative numbers, can be used as the input matrix for depthwise expandable convolutional non-negative matrix decomposition. It also has a physical meaning that reflects the time-frequency amplitude characteristics of the signal, and subsequent processing can decompose the signal based on these characteristics.

[0058] The decomposition result is a non-negative weight matrix as follows: Figure 5 In (a) and the two-dimensional time-frequency convolution template, as shown in Figure 5 (b) in the middle Figure 5 (c) and Figure 5 (d) in Figure 5 In (a), the first component hs1, the second component hs2, and the third component hs3 are the two-dimensional time-frequency convolution templates corresponding to the three components. Figure 5 (b) Figure 5 (c) and Figure 5 (d) shows the three two-dimensional time-frequency convolution templates for the orbital. Figure 5 The larger amplitude of the first component hs1 in (b) only appears within the study area, and the energy is also mainly distributed within the study area. This component should be a local component affected by the earthquake. Figure 5 The amplitude of the second component hs2 in (c) is not fluctuating much, and its energy is evenly distributed throughout the entire orbit, which may be the background component. Figure 5The maximum amplitude of the third component hs3 in (d) appears outside the study area, and there are also values ​​with large amplitudes within the study area. This component may be a global component influenced by solar or geomagnetic activity. To show its separation effect in the time domain, the decomposition results are reconstructed to the time domain results as follows: Figure 6 As shown. Among them. Figure 6 (a) in the figure represents the time-domain component corresponding to the first component hs1, which has obvious singular values ​​only within the study region; Figure 6 (b) is the time-domain component corresponding to the second component hs2, which has no obvious singular values ​​and is uniformly distributed throughout the entire orbit; Figure 6 In the diagram, (c) represents the time-domain component corresponding to the third component hs3. The largest singular value appears outside the study area, while singular values ​​with smaller amplitudes also exist within the study area. The decomposition results show that using depthwise expandable convolutional nonnegative matrix decomposition to decompose the nonnegative time-frequency amplitude matrix of the data can yield its frequency distribution characteristics and temporal variation amplitude. This separates the local impact components caused by earthquakes from the global impact components caused by solar and geomagnetic activity, and utilizes the components more relevant to seismic activity for earthquake anomaly detection.

[0059] The normalized result of the daily cumulative curve of abnormal orbits and the fitted line to the background is as follows: Figure 7 As shown, there are three components in total, namely the first component background fitting line, the second component background fitting line, and the third component background fitting line. Figure 7 Different colored dots represent the cumulative results of the first, second, and third components, which are distributed along their respective background fitted lines. The sigmoid function fitting result of the first component includes three segments of accelerated growth in anomalous cumulative growth. The resulting cumulative curve initially deviates abnormally from the background fitted line, reaching its maximum deviation near the earthquake site, then begins to recover a few days after the earthquake, and finally returns to normal. These results demonstrate that the method effectively detects anomalies generated by the target earthquake and their changes with seismic activity, including the onset of pre-earthquake anomalies and their recovery after the earthquake.

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

Claims

1. A method for detecting seismic anomalies in satellite magnetic field data, characterized in that, include: The satellite magnetic field observation data is preprocessed and converted into a time-frequency representation to obtain a non-negative time-frequency amplitude matrix; The non-negative time-frequency amplitude matrix is ​​used as the decomposition object. Using depthwise unfoldable convolutional non-negative matrix decomposition, the non-negative time-frequency amplitude matrix is ​​decomposed into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template. The two-dimensional time-frequency convolution template is used to characterize several basis features of the non-negative time-frequency amplitude matrix, and the non-negative weight matrix is ​​used to reflect the amplitude of each basis feature as time changes. Based on the decomposition results, earthquake-related local impact components are extracted, and abnormal trajectories are determined for these components. The number of abnormal orbits is counted daily and accumulated. The degree of deviation of the accumulated results from the fitted straight line relative to the background is used to determine the earthquake anomalies.

2. The method for detecting seismic anomalies in satellite magnetic field data according to claim 1, characterized in that, The satellite magnetic field observation data is preprocessed and represented in a time-frequency format to obtain a non-negative time-frequency amplitude matrix, including: Subtracting the CHAOS-6 magnetic field model data from the satellite magnetic field observation data yields the residual data, and the first difference of the residual data is then calculated to obtain the difference data. A short-time Fourier transform is performed on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the satellite magnetic field observation data.

3. The method for detecting seismic anomalies in satellite magnetic field data according to claim 2, characterized in that, Perform a short-time Fourier transform on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the satellite magnetic field observation data, including: The time-frequency transformation of the differential data is performed by short-time Fourier transform to obtain a time-frequency matrix that includes amplitude and phase information; The amplitude information of the time-frequency matrix is ​​extracted to construct a non-negative time-frequency amplitude matrix.

4. The method for detecting seismic anomalies in satellite magnetic field data according to claim 1, characterized in that, Based on the decomposition results, earthquake-related local impact components are extracted, and abnormal trajectories are determined for these components. This includes using the energy ratio criterion to screen the earthquake-related local impact components from the decomposition results, and determining abnormal trajectories for these components based on the exceedance rate.

5. The method for detecting seismic anomalies in satellite magnetic field data according to claim 1, characterized in that, Using depthwise expandable convolutional nonnegative matrix factorization, the nonnegative time-frequency magnitude matrix is ​​decomposed into a combination of a nonnegative weight matrix and a two-dimensional time-frequency convolutional template, including: Construct an objective function that minimizes the reconstruction error and the weighted sum of the L1 norm of the multilayer activation coefficients based on the β-divergence. Based on the objective function, the iterative update process is unfolded into a deep network containing multiple cascaded network layers. Each layer corresponds to an update operation along the negative gradient direction of the objective function, transforming the iterative process into an end-to-end trainable network structure. The forward propagation of each layer performs reconstruction calculations, updates the non-negative weight matrix, and updates the two-dimensional time-frequency convolution template. After all network layers, the final two-dimensional time-frequency convolution template and non-negative weight matrix are output.

6. The method for detecting seismic anomalies in satellite magnetic field data according to claim 5, characterized in that, The forward propagation of each layer performs reconstruction computation, updates the non-negative weight matrix, and updates the 2D time-frequency convolution template, including: The time-frequency matrix of the current layer is reconstructed by performing a one-dimensional convolution with the two-dimensional time-frequency convolution template of the previous layer and the non-negative weight matrix. Based on the gradient of the reconstruction error with respect to the non-negative weight matrix, combined with the learnable step size parameter and the non-negative activation function, the non-negative weight matrix is ​​updated. Based on the gradient of the reconstruction error with respect to the 2D time-frequency convolution template, and combined with the learnable stride parameter and non-negative activation function, the 2D time-frequency convolution template is updated.

7. The method for detecting seismic anomalies in satellite magnetic field data according to claim 4, characterized in that, The energy ratio criterion is used to screen the earthquake-related local influence components from the decomposition results, and the abnormal trajectories of the earthquake-related local influence components are determined based on the exceedance rate, including: By calculating the ratio of energy in the study area to the energy of the entire orbit in each row vector of the non-negative weight matrix obtained by non-negative time-frequency amplitude matrix decomposition, the component with the largest energy ratio is selected as the earthquake-related local influence component. The root mean square of the average energy level of the reaction data is used as the threshold parameter. Anomalies are judged based on whether the value of the largest component is greater than k times the root mean square of the entire orbit. Anomalies are defined as orbits that only appear in the study area and whose corresponding flags indicate that the anomalies are not caused by known ionospheric activity.

8. The method for detecting seismic anomalies in satellite magnetic field data according to claim 1, characterized in that, The number of anomalous orbits is counted daily and accumulated. Seismic anomalies are determined based on the degree of deviation of the accumulated results from the fitted straight line relative to the background, including: The daily number of abnormal orbits and the daily total number of orbits were standardized by deviation to obtain the cumulative normalized result of the daily number of abnormal orbits and the cumulative normalized result of the daily total number of orbits. The cumulative normalized results of the total number of orbits for each day are subjected to least squares linear fitting to obtain the background fitting line; The Sigmoid function is used to fit the cumulative normalized result of the number of daily abnormal orbits and the deviation of the fitted line from the background.

9. A satellite magnetic field data seismic anomaly detection system, used to execute the satellite magnetic field data seismic anomaly detection method according to any one of claims 1-8, characterized in that, include: The data preprocessing and time-frequency conversion module is used to preprocess satellite magnetic field observation data and convert it into a time-frequency representation to obtain a non-negative time-frequency amplitude matrix. The feature decomposition module is used to decompose the non-negative time-frequency amplitude matrix as the decomposition object, and to decompose the non-negative time-frequency amplitude matrix into a combination of a non-negative weight matrix and a two-dimensional time-frequency convolution template using depthwise unfoldable convolution non-negative matrix decomposition; wherein, the two-dimensional time-frequency convolution template is used to characterize several basis features of the non-negative time-frequency amplitude matrix, and the non-negative weight matrix is ​​used to reflect the amplitude of each basis feature changing over time. The abnormal trajectory determination module is used to extract earthquake-related local influence components based on the decomposition results and determine the abnormal trajectories of the earthquake-related local influence components. The anomaly identification module is used to count the number of abnormal orbits on a daily basis and accumulate them. Based on the degree of deviation of the accumulated results from the background fitted straight line, the seismic anomaly is determined.

Images