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

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

Technical Field

[0001] The present invention belongs to the field of ionospheric satellite magnetic field seismic anomaly extraction, and specifically provides a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition. Background Art

[0002] Seismic electromagnetics uses non-mechanical electromagnetic methods to study the development and occurrence of earthquakes, becoming an effective way to deepen our understanding of earthquake development. With the rapid development of space exploration technology, ionospheric exploration, especially satellite exploration, is increasingly being used in earthquake precursor research. The study of satellite electromagnetic seismic anomalies has become a major focus in international geoscience research.

[0003] Satellite magnetic field observation data is subject not only to interference from Earth but also to strong interference from non-seismic factors such as solar activity and magnetic storms. Traditional research methods, to avoid non-seismic interference, use nighttime data while removing data with high solar activity and geomagnetic indices. To improve the utilization of seismic data and extract seismic anomalies from satellite magnetic field observation data, which exhibit "strong background and weak information," blind source separation methods such as non-negative matrix decomposition are used to separate seismic signals, background, and strong geomagnetic interference from the data. However, this method currently has certain limitations: magnetic field signals contain not only waveform and amplitude information, but also other important information such as phase. These methods often perform time-frequency analysis on the magnetic field signal and then decompose the amplitude spectrum matrix. The resulting components are then reconstructed back into the time domain using the original phase matrix. This method ignores the influence of the original magnetic field signal phase information on anomaly extraction during reconstruction, and thus fails to fully utilize the complex frequency domain information to improve the accuracy of seismic magnetic field anomaly identification and extraction.

[0004] Therefore, in order to effectively utilize phase information, it is necessary to conduct in-depth research on the current complex domain matrix decomposition method and expand the satellite magnetic field seismic anomaly extraction from the real domain to the complex domain, so as to deepen the understanding of seismic anomaly phenomena. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition, which makes full use of the amplitude and phase information contained in satellite observation data in the time and frequency domain, thereby identifying earthquake-related components and making seismic anomaly extraction more accurate and reliable.

[0006] The present invention is implemented as follows: a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition comprises the following steps:

[0007] Step 1: Read the Swarm A satellite magnetic field data and split the daily magnetic field data into multiple tracks;

[0008] Step 2: Subtract the CHAOS-7 model internal source field from the raw Y-component magnetic field data of each satellite orbit to remove the main magnetic field of the core and the lithospheric magnetic field of the crust to obtain the residual magnetic field; perform differential analysis on the residual magnetic field to remove the low-frequency background and obtain the magnetic field preprocessing data;

[0009] Step 3: Perform short-time Fourier transform on the preprocessed data to obtain a complex time-frequency matrix in the rectangular coordinate representation of the complex domain;

[0010] Step 4: Decompose the complex time-frequency matrix using complex non-negative matrix decomposition to obtain R eigencomponents, each of which contains a basis vector, a coefficient vector, and a phase matrix.

[0011] Step 5: Obtain the time domain feature components by inverse short-time Fourier transform of each feature component;

[0012] Step 6: Based on the characteristics of earthquake energy and spatiotemporal information being concentrated in the earthquake study area, the energy-entropy ratio of each characteristic component is calculated using the coefficient vector of each characteristic component of the orbit and the time domain characteristic component. The components are sorted in descending order and the characteristic component with the largest energy-entropy ratio is selected as the earthquake component.

[0013] Step 7: By calculating the root mean square of the time domain seismic component of each track, using the exceedance rate method, setting the threshold, and extracting abnormal points in the track seismic component.

[0014] Furthermore, the step 1 includes:

[0015] Read the Swarm A satellite magnetic field data and split the daily magnetic field data into multiple tracks. Specifically, the magnetic field data is stored in daily units and is split with 50°S and 50°N as the orbit endpoints, resulting in 32 satellite tracks per day.

[0016] Furthermore, the step 2 includes:

[0017] The original Y-component magnetic field data of each satellite orbit minus the source field of the CHAOS-7 model is calculated as:

[0018]

[0019] Among them, B y0 is the original Y component magnetic field, The Y component data of the internal source field calculated by the CHAOS-7 model, B y is the residual Y component magnetic field;

[0020] The residual magnetic field is differentiated and calculated as:

[0021] x(n)=B y (n+1)-B y (n) (2)

[0022] Among them, B y (n) is the nth data point of the residual Y component magnetic field, B y The data length is set to N, and x(n) is the differential Y component magnetic field.

[0023] Furthermore, the step 3 includes:

[0024] Perform short-time Fourier transform on the differential Y-component magnetic field x(n). The discrete short-time Fourier transform is:

[0025]

[0026] Where m represents the time index of the input signal x(n), g(n) is the selected window function, g(nm) 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 discrete short-time Fourier transform.

[0027] Furthermore, the step 4 includes:

[0028] The complex time-frequency matrix is ​​decomposed using complex non-negative matrix decomposition. Specifically, the complex time-frequency matrix V(n,ω) obtained by short-time Fourier transform is expressed as V with k rows and l columns. k×l , given the number of decomposition eigencomponents R, R is a positive integer and satisfies R<<min(k,l), the complex time-frequency matrix V k×l Decompose into matrix W k×R , matrix H R×l and The product of W k×R It is called the basis matrix, which represents the frequency distribution characteristics of the data. R×l It is called the coefficient matrix, which represents the weight of different frequency distribution features in time / space. Refers to the time-varying phase spectrum, which contains the phase matrix of k rows and l columns corresponding to R characteristic components, hereinafter abbreviated as V, W, H, The mathematical model of complex non-negative matrix decomposition is:

[0029]

[0030] in, Represents the Hadamard product of two matrices, whose elements are defined as the product of the corresponding elements of the two matrices, W k×r represents the rth column of the basis matrix W, H r×l represents the rth row of the coefficient matrix H, Represents the time-varying phase spectrum The r-th phase matrix of , 1≤r≤R.

[0031] In order to make the decomposition result as close as possible to the complex matrix V, the KL divergence is used to measure the difference between the original matrix V and the reconstructed matrix The error between them, its objective function is:

[0032]

[0033] in, is the objective function, X k,r,l is the estimated value of the reconstructed complex time-frequency matrix of the rth eigencomponent, R(H r,l ) is a penalty term based on the p-norm. When 0<p<2, H r,l The sparsity of can be effectively controlled. The penalty term is:

[0034] R(H r,l )=2λ|H r,l | p ,λ>0 (5)

[0035] λ is the weight coefficient of the penalty term;

[0036] Auxiliary function Z k,r,l , d k,r,l , A k,r,l and B k,r,l The calculation formula is as follows:

[0037] Z k,r,l =|X k,r,l | (6)

[0038]

[0039] Finally, the iterative update rule when 0<p<1 is:

[0040]

[0041] Z k,r,l =|X k,r,l | (14)

[0042] Y k,r,l =|X k,r,l | (15)

[0043] U r,l =H r,l (16)

[0044] The iterative algorithm steps are as follows:

[0045] Input: complex time-frequency matrix V, number of decomposed eigenvectors R, upper limit of iteration number N iter

[0046] (1) Initialize W, H, X k,r,l ;

[0047] (2) Update the objective function variables using equations (11)-(17) and continue iterating until the value of the objective function equation (5) converges or reaches the set number of iterations N. iter , stop iteration;

[0048] Output: W, H,

[0049] After decomposition, R eigencomponents are obtained, and the rth eigencomponent contains the rth column vector W of the basis matrix W r , the rth row vector H of the coefficient matrix H r and the corresponding phase matrix H r Coefficient vector, W r Basis vectors.

[0050] Furthermore, the step 5 includes:

[0051] Each characteristic component is transformed by inverse short-time Fourier transform to obtain the time domain characteristic component, specifically:

[0052]

[0053] ISTFT() represents the inverse short-time Fourier transform function.

[0054] Furthermore, the step 6 includes:

[0055] The energy-entropy ratio of each characteristic component of the orbit consists of two parts: the energy ratio and the entropy ratio. The energy ratio refers to the ratio of the energy in the seismic study area of ​​the r component of the orbit to the entire orbit; the entropy ratio refers to the ratio of the Shannon entropy in the seismic study area of ​​the time domain characteristic component to the entire orbit.

[0056] Shannon entropy Z in the earthquake study area i (y r ) and the Shannon entropy Z(y r ) is calculated as follows:

[0057]

[0058] Among them, the points between s1 and s2 are within the earthquake study area, and N is the time domain characteristic component y r The length of p(y r ) indicates that the random event Y is y r , s1 and s2 are the minimum latitude and maximum latitude of the earthquake study area respectively.

[0059] The energy-entropy ratio is then calculated as follows:

[0060]

[0061] Among them, the points between Ps and Pe are within the earthquake study area, and L is H r,l Ps and Pe are the minimum latitude and maximum latitude of the corresponding earthquake study area obtained by short-time Fourier transform.

[0062] Sort the eigenvalues ​​from large to small according to the energy-entropy ratio, and denote the basis vectors of the sorted eigenvalues ​​as W s1 、W s2 ,...,W sR , the coefficient vectors are recorded as H s1 、H s2 ,…,H sR , the phase matrices are recorded as The reconstructed time domain data are recorded as y s1 、y s2 ,…,y sR , and the characteristic component with the largest energy-entropy ratio is selected as the seismic component.

[0063] Furthermore, the step 7 includes:

[0064] Sub-step ① The root mean square of the time domain seismic component of each track is calculated as:

[0065]

[0066] The threshold value in sub-step ② is set as:

[0067] Thre=P×RMS (22)

[0068] P is an empirical parameter, generally set to 3; Thre is the anomaly extraction threshold

[0069] When the amplitude of a time-domain seismic component data point in the earthquake study area is greater than Thre, the point is considered to be a seismic anomaly point, and the track is marked as an abnormal track.

[0070] Compared with the prior art, the present invention has the following beneficial effects:

[0071] The present invention is a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition. The method removes the background field of the satellite magnetic field data by subtracting the CHAOS-7 model to remove the endogenous field and then performing first-order difference. The short-time Fourier transform is used to perform time-frequency transformation on the pre-processed data after the background field is removed to obtain a complex time-frequency matrix. According to the characteristics of the influence of non-seismic factors such as solar activity and magnetic storms on the global magnetic field and the influence of earthquakes on the local magnetic field, the advantages of local feature extraction of complex non-negative matrix decomposition are used to obtain multiple feature components. Each feature component includes a basis vector, a coefficient vector and a corresponding phase matrix, which fully utilizes the amplitude and phase information contained in the satellite magnetic field signal in the time-frequency domain. The reconstructed time domain feature components are obtained by inverse short-time Fourier transform. According to the characteristics of the distribution of seismic energy and spatiotemporal information, the energy-entropy ratio is used to identify the seismic component. The root mean square of the time domain seismic component is calculated, and the threshold is set by the exceedance rate method to extract the seismic anomaly. Based on the characteristic that the distribution of seismic magnetic field anomalies is concentrated in the earthquake-affected area, the present invention utilizes the advantages of local feature extraction of complex non-negative matrix decomposition and the full utilization of complex domain amplitude information and phase information to separate the anomaly features caused by seismic and non-seismic factors, making seismic anomaly extraction more accurate and reliable, thereby deepening the understanding of ionospheric seismic magnetic field anomaly phenomena.

[0072] Based on the advantages of local feature extraction, the method of the present invention can decompose the original complex matrix V into two non-negative matrices W∈R k×r , H∈R r×l Its phase spectrum The product form of W is a matrix representing the frequency distribution characteristics of the data, and H is a coefficient matrix representing the time and space weights of different frequency distribution characteristics. After removing the background from the satellite magnetic field data, a short-time Fourier transform is performed to obtain a time-frequency spectrum matrix in the complex domain. This matrix is ​​then subjected to complex non-negative matrix decomposition to obtain seismic components, which are then reconstructed back into the time domain using an inverse short-time Fourier transform. Seismic anomalies are extracted using the excess rate method and the abnormal orbits are determined. This method fully utilizes the amplitude and phase information contained in the magnetic field signal, increasing the reliability of the seismic anomaly extraction results. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 This is a flow chart of the satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition;

[0074] Figure 2 It is the error value iteration curve of the objective function of the complex non-negative matrix decomposition algorithm;

[0075] Figure 3 is the basis matrix and coefficient matrix obtained by complex non-negative matrix decomposition;

[0076] Figure 4The basis matrix and coefficient matrix are sorted from large to small according to the energy-entropy ratio;

[0077] Figure 5 Reconstruct the time domain feature components after sorting. DETAILED DESCRIPTION

[0078] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0079] See also Figure 1 As shown, a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition includes:

[0080] Step 1: Read the Swarm A satellite magnetic field data and split the daily data into multiple tracks with 50°S-50°N as the endpoints;

[0081] Step 2: Select the Y component of the magnetic field data as the original magnetic field data. Subtract the CHAOS-7 model internal source field from the original magnetic field data of each track to remove the main magnetic field and the lithospheric magnetic field to obtain the residual magnetic field. In order to remove the low-frequency background, perform a first-order difference on the residual magnetic field to obtain the magnetic field preprocessing data. The CHAOS-7 model internal source field is obtained using existing technology.

[0082] Step 3: Use short-time Fourier transform to perform time-frequency transformation on the preprocessed data to obtain a complex time-frequency matrix;

[0083] Step 4: Decompose the complex time-frequency matrix using complex non-negative matrix decomposition to obtain R eigencomponents, each of which contains a basis vector, a coefficient vector, and a phase matrix.

[0084] Step 5: Obtain the time domain feature components by performing inverse short-time Fourier transform on each feature component;

[0085] Step 6: Based on the characteristics that earthquake energy and spatiotemporal information are concentrated in the earthquake study area, the energy-entropy ratio of each characteristic component is calculated using the coefficient vector of each characteristic component of the orbit and the time domain characteristic component, and they are sorted in order from large to small, and the characteristic component with the largest energy-entropy ratio is selected as the earthquake component.

[0086] Step 7: Calculate the root mean square of the time domain seismic component of each track and set a threshold. The data points with time domain amplitude greater than the threshold in the seismic study area are considered outliers.

[0087] The step 1 comprises:

[0088] The magnetic field data of the Swarm A satellite is stored in daily units and is divided into 32 satellite orbits per day by using 50°S and 50°N as the orbit endpoints.

[0089] The step 2 includes:

[0090] The Y component of the satellite magnetic field is selected as the original magnetic field data. In order to remove the internal source fields including the main magnetic field generated by the core and the lithospheric magnetic field generated by the crust, the original magnetic field data of each satellite orbit is subtracted from the CHAOS-7 model to obtain the residual magnetic field, which is calculated as:

[0091]

[0092] Among them, B y0 is the original Y component magnetic field, The Y component data of the internal source field calculated by the CHAOS-7 model, B y is the residual magnetic field.

[0093] In order to remove the low-frequency background, the residual magnetic field is subjected to first-order difference, and the preprocessed data is calculated as:

[0094] x(n)=B y (n+1)-B y (n) (2)

[0095] Among them, B y (n) is the nth data point of the residual magnetic field, B y The data length is set to N, and x(n) is the differential Y component magnetic field result.

[0096] The step 3 comprises:

[0097] Perform discrete short-time Fourier transform on the preprocessed data x(n) to obtain a complex time-frequency matrix. The discrete short-time Fourier transform formula is as follows:

[0098]

[0099] Where m represents the time index of the input signal x(n), g(n) is the selected window function, g(nm) 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 discrete short-time Fourier transform.

[0100] The step 4 comprises:

[0101] The complex time-frequency matrix is ​​decomposed using complex non-negative matrix decomposition. Specifically, the complex time-frequency matrix V(n,ω) obtained by short-time Fourier transform is expressed as V with k rows and l columns. k×l, given the number of decomposition eigencomponents R, R is a positive integer and satisfies R<<min(k,l), the complex time-frequency matrix V k×l Decompose into matrix W k×R , matrix H R×l and The product of W k×R It is called the basis matrix, which represents the frequency distribution characteristics of the data. R×l It is called the coefficient matrix, which represents the weight of different frequency distribution features in time / space. Refers to the time-varying phase spectrum, which contains the phase matrix of k rows and l columns corresponding to R characteristic components, hereinafter abbreviated as V, W, H, The mathematical model of complex non-negative matrix decomposition is:

[0102]

[0103] in, Represents the Hadamard product of two matrices, defined as the product of corresponding elements of the two matrices.

[0104] In order to make the decomposition result as close as possible to the complex matrix V, the KL divergence is used to measure the error between the original matrix and the reconstructed matrix, so that the decomposition result is as close as possible to the complex matrix V. The objective function is:

[0105]

[0106] in, is the objective function, X k,r,l is the estimated value of the reconstructed complex time-frequency matrix of the rth eigencomponent, R(H r,l ) is a penalty term based on the p-norm. When 0<p<2, H r,l The sparsity of can be effectively controlled. The penalty term is:

[0107] R(H r,l )=2λ|H r,l | p ,λ>0 (6)

[0108] λ is the weight coefficient of the penalty term;

[0109] Auxiliary function Z k,r,l , d k,r,l , A k,r,l and B k,r,l The calculation formula is as follows:

[0110] Z k,r,l =|X k,r,l | (7)

[0111]

[0112] Finally, the iterative update rule when 0<p<1 is:

[0113]

[0114] Z k,r,l =|X k,r,l | (15)

[0115] Y k,r,l =|X k,r,l | (16)

[0116] U r,l =H r,l (17)

[0117] The iterative algorithm steps are as follows:

[0118] Input: complex time-frequency matrix V, number of decomposed eigenvectors R, upper limit of iteration number N iter

[0119] (1) Initialize W, H, X k,r,l ;

[0120] (2) Update the objective function variables using equations (11)-(17) and continue iterating until the value of the objective function equation (5) converges or reaches the set number of iterations N. iter , stop iteration;

[0121] Output: W, H,

[0122] After decomposition, R eigencomponents are obtained, and the rth eigencomponent contains the rth column vector W of the basis matrix W r , the rth row vector H of the coefficient matrix H r and the corresponding phase matrix H r Coefficient vector, W r Basis vectors.

[0123] The step 5 comprises:

[0124] Through the inverse short-time Fourier transform, each characteristic component is reconstructed back to the time domain to obtain each time domain characteristic component. The formula of the inverse short-time Fourier transform is as follows:

[0125]

[0126] ISTFT() represents the inverse short-time Fourier transform function.

[0127] The step 6 comprises:

[0128] According to the characteristics that earthquake energy and spatiotemporal information are concentrated in the earthquake study area, the energy-entropy ratio of each characteristic component of the orbit is calculated. The energy-entropy ratio consists of two parts: the energy ratio and the entropy ratio. The energy ratio refers to the ratio of the energy of the r component of the orbit in the earthquake study area to the entire orbit; the entropy ratio refers to the ratio of the Shannon entropy of the time domain characteristic component in the earthquake study area to the entire orbit.

[0129] Shannon entropy Z in the earthquake study area i (y r ) and the Shannon entropy Z(y r ) is calculated as follows:

[0130]

[0131] Among them, the points between s1 and s2 are within the earthquake study area, and N is the time domain characteristic component y r The length of p(y r ) indicates that the random event Y is y r , s1 and s2 are the minimum latitude and maximum latitude of the earthquake study area respectively.

[0132] The energy-entropy ratio is then calculated as follows:

[0133]

[0134] Among them, the points between Ps and Pe are within the earthquake study area, and L is H r,l Ps and Pe are the minimum latitude and maximum latitude of the corresponding earthquake study area obtained by short-time Fourier transform.

[0135] Sort the eigenvalues ​​from large to small according to the energy-entropy ratio, and denote the basis vectors of the sorted eigenvalues ​​as W s1 、W s2 ,…,W sR , the coefficient vectors are recorded as H s1 、H s2 ,…,H sR , the phase matrices are recorded as The reconstructed time domain data are recorded as y s1 、y s2 ,…,y sR , and the characteristic component with the largest energy-entropy ratio is selected as the seismic component.

[0136] The step 7 comprises:

[0137] The root mean square of the time-domain seismic component for each track is calculated as:

[0138]

[0139] The threshold is set to:

[0140] Thre=P×RMS (23)

[0141] If the time domain seismic component of a data point in the seismic study area is greater than Thre, the point is considered to be a seismic anomaly point, and the track is marked as an abnormal track.

[0142] Take the Y component data of the magnetic field of Swarm A within the earthquake impact area and the study time range of the 7.8 magnitude earthquake that occurred in Ecuador on April 16, 2016 as an example. The earthquake impact area, that is, the earthquake study area, is calculated based on the Dobrovolsky formula R = 10 0.43M The determined research area is a square with the epicenter as the center and R as half the length of the side; the research time range is 60 days before the earthquake to 30 days after the earthquake.

[0143] include:

[0144] Step 1: Download and read Swarm A's magnetic field data from February 16 to May 16, 2016. Since the data is stored on a daily basis, it needs to be segmented into orbits centered at 50°S and 50°N. Swarm A orbits the Earth approximately every 90 minutes, meaning one orbit is obtained every 45 minutes, or 32 orbits per day.

[0145] Step 2: Taking orbit 5 on April 9, 2016 as an example, the vector magnetic field includes X (north), Y (east), and Z (vertical) components. Since the Y component of the magnetic field may be affected by rock layer activity but is less affected by external magnetic field interference, the Y component magnetic field data is used as the raw magnetic field data. The CHAOS-7 model is subtracted from the raw magnetic field data of each satellite orbit to obtain the residual magnetic field, which is calculated as:

[0146]

[0147] Among them, B y0 is the original Y component magnetic field, The Y component data of the internal source field calculated by the CHAOS-7 model, B y is the residual magnetic field.

[0148] In order to remove the low-frequency background, the residual magnetic field is subjected to first-order difference, and the preprocessed data is calculated as:

[0149] x(n)=B y (n+1)-B y (n)(2)

[0150] Among them, B y (n) is the nth data point of the residual magnetic field, B yThe data length is set to N, and x(n) is the differential Y component magnetic field result.

[0151] Step 3: In order to obtain the characteristics of the frequency change of the magnetic field signal over time, the preprocessed data x(n) is subjected to discrete short-time Fourier transform to obtain a complex time-frequency matrix. The discrete short-time Fourier transform formula is as follows:

[0152]

[0153] Where n represents the current sliding position, g(n) is the selected window function, and g(nm) is the sliding window.

[0154] Step 4: Since the impact of non-seismic factors such as solar activity and magnetic storms on magnetic field signals is global, while the impact of earthquakes is more likely to be local, the advantage of local feature extraction of complex non-negative matrix decomposition is used to decompose the complex time-frequency matrix to separate the features caused by earthquakes and non-seismic factors. Complex non-negative matrix decomposition considers the following problem: Given a complex matrix V∈C k×l and a positive integer R, satisfying R<<min(k,l), find the matrix W∈R k×r , H∈R r×l as well as Its mathematical model is:

[0155]

[0156] in, Represents the Hadamard product of two matrices, defined as the product of corresponding elements of the two matrices.

[0157] The KL divergence is used to measure the error between the original matrix and the reconstructed matrix, so that the decomposition result is as close as possible to the complex matrix V. The objective function is:

[0158]

[0159] Among them, R(H r,l ) is a penalty term based on the p-norm. When 0<p<2, H r,l The sparsity of can be effectively controlled. The penalty term is:

[0160] R(H r,l )=2λ|H r,l | p ,λ>0 (6)

[0161] Auxiliary function Z k,r,l , d k,r,l , A k,r,l and B k,r,l The calculation formula is as follows:

[0162] Z k,r,l =|Xk,r,l | (7)

[0163]

[0164] Finally, the iterative update rule when 0<p<1 is:

[0165]

[0166] Z k,r,l =|X k,r,l | (15)

[0167] Y k,r,l =|X k,r,l | (16)

[0168] U r,l =H r,l (17)

[0169] The iterative algorithm steps are as follows:

[0170] (1) Initialize W k,r , H r,l , X k,r,l ;

[0171] (2) Update the objective function variables using equations (10)-(16) and continue iterating until the value of the objective function equation (4) converges or reaches the set number of iterations N. iter =250, stop iteration.

[0172] Taking track 5 on April 9, 2016 as an example, the transformation curve of the error value of the complex non-negative matrix decomposition objective function with the number of iterations is as follows: Figure 2 As shown in Figure 2, it can be seen that the final objective function error value converges. After decomposition, three characteristic components are obtained, such as Figure 3 As shown, the rth eigencomponent contains the rth (r=1,2,3) column vector W of the basis matrix r (basis vector), the rth row vector of the coefficient matrix H r (coefficient vector) and the corresponding phase matrix

[0173] Step 5: Reconstruct each characteristic component back to the time domain through inverse short-time Fourier transform to obtain each time domain characteristic component. The formula of inverse short-time Fourier transform is as follows:

[0174]

[0175] Step 6: Based on the characteristics that earthquake energy and spatiotemporal information are concentrated in the earthquake study area, calculate the energy-entropy ratio of each characteristic component of the orbit. The energy-entropy ratio consists of two parts: the energy ratio and the entropy ratio. The energy ratio refers to the ratio of the energy in the earthquake study area of ​​the r component of the orbit to the energy of the entire orbit, which represents the distribution of earthquake energy. The entropy ratio refers to the ratio of the Shannon entropy in the earthquake study area of ​​the time domain characteristic component to the entire orbit, which represents the distribution of earthquake information.

[0176] Shannon entropy Z in the earthquake study area i (y r ), the Shannon entropy Z(y r ) and the energy-entropy ratio are calculated as follows:

[0177]

[0178] Among them, the points between s1 and s2 and the points between Ps and Pe are in the earthquake study area, and N is the time domain characteristic component y r The length of H r,l The number of columns.

[0179] Sort the eigenvalues ​​from large to small according to the energy-entropy ratio, and denote the basis vectors of the sorted eigenvalues ​​as W s1 、W s2 、W s3 , the coefficient vectors are recorded as H s1 、H s2 、H s3 , the phase matrices are recorded as The reconstructed time domain data are recorded as y s1 、y s2 、y s3 , and select the characteristic component with the largest energy-entropy ratio as the earthquake component. The sorting results of track 5 on April 9, 2016 are as follows Figure 4 、 Figure 5 As shown. The basis vector frequency W of the seismic component s1 It is mainly distributed around 0.04Hz, which is consistent with the ULF wave observed in low earth orbit, which is usually in the frequency range of Pc3 (20-100mHz). At the same time, the coefficient vector H of the seismic component s1 Only in the study area there are obvious anomalies, and the coefficient vectors H of the other characteristic components s2 Distributed throughout the track, H s3 The anomaly exists outside the study area, which is consistent with the law that the energy of earthquake components is mainly distributed within the earthquake-affected area.

[0180] Step 7: The RMS of the time-domain seismic component of each track is calculated as:

[0181]

[0182] The threshold is set to:

[0183] Thre=P×RMS (23)

[0184] Taking track 5 on April 9, 2016 as an example, the root mean square of the time domain earthquake component of this track is 0.0142, and the threshold is set to 3 times the root mean square, that is, Thre = 0.0426. Figure 5 seismic component y s1 As shown by the horizontal dotted line, the four data points whose time domain seismic components of the track are greater than Thre in the earthquake study area are marked as seismic anomaly points, and the track is also marked as an abnormal track.

[0185] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition, characterized in that: The following steps are involved: Step 1: Read the Swarm A satellite magnetic field data and split the daily magnetic field data into multiple tracks; Step 2: Subtract the CHAOS-7 model internal source field from the raw Y-component magnetic field data of each satellite orbit to remove the main magnetic field of the core and the lithospheric magnetic field of the crust to obtain the residual magnetic field; perform differential analysis on the residual magnetic field to remove the low-frequency background and obtain the magnetic field preprocessing data; Step 3: Perform short-time Fourier transform on the preprocessed data to obtain a complex time-frequency matrix in the rectangular coordinate representation of the complex domain; Step 4: Decompose the complex time-frequency matrix using complex non-negative matrix decomposition to obtain R eigencomponents, each of which contains a basis vector, a coefficient vector, and a phase matrix. Step 5: Obtain the time domain feature components by inverse short-time Fourier transform of each feature component; Step 6: Based on the characteristics of earthquake energy and spatiotemporal information being concentrated in the earthquake study area, the energy-entropy ratio of each characteristic component is calculated using the coefficient vector of each characteristic component of the orbit and the time domain characteristic component. The components are sorted in descending order and the characteristic component with the largest energy-entropy ratio is selected as the earthquake component. Step 7: By calculating the root mean square of the time domain seismic component of each track, using the exceedance rate method, setting the threshold, and extracting abnormal points in the track seismic component.

2. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition according to claim 1 is characterized in that: The step 1 comprises: Read the Swarm A satellite magnetic field data and split the daily magnetic field data into multiple tracks. Specifically, the magnetic field data is stored in daily units and is split with 50°S and 50°N as the orbit endpoints, resulting in 32 satellite tracks per day.

3. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition according to claim 1 is characterized in that: The step 2 includes: The original Y-component magnetic field data of each satellite orbit minus the source field of the CHAOS-7 model is calculated as: Among them, B y0 is the original Y component magnetic field, The Y component data of the internal source field calculated by the CHAOS-7 model, B y is the residual Y component magnetic field; The residual magnetic field is differentiated and calculated as: x(n)=B y (n+1)-B y (n) (2) Among them, B y (n) is the nth data point of the residual Y component magnetic field, B y The data length 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 decomposition according to claim 1 is characterized in that: The step 3 comprises: Perform short-time Fourier transform on the differential Y-component magnetic field x(n). The discrete short-time Fourier transform is: Where m represents the time index of the input signal x(n), g(n) is the selected window function, g(nm) 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 discrete short-time Fourier transform.

5. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition according to claim 4 is characterized in that: The step 4 comprises: The complex time-frequency matrix is ​​decomposed using complex non-negative matrix decomposition. Specifically, the complex time-frequency matrix V(n,ω) obtained by short-time Fourier transform is expressed as V with k rows and l columns. k×l , given the number of decomposition eigencomponents R, R is a positive integer and satisfies R<<min(k,l), the complex time-frequency matrix V k×l Decompose into matrix W k×R , matrix H R×l and The product of W k×R It is called the basis matrix, which represents the frequency distribution characteristics of the data. R×l It is called the coefficient matrix, which represents the weight of different frequency distribution features in time / space. Refers to the time-varying phase spectrum, which contains the phase matrix of k rows and l columns corresponding to R characteristic components, hereinafter abbreviated as V, W, H, The mathematical model of complex non-negative matrix decomposition is: in, Represents the Hadamard product of two matrices, whose elements are defined as the product of the corresponding elements of the two matrices, W k×r represents the rth column of the basis matrix W, H r×l represents the rth row of the coefficient matrix H, Represents the time-varying phase spectrum The r-th phase matrix of , 1≤r≤R.

6. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition according to claim 5 is characterized in that: Use KL divergence to measure the original matrix V and the reconstructed matrix The error between them, its objective function is: in, is the objective function, X k,r,l is the estimated value of the reconstructed complex time-frequency matrix of the rth eigencomponent, R(H r,l ) is a penalty term based on the p-norm. When 0<p<2, H r,l The sparsity of can be effectively controlled, and the penalty term is: RH r,l )=2λ|H r,l | p ,λ>0 (5) λ is the weight coefficient of the penalty term; Auxiliary function Z k,r,l , d k,r,l , A k,r,l and B k,r,l The calculation formula is as follows: Z k,r,l =|X k,r,l | (6) Finally, the iterative update rule when 0<p<1 is: Z k,r,l =|X k,r,l | (14) AND 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 number N iter (1) Initialize W, H, X k,r,l ; (2) Update the objective function variables using equations (11)-(17) and continue iterating until the value of the objective function equation (5) converges or reaches the set number of iterations N. iter , stop iteration; Output: W, H, After decomposition, R eigencomponents are obtained, and the rth eigencomponent contains the rth column vector W of the basis matrix W r , the rth row vector H of the coefficient matrix H r and the corresponding phase matrix H r Coefficient vector, W r Basis vectors.

7. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition according to claim 6 is characterized in that: The step 5 comprises: Each characteristic component is transformed by inverse short-time Fourier transform to obtain the time domain characteristic component, specifically: ISTFT() represents the inverse short-time Fourier transform function.

8. The satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition according to claim 1 is characterized in that: The step 6 comprises: The energy-entropy ratio of each characteristic component of the orbit consists of two parts: the energy ratio and the entropy ratio. The energy ratio refers to the ratio of the energy in the seismic study area of ​​the r component of the orbit to the entire orbit; the entropy ratio refers to the ratio of the Shannon entropy in the seismic study area of ​​the time domain characteristic component to the entire orbit. Shannon entropy Z in the earthquake study area i (y r ) and the Shannon entropy Z(y r ) is calculated as follows: Among them, the points between s1 and s2 are within the earthquake study area, and N is the time domain characteristic component y r The length of p(y r ) indicates that the random event Y is y r The probability of , s1, s2 are the minimum latitude and maximum latitude of the earthquake study area respectively; The energy-entropy ratio is then calculated as follows: Among them, the points between Ps and Pe are within the earthquake study area, and L is H r,l The number of columns, Ps and Pe are the minimum latitude and maximum latitude of the corresponding earthquake study area obtained by short-time Fourier transform; Sort the eigenvalues ​​from large to small according to the energy-entropy ratio, and denote the basis vectors of the sorted eigenvalues ​​as W s1 、W s2 ,…,W sR , the coefficient vectors are recorded as H s1 、H s2 ,…,H sR , the phase matrices are recorded as The reconstructed time domain data are recorded as y s1 、y s2 ,…,y sR , and the characteristic 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 decomposition according to claim 1 is characterized in that: The step 7 comprises: The root mean square of the time domain seismic component of each track is calculated as: The threshold is set to: Thre=P×RMS (22) P is the empirical parameter, Thre is the anomaly extraction threshold; When the amplitude of a time-domain seismic component data point in the earthquake study area is greater than Thre, the point is considered to be a seismic anomaly point, and the track is marked as an abnormal track.

Images