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

Abstract

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, which comprises the steps of reading Swarm A satellite magnetic field data, and removing a main magnetic field of a ground core and a lithosphere magnetic field of an earth crust to obtain a residual magnetic field; performing 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 using complex non-negative matrix decomposition to obtain R feature components; performing inverse short-time Fourier transform on each feature component to obtain a time domain feature component; a mean square root of a time domain seismic component of each orbit is calculated, a threshold value is set by using an over-limit rate method, and abnormal points are extracted from the seismic components of the orbits. The amplitude and phase information contained in the time-frequency domain of the satellite observation data is fully utilized, so that the component related to the earthquake is identified, and the earthquake anomaly extraction is 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 is 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 process of earthquake development and occurrence, and has become one of the effective ways to deepen our understanding of the earthquake development process. With the rapid development of space exploration technology, ionospheric detection technology, especially satellite detection, is increasingly used in earthquake precursor research, and satellite electromagnetic seismic anomaly research has become a hot topic in international geoscience research.

[0003] Satellite magnetic field observation data is not only subject to interference from the earth, but also strongly interfered by non-seismic factors such as solar activity and magnetic storms. In order to avoid non-seismic interference, traditional research methods use night data and remove data with high solar activity index and geomagnetic index. In order to improve the utilization rate of seismic data and extract seismic anomalies from satellite magnetic field observation data with the characteristics of "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 activity interference in the data. However, this method currently has certain limitations: the magnetic field signal contains not only waveform and amplitude information, but also other important information such as phase. This method is mostly to decompose the amplitude spectrum matrix after time-frequency analysis of the magnetic field signal, and reconstruct the components obtained by combining the original phase matrix back to the time domain, while ignoring the influence of the phase information of the original magnetic field signal on anomaly extraction during reconstruction, and cannot make full use of the complex frequency domain information to improve the accuracy of seismic magnetic field anomaly identification and extraction.

[0004] Therefore, in order to effectively utilize the 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 the satellite observation data in the time and frequency domain, thereby identifying the components related to earthquakes and making the 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 divide the daily magnetic field data into multiple tracks;

[0008] Step 2: The original Y-component magnetic field data of each satellite orbit is subtracted from the CHAOS-7 model internal source field, and the main magnetic field of the core and the lithospheric magnetic field of the crust are removed to obtain the residual magnetic field; the residual magnetic field is differentiated 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 form of rectangular coordinates in 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: 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 is calculated using the coefficient vector of each characteristic component of the orbit and the characteristic component in the time domain, 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;

[0013] Step 7: By calculating the RMS 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 comprises:

[0015] Read the Swarm A satellite magnetic field data and divide the daily magnetic field data into multiple tracks. Specifically, the magnetic field data is stored in units of days. It is divided with 50°S and 50°N as the orbit endpoints, and 32 satellite tracks are obtained every day.

[0016] Furthermore, the step 2 comprises:

[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 comprises:

[0024] The differential Y component magnetic field x(n) is subjected to short-time Fourier transform, and the discrete short-time Fourier transform is specifically:

[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 comprises:

[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 Decomposed into matrix W k×R , matrix H R×l and The product of 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 weights 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, which are abbreviated as V, W, H below. 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. 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 comprises:

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

[0052]

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

[0054] Furthermore, the step 6 comprises:

[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 denoted as H s1 , H s2 ,…,H sR , the phase matrices are recorded as The reconstructed time domain data are denoted as y s1 ,y s2 , …, y sR , and select the characteristic component with the largest energy-entropy ratio as the seismic component.

[0063] Furthermore, the step 7 comprises:

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

[0065]

[0066] The threshold value of 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 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.

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

[0071] The invention discloses a satellite magnetic field seismic anomaly extraction method based on complex non-negative matrix decomposition. The method removes the background field of satellite magnetic field data by subtracting a CHAOS-7 model to remove the endogenous field and then performing first-order difference. The method uses short-time Fourier transform to perform time-frequency transformation on preprocessed data with the background field removed to obtain a complex time-frequency matrix. According to the characteristics of the impact of non-seismic factors such as solar activity and magnetic storms on the global magnetic field and the impact of earthquakes on the local magnetic field, the method uses the advantages of local feature extraction of complex non-negative matrix decomposition to obtain multiple feature components, each of which includes a basis vector, a coefficient vector and a corresponding phase matrix, and makes full use of the amplitude and phase information contained in the satellite magnetic field signal in the time-frequency domain. The method uses inverse short-time Fourier transform to obtain reconstructed time domain feature components. According to the characteristics of seismic energy and spatiotemporal information distribution, the method uses energy-entropy ratio to identify seismic components. The method uses root mean square of time domain seismic components and uses an excess rate method to set a threshold value to extract seismic anomalies. Based on the characteristic that the anomaly distribution of seismic magnetic field 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 amplitude information and phase information in the complex domain to separate the anomaly features caused by seismic and non-seismic factors, making the extraction of seismic anomaly more accurate and reliable, thereby deepening the understanding of the anomaly phenomenon of ionospheric seismic magnetic field.

[0072] Based on the advantages of local feature extraction, the complex non-negative matrix decomposition method of the present invention can approximately 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, where W is called the basis matrix, which represents the frequency distribution characteristics of the data, and H is called the coefficient matrix, which represents the weights of different frequency distribution characteristics in time / space. After removing the background of the satellite magnetic field data, a short-time Fourier transform is performed to obtain the time-frequency spectrum matrix in the complex domain. The complex non-negative matrix decomposition is used to obtain the seismic component and reconstruct it back to the time domain using the inverse short-time Fourier transform. The seismic anomaly is extracted by the excess rate method and the abnormal orbit is determined. This method makes full use of 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 decomposing the complex non-negative matrix;

[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 solution 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 used 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 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, thereby removing 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 by 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 component by performing inverse short-time Fourier transform on each feature component;

[0085] Step 6: 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 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 amplitudes greater than the threshold in the seismic study area are considered abnormal points.

[0087] The step 1 comprises:

[0088] The magnetic field data of the Swarm A satellite is stored in units of days. It is divided with 50°S and 50°N as the orbit endpoints to obtain 32 satellite orbits per day.

[0089] The step 2 comprises:

[0090] The Y component of the satellite magnetic field is selected as the original magnetic field data. In order to remove the endogenous 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 first-order differentiated, 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 Decomposed into matrix W k×R , matrix H R×l and The product of 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 weights 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, which are abbreviated as V, W, H below. 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 to the complex matrix V as possible, 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 to the complex matrix V as possible. 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. 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 in the earthquake 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 earthquake study area of ​​the time domain characteristic component 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 denoted as H s1 , H s2 ,…,H sR , the phase matrices are recorded as The reconstructed time domain data are denoted as y s1 ,y s2 , …, y sR , and select the characteristic component with the largest energy-entropy ratio as the seismic component.

[0136] The step 7 comprises:

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

[0138]

[0139] The threshold is then 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 research 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 research area, is based on the Dobrovolsky formula R = 10 0.43M The determined square research area is centered on the epicenter and R is 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 the magnetic field data of Swarm A from February 16 to May 16, 2016. Since the data is stored in days, it needs to be divided into orbits with the geomagnetic latitudes 50°S and 50°N as endpoints. The period of Swarm A orbiting the earth is about 90 minutes, that is, one orbit can be obtained every 45 minutes, and 32 orbits can be obtained every day.

[0145] Step 2: Taking orbit 5 on April 9, 2016 as an example, the vector magnetic field includes X (north) component, Y (east) component, and Z (vertical) component. Since the magnetic field of the Y component may be affected by the activity of the rock layer, but is less affected by the external magnetic field interference, the Y component magnetic field data is used as the original magnetic field data. 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:

[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 first-order differentiated, 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] Among them, 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 complex time-frequency matrix is ​​decomposed using the advantages of local feature extraction of complex non-negative matrix decomposition 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 to the complex matrix V as possible. 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 orbit 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, 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 feature component back to the time domain through inverse short-time Fourier transform to obtain each time domain feature component. The formula of inverse short-time Fourier transform is as follows:

[0174]

[0175] Step 6: According to 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 in the earthquake study area are N, which is the time domain characteristic component y r length, L is 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 denoted as H s1 , H s2 , H s3 , the phase matrices are recorded as The reconstructed time domain data are denoted 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 in the earthquake-affected area.

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

[0181]

[0182] The threshold is then set to:

[0183] Thre=P×RMS (23)

[0184] Taking track 5 on April 9, 2016 as an example, the RMS of the time domain seismic component of this track is 0.0142, and the threshold is set to 3 times the RMS, 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 seismic study area are marked as seismic anomaly points, and the track is 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 protection scope 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 divide the daily magnetic field data into multiple tracks; Step 2: The original Y-component magnetic field data of each satellite orbit is subtracted from the CHAOS-7 model internal source field, and the main magnetic field of the core and the lithospheric magnetic field of the crust are removed to obtain the residual magnetic field; the residual magnetic field is differentiated 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 form of rectangular coordinates in 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: 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 is calculated using the coefficient vector of each characteristic component of the orbit and the characteristic component in the time domain, 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; Step 7: By calculating the RMS 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 divide the daily magnetic field data into multiple tracks. Specifically, the magnetic field data is stored in units of days. It is divided with 50°S and 50°N as the orbit endpoints, and 32 satellite tracks are obtained every 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 comprises: 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: The differential Y component magnetic field x(n) is subjected to short-time Fourier transform, and the discrete short-time Fourier transform is specifically: 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 Decomposed into matrix W k×R , matrix H R×l and The product of 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 weights 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, which are abbreviated as V, W, H below. 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. 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. 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, which is: 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: energy ratio and 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 whole 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 whole 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 denoted as H s1 , H s2 ,…,H sR , the phase matrices are recorded as The reconstructed time domain data are denoted as y s1 ,y s2 , …, y sR , and select the characteristic component with the largest energy-entropy ratio 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 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.

Images