Satellite-earth multi-layered seismic electromagnetic anomaly fusion extraction method and system

Abstract

The application belongs to the technical field of earthquake precursor detection, and particularly relates to a satellite-ground multi-layer earthquake electromagnetic anomaly fusion extraction method and system, which comprises the following steps: time alignment and truncation are performed on the preprocessed satellite-ground magnetic field data to form a synchronous time sequence; time-frequency transformation is performed on satellite magnetic field data and station magnetic field data in the synchronous time sequence to obtain corresponding time-frequency amplitude spectra; the two time-frequency amplitude spectra are superimposed along a third dimension to construct a third-order non-negative tensor; non-negative tensor decomposition is performed on the third-order non-negative tensor to obtain a basis matrix, a weight matrix and a contribution matrix; based on the contribution difference between the satellite magnetic field data and the station magnetic field data, earthquake-related components are screened out; and according to the earthquake-related components, the abnormal extraction range is dynamically determined and abnormal points are extracted from the weight matrix. The application realizes the fusion of satellite-station heterogeneous electromagnetic data, and breaks through the limitation of the prior art which can only process satellite homogenous data.

Description

Technical Field

[0001] This application belongs to the field of earthquake precursor detection technology, specifically relating to a method and system for fusion extraction of earthquake electromagnetic anomalies across multiple satellite and ground layers. Background Technology

[0002] Earthquake precursor research is one of the keys to improving earthquake prediction and forecasting. During the earthquake gestation process, precursor anomalies often appear in the lithosphere, atmosphere, and ionosphere. Among these, earthquake electromagnetic precursor anomalies are considered one of the research directions most likely to achieve breakthroughs in short-term earthquake prediction.

[0003] Chinese invention patent CN113435259A discloses a method for seismic anomaly extraction based on tensor decomposition of satellite magnetic field data fusion. This scheme targets the same-source magnetic field data from multiple satellites, achieving multi-satellite magnetic field data fusion and seismic anomaly extraction through tensor decomposition, providing a feasible approach for seismic anomaly identification from satellite electromagnetic data. However, this scheme suffers from the following insurmountable core defects: First, the data source is limited to satellite magnetic field data, only acquiring macroscopic trend information of the upper-air ionospheric electromagnetic field, completely ignoring high-precision local lithosphere observation data from stations, resulting in a severe information dimensional gap and an inability to capture the microscopic details of earthquake precursors; Second, it lacks heterogeneous data compatibility design, completely neglecting the differences in scale, time duration, and spatial coverage between satellite and ground station magnetic field data, thus failing to achieve the fusion of multi-sphere electromagnetic data from the lithosphere and ionosphere; Third, the seismic component screening rules are designed only for single-source satellite magnetic field data, employing the energy-entropy ratio criterion, which cannot identify seismic coupling anomalies occurring simultaneously in multiple spheres, easily leading to misjudgments of non-seismic anomalies and limiting the accuracy of anomaly extraction.

[0004] Furthermore, existing methods for extracting seismic electromagnetic anomalies are typically limited to extracting anomaly information within a single sphere. For example, fusing dual-satellite electromagnetic data to separate weak seismic electromagnetic signals from the third-order time-frequency tensor of the dual-satellite electromagnetic data, thereby extracting seismic magnetic field fusion anomalies, focuses on "satellite-to-satellite" dual-data fusion. This method also suffers from the aforementioned problems of incompatibility between single-source data dimensional faults and heterogeneous data, making it difficult to accurately extract seismic electromagnetic anomalies across multiple spheres.

[0005] Therefore, there is an urgent need for a seismic electromagnetic anomaly extraction method that can integrate heterogeneous electromagnetic data from satellites and stations and identify synchronous anomalies across multiple spheres of the lithosphere and ionosphere, in order to solve the problems of limited information dimensions, incompatibility of heterogeneous data, and low reliability of seismic coupling anomaly identification in existing technologies. Summary of the Invention

[0006] This application provides a method and system for fusion extraction of seismic electromagnetic anomalies across multiple satellite and ground layers, addressing the problems of limited information dimensions, incompatibility of heterogeneous data, and low reliability of seismic coupling anomaly identification.

[0007] The first aspect of this application provides a method for fusing and extracting multi-sphere seismic electromagnetic anomalies from space and Earth, including: The satellite magnetic field data and the station magnetic field data were preprocessed to remove internal source fields and noise, resulting in preprocessed satellite-to-ground magnetic field data. The preprocessed satellite-Ground magnetic field data is time-aligned and truncated to form a synchronized time series; Time-frequency transformations were performed on the satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectra. The two time-frequency amplitude spectra are superimposed along the third dimension to construct a third-order non-negative tensor; The third-order nonnegative tensor is decomposed into a nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix. Based on the difference in contribution between satellite magnetic field data and station magnetic field data in the contribution matrix, earthquake-related components are selected. Based on the earthquake correlation components, the anomaly extraction range is dynamically determined and anomaly points are extracted from the weight matrix.

[0008] Furthermore, the satellite magnetic field data and the station magnetic field data are preprocessed separately to remove internal source fields and noise, resulting in preprocessed satellite-to-ground magnetic field data, including: The db4 wavelet is selected as the wavelet basis function, and the total number of decomposition levels L is preset; the magnetic field signal is subjected to L-level multi-scale wavelet decomposition, and the magnetic field signal includes satellite magnetic field data and station magnetic field data; In the first-level decomposition, the magnetic field signal is convolved using the low-pass filter coefficients and the high-pass filter coefficients respectively, and then downsampled to obtain the approximate components and detail components of the first level. From the second layer to the Lth layer, each layer repeats the convolution operation and downsampling above on the approximate components obtained from the previous layer, and obtains the approximate components and detail components of the current layer in turn. After decomposition, all effective components are retained and inverse wavelet transform is performed: each layer of components is upsampled layer by layer, convolved with the reconstruction filter respectively, and then the convolution results are superimposed to obtain the preprocessed star-ground magnetic field data.

[0009] Furthermore, the preprocessed satellite-to-ground magnetic field data is time-aligned and truncated to form a synchronized time series, including: based on the time information of the satellite magnetic field data, selecting station magnetic field data with the same time point as the satellite magnetic field data; and performing alignment and truncation processing on the station magnetic field data to obtain two electromagnetic signals with the same time length. In this context, the length of the satellite Y-component magnetic field data and its corresponding time series in the satellite magnetic field data is denoted as N; the length of the station magnetic field data and its corresponding time series is denoted as M, and N is less than M; the truncated station magnetic field data is denoted as a new time series, which is taken from the data segment in the station magnetic field data before truncation, within the time interval corresponding to the first satellite time point to the Nth satellite time point.

[0010] Furthermore, time-frequency transformations are performed on the satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectra, including: Synchronous compressed wavelet transforms were performed on satellite magnetic field data and station magnetic field data in the synchronous time series, respectively. The instantaneous frequency at each moment is obtained from the wavelet coefficients obtained by wavelet transform; Based on a series of pre-set discrete compression target frequencies, for each calculated instantaneous frequency, it is determined which specified frequency interval the instantaneous frequency falls within, and all wavelet coefficient energy within the specified frequency interval is compressed to the corresponding compression target frequency.

[0011] Furthermore, the third-order nonnegative tensor is decomposed into a nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix, including: The third-order nonnegative tensor is decomposed into the sum of R rank tensors using the tensor decomposition method, resulting in R components. Based on KL divergence to measure the difference between a third-order nonnegative tensor and an approximate tensor, an objective function is established, in which each element of the third-order nonnegative tensor and the corresponding element of the approximate third-order tensor participate in the KL divergence calculation, and the approximate third-order tensor is obtained by summing R rank tensors. The third-order nonnegative tensor is expanded along three different moduli directions using n-modulus expansions to obtain the matrix representation of the original tensor. The approximate third-order tensor is also expanded using n-modulus expansions and is represented as the transpose of the Khatri-Rao product of the decomposition matrix and the remaining factor matrices. Let the intermediate variable be equal to the transpose of the Khatri-Rao product of the remaining factor matrices, and rewrite the objective function as a function of the decomposition matrix of the current modulus direction; The decomposition matrices of each modulus direction are optimized and updated sequentially using an alternating iterative method to obtain the basis matrix, weight matrix, and contribution matrix.

[0012] Furthermore, based on the difference in contribution rates between satellite magnetic field data and station magnetic field data in the contribution matrix, earthquake-related components are selected, including: calculating the absolute value of the difference between the contribution rate of satellite magnetic field data and the contribution rate of station magnetic field data in each component of the contribution matrix, wherein the absolute value reflects the degree of consistency between the contributions of the two data sources in the corresponding component. Arrange all R components in ascending order of their absolute values; The component ranked first is selected as the earthquake-related component.

[0013] Furthermore, based on the earthquake correlation components, the anomaly extraction range is dynamically determined, including: Determine whether the absolute value of the difference between the satellite contribution rate and the station contribution rate in the contribution matrix corresponding to the earthquake correlation component is greater than a preset threshold; If the absolute value is greater than the preset threshold, it will not be considered a valid seismic component. If the absolute value is less than or equal to a preset threshold, determine the main contributing data source of the earthquake-related component. If the main contributing data source is the station's magnetic field data, the research scope will be set to the complete time interval of the station's magnetic field data. If the primary source of contributing data is satellite magnetic field data, the research scope will be set to the corresponding time interval within the earthquake research area traversed by the satellite.

[0014] Furthermore, outliers are extracted from the weight matrix, including: Calculate the root mean square of all data for the seismic component in the corresponding weight matrix within the study area; Multiply the preset empirical parameter by the root mean square value, and set the product result as the anomaly extraction threshold; When the value of a data point within the study area corresponding to the seismic component in the weight matrix is ​​greater than the anomaly extraction threshold, the data point is determined to be a seismic anomaly.

[0015] A second aspect of this application provides a multi-sphere seismic electromagnetic anomaly fusion extraction system, comprising: The preprocessing module is used to preprocess satellite magnetic field data and station magnetic field data respectively, remove internal source fields and noise, and obtain preprocessed satellite-to-ground magnetic field data; The time alignment module is used to perform time alignment and truncation on the preprocessed satellite-Ground magnetic field data to form a synchronized time series. The time-frequency conversion module is used to perform time-frequency conversion on satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectrum. The tensor construction module is used to superimpose two time-frequency amplitude spectra along the third dimension to construct a third-order non-negative tensor. The tensor decomposition module is used to perform non-negative tensor decomposition on the third-order non-negative tensor to obtain the basis matrix, weight matrix and contribution matrix. The component filtering module is used to filter out earthquake-related components based on the difference in contribution between satellite magnetic field data and station magnetic field data in the contribution matrix. The anomaly extraction module is used to dynamically determine the anomaly extraction range based on the earthquake correlation components and extract anomaly points from the weight matrix.

[0016] Compared with the prior art, the advantages of this application are as follows: This application realizes the fusion of heterogeneous electromagnetic data from satellites and stations, breaks through the limitation of existing technologies that can only process satellite data of the same source, realizes the joint analysis of electromagnetic data from multiple spheres of the lithosphere and ionosphere, and fills the information dimension gap of a single data source. This application solves the problem of mismatch in time, length and magnitude of heterogeneous satellite-ground data by time alignment and truncation, providing a foundation for multi-concentric data fusion and solving the pain point of existing technologies being incompatible with heterogeneous data. Based on the multi-sphere coupling mechanism of earthquake gestation, this application designs a screening rule for the seismic related components of the star-ground channel contribution difference, which can accurately identify seismic coupling anomalies that occur synchronously in multiple spheres, greatly reduce the misjudgment rate of non-seismic interference, and effectively improve the reliability of anomaly identification. The design of dynamic anomaly extraction range is adopted, which adaptively adjusts the extraction range according to the contribution ratio of the components, thus solving the problem of easy omission / false judgment in the fixed range of the existing technology.

[0017] In achieving the aforementioned beneficial effects, this application overcomes the following technical challenges: First, satellite magnetic field data and station magnetic field data differ significantly in physical scale, magnitude, and time sampling frequency, making direct fusion impossible. This application solves the problem of heterogeneous data compatibility by constructing a unified time-frequency tensor through time-aligned truncation and synchronous compressed wavelet transform. Second, multi-sphere seismic electromagnetic anomalies exhibit coupling characteristics, but existing methods cannot effectively separate common anomalies from non-seismic interference. Based on the multi-sphere coupling mechanism of seismic origination, this application designs a satellite-ground contribution difference screening criterion, ensuring that the selected components have clear geophysical significance rather than being mere mathematical statistical results. Attached Figure Description

[0018] Figure 1 A flowchart of the multi-sphere seismic electromagnetic anomaly fusion extraction method provided in the embodiments of this application; Figure 2 Example diagrams of background removal for satellite magnetic field data stabilization provided in this application embodiment, wherein (a) is the original satellite magnetic field data and (b) is the satellite magnetic field data after background removal; Figure 3The following is an example diagram of wavelet decomposition for removing low-frequency trends and high-frequency noise from station magnetic field data, provided in an embodiment of this application. (a) shows the original station magnetic field data, and (b) shows the denoised station magnetic field data. Figure 4 (a) Satellite magnetic field data and (b) Example diagram of aligned truncation of station electromagnetic signals provided for embodiments of this application; Figure 5 Examples of time-frequency amplitude spectra of satellite magnetic field data synchronously compressed wavelet transform provided in the embodiments of this application: (a) time-frequency amplitude spectrum of satellite magnetic field data, (b) time-frequency amplitude spectrum of station magnetic field data; Figure 6 Example diagram of constructing a third-order nonnegative tensor by spectral superposition provided in the embodiments of this application; Figure 7 An example diagram of the earthquake-related component anomaly extraction results provided in the embodiments of this application. Detailed Implementation

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

[0020] This application employs the fusion of two types of data: satellite magnetic field data and station magnetic field data. Satellite magnetic field data refers to the Earth's magnetic field observation data acquired by the magnetometer carried by a low-orbit geophysical exploration satellite, which includes three components of the magnetic field vector: the northward component X, the eastward component Y, and the vertical component Z. Station magnetic field data refers to the geomagnetic field variation data continuously recorded by fixed geomagnetic observation stations deployed on the Earth's surface.

[0021] This application integrates the two types of magnetic field data mentioned above to achieve the fusion extraction of seismic electromagnetic anomalies across multiple spheres of the Earth and space.

[0022] See Figure 1 As shown in the figure, an embodiment of this application provides a method for fusing and extracting multi-sphere seismic electromagnetic anomalies, including: S101 preprocesses both satellite magnetic field data and station magnetic field data to remove internal source fields and noise, resulting in preprocessed satellite-to-Earth magnetic field data. For satellite magnetic field data, firstly, an international reference geomagnetic field model or a high-precision geomagnetic field model (such as the CHAOS model or IGRF model) is used to calculate and subtract internal source field components such as the core field and crustal field. Then, large-scale external source field interference generated by the magnetospheric current system and the ionospheric current system is filtered out. Finally, anomalous jumps and high-frequency noise caused by solar high-energy particle events, satellite platform attitude disturbances, and instrument noise are removed. The preprocessed satellite magnetic field data mainly retains magnetic field anomaly information related to medium- and small-scale geological activities.

[0023] For station magnetic field data, the geomagnetic field model is also used to subtract the internal source field components. Considering the station's observation environment, it is also necessary to correct for instrument baseline drift and temperature effects; pulse interference and background noise caused by thunderstorms, human activities (such as rail traffic and industrial equipment) are removed; and short-term missing data are appropriately interpolated or labeled. The preprocessed station magnetic field data reflects the continuous change of the magnetic field at a fixed geographical location over time.

[0024] S102, the preprocessed satellite-to-ground magnetic field data is time-aligned and truncated to form a synchronized time series. Because the satellite flies along its orbit, its data sampling interval is relatively large, while the station magnetic field data is a continuous high-sampling-rate record; the time length and number of sampling points of the two are usually inconsistent. Using the time interval of the satellite magnetic field data as a reference, a data segment that completely corresponds to the satellite's time interval is extracted from the station magnetic field data to ensure strict time alignment between the two sets of data. In the aligned synchronized time series, the number of points in the satellite magnetic field data is N, and the number of points in the station magnetic field data is M, where N is less than M.

[0025] S103, time-frequency transformations are performed on satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectra; synchronous compressed wavelet transform is used as the time-frequency transformation method. First, a continuous wavelet transform is performed on the magnetic field data to obtain wavelet coefficients that depend on the scale factor and translation factor; then, the instantaneous frequency at each moment is calculated based on the wavelet coefficients; finally, the values ​​around any frequency are compressed to the current frequency, thereby obtaining a high-resolution time-frequency amplitude spectrum. Synchronous compressed wavelet transform can effectively improve the time-frequency clustering of transient anomalies in the magnetic field signal.

[0026] S104, two time-frequency amplitude spectra are superimposed along the third dimension to construct a third-order non-negative tensor; the time-frequency amplitude spectrum of satellite magnetic field data is used as the first channel of the tensor, and the time-frequency amplitude spectrum of station magnetic field data is used as the second channel of the tensor, stacked along the third dimension perpendicular to the time-frequency plane to form a third-order non-negative tensor. This third-order non-negative tensor simultaneously integrates the spatial distribution characteristics of satellite magnetic field data and the continuous evolution characteristics of station magnetic field data in the temporal dimension.

[0027] S105, perform nonnegative tensor decomposition on the third-order nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix; using a nonnegative tensor decomposition method based on KL divergence, the original third-order nonnegative tensor is approximated as the sum of R rank tensors, corresponding to R components. During the decomposition process, an alternating iterative optimization algorithm is used to update the decomposition matrix along the three modulus directions respectively. After each iteration, convergence is checked, and iteration stops when the convergence tolerance is met or the maximum number of iterations is reached.

[0028] The decomposition yields three matrices: the first matrix, called the basis matrix, denoted as W, is related to the frequency dimension and is used to characterize the basic features of the magnetic field signal at different frequencies; the second matrix, called the weight matrix, denoted as H, is related to the time or geographical location dimension and is used to characterize the intensity of the anomalous signal changes over time or space; the third matrix, called the contribution matrix, denoted as C, is 2 times R, where R represents the number of rank tensors, and each column corresponds to the contribution rate of satellite data and the contribution rate of station data in a decomposition component.

[0029] S106, Based on the difference in contribution rates between satellite magnetic field data and station magnetic field data in the contribution matrix, earthquake-related components are selected. For each decomposed component, the absolute value of the difference between the contribution rate of satellite data and the contribution rate of station data in that component is calculated. The smaller the difference, the more consistent the characteristics of satellite magnetic field data and station magnetic field data are in that decomposed component, and the more likely it is to reflect a common anomaly in multiple spheres caused by earthquake-related geological activities. All R components are sorted in ascending order of absolute value, and the component with the smallest difference is selected as the earthquake-related component.

[0030] For the selected seismic correlation components, it is further determined whether their absolute difference is less than or equal to a preset threshold. If it is greater than the threshold, the component is considered to lack common satellite-ground characteristics and is not considered a valid seismic component. When the absolute value meets the condition, the main contribution source of the decomposition component is determined based on the contribution matrix: if it is mainly contributed by station data, the corresponding column of the decomposition component in the weight matrix mainly reflects the temporal characteristics; if it is mainly contributed by satellite magnetic field data, the corresponding column of the decomposition component in the weight matrix mainly reflects the geographical location characteristics.

[0031] S107, Based on the earthquake correlation components, dynamically determine the anomaly extraction range and extract anomaly points from the weight matrix.

[0032] If the earthquake-related components are mainly contributed by station magnetic field data, the time range for anomaly extraction is set to the complete time interval of the entire synchronous time series; if they are mainly contributed by satellite magnetic field data, the research scope is set to the time interval during which the satellite passes through the earthquake research area.

[0033] Within the defined study area, all data points corresponding to the seismic components in the weight matrix are extracted, and the root mean square (RMS) values ​​of these data points are calculated. A preset empirical parameter is multiplied by this RMS value to obtain the anomaly extraction threshold. Each data point within the study area is iterated over; when the value of a data point is greater than the threshold, the point is determined to be a seismic anomaly. Simultaneously, the satellite orbit data and station data corresponding to this data point are marked as seismic-related data, and the detected anomalies are cumulatively statistically analyzed.

[0034] Through the above steps, the fusion processing of satellite magnetic field data and station magnetic field data was realized, and the automated extraction of earthquake-related anomalies was completed, providing quantitative data support for earthquake monitoring and earthquake precursor research.

[0035] In step S101, Swarm A satellite magnetic field data is read, daily magnetic field data is divided into multiple orbits, orbits passing through the seismic study area are selected, and these orbits are converted into time series. The time series of station magnetic field data is then read. The original Y-component magnetic field data of each satellite orbit is subtracted from the source field of the CHAOS-7 model to remove the main magnetic field of the Earth's core and the lithosphere magnetic field of the crust, resulting in the remaining Y-component magnetic field. That is: the original Y-component magnetic field data of each satellite orbit minus the source field of the CHAOS-7 model is calculated as follows: , in, For the original Y component magnetic field, Y-component data of the internal source field calculated for the CHAOS-7 model. This represents the magnetic field of the remaining Y component.

[0036] In one example, the daily satellite magnetic field data is divided into multiple orbits as follows: the satellite magnetic field data is stored on a daily basis, and it is divided with 50°S and 50°N as orbit endpoints, resulting in 32 satellite orbits per day.

[0037] The orbits passing through the earthquake study area were selected, and the radius of the study area was calculated using the formula studied by Dobrovolsky: , in, To study the magnitude of earthquakes, The unit is kilometers.

[0038] The magnetic field-latitude signal was then converted into a magnetic field-time series, and station magnetic field data for the same day were read.

[0039] Preprocessing was performed on both satellite magnetic field data and station magnetic field data to remove intrinsic fields and noise, resulting in preprocessed satellite-to-Earth magnetic field data, including: The db4 wavelet is selected as the wavelet basis function, and the total number of decomposition levels L is preset; the magnetic field signal is subjected to L-level multi-scale wavelet decomposition, and the magnetic field signal includes satellite magnetic field data and station magnetic field data; In the first-level decomposition, the magnetic field signal is convolved using the low-pass filter coefficients and the high-pass filter coefficients respectively, and then downsampled to obtain the approximate components and detail components of the first level. From the second layer to the Lth layer, each layer repeats the convolution operation and downsampling above on the approximate components obtained from the previous layer to obtain the approximate components and detail components of the current layer in turn. After decomposition, all effective components are retained and inverse wavelet transform is performed: each layer of components is upsampled layer by layer, convolved with the reconstruction filter respectively, and then the convolution results are superimposed to obtain the preprocessed star-ground magnetic field data.

[0040] In one example, the Daubechies (db4) wavelet is selected as the wavelet basis function, the number of decomposition levels L is determined, and the magnetic field signal B (including satellite magnetic field data and station magnetic field signals) is decomposed into L levels of multi-scale wavelet decomposition using the Mallat fast algorithm: Level 1 decomposition: , , in, These are the coefficients of the low-pass filter. These are the coefficients of the high-pass filter. For the first-level approximate component, For the first layer of detail components, For the decomposition of the first Data points, For the summation variable.

[0041] No. Layer decomposition ( ): For the front +1 layer approximation components Repeat the above operation to obtain the first... Layer approximate components and the Layer detail components ; The effective components are retained, and inverse wavelet transform is performed using the inverse Mallat algorithm. After upsampling (zero-placing) at each level, the signals are convolved with the reconstruction filter and then superimposed to obtain the denoised target signal. This refers to the preprocessed satellite-Ground magnetic field data. The reconstruction filter is a filter bank used in the inverse wavelet transform process. It is convolved with the upsampled approximate components and detail components respectively and then superimposed to achieve signal reconstruction and recovery.

[0042] In step S102, since satellite magnetic field data and station magnetic field data originate from different observation platforms, they have different sampling frequencies, time start points, and data lengths. Satellites fly along their orbits, and their data sampling intervals are relatively large and the time intervals are relatively fixed; stations conduct continuous observations from fixed locations, with high sampling rates and continuous time, but the data length is usually much longer than the time coverage of satellite data. To achieve the fusion analysis of the two types of data, the station magnetic field data must be time-aligned and truncated to ensure strict time synchronization with the satellite magnetic field data.

[0043] The preprocessed satellite-to-ground magnetic field data is time-aligned and truncated to form a synchronized time series, including: selecting station magnetic field data with the same time point as the satellite magnetic field data based on the time information of the satellite magnetic field data; and aligning and truncating the station magnetic field data to obtain two electromagnetic signals with the same time length. In this context, the length of the satellite Y-component magnetic field data and its corresponding time series in the satellite magnetic field data is denoted as N; the length of the station magnetic field data and its corresponding time series is denoted as M, and N is less than M; the truncated station magnetic field data is denoted as a new time series, which is taken from the data segment in the station magnetic field data before truncation, within the time interval corresponding to the first satellite time point to the Nth satellite time point.

[0044] The formula is explained below: The electromagnetic data from the stations is read, and satellite electromagnetic data from stations with the same time are filtered based on the time of the satellite magnetic field data. The station data is then aligned and truncated to obtain two electromagnetic signals with the same time: ; ; .

[0045] In the formula This refers to the satellite Y-component magnetic field data in the satellite magnetic field data, and the corresponding time series is... The length is N. The data consists of the station's magnetic field data, and its corresponding time series is as follows: The length is M, where , This is the truncated station magnetic field data. The first component of the satellite's Y-component magnetic field data represents the... Data points, This indicates the first [item] in the station's magnetic field data. Data points, The sequence in the satellite magnetic field data is The point in time, For the station magnetic field data sequence is The point in time, This represents the starting time point of the truncated station magnetic field data. This represents the end time of the truncated station magnetic field data.

[0046] Time series of satellite magnetic field data extraction Obtain the satellite's start time point and satellite end time Satellite start time and satellite end time These two time points constitute the satellite time interval covered by satellite data. .

[0047] Time series of station magnetic field data In the middle, find the time interval with the satellite. Corresponding start and end indices. Due to the high sampling frequency of station electromagnetic data, its time points usually cannot be precisely matched one-to-one with satellite time points. Therefore, the nearest neighbor matching principle is adopted: the station time series with a time point greater than or equal to the satellite start time point is selected as the truncation start point, and the time series with a time point less than or equal to the satellite end time point is selected as the truncation end point.

[0048] Based on the determined cutoff start and end points, the original station magnetic field data were analyzed. The data segments within the corresponding intervals are extracted to form the truncated station magnetic field data. These are data points of the truncated station magnetic field data. The time points for the truncated station magnetic field data. This represents the starting time point of the truncated station magnetic field data. This represents the end time of the truncated station magnetic field data.

[0049] In S103, time-frequency transformations are performed on the satellite magnetic field data and ground station magnetic field data in the synchronized time series after time alignment and truncation processing to obtain the energy distribution representation of the two types of data in the time-frequency domain. This application employs the Synchrosqueezing Wavelet Transform (SWT) method, which can significantly improve time-frequency clustering while maintaining reversibility, thereby more accurately identifying transient seismic anomaly features in the magnetic field signal.

[0050] Time-frequency transformations are performed on satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectra. This includes performing synchronous compressed wavelet transforms on satellite magnetic field data and station magnetic field data in the synchronous time series. The instantaneous frequency at each moment is obtained from the wavelet coefficients obtained by wavelet transform; Based on a series of pre-set discrete compression target frequencies, for each calculated instantaneous frequency, it is determined which specified frequency interval it falls within near the compression target frequency, and all wavelet coefficient energy within the specified frequency interval is compressed to the corresponding compression target frequency.

[0051] The formula is expressed as follows: , , , in, This is the compressed time-frequency representation. For frequency intervals, These are wavelet coefficients. As a scale factor, The translation factor is... For the mother wavelet function, Selected scale factor, To compress the target frequency, The imaginary unit, This refers to the instantaneous frequency.

[0052] In step S104, after synchronous compressed wavelet transform, the time-frequency amplitude spectra of satellite magnetic field data and ground station magnetic field data are obtained respectively. Each of the two time-frequency amplitude spectra independently reflects the energy distribution characteristics of a single data source in the time-frequency domain. However, earthquake-related magnetic field anomalies usually exhibit common anomaly characteristics in both satellite and ground station magnetic field data, making it difficult for a single data source to effectively distinguish earthquake precursor anomalies from other interference signals.

[0053] To integrate the complementary information from the two types of data, this embodiment superimposes the two time-frequency amplitude spectra along the third dimension to construct a third-order non-negative tensor. This third-order non-negative tensor simultaneously incorporates the spatial distribution characteristics of the satellite data and the temporal evolution characteristics of the station data.

[0054] The two obtained after synchronous compressed wavelet transform The time-frequency amplitude spectra of the magnitudes are superimposed along a direction perpendicular to their third dimension to obtain... A third-order nonnegative tensor of size is called a third-order time-spectrum tensor.

[0055] In step S105, the third-order nonnegative tensor is decomposed into a nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix, including: The third-order nonnegative tensor is decomposed into the sum of R rank tensors using the tensor decomposition method; Based on KL divergence to measure the difference between a third-order nonnegative tensor and an approximate tensor, an objective function is established, in which each element of the third-order nonnegative tensor and the corresponding element of the approximate third-order tensor participate in the KL divergence calculation, and the approximate third-order tensor is obtained by summing R rank tensors. The third-order nonnegative tensor is expanded along three different moduli directions using n-modulus expansions to obtain the matrix representation of the original tensor. The approximate third-order tensor is also expanded using n-modulus expansions and is represented as the transpose of the Khatri-Rao product of the decomposition matrix and the remaining factor matrices. Let the intermediate variable be equal to the transpose of the Khatri-Rao product of the remaining factor matrices, and rewrite the objective function as a function of the decomposition matrix of the current modulus direction; The decomposition matrices of each modulus direction are optimized and updated sequentially using an alternating iterative method to obtain the basis matrix, weight matrix, and contribution matrix.

[0056] The formula is described in detail as follows: Tensor decomposition is used to decompose a third-order nonnegative tensor into the sum of R rank-one tensors. First, the objective function is established by measuring the difference between the original third-order nonnegative tensor and the sum of the R rank-one tensors using the KL divergence. , In the formula, For a third-order nonnegative tensor The One element, The approximate third-order tensor obtained by summing R rank tensors. It is the first of the approximate third-order tensors Each element.

[0057] conduct The module is expanded to form a third-order nonnegative tensor. It can be represented as Approximate third-order tensor It can be expressed as the following formula: , , in, For the first The factor matrix of each modulus, For the first The factor matrix of each modulus, For the first The factor matrix of each modulus, This indicates transpose.

[0058] make The objective function can be further written as: , , in, Let be the objective function. This is an intermediate variable matrix, representing the matrix excluding... Khatri-Rao product of all factor matrices outside the matrix. Subsequently, an alternating iterative method was used to... Optimize, The basic iterative formula is: ; ; In the formula, % represents element division. For the first The factor matrix at the next iteration For the first The factor matrix at the next iteration This indicates updating the factor matrix. and They have the same meaning, representing an intermediate variable matrix.

[0059] To obtain the decomposition matrix, the Karush-Kuhn-Tucker criterion is used after each iteration to determine whether the result has converged. The convergence condition is as follows: , In the formula, Indicates the convergence tolerance. Is with identity matrices of the same size For the first The update factor of each factor matrix. This occurs when the result satisfies the convergence condition or the maximum number of iterations is reached. When the iteration stops, the result is obtained. , , This is a decomposition matrix. , , Representing the basis matrix, weight matrix, and contribution matrix respectively, using... , , express The iterative algorithm steps are as follows: (1) Initialization , , ; (2) Utilization The basic iterative formula updates the objective function variable, iterating continuously until the convergence condition is met or the number of iterations reaches the upper limit. The iteration is terminated, and the decomposition result is obtained.

[0060] The basis matrix, weight matrix, and contribution matrix respectively characterize the main features of the third-order time-spectrum tensor of the magnetic field data in the dimensions of frequency, time (geographical location), and satellite observation channel.

[0061] In step S106, after nonnegative tensor decomposition, the original third-order time-spectrum tensor is decomposed into the sum of R rank tensors, each component corresponding to a set of characteristic combinations. However, not all components are related to seismic activity. Magnetic field anomalies generated during earthquake gestation often appear simultaneously in satellite data and ground station data, exhibiting common satellite-ground anomaly characteristics.

[0062] The purpose of step S106 is to select the component that best reflects the common characteristics of the satellite and the ground from the R components based on the difference in contribution between the satellite magnetic field data and the station magnetic field data in the contribution matrix, and use it as the target component (called the earthquake correlation component) for subsequent seismic anomaly extraction.

[0063] The contribution matrix reflects the contribution of satellite magnetic field data and station magnetic field data to the feature. Its size is 2×R, and each column corresponds to the contribution rate of satellite and station in each of the R components obtained by decomposition.

[0064] Based on the difference in contribution rates between satellite magnetic field data and station magnetic field data in the contribution matrix, earthquake-related components are selected, including: calculating the absolute value of the difference between the contribution rate of satellite magnetic field data and the contribution rate of station magnetic field data in each component of the contribution matrix. This absolute value reflects the degree of consistency between the contributions of the two data sources in the corresponding component. The formula is as follows: , It is the absolute value of the difference. Contribution rate to satellite magnetic field data Contribution rate to station magnetic field data; Arrange all R components in ascending order of their absolute values; The component ranked first is selected as the earthquake-related component.

[0065] The above screening criteria are based on the physical understanding that earthquake anomalies should appear simultaneously in both satellite magnetic field data and station magnetic field data. This can effectively eliminate unilateral noise interference and achieve automatic identification of earthquake-related features.

[0066] Based on the seismic correlation components obtained in step S106, in step S107, the anomaly extraction range is dynamically determined and anomaly points are extracted from the weight matrix. After filtering out the seismic correlation components through the contribution matrix, the effectiveness and reliability of the seismic correlation components need to be further verified. Not all seismic correlation components with the smallest absolute value truly reflect the common satellite-ground anomalies related to earthquakes.

[0067] Therefore, after screening out the earthquake-related components, this application further introduces a preset threshold judgment mechanism to verify the effectiveness of the earthquake-related components and dynamically sets the research scope based on the verification results.

[0068] First, calculate the absolute value of the difference between the contribution rate of satellite magnetic field data and the contribution rate of station magnetic field data corresponding to the seismic correlation component in the contribution matrix. Set a preset threshold, which can usually be adjusted according to the actual data characteristics. Then compare the calculated absolute value with the preset threshold.

[0069] If the absolute value is greater than the preset threshold, the earthquake-related component is determined to be an invalid earthquake component, indicating that the contribution rate of satellite magnetic field data and station magnetic field data to the earthquake-related component is too different. The earthquake-related component may mainly reflect noise or interference from a single data source, rather than a common geophysical process between satellite and ground. Therefore, the processing is terminated and it will not enter the subsequent anomaly extraction process.

[0070] If the absolute value is less than or equal to the preset threshold, the earthquake-related component is determined to be a valid earthquake component, and the process proceeds to the second level of judgment.

[0071] The second layer is determining the primary source of contribution: For earthquake-related components deemed valid, their primary contributing data sources are further determined. The contribution rates of satellite magnetic field data and station magnetic field data are compared.

[0072] If the contribution rate of station magnetic field data is greater than that of satellite magnetic field data, then the main contributing data source is determined to be the station magnetic field data. This indicates that the seismic correlation component primarily reflects the temporal evolution characteristics at a fixed station.

[0073] If the contribution rate of satellite magnetic field data is greater than that of station magnetic field data, then the main contributing data source is determined to be satellite magnetic field data. This indicates that the seismic correlation component primarily reflects the spatial distribution characteristics of the satellite along its orbit.

[0074] If the two are equal, they can be treated as a joint contribution from satellite and ground, and are usually classified as contributions from station magnetic field data.

[0075] The research scope for anomaly extraction is dynamically set based on the determination of the main sources of contribution.

[0076] The first scenario: the magnetic field data from the station is the primary contributor; When the primary contributing data source is determined to be station magnetic field data, the research scope is set to the complete time interval of the station magnetic field data. Specifically, the research scope starts from the beginning time of the truncated station magnetic field data and ends at the end time of the truncated station magnetic field data, which corresponds to the last time point of the satellite magnetic field data.

[0077] In this scenario, all time points corresponding to the earthquake-related components in the weight matrix are included in the anomaly detection range. The physical basis for this setting is that the station's magnetic field data are continuously observed from fixed locations, and earthquake anomalies may occur at any time before or after an earthquake; therefore, it is necessary to search for anomalies within a complete time interval.

[0078] The second scenario: satellite magnetic field data contributes the main contribution; When the primary contributing data source is determined to be satellite magnetic field data, the research scope is defined as the corresponding time interval during which the satellite passes through the seismic research area. Specifically, the research scope begins from the moment the satellite enters the seismic research area and ends at the moment the satellite leaves the seismic research area.

[0079] The spatial extent of the earthquake study area is determined using the Dobrovolsky formula for the radius of the seismogenic zone. The formula is: , To study the magnitude of earthquakes, The unit is kilometers. The circular area centered on the future epicenter and with the calculated radius of the earthquake-prone zone as its radius is defined as the earthquake research area.

[0080] Based on satellite orbit data, which includes the latitude and longitude information of each sampling point, the times when the satellite enters and leaves the earthquake study area are determined. The entry time is the moment when the satellite first meets the condition of being less than or equal to the radius of the seismogenic zone from the epicenter, and the departure time is the moment when the satellite last meets the condition of being less than or equal to the radius of the seismogenic zone from the epicenter.

[0081] In this scenario, only the data points corresponding to the earthquake-related components in the weighting matrix that fall within the time interval of the satellite's transit study area are included in the anomaly detection range. The physical basis for this setting is that as the satellite flies along its orbit, earthquake anomalies mainly occur during the transit times corresponding to the area near the epicenter. Limiting the study range can effectively exclude irrelevant data far from the epicenter.

[0082] Extracting outliers from the weight matrix includes: Calculate the root mean square of all data for the seismic component in the corresponding weight matrix within the study area: , in, Represents the weight matrix of the first element. Each earthquake component To study the root mean square of the data within the region, The total number of data points within the research scope; Multiply the preset empirical parameter by the root mean square value, and set the product as the anomaly extraction threshold. The formula is as follows: ,in, These are empirical parameters; When the value of a data point within the study area corresponding to the seismic component in the weight matrix exceeds the anomaly extraction threshold, the data point is determined to be a seismic anomaly. Simultaneously, the corresponding orbital and station magnetic field data are marked as seismic-related data, and anomalies are accumulated.

[0083] On the other hand, this application provides a multi-sphere seismic electromagnetic anomaly fusion and extraction system. This system, through the collaborative work of multiple functional modules, achieves joint processing of satellite magnetic field data and ground station magnetic field data, automatically identifying magnetic field anomaly signals related to seismic activity. It includes: a preprocessing module, used to preprocess both satellite and ground station magnetic field data, removing internal source fields and noise to obtain preprocessed satellite-Ground magnetic field data; for satellite magnetic field data, the preprocessing module first calculates and subtracts internal source field components such as the core field and crustal field using an international reference geomagnetic field model or a high-precision geomagnetic field model. Commonly used geomagnetic field models include the CHAOS model and the IGRF model. After subtracting internal source fields, the preprocessing module further filters out large-scale external source field interference generated by the magnetosphere current system and the ionosphere current system. Finally, the preprocessing module removes anomalous jumps and high-frequency noise caused by solar high-energy particle events, satellite platform attitude disturbances, and instrument noise. After the above processing, the preprocessed satellite magnetic field data mainly retains magnetic field anomaly information related to medium- and small-scale geological activities. For ground-based station magnetic field data, the preprocessing module also uses a geomagnetic field model to subtract internal source field components. Considering the unique characteristics of the station's observation environment, the preprocessing module also needs to correct for instrument baseline drift and temperature effects, remove pulse interference and background noise caused by thunderstorms, human activities (such as rail traffic and industrial equipment), and appropriately interpolate or label short-term missing data. The preprocessed station magnetic field data reflects the continuous change of the magnetic field at a fixed geographical location over time.

[0084] The time alignment module is used to perform time alignment and truncation on the preprocessed satellite-Ground magnetic field data to form a synchronized time series. Because satellites fly along their orbits with relatively large data sampling intervals, while station magnetic field data is a continuous high-sampling-rate record, the time lengths and number of sampling points of the two are usually inconsistent. The time alignment module uses the time interval of the satellite magnetic field data as a reference and extracts a data segment from the station magnetic field data that completely corresponds to the satellite's time interval, ensuring strict time alignment between the two sets of data.

[0085] Let the preprocessed satellite magnetic field data be the Y-component satellite magnetic field data, with a length of N, and the corresponding time series extending from the start time to the end time. Let the preprocessed station magnetic field data have a length of M, and the corresponding time series extending from the start time to the end time. Typically, stations conduct continuous high-sampling-rate observations, and their data length M is much larger than the satellite data length N.

[0086] The alignment and truncation processing of the time alignment module includes the following operations. First, the time boundaries of the satellite magnetic field data are extracted to obtain the start and end time points, which constitute the complete time interval covered by the satellite data. Then, in the time series of the station magnetic field data, the start and end indices corresponding to the satellite time interval are found. Since the sampling frequency of station data is high, its time points usually cannot correspond precisely one-to-one with the satellite time points. Therefore, the nearest neighbor matching principle is adopted, selecting the first time point in the station time series that is greater than or equal to the satellite start time as the truncation start point, and selecting the last time point that is less than or equal to the satellite end time as the truncation end point. Finally, based on the determined start and end indices, the data segments within the corresponding intervals are extracted from the original station magnetic field data to form the truncated station magnetic field data.

[0087] The output of the time alignment module is a synchronized time series, in which the satellite magnetic field data and the truncated station magnetic field data are strictly aligned in time, providing aligned input data for subsequent time-frequency transformation. The time-frequency transformation module performs time-frequency transformation on the satellite magnetic field data and station magnetic field data in the synchronized time series, respectively, to obtain the corresponding time-frequency amplitude spectra. The time-frequency transformation module uses synchronous compressed wavelet transform as the time-frequency transformation method. Synchronous compressed wavelet transform is a post-processing enhanced time-frequency analysis method that can significantly improve time-frequency clustering while maintaining reversibility, thereby more accurately identifying transient seismic anomaly features in the magnetic field signal.

[0088] The time-frequency transformation module first performs a continuous wavelet transform on the magnetic field data (satellite magnetic field data and station magnetic field data) to obtain complex wavelet coefficients that depend on the scaling factor and translation factor. Specifically, for a given magnetic field signal, the time-frequency transformation module selects a mother wavelet function and performs a continuous wavelet transform on the signal to obtain a wavelet coefficient matrix.

[0089] Then, the time-frequency transformation module calculates the instantaneous frequency at each moment based on the wavelet coefficients. For each point on the wavelet coefficient plane, the instantaneous frequency is calculated by the ratio of the partial derivative of the wavelet coefficient with respect to the translation factor to the wavelet coefficient itself. This calculation essentially tracks the phase change rate of the wavelet coefficient.

[0090] Finally, the time-frequency transformation module performs energy compression and redistribution. The module pre-sets a series of discrete compression target frequencies, which constitute the compressed frequency axis. For each calculated instantaneous frequency, the module determines which specified frequency interval it falls within near the compression target frequency and compresses all wavelet coefficient energy within that specified frequency interval to the corresponding compression target frequency. After synchronous compression, the original wavelet coefficients are converted into a two-dimensional array defined on the time-frequency plane, i.e., the time-frequency amplitude spectrum.

[0091] The time-frequency transformation module performs the above transformation on both the satellite magnetic field data and the station magnetic field data, resulting in two two-dimensional time-frequency amplitude spectra of the same size. Each time-frequency amplitude spectrum has a size of I1 multiplied by I2, where I1 corresponds to the number of points in the frequency dimension and I2 corresponds to the number of points in the time dimension. The two time-frequency amplitude spectra have identical sizes, which is ensured by using the same transformation parameters, the same mother wavelet function, the same scale factor range, and the same set of compressed target frequencies.

[0092] The tensor construction module is used to superimpose two time-frequency amplitude spectra along the third dimension to construct a third-order non-negative tensor. The module uses the time-frequency amplitude spectrum of satellite magnetic field data as the first channel of the third-order non-negative tensor and the time-frequency amplitude spectrum of station magnetic field data as the second channel, stacking them along the third dimension perpendicular to the time-frequency plane. Specifically, for any frequency and time index, the third-order non-negative tensor element takes the corresponding element of the satellite time-frequency amplitude spectrum when the third-dimensional index is one, and takes the corresponding element of the station time-frequency amplitude spectrum when the third-dimensional index is two.

[0093] The size of the constructed third-order nonnegative tensor is I1 multiplied by I2 multiplied by 2, where I1 is the number of points in the frequency dimension, I2 is the number of points in the time dimension, and the third dimension is 2, corresponding to the satellite magnetic field data channel and the station magnetic field data channel, respectively.

[0094] This third-order nonnegative tensor possesses nonnegativity, meaning all elements are greater than or equal to zero. This property stems from the nonnegativity of the time-frequency amplitude spectrum and the fact that no negative value operations were introduced during tensor construction. Nonnegativity ensures that subsequent nonnegative tensor decomposition algorithms can be applied correctly, while also guaranteeing the physical interpretability of the decomposition results.

[0095] The tensor decomposition module is used to perform nonnegative tensor decomposition on the third-order nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix. This application adopts a nonnegative tensor decomposition method based on KL divergence, which approximates the original third-order nonnegative tensor as the sum of R rank tensors, corresponding to R components. During the decomposition process, the decomposition matrix is ​​updated along the three modulus directions through an alternating iterative optimization algorithm.

[0096] Specifically, the tensor decomposition module first randomly initializes three decomposition matrices, denoted as the base matrix, weight matrix, and contribution matrix, respectively. Then, the tensor decomposition module iteratively optimizes the three decomposition matrices according to an alternating update rule. In each iteration, the tensor decomposition module fixes the other two decomposition matrices and updates the current decomposition matrix using a multiplicative update rule. After each iteration, the tensor decomposition module uses the Karush-Kuhn-Tucker criterion to determine whether the result has converged. Iteration stops when the convergence condition is met or the preset maximum number of iterations is reached.

[0097] After decomposition, the tensor decomposition module outputs three matrices: the first matrix is ​​called the basis matrix, denoted as W, whose size is related to the frequency dimension and is used to characterize the basic features of the magnetic field signal at different frequencies; the second matrix is ​​called the weight matrix, denoted as H, whose size is related to the time dimension and is used to characterize the intensity of the anomalous signal change over time; the third matrix is ​​called the contribution matrix, denoted as C, whose size is 2 times R, and each column corresponds to the contribution rate of satellite data and the contribution rate of station data in a decomposed component, respectively, where R represents the number of rank tensors.

[0098] The component filtering module is used to filter out earthquake-related components based on the difference in contribution rates between satellite magnetic field data and station magnetic field data in the contribution matrix. The module first calculates the absolute value of the difference between the satellite contribution rate and the station contribution rate for each component in the contribution matrix. The module then arranges all R components in ascending order of their absolute values, resulting in a sorted component sequence.

[0099] The component filtering module selects the component that is ranked first, i.e. the component with the smallest absolute value, as the earthquake-related component.

[0100] The anomaly extraction module dynamically determines the anomaly extraction range based on the seismic correlation components and extracts outliers from the weight matrix. The anomaly extraction module first determines the validity of the seismic correlation components. It calculates the absolute value of the difference in contribution rates between the seismic correlation components and compares this absolute value with a preset threshold. If the absolute value is greater than the preset threshold, it indicates that the contribution rate difference between satellite magnetic field data and station magnetic field data on this component is too large, and the component is determined to be an invalid seismic component, terminating the processing. If the absolute difference is less than or equal to the preset threshold, the component is determined to be a valid seismic component, and processing continues.

[0101] For valid seismic components, the anomaly extraction module further determines the main contributing data source. If the station contribution rate is greater than the satellite contribution rate, the main contributing source is determined to be station magnetic field data; if the satellite magnetic field data contribution rate is greater than the station magnetic field data contribution rate, the main contributing source is determined to be satellite magnetic field data.

[0102] Based on the determination of the primary source of contribution, the anomaly extraction module dynamically sets the anomaly extraction range. When the primary source of contribution is station magnetic field data, the research range is set to the complete time interval of the station magnetic field data, i.e., from the start time to the end time after the station data is truncated. When the primary source of contribution is satellite magnetic field data, the research range is set to the corresponding time interval within the earthquake research area traversed by the satellite.

[0103] After defining the study scope, the anomaly extraction module extracts the data sequences corresponding to the earthquake-related components within the study scope from the weight matrix. The anomaly extraction module then calculates the root mean square (RMS) value of this data sequence, which is obtained by summing the squares of all data points, dividing by the total number of data points, and then taking the square root. Finally, a preset empirical parameter is multiplied by this RMS value to obtain the anomaly extraction threshold.

[0104] The anomaly extraction module iterates through each data point within the study area, comparing the value of each data point with an anomaly extraction threshold. When the value of a data point exceeds the threshold, it is determined to be an earthquake anomaly. Simultaneously, the corresponding satellite orbit data and station data are marked as earthquake-related data, and the detected anomalies are accumulated and statistically analyzed.

[0105] The output of the anomaly extraction module includes the location information of seismic anomalies, anomaly amplitude, cumulative number of anomalies, and marked satellite-ground data, providing quantitative data support for earthquake monitoring and earthquake precursor research.

[0106] In one embodiment, taking a magnitude 7.3 earthquake in a certain area on May 22, 2021 as an example, Swarm A satellite data is used. This data contains all orbital data for a single day, with a sampling rate of 1Hz. It needs to be segmented with 50°S and 50°N as orbital endpoints, dividing the entire day's data into 32 satellite orbits. The radius of the study area is calculated using the formula studied by Dobrovolsky: Converted to latitude and longitude, this is approximately 12.4 degrees; therefore, the research scope is... , The orbits passing through the earthquake research area were selected, and the magnetic field-latitude signal was converted into a magnetic field-time series. The magnetic field data of the stations on the same day were read. The sampling rate of this data was 1Hz, and it included all 24 hours of data within a day.

[0107] Taking orbit 23 on April 22, 2021 as an example, the satellite magnetic field data for each satellite orbit is subtracted from the source field within the CHAOS-7 model; for example... Figure 2 As shown, where Figure 2 (a) in the image represents satellite magnetic field data. Figure 2 (b) in the figure shows the data after removing the stable background. It can be seen that the long-period trend in the original satellite magnetic field data has been effectively removed, and the data curve fluctuates around the zero line. Figure 2 In (b), obvious short-cycle fluctuations and local abnormal peaks can be observed.

[0108] Taking the station data from April 22, 2021 as an example, with a decomposition level L=5, the result after wavelet decomposition and reconstruction is as follows: Figure 3 As shown, where Figure 3 (a) in the figure represents the magnetic field data of the station. Figure 3 (b) shows the station magnetic field data after wavelet decomposition to remove low-frequency trends and high-frequency noise. It can be seen that wavelet decomposition effectively separates different frequency components in the station magnetic field signal. Background noise and interference in the original station magnetic field data are effectively suppressed, and mid-frequency anomalous signals related to geological activity become more prominent.

[0109] Taking data from April 22, 2021 as an example, the length of the station's magnetic field data on that day was M=86400, and the length of the satellite's magnetic field data was N=1558. After alignment and truncation, the result is as follows: Figure 4 The two magnetic field data with the same time length are shown, among which Figure 4 (a) in the image represents satellite magnetic field data. Figure 4 (b) in the figure represents the magnetic field data of the station, and the time lengths of the two are the same.

[0110] Taking the magnetic field data from April 22, 2021 as an example, the synchronous compressed wavelet transform result is as follows: Figure 5 As shown, where Figure 5 (a) in the figure represents the time-frequency amplitude spectrum of the satellite magnetic field data. Figure 5 (b) shows the time-frequency amplitude spectrum of the station's magnetic field data. The time-frequency amplitude spectrum of the satellite magnetic field data shows anomalies around 16:58 and 17:08. The anomalous energy around 16:58 is mainly concentrated between 0.01Hz and 0.04Hz, while the anomalous energy around 17:08 is mainly concentrated between 0.02Hz and 0.09Hz, which are relatively high frequencies. The time-frequency amplitude spectrum of the station's magnetic field data shows that the geomagnetic field of the lithosphere is disturbed around 16:59 and 17:12. The anomalous energy around 16:59 is mainly concentrated between 0.01Hz and 0.05Hz, while the anomalous energy around 17:12 is mainly concentrated between 0.01Hz and 0.06Hz. There is a certain quasi-synchronous or lagging correlation between the period of energy enhancement in the ionospheric time-frequency amplitude spectrum and the anomalies in the lithosphere.

[0111] The two obtained after synchronous compressed wavelet transform The time-frequency amplitude spectra of the magnitudes are superimposed along a direction perpendicular to their third dimension to obtain... A third-order nonnegative tensor X of size X is called a third-order time-spectrum tensor. Taking April 22, 2021 as an example, the two values ​​obtained by synchronous compressed wavelet transform... The time-frequency amplitude spectra of the magnitudes are superimposed. The third-order nonnegative tensor is as follows Figure 6 As shown.

[0112] Figure 7 The example image shown is of the earthquake-related component anomaly extraction results. The data points above the dashed line are those exceeding the anomaly extraction threshold and are marked as anomalies.

[0113] 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 fusing and extracting seismic electromagnetic anomalies across multiple spheres of space and Earth, characterized in that, include: The satellite magnetic field data and the station magnetic field data were preprocessed to remove internal source fields and noise, resulting in preprocessed satellite-to-ground magnetic field data. The preprocessed satellite-Ground magnetic field data is time-aligned and truncated to form a synchronized time series; Time-frequency transformations were performed on the satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectra. The two time-frequency amplitude spectra are superimposed along the third dimension to construct a third-order non-negative tensor; The third-order nonnegative tensor is decomposed into a nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix. Based on the difference in contribution between satellite magnetic field data and station magnetic field data in the contribution matrix, earthquake-related components are selected. Based on the earthquake correlation components, the anomaly extraction range is dynamically determined and anomaly points are extracted from the weight matrix.

2. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 1, characterized in that, Preprocessing was performed on both satellite magnetic field data and station magnetic field data to remove intrinsic fields and noise, resulting in preprocessed satellite-to-Earth magnetic field data, including: The db4 wavelet is selected as the wavelet basis function, and the total number of decomposition levels L is preset; the magnetic field signal is subjected to L-level multi-scale wavelet decomposition, and the magnetic field signal includes satellite magnetic field data and station magnetic field data; In the first-level decomposition, the magnetic field signal is convolved using the low-pass filter coefficients and the high-pass filter coefficients respectively, and then downsampled to obtain the approximate components and detail components of the first level. From the second layer to the Lth layer, each layer repeats the convolution operation and downsampling above on the approximate components obtained from the previous layer, and obtains the approximate components and detail components of the current layer in turn. After decomposition, all effective components are retained and inverse wavelet transform is performed: each layer of components is upsampled layer by layer, convolved with the reconstruction filter respectively, and then the convolution results are superimposed to obtain the preprocessed star-ground magnetic field data.

3. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 1, characterized in that, The preprocessed satellite-to-ground magnetic field data is time-aligned and truncated to form a synchronized time series, including: selecting station magnetic field data with the same time point as the satellite magnetic field data based on the time information of the satellite magnetic field data; and aligning and truncating the station magnetic field data to obtain two electromagnetic signals with the same time length. In this context, the length of the satellite Y-component magnetic field data and its corresponding time series in the satellite magnetic field data is denoted as N; the length of the station magnetic field data and its corresponding time series is denoted as M, and N is less than M; the truncated station magnetic field data is denoted as a new time series, which is taken from the data segment in the station magnetic field data before truncation, within the time interval corresponding to the first satellite time point to the Nth satellite time point.

4. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 1, characterized in that, Time-frequency transformations were performed on the satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectra, including: Synchronous compressed wavelet transforms were performed on satellite magnetic field data and station magnetic field data in the synchronous time series, respectively. The instantaneous frequency at each moment is obtained from the wavelet coefficients obtained by wavelet transform; Based on a series of pre-set discrete compression target frequencies, for each calculated instantaneous frequency, it is determined which specified frequency interval the instantaneous frequency falls within, and all wavelet coefficient energy within the specified frequency interval is compressed to the corresponding compression target frequency.

5. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 1, characterized in that, The third-order nonnegative tensor is decomposed into a nonnegative tensor to obtain the basis matrix, weight matrix, and contribution matrix, including: The third-order nonnegative tensor is decomposed into the sum of R rank tensors using the tensor decomposition method, resulting in R components. Based on KL divergence to measure the difference between a third-order nonnegative tensor and an approximate tensor, an objective function is established, in which each element of the third-order nonnegative tensor and the corresponding element of the approximate third-order tensor participate in the KL divergence calculation, and the approximate third-order tensor is obtained by summing R rank tensors. The third-order nonnegative tensor is expanded along three different moduli directions using n-modulus expansions to obtain the matrix representation of the original tensor. The approximate third-order tensor is also expanded using n-modulus expansions and is represented as the transpose of the Khatri-Rao product of the decomposition matrix and the remaining factor matrices. Let the intermediate variable be equal to the transpose of the Khatri-Rao product of the remaining factor matrices, and rewrite the objective function as a function of the decomposition matrix of the current modulus direction; The decomposition matrices of each modulus direction are optimized and updated sequentially using an alternating iterative method to obtain the basis matrix, weight matrix, and contribution matrix.

6. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 1, characterized in that, Based on the difference in contribution rates between satellite magnetic field data and station magnetic field data in the contribution matrix, earthquake-related components are selected, including: calculating the absolute value of the difference between the contribution rate of satellite magnetic field data and the contribution rate of station magnetic field data in each component of the contribution matrix, wherein the absolute value reflects the degree of consistency between the contributions of the two data sources in the corresponding component. Arrange all R components in ascending order of their absolute values; The component ranked first is selected as the earthquake-related component.

7. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 6, characterized in that, Based on the earthquake correlation components, the anomaly extraction range is dynamically determined, including: Determine whether the absolute value of the difference between the satellite contribution rate and the station contribution rate in the contribution matrix corresponding to the earthquake correlation component is greater than a preset threshold; If the absolute value is greater than the preset threshold, it will not be considered a valid seismic component. If the absolute value is less than or equal to a preset threshold, determine the main contributing data source of the earthquake-related component. If the main contributing data source is the station's magnetic field data, the research scope will be set to the complete time interval of the station's magnetic field data. If the primary source of contributing data is satellite magnetic field data, the research scope will be set to the corresponding time interval within the earthquake research area traversed by the satellite.

8. The method for fusing and extracting multi-sphere seismic electromagnetic anomalies according to claim 7, characterized in that, Extracting outliers from the weight matrix includes: Calculate the root mean square of all data for the seismic component in the corresponding weight matrix within the study area; Multiply the preset empirical parameter by the root mean square value, and set the product result as the anomaly extraction threshold; When the value of a data point within the study area corresponding to the seismic component in the weight matrix is ​​greater than the anomaly extraction threshold, the data point is determined to be a seismic anomaly.

9. A multi-sphere seismic electromagnetic anomaly fusion extraction system, characterized in that, include: The preprocessing module is used to preprocess satellite magnetic field data and station magnetic field data respectively, remove internal source fields and noise, and obtain preprocessed satellite-to-ground magnetic field data; The time alignment module is used to perform time alignment and truncation on the preprocessed satellite-Ground magnetic field data to form a synchronized time series. The time-frequency conversion module is used to perform time-frequency conversion on satellite magnetic field data and station magnetic field data in the synchronous time series to obtain the corresponding time-frequency amplitude spectrum. The tensor construction module is used to superimpose two time-frequency amplitude spectra along the third dimension to construct a third-order non-negative tensor. The tensor decomposition module is used to perform non-negative tensor decomposition on the third-order non-negative tensor to obtain the basis matrix, weight matrix and contribution matrix. The component filtering module is used to filter out earthquake-related components based on the difference in contribution between satellite magnetic field data and station magnetic field data in the contribution matrix. The anomaly extraction module is used to dynamically determine the anomaly extraction range based on the earthquake correlation components and extract anomaly points from the weight matrix.

Images