CN119414452B - A method for extracting seismic anomalies from satellite magnetic fields based on complex non-negative matrix decomposition - Google Patents

Abstract

The present invention belongs to the field of ionospheric satellite magnetic field seismic anomaly extraction. Specifically, it is a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition. The method includes reading the magnetic field data of the Swarm A satellite, removing the main magnetic field of the Earth's core and the lithospheric magnetic field of the Earth's crust to obtain the residual magnetic field; performing a short-time Fourier transform on the preprocessed data to obtain a complex time-frequency matrix in the form of rectangular coordinates in the complex domain; decomposing the complex time-frequency matrix using complex non-negative matrix decomposition to obtain R characteristic components; performing an inverse short-time Fourier transform on each characteristic component to obtain a time-domain characteristic component; and calculating the root mean square of the time-domain seismic component of each orbit, using the exceedance rate method to set a threshold, and extracting anomalies from the orbital seismic components. This method fully utilizes the amplitude and phase information contained in the satellite observation data in the time-frequency domain to identify components related to earthquakes, making seismic anomaly extraction more accurate and reliable.

Description

Satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization

Technical Field

The invention belongs to the field of ionosphere satellite magnetic field seismic anomaly extraction, and particularly relates to a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization.

Background

Seismology utilizes non-mechanical electromagnetic methods to study the processes of inoculation and occurrence of earthquakes, and becomes one of the effective ways to deepen the understanding of the earthquake inoculation process. With the rapid development of space detection technology, ionosphere detection technology, especially satellite detection, is increasingly applied to seismic precursor research, and satellite electromagnetic seismic anomaly research becomes a great hotspot of international research.

The satellite magnetic field observation data are not only interfered by the earth, but also are strongly interfered by non-seismic factors such as solar activity, magnetic storm and the like, and the traditional research method selects night data for avoiding the non-seismic interference and simultaneously removes data with higher solar activity index and geomagnetic index. In order to improve the utilization rate of seismic data, seismic anomalies are extracted from satellite magnetic field observation data with the characteristics of strong background and weak information, and seismic signals, background and strong geomagnetic activity interference in the data are separated by adopting blind source separation modes such as nonnegative matrix factorization. However, the method has certain limitation at present that the magnetic field signal not only contains waveform and amplitude information, but also contains other important information such as phase and the like. In the method, the amplitude spectrum matrix is decomposed after the time-frequency analysis is carried out on the magnetic field signals, all the obtained components are combined with the original phase matrix to reconstruct into a time domain, the influence of the original magnetic field signal phase information on the abnormal extraction during reconstruction is ignored, and complex frequency domain information cannot be fully utilized to improve the accuracy of the seismic magnetic field abnormal identification extraction.

Therefore, in order to effectively use the phase information, the current complex domain matrix decomposition method needs to be studied in depth, and satellite magnetic field seismic anomaly extraction is expanded from a real number domain to a complex number domain, so that the understanding of seismic anomaly phenomenon is deepened.

Disclosure of Invention

The invention aims to solve the technical problem of providing a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization, which fully utilizes amplitude and phase information contained in satellite observation data in a time-frequency domain, so as to identify components related to a seismic, and enable seismic anomaly extraction to be more accurate and reliable.

The invention discloses a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization, which comprises the following steps:

step 1, reading Swarm A satellite magnetic field data, and dividing the daily magnetic field data into a plurality of orbits;

Step 2, subtracting the CHASS-7 model endogenous field from the original Y component magnetic field data of each satellite orbit, removing the main magnetic field of the earth core and the rock ring magnetic field of the earth crust to obtain a residual magnetic field;

Step 3, carrying out short-time Fourier transform on the preprocessed data to obtain a complex time-frequency matrix in a rectangular coordinate representation form of a complex domain;

decomposing the complex time-frequency matrix by utilizing complex non-negative matrix decomposition to obtain R characteristic components, wherein each characteristic component comprises a base vector, a coefficient vector and a phase matrix;

step 5, obtaining time domain characteristic components from the characteristic components through inverse short-time Fourier transform;

Step 6, according to the characteristic that the seismic energy and the space-time information are concentrated in the seismic research area, calculating the energy-entropy ratio of each characteristic component by utilizing the coefficient vector of each characteristic component of the track and the time domain characteristic component, and selecting the characteristic component with the largest energy-entropy ratio as the seismic component according to the sequence from big to small;

and 7, calculating the root mean square of each track time domain seismic component, setting a threshold value by using an overrun method, and extracting abnormal points from the track seismic components.

Further, the step 1 includes:

The method comprises the steps of reading Swarm A satellite magnetic field data, dividing daily magnetic field data into a plurality of orbits, specifically, storing the magnetic field data in a daily unit, dividing the magnetic field data by using 50 DEG S and 50 DEG N as orbit endpoints, and obtaining 32 satellite orbits daily.

Further, the step 2 includes:

The raw Y-component magnetic field data for each satellite orbit minus the CHAOS-7 model endogenous field is calculated as:

Wherein B _y0 is the original Y component magnetic field, Endogenous field Y component data calculated for the CHASS-7 model, B _y being the residual Y component magnetic field;

the residual magnetic field is differentiated, and the calculation is as follows:

x(n)=B_y(n+1)-B_y(n) (2)

Where B _y (N) is the nth data point of the remaining Y component magnetic field, the data length of B _y is set to N, and x (N) is the differential Y component magnetic field.

Further, the step 3 includes:

short-time Fourier transform is carried out on the differential Y component magnetic field x (n), and the discrete short-time Fourier transform is specifically as follows:

Where m represents the time index of the input signal x (n), g (n) is the selected window function, g (n-m) is the sliding window, n determines the current sliding position, j represents the imaginary unit, ω is the frequency index, and V (n, ω) represents the complex time-frequency matrix obtained by the discrete short-time fourier transform.

Further, the step4 includes:

Decomposing the complex time-frequency matrix by complex non-negative matrix decomposition, specifically, representing the complex time-frequency matrix V (n, omega) obtained by short-time Fourier transformation as V _k×l of k rows and l columns, giving the number R of decomposition characteristic components, wherein R is a positive integer and satisfies R < < min (k, l), decomposing the complex time-frequency matrix V _k×l into a matrix W _k×R, a matrix H _R×l and a matrix Wherein W _k×R is referred to as a basis matrix, representing the frequency distribution characteristics of the data, H _R×l is referred to as a coefficient matrix, representing the weights of the different frequency distribution characteristics in time/space,Refers to a time-varying phase spectrum, comprising a phase matrix of k rows and l columns corresponding to R characteristic components, hereinafter abbreviated as V, W, H,The complex non-negative matrix factorization mathematical model is:

Wherein, the A Hadamard product representing two matrices, the elements of which are defined as the product of the corresponding elements of the two matrices, W _k×r representing the r-th column of the base matrix W, H _r×l representing the r-th row of the coefficient matrix H,Representing time-varying phase spectraR is more than or equal to 1 and less than or equal to R.

To approximate the decomposition result to complex matrix V as much as possible, KL divergence is used to measure the original and reconstructed matricesThe error between them, the objective function is:

Wherein, the As an objective function, X _k,r,l is an estimated value of the reconstructed complex time-frequency matrix of the R-th feature component, R (H _r,l) is a penalty term based on p-norm, and when 0< p <2, the sparsity of H _r,l can be effectively controlled. The penalty term is:

R(H_r,l)＝2λ|H_r,l|^p,λ>0 (5)

λ is the weight coefficient of the penalty term;

The formulas for the auxiliary functions Z _k,r,l,d_k,r,l,A_k,r,l and B _k,r,l are as follows:

Z_k,r,l＝|X_k,r,l| (6)

finally, the iterative updating rule when 0< p <1 is obtained is as follows:

Z_k,r,l＝|X_k,r,l| (14)

Y_k,r,l＝|X_k,r,l| (15)

U_r,l＝H_r,l (16)

the iterative algorithm steps are as follows:

Input complex time-frequency matrix V, number of decomposed eigenvectors R, upper limit of iteration times N _iter

(1) Initializing W, H and X _k,r,l;

(2) Updating the objective function variable by using the formulas (11) - (17), iterating continuously until the value of the objective function formula (5) converges or reaches the set iteration times N _iter, and stopping iterating;

the output is W, H,

R characteristic components are obtained through decomposition, wherein the R characteristic components comprise an R column vector W _r of a base matrix W, an R row vector H _r of a coefficient matrix H and a corresponding phase matrixH _r coefficient vector, W _r base vector.

Further, the step 5 includes:

the time domain characteristic components are obtained by carrying out inverse short-time Fourier transform on the characteristic components, and the time domain characteristic components are specifically:

ISTFT () represents a function of the inverse short-time Fourier transform.

Further, the step 6 includes:

The energy-entropy ratio of each characteristic component of the track consists of an energy ratio and an entropy ratio, wherein the energy ratio refers to the ratio of the energy to the whole track in the seismic investigation region of the r component of the track, and the entropy ratio refers to the ratio of shannon entropy to the whole track in the seismic investigation region of the time domain characteristic component.

Shannon entropy Z _i(y_r) in the seismic investigation region and shannon entropy Z (y _r) for the entire trajectory are calculated as follows:

wherein, the point between s ₁ and s ₂ is in the seismic investigation region, N is the length of the time domain feature component Y _r, p (Y _r) represents the probability that the random event Y is Y _r, and s ₁、s₂ is the minimum latitude and the maximum latitude of the seismic investigation region respectively.

The energy-entropy ratio is calculated as follows:

the points between Ps and Pe are in the seismic research area, L is the number of columns of H _r,l, and Ps and Pe are the minimum latitude and the maximum latitude of the corresponding seismic research area obtained by short-time Fourier transform respectively.

Sorting the characteristic components according to the energy-entropy ratio from large to small, respectively marking the basis vectors of the sorted characteristic components as W _s1、W_s2、...、W_sR, respectively marking the coefficient vectors as H _s1、H_s2、…、H_sR, respectively marking the phase matrixes asThe reconstructed time domain data are denoted as y _s1、y_s2、…、y_sR, respectively, and the feature component with the largest energy-entropy ratio is selected as the seismic component.

Further, the step 7 includes:

In substep ①, the root mean square of the time domain seismic components of each track is calculated as:

The substep ② threshold is set as:

Thre=P×RMS (22)

p is an empirical parameter, typically set to 3, and Thre is an anomaly extraction threshold

When the amplitude of the time domain seismic component data point in the seismic investigation region is greater than Thre, the point is considered to be a seismic outlier while the trajectory is marked as an outlier trajectory.

Compared with the prior art, the invention has the beneficial effects that:

The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization comprises the steps of removing an endogenous field through subtracting a CHAOS-7 model, performing first-order difference, removing a background field of satellite magnetic field data, performing time-frequency conversion on preprocessed data with the background field removed through short-time Fourier transform to obtain a complex time-frequency matrix, according to the characteristics of global magnetic field influence and seismic influence on a local area magnetic field caused by non-seismic factors such as solar activity and magnetic storm, obtaining a plurality of characteristic components through utilizing the advantages of complex non-negative matrix factorization local feature extraction, wherein each characteristic component comprises a base vector, a coefficient vector and a corresponding phase matrix, fully utilizing amplitude and phase information contained in a time-frequency domain of satellite magnetic field signals, obtaining a reconstructed time domain characteristic component through inverse short-time Fourier transform, identifying a seismic component through utilizing an energy-entropy ratio according to the characteristics of seismic energy and space-time information distribution, setting a threshold value through calculating root mean square of the time domain seismic component, and extracting seismic anomalies through adopting an overrun method. According to the characteristic that the abnormal distribution of the seismic magnetic field is concentrated in the seismic influence area, the abnormal characteristics generated by the earthquake and the non-earthquake factors are separated by utilizing the advantages of complex non-negative matrix factorization local characteristic extraction and the characteristic of fully utilizing complex domain amplitude information and phase information, so that the abnormal seismic magnetic field is extracted more accurately and reliably, and the recognition of the abnormal phenomenon of the ionosphere seismic magnetic field is deepened.

On the basis of having the advantage of extracting local features, the complex non-negative matrix decomposition method can approximately decompose the original complex matrix V into two non-negative matrices W E R ^k×r、H∈R^r×l and a phase spectrum thereofWhere W is called a basis matrix, representing the frequency distribution characteristics of the data, and H is called a coefficient matrix, representing the weights of the different frequency distribution characteristics in time/space. After removing the background of satellite magnetic field data, carrying out short-time Fourier transform to obtain a time spectrum matrix of a complex domain, decomposing the time spectrum matrix by adopting a complex non-negative matrix to obtain seismic components, reconstructing back to a time domain by utilizing inverse short-time Fourier transform, extracting seismic anomalies by an overrun method, and determining an anomaly orbit. The method fully utilizes the amplitude and phase information contained in the magnetic field signal, and increases the reliability of the seismic anomaly extraction result.

Drawings

FIG. 1 is a flow chart of a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization;

FIG. 2 is an iterative curve of the objective function error values of the complex non-negative matrix factorization algorithm;

FIG. 3 is a base matrix and coefficient matrix obtained by complex non-negative matrix decomposition;

FIG. 4 is a base matrix and coefficient matrix ordered from large to small in energy-to-entropy ratio;

fig. 5 is a view of the ordered reconstructed time domain feature components.

Detailed Description

The present invention will be described in further detail with reference to the following examples in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.

Referring to fig. 1, a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization includes:

Step1, reading Swarm A satellite magnetic field data, and dividing daily data into a plurality of orbits taking 50 DEG S-50 DEG N as endpoints;

and 2, selecting a Y component of magnetic field data as original magnetic field data, subtracting an endogenous field of the CHASS-7 model from the original magnetic field data of each track, thereby removing a main magnetic field and a rock ring magnetic field to obtain a residual magnetic field, and performing first-order difference on the residual magnetic field to obtain magnetic field pretreatment data in order to remove a low-frequency background, wherein the endogenous field of the CHASS-7 model belongs to the prior art.

Step 3, performing time-frequency transformation on the preprocessed data by adopting short-time Fourier transformation to obtain a complex time-frequency matrix;

decomposing the complex time-frequency matrix by utilizing complex non-negative matrix decomposition to obtain R characteristic components, wherein each characteristic component comprises a base vector, a coefficient vector and a phase matrix;

Step 5, obtaining time domain characteristic components by performing inverse short-time Fourier transform on the characteristic components;

and 6, calculating the energy-entropy ratio of each characteristic component by utilizing the coefficient vector of each characteristic component of the track and the time domain characteristic component according to the characteristic that the seismic energy and the time-space information are concentrated in the seismic research area, and selecting the characteristic component with the largest energy-entropy ratio as the seismic component according to the sequence from large to small.

And 7, calculating the root mean square of the time domain seismic components of each track, setting a threshold value, and taking the data point with the time domain amplitude larger than the threshold value in the seismic research area as an abnormal point.

The step 1 comprises the following steps:

the read Swarm A satellite magnetic field data are stored in daily units, and are divided into the track endpoints of 50 DEG S and 50 DEG N to obtain 32 satellite tracks per day.

The step 2 comprises the following steps:

Selecting a satellite magnetic field Y component as original magnetic field data, subtracting a CHASS-7 model from the original magnetic field data of each satellite orbit to obtain a residual magnetic field, and calculating as:

Wherein B _y0 is the original Y component magnetic field, Endogenous field Y component data calculated for the CHASS-7 model, B _y is the residual magnetic field.

In order to remove the low-frequency background, the residual magnetic field is subjected to first order difference, and the obtained preprocessing data is calculated as:

x(n)=B_y(n+1)-B_y(n) (2)

Where B _y (N) is the nth data point of the residual magnetic field, the data length of B _y is set to N, and x (N) is the differential Y component magnetic field result.

The step 3 comprises the following steps:

performing discrete short-time Fourier transform on the preprocessed data x (n) to obtain a complex time-frequency matrix, wherein the discrete short-time Fourier transform formula is as follows:

Where m represents the time index of the input signal x (n), g (n) is the selected window function, g (n-m) is the sliding window, n determines the current sliding position, j represents the imaginary unit, ω is the frequency index, and V (n, ω) represents the complex time-frequency matrix obtained by the discrete short-time fourier transform.

The step 4 includes:

Decomposing the complex time-frequency matrix by complex non-negative matrix decomposition, specifically, representing the complex time-frequency matrix V (n, omega) obtained by short-time Fourier transformation as V _k×l of k rows and l columns, giving the number R of decomposition characteristic components, wherein R is a positive integer and satisfies R < < min (k, l), decomposing the complex time-frequency matrix V _k×l into a matrix W _k×R, a matrix H _R×l and a matrix Wherein W _k×R is referred to as a basis matrix, representing the frequency distribution characteristics of the data, H _R×l is referred to as a coefficient matrix, representing the weights of the different frequency distribution characteristics in time/space,Refers to a time-varying phase spectrum, comprising a phase matrix of k rows and l columns corresponding to R characteristic components, hereinafter abbreviated as V, W, H,The complex non-negative matrix factorization mathematical model is:

Wherein, the The Hadamard product, representing two matrices, is defined as the product of the corresponding elements of the two matrices.

In order to make the decomposition result approach the complex matrix V as much as possible, the KL divergence is used to measure the error between the original matrix and the reconstructed matrix, so that the decomposition result approaches the complex matrix V as much as possible, and the objective function is:

Wherein, the As an objective function, X _k,r,l is an estimated value of the reconstructed complex time-frequency matrix of the R-th feature component, R (H _r,l) is a penalty term based on p-norm, and when 0< p <2, the sparsity of H _r,l can be effectively controlled. The penalty term is:

R(H_r,l)＝2λ|H_r,l|^p,λ>0 (6)

λ is the weight coefficient of the penalty term;

The formulas for the auxiliary functions Z _k,r,l,d_k,r,l,A_k,r,l and B _k,r,l are as follows:

Z_k,r,l＝|X_k,r,l| (7)

finally, the iterative updating rule when 0< p <1 is obtained is as follows:

Z_k,r,l＝|X_k,r,l| (15)

Y_k,r,l＝|X_k,r,l| (16)

U_r,l＝H_r,l (17)

the iterative algorithm steps are as follows:

Input complex time-frequency matrix V, number of decomposed eigenvectors R, upper limit of iteration times N _iter

(1) Initializing W, H and X _k,r,l;

(2) Updating the objective function variable by using the formulas (11) - (17), iterating continuously until the value of the objective function formula (5) converges or reaches the set iteration times N _iter, and stopping iterating;

the output is W, H,

R characteristic components are obtained through decomposition, wherein the R characteristic components comprise an R column vector W _r of a base matrix W, an R row vector H _r of a coefficient matrix H and a corresponding phase matrixH _r coefficient vector, W _r base vector.

The step 5 includes:

Reconstructing each characteristic component back to the time domain through the inverse short-time Fourier transform to obtain each time domain characteristic component, wherein the formula of the inverse short-time Fourier transform is as follows:

ISTFT () represents a function of the inverse short-time Fourier transform.

The step 6 includes:

According to the characteristic that the seismic energy and the space-time information are concentrated in the seismic research area, the energy-entropy ratio of each characteristic component of the track is calculated, wherein the energy-entropy ratio consists of an energy ratio and an entropy ratio, the energy ratio refers to the ratio of the energy in the seismic research area of the r component of the track to the whole track, and the entropy ratio refers to the ratio of shannon entropy in the seismic research area of the time domain characteristic component to the whole track.

Shannon entropy Z _i(y_r) in the seismic investigation region and shannon entropy Z (y _r) for the entire trajectory are calculated as follows:

wherein, the point between s ₁ and s ₂ is in the seismic investigation region, N is the length of the time domain feature component Y _r, p (Y _r) represents the probability that the random event Y is Y _r, and s ₁、s₂ is the minimum latitude and the maximum latitude of the seismic investigation region respectively.

The energy-entropy ratio is calculated as follows:

the points between Ps and Pe are in the seismic research area, L is the number of columns of H _r,l, and Ps and Pe are the minimum latitude and the maximum latitude of the corresponding seismic research area obtained by short-time Fourier transform respectively.

Sorting the characteristic components according to the energy-entropy ratio from large to small, respectively marking the basis vectors of the sorted characteristic components as W _s1、W_s2、…、W_sR, respectively marking the coefficient vectors as H _s1、H_s2、…、H_sR, respectively marking the phase matrixes asThe reconstructed time domain data are denoted as y _s1、y_s2、…、y_sR, respectively, and the feature component with the largest energy-entropy ratio is selected as the seismic component.

The step 7 includes:

The root mean square of the time domain seismic components for each track is calculated as:

the threshold is set to:

Thre=P×RMS (23)

If the data point of the time domain seismic component in the seismic investigation region is greater than Thre, the point is considered to be a seismic outlier, while the trajectory is marked as an outlier trajectory.

Taking the example of the Swarm A star magnetic field Y component data in the area of seismic influence and the study time range of 7.8 grade earthquake occurring in early Cyclodor at 4 months and 16 days of 2016. The seismic influence area, namely the seismic research area, is a square research area which is determined according to Dobrovolsky formula R=10 ^0.43M and takes the center of the earthquake as the center, wherein R is half side length, and the research time range is 60 days before the earthquake to 30 days after the earthquake.

Comprising the following steps:

step 1, downloading and reading magnetic field data of the Swarm A star 2016 from 2 months to 16 months to 5 months, wherein the magnetic field data are stored in a unit of day, and the magnetic field data are required to be divided into tracks with geomagnetic latitude of 50 DEG S and 50 DEG N as endpoints. The period of travel of the Swarm A star around the earth is about 90 minutes, i.e. one orbit can be obtained every 45 minutes, and 32 orbits can be obtained per day.

Taking track 5 of 2016, 4 and 9 days as an example, the vector magnetic field comprises an X (north) component, a Y (east) component and a Z (vertical) component, and the magnetic field of the Y component is possibly influenced by the movement of a rock layer and is less influenced by the disturbance of an external magnetic field, so that the magnetic field data of the Y component is used as the original magnetic field data. Subtracting the CHASS-7 model from the original magnetic field data of each satellite orbit to obtain a residual magnetic field, and calculating as:

Wherein B _y0 is the original Y component magnetic field, Endogenous field Y component data calculated for the CHASS-7 model, B _y is the residual magnetic field.

In order to remove the low-frequency background, the residual magnetic field is subjected to first order difference, and the obtained preprocessing data is calculated as:

x(n)=B_y(n+1)-B_y(n)(2)

Where B _y (N) is the nth data point of the residual magnetic field, the data length of B _y is set to N, and x (N) is the differential Y component magnetic field result.

Step 3, in order to obtain the characteristic of the frequency change of the magnetic field signal along with time, performing discrete short-time Fourier transform on the preprocessed data x (n) to obtain a complex time-frequency matrix, wherein the discrete short-time Fourier transform formula is as follows:

where n represents the current sliding position, g (n) is the selected window function, and g (n-m) is the sliding window.

And 4, decomposing the complex time-frequency matrix by utilizing the advantages of the local feature extraction of complex non-negative matrix decomposition to separate the features caused by the earthquake and the non-earthquake factors because the influence of the non-earthquake factors such as the sun activity, the magnetic storm and the like on the magnetic field signals is global and the influence of the earthquake is more likely to be local. Complex non-negative matrix factorization considers the problem that given a complex matrix V ε C ^k×l and a positive integer R, R < < min (k, l) is satisfied, and a matrix W ε R ^k×r、H∈R^r×l andThe mathematical model is as follows:

Wherein, the The Hadamard product, representing two matrices, is defined as the product of the corresponding elements of the two matrices.

The KL divergence is used for measuring the error between the original matrix and the reconstruction matrix, so that the decomposition result approaches the complex matrix V as much as possible, and the objective function is as follows:

Wherein R (H _r,l) is a penalty term based on p-norm, and when 0< p <2, the sparsity of H _r,l can be effectively controlled. The penalty term is:

R(H_r,l)＝2λ|H_r,l|^p,λ>0 (6)

The formulas for the auxiliary functions Z _k,r,l,d_k,r,l,A_k,r,l and B _k,r,l are as follows:

Z_k,r,l＝|X_k,r,l| (7)

finally, the iterative updating rule when 0< p <1 is obtained is as follows:

Z_k,r,l＝|X_k,r,l| (15)

Y_k,r,l＝|X_k,r,l| (16)

U_r,l＝H_r,l (17)

the iterative algorithm steps are as follows:

(1) Initializing W _k,r,H_r,l,X_k,r,l;

(2) Updating the objective function variable by using the formulas (10) - (16), and continuously iterating until the value of the objective function formula (4) converges or reaches the set iteration number N _iter =250, and stopping iterating.

Taking 2016, 4, 9, and 5 as an example, a transformation curve of the complex non-negative matrix factorization objective function error value with the number of iterations is shown in fig. 2, and it can be seen that the final objective function error value converges. The 3 eigenvectors are obtained by decomposition, and as shown in fig. 3, the r eigenvector includes the r (r=1, 2, 3) th column vector W _r (base vector) of the base matrix, the r row vector H _r (coefficient vector) of the coefficient matrix, and the corresponding phase matrix

And 5, reconstructing each characteristic component back to the time domain through inverse short-time Fourier transform to obtain each time domain characteristic component, wherein the formula of the inverse short-time Fourier transform is as follows:

And 6, calculating the energy-entropy ratio of each characteristic component of the track according to the characteristic that the seismic energy and the space-time information are concentrated in the seismic research area, wherein the energy-entropy ratio consists of an energy ratio and an entropy ratio, the energy ratio is the ratio of the energy in the seismic research area of the r component of the track to the whole track and represents the seismic energy distribution, and the entropy ratio is the ratio of the Shannon entropy in the seismic research area of the time domain characteristic component to the whole track and represents the seismic information distribution.

Shannon entropy Z _i(y_r in the seismic investigation region), shannon entropy Z (y _r) for the entire trajectory, and the energy-to-entropy ratio is calculated as follows:

Where N is the length of the time domain feature component y _r and L is the number of columns of H _r,l for points between s ₁ to s ₂ and points between Ps to Pe within the seismic region.

Sorting the characteristic components according to the energy-entropy ratio from large to small, respectively marking the basis vectors of the sorted characteristic components as W _s1、W_s2、W_s3, respectively marking the coefficient vectors as H _s1、H_s2、H_s3, respectively marking the phase matrixes asThe reconstructed time domain data are respectively recorded as y _s1、y_s2、y_s3, and the characteristic component with the largest energy-entropy ratio is selected as the seismic component, and the sequencing result of the track 5 of 2016, 4 and 9 days is shown in fig. 4 and 5. The fundamental vector frequency W _s1 of the seismic component is mainly distributed around 0.04Hz, and ULF waves conforming to low earth orbit observations are typically in the frequency range of Pc3 (20-100 mHz). Meanwhile, the coefficient vector H _s1 of the seismic component only has obvious abnormality in a research area, the coefficient vector H _s2 of the other characteristic components is distributed on the whole track, the abnormality of H _s3 exists outside the research area, and the rule that the energy of the seismic component is mainly distributed in the seismic influence area is met.

Step 7, calculating the root mean square of the time domain seismic components of each track as follows:

the threshold is set to:

Thre=P×RMS (23)

Taking the example of the track 5 of the year 2016, month 4 and day 9, the root mean square of the time domain seismic component of the track is 0.0142, the threshold is set to be 3 times root mean square, namely thre= 0.0426, as shown by the horizontal dashed line of the seismic component y _s1 in fig. 5, 4 data points of the time domain seismic component of the track, which are larger than Thre, in the seismic research area are marked as seismic outliers, and the track is marked as an outlier.

The foregoing description of the preferred embodiments of the invention is not intended to be limiting, but rather is intended to cover all modifications, equivalents, and alternatives falling within the spirit and principles of the invention.

Claims (9)

1. A satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization is characterized by comprising the following steps:

step 1, reading Swarm A satellite magnetic field data, and dividing the daily magnetic field data into a plurality of orbits;

Step 2, subtracting the CHASS-7 model endogenous field from the original Y component magnetic field data of each satellite orbit, removing the main magnetic field of the earth core and the rock ring magnetic field of the earth crust to obtain a residual magnetic field;

Step 3, carrying out short-time Fourier transform on the preprocessed data to obtain a complex time-frequency matrix in a rectangular coordinate representation form of a complex domain;

decomposing the complex time-frequency matrix by utilizing complex non-negative matrix decomposition to obtain R characteristic components, wherein each characteristic component comprises a base vector, a coefficient vector and a phase matrix;

step 5, obtaining time domain characteristic components from the characteristic components through inverse short-time Fourier transform;

Step 6, according to the characteristic that the seismic energy and the space-time information are concentrated in the seismic research area, calculating the energy-entropy ratio of each characteristic component by utilizing the coefficient vector of each characteristic component of the track and the time domain characteristic component, and selecting the characteristic component with the largest energy-entropy ratio as the seismic component according to the sequence from big to small;

and 7, calculating the root mean square of each track time domain seismic component, setting a threshold value by using an overrun method, and extracting abnormal points from the track seismic components.

2. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 1, wherein,

The step 1 comprises the following steps:

The method comprises the steps of reading Swarm A satellite magnetic field data, dividing daily magnetic field data into a plurality of orbits, specifically, storing the magnetic field data in a daily unit, dividing the magnetic field data by using 50 DEG S and 50 DEG N as orbit endpoints, and obtaining 32 satellite orbits daily.

3. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 1, wherein,

The step 2 comprises the following steps:

The raw Y-component magnetic field data for each satellite orbit minus the CHAOS-7 model endogenous field is calculated as:

Wherein B _y0 is the original Y component magnetic field, Endogenous field Y component data calculated for the CHASS-7 model, B _y being the residual Y component magnetic field;

the residual magnetic field is differentiated, and the calculation is as follows:

x(n)=B_y(n+1)-B_y(n) (2)

Where B _y (N) is the nth data point of the remaining Y component magnetic field, the data length of B _y is set to N, and x (N) is the differential Y component magnetic field.

4. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 1, wherein,

The step 3 comprises the following steps:

short-time Fourier transform is carried out on the differential Y component magnetic field x (n), and the discrete short-time Fourier transform is specifically as follows:

Where m represents the time index of the input signal x (n), g (n) is the selected window function, g (n-m) is the sliding window, n determines the current sliding position, j represents the imaginary unit, ω is the frequency index, and V (n, ω) represents the complex time-frequency matrix obtained by the discrete short-time fourier transform.

5. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 4, wherein,

The step 4 includes:

Decomposing the complex time-frequency matrix by complex non-negative matrix decomposition, specifically, representing the complex time-frequency matrix V (n, omega) obtained by short-time Fourier transformation as V _k×l of k rows and l columns, giving the number R of decomposition characteristic components, wherein R is a positive integer and satisfies R < < min (k, l), decomposing the complex time-frequency matrix V _k×l into a matrix W _k×R, a matrix H _R×l and a matrix Wherein W _k×R is referred to as a basis matrix, representing the frequency distribution characteristics of the data, H _R×l is referred to as a coefficient matrix, representing the weights of the different frequency distribution characteristics in time/space,Refers to a time-varying phase spectrum, comprising a phase matrix of k rows and l columns corresponding to R characteristic components, hereinafter abbreviated as V, W, H,The complex non-negative matrix factorization mathematical model is:

Wherein, the A Hadamard product representing two matrices, the elements of which are defined as the product of the corresponding elements of the two matrices, W _k×r representing the r-th column of the base matrix W, H _r×l representing the r-th row of the coefficient matrix H,Representing time-varying phase spectraR is more than or equal to 1 and less than or equal to R.

6. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 5, wherein,

Measuring original matrix V and reconstructed matrix by KL divergenceThe error between them, the objective function is:

Wherein, the As an objective function, X _k,r,l is an estimated value of a reconstructed complex time-frequency matrix of the R-th feature component, R (H _r,l) is a penalty term based on a p-norm, and when 0< p <2, the sparsity of H _r,l can be effectively controlled, and the penalty term is:

R(H_r,l)＝2λ|H_r,l|^p,λ>0 (5)

λ is the weight coefficient of the penalty term;

The formulas for the auxiliary functions Z _k,r,l,d_k,r,l,A_k,r,l and B _k,r,l are as follows:

Z_k,r,l＝|X_k,r,l| (6)

finally, the iterative updating rule when 0< p <1 is obtained is as follows:

Z_k,r,l＝|X_k,r,l| (14)

Y_k,r,l＝|X_k,r,l| (15)

U_r,l＝H_r,l (16)

the iterative algorithm steps are as follows:

Input complex time-frequency matrix V, number of decomposed eigenvectors R, upper limit of iteration times N _iter

(1) Initializing W, H and X _k,r,l;

(2) Updating the objective function variable by using the formulas (11) - (17), iterating continuously until the value of the objective function formula (5) converges or reaches the set iteration times N _iter, and stopping iterating;

the output is W, H,

R characteristic components are obtained through decomposition, wherein the R characteristic components comprise an R column vector W _r of a base matrix W, an R row vector H _r of a coefficient matrix H and a corresponding phase matrixH _r coefficient vector, W _r base vector.

7. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 6, wherein,

The step 5 includes:

the time domain characteristic components are obtained by carrying out inverse short-time Fourier transform on the characteristic components, and the time domain characteristic components are specifically:

ISTFT () represents a function of the inverse short-time Fourier transform.

8. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 1, wherein,

The step 6 includes:

the energy-entropy ratio of each characteristic component of the track consists of an energy ratio and an entropy ratio, wherein the energy ratio refers to the ratio of the energy to the whole track in the seismic research area of the r component of the track, and the entropy ratio refers to the ratio of shannon entropy to the whole track in the seismic research area of the time domain characteristic component;

Shannon entropy Z _i(y_r) in the seismic investigation region and shannon entropy Z (y _r) for the entire trajectory are calculated as follows:

Wherein, the point between s ₁ and s ₂ is in the seismic research area, N is the length of the time domain feature component Y _r, p (Y _r) represents the probability that the random event Y is Y _r, and s ₁、s₂ is the minimum latitude and the maximum latitude of the seismic research area respectively;

The energy-entropy ratio is calculated as follows:

Wherein, the points between Ps and Pe are in the seismic research area, L is the column number of H _r,l, and Ps and Pe are the minimum latitude and the maximum latitude of the corresponding seismic research area obtained by short-time Fourier transform respectively;

Sorting the characteristic components according to the energy-entropy ratio from large to small, respectively marking the basis vectors of the sorted characteristic components as W _s1、W_s2、…、W_sR, respectively marking the coefficient vectors as H _s1、H_s2、…、H_sR, respectively marking the phase matrixes as The reconstructed time domain data are denoted as y _s1、y_s2、…、y_sR, respectively, and the feature component with the largest energy-entropy ratio is selected as the seismic component.

9. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix factorization of claim 1, wherein,

The step 7 includes:

the root mean square of the time domain seismic components of each track is calculated as:

the threshold is set as:

Thre=P×RMS (22)

p is an empirical parameter, and Thre is an abnormality extraction threshold;

when the amplitude of the time domain seismic component data point in the seismic investigation region is greater than Thre, the point is considered to be a seismic outlier while the trajectory is marked as an outlier trajectory.