CN110673206A - A seismic anomaly detection method for satellite magnetic field data based on non-negative matrix decomposition - Google Patents

Abstract

The present invention is a method for detecting anomalies of satellite magnetic field data based on non-negative matrix decomposition, comprising: removing invalid data and subtracting CHAOS-6 magnetic field model data according to a flag bit, and obtaining a first-order difference of the results to obtain differential data; The short-time Fourier transform is used to construct the non-negative time-frequency amplitude matrix corresponding to the magnetic field data; the non-negative time-frequency amplitude matrix is decomposed by the non-negative matrix decomposition method, and the local influence components caused by earthquakes and the components caused by solar activity and geomagnetic The global influence components generated by the activity are separated; the local influence components caused by the earthquake are selected according to the energy ratio, and the abnormal orbits of the components are judged by the method of the overrun rate; the daily abnormal orbits are accumulated, and the deviation from the background The degree of fitting a straight line detects seismic anomalies. The present invention can retain and utilize all the data obtained by measurement to study the earthquake, and at the same time obtain the components more relevant to the seismic activity to effectively perform seismic anomaly detection.

Description

A seismic anomaly detection method for satellite magnetic field data based on non-negative matrix decomposition

technical field

The invention belongs to the field of earthquake monitoring, in particular to a method for detecting anomalies of satellite magnetic field data based on non-negative matrix decomposition.

Background technique

Seismic anomalies refer to abnormal phenomena that occur before and after an earthquake, including geophysical anomalies such as geomagnetism, geoelectricity, and gravity near the epicenter. The researchers often realize the detection of earthquake anomalies by measuring the above-mentioned physical quantities and studying their anomalies. Measurement methods mainly include ground stations and satellite monitoring. Satellite monitoring methods have been widely studied because of their wide measurement range and proven ability to detect seismic anomalies. Among them, the magnetic field data measured by satellites has been extensively studied. At present, the influence of solar and geomagnetic activities on the data is excluded by selecting the nighttime data and simultaneously removing the data with higher solar activity index and geomagnetic index, and directly using the unremoved data for abnormal detection.

A method for seismic anomaly detection of satellite magnetic field data based on non-negative matrix decomposition is proposed. The non-negative matrix factorization method can decompose any high-dimensional non-negative matrix into two low-dimensional non-negative matrices W and H, so that the product of W and H approximately approximates the original high-dimensional matrix V, where W is called the basis matrix, and H is called is the weight matrix. The basis matrix W can be regarded as a set of basis of the original matrix V, and its column vector represents the characteristics of the data; and the weight matrix H is the projection of the matrix V in the basis vector matrix space, and its row vector represents the corresponding column vector of the basis matrix W with time. magnitude of change. The elements in the original matrix V are all non-negative numbers, and the elements in the decomposition matrices W and H are also non-negative numbers, which makes the decomposition result have real physical meaning. Based on this characteristic, the non-negative time-frequency matrix obtained by short-time Fourier transform after preprocessing the satellite Y-component magnetic field data is decomposed as the original matrix, and several time-frequency distribution characteristics of the data and their amplitudes in the time domain can be obtained. value. Separating the local influence components caused by earthquakes and the global influence components caused by solar activity and geomagnetic activity does not need to remove the data of days with strong solar and geomagnetic activity, and at the same time, the components more related to seismic activity can be obtained. Perform seismic anomaly detection.

CN101793972A discloses a technology for predicting satellite thermal infrared bright temperature anomalies in the short-term of strong earthquakes. It mainly uses polar orbiting satellites, stationary meteorological satellites and satellite constellation microwaves for observation to obtain satellite thermal infrared cloud images, and corrects through atmospheric models to obtain actual temperature values. , according to the topography and air mass shape, the interference information caused by topography and weather factors is excluded, and finally the relationship between the thermal infrared temperature anomaly characteristics and evolution law of earthquake precursors and the short-term forecast of the three elements of earthquakes is obtained.

Zeren Zhima et al. used the changing magnetic field data recorded by the DEMETER satellite to statistically study the disturbance characteristics of the space magnetic field before and after the strong earthquake of magnitude 7 or above in the northern hemisphere from 2005 to 2009. According to the research space and time range selected according to the earthquake magnitude and occurrence time, the nighttime data was selected to remove the influence of the sun, the satellite observation data with strong geomagnetic activity was excluded, and the remaining data were used to study the earthquake. The background field is constructed by using the observation data of the same period, and the seismic anomaly is observed by the disturbance amplitude of the magnetic field relative to the background field during the earthquake period.

Mehdi Akhoondzadeh et al. studied the Ecuador earthquake on April 16, 2016, studied the three-component data of the magnetic field of the Swarm satellite group Alpha star from November 1, 2015 to April 30, 2016, and calculated the median value of the daily data , through the quartiles of all median values to judge whether the day is an abnormal day, and remove the days with greater influence of geomagnetic activity in the results, and think that the remaining abnormal days may be caused by earthquakes.

The above technologies use traditional methods to remove the influence of solar and geomagnetic activities on the data, and need to remove some data and directly perform anomaly detection on the remaining data.

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is to provide a seismic anomaly detection method of satellite magnetic field data based on non-negative matrix decomposition, which can retain and use all the measured data to study earthquakes, and simultaneously obtain components more relevant to seismic activity. for seismic anomaly detection.

The present invention is realized in this way,

A method for seismic anomaly detection of satellite magnetic field data based on non-negative matrix decomposition, the method comprising:

a. Read the magnetic field data of the satellite, and remove the invalid data according to the flag bit as the original data;

b. Subtract the CHAOS-6 magnetic field model data from the original data, and obtain the differential data by taking the first-order difference of the results;

c. Perform short-time Fourier transform on the differential data of step b to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data;

d. Use the non-negative matrix decomposition method to decompose the non-negative time-frequency amplitude matrix of the magnetic field data, and separate the local influence components caused by earthquakes and the global influence components caused by solar activity and geomagnetic activity;

e. Select the local influence component caused by the earthquake according to the energy ratio, and judge the abnormal orbit of this component by the method of overrun rate;

f. Accumulate the number of abnormal orbits per day, and detect seismic anomalies by their deviation from the background fitting straight line;

g. Output deviation degree index curve result graph.

Further, in step a, the magnetic field data of the read satellite includes data of three components of the north X component B _x , the east Y component By and the vertical _{Z component B z .} Further, step b, specifically: subtracting the background field calculated by the CHAOS-6 magnetic field model from the Y-component magnetic field data of the original data to obtain residual data, and according to the research, the Y-component magnetic field data is most easily affected by changes in the lithosphere and reflects earthquakes. abnormal. Here, the Y-component magnetic field data of the raw data is selected.

Then take the first-order difference of the residual data to remove the remaining small background value to obtain the variation of the magnetic field data, which is calculated by formula (1):

i=1,2,3,...L-1

为减模型求一阶差分后的预处理结果。where By is the original _Y -component magnetic field data, B _m is the background field calculated by the CHAOS-6 magnetic field model, t _i is the time corresponding to the ith data, L is the length of each track data,

The preprocessing result after the first-order difference is calculated for the subtractive model.

Further, step c, performing short-time Fourier transform on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data includes:

Perform time-frequency transformation on the preprocessed magnetic field data through short-time Fourier transform to obtain a time-frequency matrix including amplitude information and phase information;

The time-frequency amplitude matrix formed by the amplitude values is taken as the input matrix of the subsequent non-negative matrix decomposition.

Further, step d specifically includes: decomposing any high-dimensional non-negative time-frequency amplitude matrix into two low-dimensional non-negative matrices W and H by a non-negative matrix decomposition method, so that the product of W and H is approximately approximating the original height. The dimension matrix V, where W is called the basis matrix, and H is called the weight matrix, and the time-frequency distribution characteristic matrix of the data and the amplitude of the time-frequency distribution characteristic matrix in the time domain are obtained;

The mathematical model of the non-negative matrix factorization method is:

V≈WH (6)

where V∈R ^≥0 is an M×N non-negative matrix, V∈R ^≥0 is an M×R non-negative matrix, V∈R ^≥0 is an R×N non-negative matrix, and the number of decomposition features The choice of R satisfies R<MN/(M+N);

The objective function of non-negative matrix factorization is formulated based on K-L divergence, where the objective function based on K-L divergence is:

i=1,2,3...M

j=1,2,3,...N

Wherein V _ij represents the i-th row and j column elements of the matrix V, and [WH] _ij represents the base matrix W and the weight matrix H multiplied to obtain the i-th row and j column elements of the matrix;

The weight matrix H is constrained by the minimum determinant criterion, and the volume of the vector stretched space of the weight matrix H is constrained to be the smallest. The solution formula for the volume of the vector stretched space of the weight matrix H is:

Then the K-L divergence non-negative matrix factorization objective function with the minimum determinant constraint added is:

i=1,2,3,...M

j=1,2,3,...N

where det represents the determinant of the matrix, and α is the balance parameter;

The objective function is optimized and solved by the multiplicative iteration method, the optimal solution is obtained by alternately updating the variables, and the update formula is calculated by the gradient descent method:

where _μw and _μh are respectively:

The final update iteration formula is:

表示矩阵元素的点乘；in

Represents the dot product of matrix elements;

The normalized formula is:

The iteration stop condition is set as the objective function is less than the set threshold Th, and the upper limit of the number of iterations is set as _Niter . If the value of the objective function is still not less than the set threshold after iteration _Niter times, the iteration is forcibly stopped;

The basis matrix W obtained by each orbital separation is:

Weight matrix H

for:

Further, step e, selects the local influence component due to the earthquake according to the energy ratio, and judges the abnormal orbit of this component by the method of overrun rate, specifically:

By calculating the ratio of the energy in the study area to the energy of the whole track in each row vector of the weight matrix H obtained by decomposing the non-negative time-frequency amplitude matrix, the component HH with the largest energy ratio is selected as the local influence component caused by the earthquake;

Then, the root mean square of the average energy level of the response data is selected as the threshold parameter, and the abnormal point is judged according to whether the value of the largest component HH is greater than k times the root mean square of the entire track. Orbits where the flag indicates that the anomaly is not due to known ionospheric activity are anomalous.

Further, in step f, the number of abnormal orbits per day is accumulated, and the seismic anomaly detected by the degree of its deviation from the background fitting straight line specifically includes:

The daily abnormal orbits and the number of all orbits per day are normalized by dispersion standardization,

The cumulative normalization result of the daily number of abnormal orbits is calculated as:

MN _ano = [MN _ano (1), MN _ano (2), MN _ano (3), ... MN _ano (d)] (26)

The cumulative normalized results of all orbits in each day are fitted with a straight line by least squares to obtain a background fitting straight line;

The cumulative normalization result of calculating the number of all tracks per day is:

MN _nom = [MN _nom (1), MN _nom (2), MN _nom (3), ... MN _nom (d)] (27)

The least squares linear fitting is performed on MN _nom , and the fitted straight line is:

y=a ₀ +a ₁ x (28)

The parameters a ₀ and a ₁ are estimated to be solved by:

和

where x is [1, 2, 3...d], d is the number of days, x _i is the i-th element in x; y is MN _nom , y _i is the i-th element in MN _nom ; x _i y _i is the i-th element in x multiplied by the i-th element in MN _nom ; after calculation, the parameters of the equation of the straight line are obtained

and

The seismic change process is displayed by the deviation of the cumulative curve of the number of abnormal orbits per day from the background fitting straight line. The calculation formula of the deviation from the linear fitting degree is:

i=1,2,3,...d

where Y _i is the i-th element in MN _ano ;

ρ is an index indicating the degree of deviation of the daily abnormal orbit cumulative curve from the background fitting straight line.

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

The invention is based on the satellite magnetic field data seismic anomaly detection method based on non-negative matrix decomposition. By subtracting the CHAOS-6 magnetic field model data and performing first-order difference, the magnetic field data background is removed to obtain the magnetic field data variation; the short-time Fourier transform is used to construct the magnetic field. The non-negative time-frequency amplitude matrix corresponding to the data; the constructed time-frequency amplitude matrix is decomposed through non-negative matrix decomposition, and the corresponding frequency distribution characteristic matrix and weight matrix are adaptively decomposed according to the time-frequency characteristics of the data. The local influence component generated is separated from the global influence component caused by solar activity and geomagnetic activity; the local influence component caused by earthquake is selected for research through the energy ratio, and the abnormal orbit of this component is judged by the method of overrun rate; finally The daily number of abnormal orbits is accumulated, and the seismic anomaly is detected by the degree of its deviation from the linear fit. The invention decomposes the non-negative time-frequency amplitude matrix of the data through non-negative matrix decomposition, and adaptively separates the interference of the sun and geomagnetic activities on the data according to the time-frequency characteristics of the magnetic field data, and obtains its frequency distribution characteristics and time variation amplitude. , and its decomposition result has real physical meaning. It is possible to retain and use all the measured data for earthquake research, use components more related to seismic activity to effectively detect seismic anomalies, make up for the above technology to remove data when the sun and geomagnetic activity are strong, reduce data usage, And use the remaining data within the insufficiency to directly perform anomaly detection.

Description of drawings

Fig. 1 is the flow chart of the seismic anomaly detection method of satellite magnetic field data based on non-negative matrix decomposition;

Figure 2 is a schematic diagram of the Ecuador earthquake epicenter, affected area and study area;

Fig. 3 is the original curve of satellite magnetic field Y component data;

Figure 4 shows the original data of the Y component of the satellite magnetic field and the CHAOS-6 model data curve, the residual curve and the first-order difference curve; the solid line in Figure 4(1) is the original data curve, and the dotted line is the model data curve. The amplitude is large and the difference between the model data is very small. Figure 4(2) enlarges the part of Figure 4(1). Figure 4(3) shows the residual data obtained by subtracting the model data from the original data. Figure 4(4) is the first-order difference data of the residual data;

Fig. 5 is a result diagram of short-time Fourier transform time-frequency amplitude matrix of satellite Y-component magnetic field differential data;

Fig. 6 is non-negative matrix decomposition frequency distribution characteristic matrix W result graph;

Fig. 7 is a non-negative matrix decomposition weight matrix H result diagram;

Figure 8 is the time-domain reconstruction result of the non-negative matrix decomposition result, Figure 8(1) is the first-order difference data, and there are obvious singular values both inside and outside the study area; Figure 8(2) is the time domain components corresponding to W1 and H1 , there are obvious singular values only in the study area; Fig. 8(3) is the time domain component corresponding to W2 and H2, there is no obvious singular value, and it is evenly distributed in the whole orbit; Fig. 8(4) is the corresponding W3 and H3 the time domain component of ;

Figure 9 is the result of abnormal orbit judgment. Figure 9(1) is the orbit 3 on April 9, and no point exceeds the threshold, so it is not judged as an abnormal orbit; Figure 9 (2) is the orbit 1 on April 10, and there is an exceeding Threshold points, which occur only in the study area, and which do not show known ionospheric activity by checking their flags, can be judged as abnormal orbits;

Figure 10 is a graph showing the result of daily abnormal orbit accumulation and background fitting straight line;

Figure 11 is the index curve of the degree of deviation of the daily abnormal orbit cumulative curve from the background fitting straight line.

Detailed ways

In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention, but not to limit the present invention.

Referring to Fig. 1, a method for detecting seismic anomalies in satellite magnetic field data based on non-negative matrix decomposition, the method includes:

a. Read the magnetic field data of the satellite, and remove the invalid data according to the flag bit;

b. The original data is the data after removing the invalid data minus the CHAOS-6 magnetic field model data, and the first-order difference is calculated for the result;

c. Perform short-time Fourier transform on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data;

d. Use the non-negative matrix decomposition method to decompose the time-frequency amplitude matrix of the magnetic field data is a non-negative time-frequency matrix, and separate the local influence components caused by earthquakes and the global influence components caused by solar activity and geomagnetic activity;

e. Select the local influence component caused by the earthquake according to the energy ratio, and judge the abnormal orbit of this component by the method of overrun rate;

f. Accumulate the number of abnormal orbits per day, and detect seismic anomalies by their deviation from the background fitting straight line;

g. Output deviation degree index curve result graph.

In step a, the magnetic field data of the satellite is read, including data of three components of the north X component B _x , the east Y component By and the vertical _{Z component B z .}

Step b, subtract the CHAOS-6 magnetic field model data from the original data, and calculate the first-order difference of the result, which is the background field calculated by subtracting the CHAOS-6 magnetic field model from the original Y-component magnetic field data. First, use the original data to subtract the CHAOS-6 For the background field calculated by the magnetic field model, remove the larger background in the magnitude of the geomagnetic field to obtain the residual data, and then take the first-order difference of the residual data to remove the remaining small background value to obtain the variation of the magnetic field data, and the result can be used For subsequent operations such as short-time Fourier transform.

i=1,2,3,...L-1

where By is the original _Y -component magnetic field data, B _m is the background field calculated by the CHAOS-6 magnetic field model, t _i is the time corresponding to the i-th data, L is the length of each track data, and dB is the first-order subtraction model The preprocessing result after the difference.

Step c, performing short-time Fourier transform on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data, including:

The time-frequency matrix is obtained by performing time-frequency transform on the preprocessed magnetic field data through short-time Fourier transform;

The time-frequency amplitude matrix formed by its amplitude is taken as the input matrix of the subsequent non-negative matrix decomposition. Non-negative matrix decomposition requires that each element in the input matrix is a non-negative value. Since the Y-component magnetic field data is single-channel data, it cannot form a matrix, and there are negative values in the preprocessing results. The data is time-frequency transformed, and a non-negative time-frequency amplitude matrix is constructed as the input matrix of the non-negative matrix decomposition. And the short-time Fourier transform has relevant physical meaning, and the result reflects the time-frequency characteristics of the data, that is, the change of the amplitude and phase distribution of the data frequency with time. Here, the change of the amplitude distribution of the data frequency with time is used. , forming a non-negative matrix for decomposition, and subsequent decomposition can be performed according to the time-frequency amplitude characteristics of the data.

Let the time series of each track data after preprocessing be:

X=[x ₁ , x ₂ , x ₃ , x ₄ , ..., x _L-1 ] (2)

Where L is the length of each track data.

The discrete form of the short-time Fourier transform is:

n=0,1,2,3,...N-1

k=0,1,2,...L-1

Where X is the time series of the data, g[l] is the window function, * is the conjugate, the length of the window function is M, the step size is S, N is the total number of windows, L is the total number of discrete frequency points, Finally, the complex number matrix of M×N is obtained as its time-frequency matrix:

相位为 Each element in the matrix is a complex number z _mn =a _mn +jb _mn , including amplitude information and phase information, and its amplitude is

Phase is

Take its amplitude information to obtain the M×N non-negative time-frequency amplitude matrix as:

Step d, using the non-negative matrix decomposition method to decompose the time-frequency amplitude matrix of the magnetic field data, in order to use the non-negative matrix decomposition method to adaptively decompose the non-negative time-frequency amplitude matrix of the data according to the time-frequency characteristics of the data to obtain the corresponding The weight matrix corresponding to the frequency distribution characteristic matrix and the frequency distribution characteristic matrix (the column vector of the W matrix corresponds to the row vector of the H matrix) separates the local influence components caused by earthquakes and the global influence components caused by solar activity and geomagnetic activity.

The non-negative matrix factorization method can decompose any high-dimensional non-negative matrix into two low-dimensional non-negative matrices W and H, so that the product of W and H approximately approximates the original high-dimensional matrix V, where W is called the basis matrix, and H is called is the weight matrix. The basis matrix W can be regarded as a set of basis of the original matrix V, and its column vector represents the characteristics of the data; and the weight matrix H is the projection of the matrix V in the basis vector matrix space, and the row vector of the weight matrix represents the column vector corresponding to the basis matrix W. Amplitude over time. The elements in the original matrix V are all non-negative numbers, and the elements in the decomposition matrices W and H are also non-negative numbers, which makes the decomposition result more physical. Based on this characteristic, the non-negative time-frequency matrix corresponding to the magnetic field data can be decomposed as the original matrix, and several time-frequency distribution characteristics of the data and their amplitudes in the time domain can be obtained. By separating the local influence components due to earthquakes and the global influence components due to solar activity and geomagnetic activity, it is not necessary to remove the data of days with strong solar and geomagnetic activity, and at the same time, components more related to seismic activity can be obtained.

The mathematical model of non-negative matrix factorization is:

V≈WH (6)

where V∈R ^≥0 is an M×N non-negative matrix, W∈R ^≥0 is an M×R non-negative matrix, H∈R ^≥0 is an R×N non-negative matrix, and the number of decomposition features The selection of R should generally satisfy R<MN/(M+N), which can be set according to the known information or empirical information of the data.

The present invention formulates the objective function of non-negative matrix decomposition based on K-L divergence, and at the same time adds determinant constraint to the weight matrix H, so that the decomposition result is relatively unique. The objective function based on generalized K-L divergence is:

i=1,2,3...M

j=1,2,3,...N

Wherein V _ij represents the i-th row and j column elements of the matrix V, and [WH] _ij represents the i-th row and j column elements of the matrix obtained by multiplying the W matrix and the H matrix.

Since V can be decomposed into an infinite number of W and H, in order to make the decomposition result relatively unique, the weight matrix H is constrained by the minimum determinant criterion, and the minimum volume of the space formed by the vector of the weight matrix H is constrained. The formula for solving the volume of the vector spanning space of the weight matrix H is:

Then the K-L divergence non-negative matrix factorization objective function with the minimum determinant constraint added is:

i=1,2,3,...M

j=1,2,3,...N

where det represents the determinant of the matrix, and α is the balance parameter.

The objective function is optimized by the multiplicative iterative method, and the optimal solution is obtained by alternately updating the variables. The update formula is calculated using the gradient descent method:

where _μw and _μh are respectively:

The final update iteration formula is:

表示矩阵元素的点乘。in

Represents the dot product of matrix elements.

The normalized formula is:

Since the addition of the minimum determinant constraint cannot guarantee the non-negativity of W and H, the negative values in the W and H matrices are set to 0 to ensure their non-negativity.

The above iteration formula satisfies that the more iterations, the more the objective function converges. The iteration stop condition is set as the objective function is less than the set threshold Th, and the upper limit of the number of iterations is set as _Niter , that is, if the value of the objective function is still not less than the set threshold after iteration _Niter times, the iteration is forcibly stopped. Since the values of the time-frequency amplitude matrix obtained from different orbits are different, the threshold Th should be set according to the values of different orbits:

where ε is the error factor.

The W matrix obtained from each orbital separation is:

The H matrix is:

Step e, select the local influence component caused by the earthquake according to the energy ratio, and judge the abnormal orbit of this component by the method of the overrun rate, which is the energy in the study area in each row vector of the H matrix obtained by calculating the non-negative matrix decomposition Compared with the energy ratio of the entire track, the component HH with the largest energy ratio is selected as the local influence component caused by the earthquake, and the root mean square that can reflect the average level of data energy is selected as the threshold parameter. The k times of the square root are selected by experience to determine the abnormal point, where k is selected by experience, and it is considered that the abnormal point only appears in the study area, and the corresponding sign shows that the abnormal point is not caused by the known ionospheric activity. The orbit is an abnormal orbit . Since the impact of seismic events on data is a local impact, it only has a certain impact on the data in the seismic research area in the orbit. At this time, the data energy is mainly concentrated in the earthquake-affected area; while solar activity and geomagnetic activity are global impact events. Therefore, it will affect the data of the entire track. At this time, the data should be evenly distributed on the entire track. Therefore, the component with the largest ratio of the energy in the study area to the energy of the entire track is selected. The data energy of this component is mainly concentrated in the seismic study area, and most likely reflects the impact of the earthquake on the local area.

The formula for calculating the energy ratio is:

i=1,2,3,...R

j=1,2,3,...N

Among them, P _s represents the starting position of the data in the study area, and P _e represents the end position of the data in the study area. R energy ratio results [η ₁ , η ₂ , η ₃ , ... η _R ] can be obtained by calculation, and the H row vector HH corresponding to the maximum energy ratio is selected.

The abnormal orbit of the selected component HH is judged by the over-limit rate method, and the root mean square that can reflect the average level of data energy is selected as the threshold parameter. The abnormal data should far exceed this value and last for a period of time. Since the data has been smoothed when the original data is subjected to short-time Fourier transform, whether the value of HH is greater than k times the root mean square of the entire track can be directly used to judge abnormal points. The formula for calculating the root mean square of a time series X _i of length N is:

Then the threshold is kX _rms , and the points larger than the threshold are abnormal points. At the same time, it is considered that abnormal points only appear in the study area, and the corresponding flag position shows that the abnormal orbit is not caused by known ionospheric activities.

In step f, the number of abnormal orbits per day is accumulated, and the seismic anomaly is detected by the degree of its deviation from the background fitting straight line. The least squares straight line fitting is used as the background, and the variation process of seismic anomalies is displayed by the degree of deviation of the cumulative results of daily anomalous orbits from the background fitting straight line. Because the abnormal orbit judgment of all the orbits that meet the research conditions cannot directly show the connection between the orbit and the earthquake, the daily number of abnormal orbits is accumulated, and the minimum number of all orbits per day is used. The square line fitting is used as the background, and the degree of deviation of the cumulative results of the daily abnormal orbits from the background line is used to reflect the impact of earthquakes on the data. The greater the deviation, the greater the impact.

The daily number of abnormal orbits is accumulated, and the result is:

_Nano = [ _Nano (1), _Nano (2), _Nano (3), ... _Nano (d)] (23)

where d is the number of days studied.

Then accumulate the number of all tracks per day, and the result is:

N _nom = [N _nom (1), N _nom (2), N _nom (3), ... N _nom (d)] (24)

Since the number of all orbits per day is quite different from the number of abnormal orbits per day, the two are normalized by dispersion normalization, where the dispersion normalization formula is:

The cumulative normalization result of the daily number of abnormal orbits is calculated as:

MN _ano = [MN _ano (1), MN _ano (2), MN _ano (3), ... MN _ano (d)] (26)

The cumulative normalization result of the number of all tracks per day is:

MN _nom = [MN _nom (1), MN _nom (2), MN _nom (3), ... MN _nom (d)] (27)

Least squares linear fitting is performed on MN _nom , and the fitting line is set as:

y=a ₀ +a ₁ x (28)

The parameters a ₀ and a ₁ are estimated to be solved by:

和

where x is [1, 2, 3...d], d is the number of days, x _i is the i-th element in x; y is MN _nom , y _i is the i-th element in MN _nom ; x _i y _i is the i-th element in x multiplied by the i-th element in MN _nom ; after calculation, the parameters of the equation of the straight line are obtained

and

The seismic change process is displayed by the deviation of the cumulative curve of the number of abnormal orbits per day from the background fitting straight line. The calculation formula of the deviation from the linear fitting degree is:

i=1,2,3,...d

where Y _i is the i-th element in MN _ano ;

ρ is an index indicating the degree of deviation of the daily abnormal orbit cumulative curve from the background fitting straight line.

Step g, outputting the result graph of the deviation degree index curve, is to output the curve representing the index change over time of the degree of deviation of the daily abnormal track cumulative result from the background fitting straight line as the final seismic abnormality detection result. Since ρ is an index indicating the degree of deviation of the daily abnormal track cumulative curve from the background fitting straight line, the larger the data, the greater the abnormality caused by the earthquake. Therefore, it can represent the abnormal change of seismic activity and output it as the seismic abnormality detection result. .

For the magnitude 7.8 Ecuador earthquake that occurred on April 16, 2016, the Y-component magnetic field data (1Hz) of the Swarm satellite group Alpha star was taken as an example. According to Dobrovilsky's formula R=10 ^0.43M (M is the magnitude of the earthquake), the study area of the Ecuador earthquake was selected, and the rectangular area with the epicenter as the center of ±20° longitude and latitude was the study area, and the study time was from 60 days before the earthquake to the earthquake The last 30 days, from February 16 to May 16, 2016. The Ecuador earthquake epicenter, earthquake-affected area, and study area are shown in Figure 2, where the five-pointed star is the epicenter, the circular area is the earthquake-affected area, and the square area is the study area.

a. Read and enter the _Y -component magnetic field data By (1Hz) of the Swarm satellite group Alpha star from January 17, 2016 to May 16, 2016, and remove the invalid data caused by the wrong measurement of the instrument itself according to the flag bit in the data set . According to the particularity of satellite data measurement, that is, the measurement value changes with time and space, and the local influence caused by earthquakes and the global components caused by the sun and geomagnetic activities must be separated in the future, so the geomagnetic latitude is used in the range of +50°～-50° The orbits within the range of ° and passing through the seismic research area are used as the research object for data processing.

Taking orbit 2 on April 8, 2016 as an example, the original curve of the Y-component magnetic field data is shown in Figure 3. It can be seen that the amplitude of the original magnetic field data is very large, and its changes cannot be observed, and it is difficult to directly decompose it. or exception extraction and other processing.

b. Use the original Y-component magnetic field data to subtract the background field calculated by the CHAOS-6 magnetic field model, remove the background with large amplitude in the geomagnetic field, and take the first-order difference of the result, remove the remaining small background value to obtain the magnetic field data amount of change.

The formula for subtracting the model magnetic field from the raw data and differencing it is:

i=1,2,3,...L-1

where By is the original _Y -component magnetic field data, B _m is the background field calculated by the CHAOS-6 magnetic field model, t _i is the time corresponding to the i-th data, L is the length of each track data, and dB is the first-order subtraction model The preprocessing result after the difference.

Taking track 2 on April 8, 2016 as an example, the original data and model data curves, the partial amplification curves of original data and model data, the residual data curves of original data minus model data, and the first-order difference data curves are shown in Figure 4. (1)(2)(3)(4). The solid line in Fig. 4(1) is the original data curve, and the dotted line is the model data curve. Since the original data has a large amplitude and the model data has a small difference with it, Fig. 4(2) compares Fig. 4(1). Partially enlarged, Fig. 4(3) is the residual data obtained by subtracting the model data from the original data, and Fig. 4(4) is the first-order difference data of the residual data. It can be seen from the results that the amplitude of the original data is very large, and the background field with a large amplitude can be removed by subtracting the magnetic field model, but the residual data still has a certain background amplitude, which needs to be further removed. Therefore, the first-order difference is calculated for the residual data, the residual background amplitude is removed, and the variation of the magnetic field data is obtained. The final result obtained is that the background amplitude is very small, and the singular value is more prominent.

c. Perform short-time Fourier transform on the preprocessed magnetic field data to obtain its time-frequency matrix, and take its amplitude to form a non-negative time-frequency amplitude matrix, which is used as the input matrix for subsequent non-negative matrix decomposition, using the data The time-frequency amplitude characteristics are decomposed.

Let the time series of each track data after preprocessing be:

X=[x ₁ , x ₂ , x ₃ , x ₄ , ..., x _L-1 ] (2)

Among them, L is the length of each track data, and the value is different according to the length of each track.

The discrete form of the short-time Fourier transform is:

n=0,1,2,3,...N-1

k=0,1,2,...L-1

Where x[l] represents X is the time series of the data, g[l] is the window function, Hanning window is used to reduce the spectral leakage, * represents the conjugate, the length of the window function is M, take M=64, the step size is S, take S=16, N is the total number of windows, and L is the total number of discrete frequency points. Finally, the complex number matrix of M×N is obtained as its time-frequency matrix:

相位为

Each element in the matrix is a complex number z _mn =a _mn +jb _mn , including amplitude information and phase information, and its amplitude is

Phase is

Take its amplitude information to obtain the M×N non-negative time-frequency amplitude matrix as:

Taking orbit 2 on April 8, 2016 as an example, the time-frequency amplitude matrix result of the short-time Fourier transform is shown in Figure 5, where the abscissa is the geographic latitude, the ordinate is the frequency, and the gray scale represents its amplitude. value size. It can be seen that this step can convert the single-channel magnetic field data into an M×N matrix. At the same time, the matrices are all non-negative numbers, which can be used as the input matrix for non-negative matrix decomposition, and has its physical meaning and reflects the time-frequency amplitude of the signal. characteristic, and subsequent processing can be decomposed according to this characteristic.

d. Use the non-negative matrix decomposition method to decompose the time-frequency amplitude matrix of the magnetic field data, and use the non-negative matrix decomposition method to adaptively decompose the non-negative time-frequency amplitude matrix of the data according to the time-frequency characteristics of the data to obtain the corresponding The frequency distribution characteristic matrix and its corresponding weight matrix separate the local influence components caused by earthquakes and the global influence components caused by solar activity and geomagnetic activity.

The non-negative matrix factorization method can decompose any high-dimensional non-negative matrix into two low-dimensional non-negative matrices W and H, so that the product of W and H approximately approximates the original high-dimensional matrix V, where W is called the basis matrix, and H is called is the weight matrix. The basis matrix W can be regarded as a set of basis of the original matrix V, and its column vector represents the characteristics of the data; and the weight matrix H is the projection of the matrix V in the basis vector matrix space, and its row vector represents the corresponding column vector of the basis matrix W with time. magnitude of change. The elements in the original matrix V are all non-negative numbers, and the elements in the decomposition matrices W and H are also non-negative numbers, which makes the decomposition result more physical. Based on this characteristic, the non-negative time-frequency matrix corresponding to the magnetic field data can be decomposed as the original matrix, and several time-frequency distribution characteristics of the data and their amplitudes in the time domain can be obtained. By separating the local influence components due to earthquakes and the global influence components due to solar activity and geomagnetic activity, it is not necessary to remove the data of days with strong solar and geomagnetic activity, and at the same time, components more related to seismic activity can be obtained.

The mathematical model of non-negative matrix factorization is:

V≈WH (6)

where V∈R ^≥0 is an M×N non-negative matrix, W∈R ^≥0 is an M×R non-negative matrix, H∈R ^≥0 is an R×N non-negative matrix, and the number of decomposition features The selection of R should generally satisfy R<MN/(M+N). According to experience and known information, take R=3.

The objective function of non-negative matrix decomposition is formulated based on K-L divergence, and the determinant constraint is added to the weight matrix H to make the separation result relatively unique. The objective function based on generalized K-L divergence is:

i=1,2,3...M

j=1,2,3,...N

Wherein V _ij represents the i-th row and j column elements of the matrix V, and [WH] _ij represents the i-th row and j column elements of the matrix obtained by multiplying the W matrix and the H matrix.

Since V can be decomposed into an infinite number of W and H, in order to make the decomposition result relatively unique, the weight matrix H is constrained by the minimum determinant criterion. Constrained by the minimum volume of the space spanned by the vectors of the matrix H. The formula for solving the volume of the vector spanned space of the matrix H is:

Then the K-L divergence non-negative matrix factorization objective function with the minimum determinant constraint added is:

i=1,2,3...M

j=1,2,3,...N

where det represents the determinant of the matrix, α is the balance parameter, and α=0.005.

The objective function is optimized by the multiplicative iterative method, and the optimal solution is obtained by alternately updating the variables. The update formula is calculated using the gradient descent method:

where _μw and _μh are respectively:

The final update iteration formula is:

in Represents the dot product of matrix elements.

The normalized formula is:

Since the addition of the minimum determinant constraint cannot guarantee the non-negativity of W and H, the negative values in the W and H matrices are set to 0 to ensure their non-negativity.

The iterative formula satisfies that the more the number of iterations, the more the objective function converges. Set the iteration stop condition as the objective function is less than the set threshold Th, and set the upper limit of the number of iterations to N _iter , take N _iter =500, that is, if the value of the objective function is still not less than the set threshold after iterative N _iter times, it will be forced to stop iterate. Since the values of the time-frequency amplitude matrix obtained from different orbits are different, the threshold Th should be set according to the values of different orbits:

Among them, ε is the error factor, which is taken as ε=0.05.

The W matrix obtained from each orbital separation is:

The H matrix is:

Taking orbit 2 on April 8, 2016 as an example, the decomposition results W matrix and H matrix are shown in Figure 6 and Figure 7, respectively, where W1, W2, and W3 show the three frequency distribution characteristics of the orbit, H1, H2 and H3 are their corresponding amplitudes, and between the dotted lines in the H result graph is the data in the study area. Among them, the larger amplitude of H1 only appears in the study area, and the energy is mainly distributed in the study area. This component should be the local component affected by the earthquake; the amplitude of H2 does not fluctuate much, and the energy is evenly distributed in the entire track , it may be the background component; the maximum amplitude of H3 occurs outside the study area, and there are also large amplitude values in the study area, which may be the global component affected by the sun or geomagnetic activity. In order to show its separation effect from the time domain, the decomposition result is reconstructed to the time domain result as shown in Figure 8. Among them, Figure 8(1) is the first-order difference data, and there are obvious singular values both inside and outside the study area; Figure 8(2) is the time domain components corresponding to W1 and H1, and there are obvious singular values only in the study area; Figure 8( 3) is the time domain component corresponding to W2 and H2, there is no obvious singular value, and is evenly distributed in the entire orbit; Figure 8(4) is the time domain component corresponding to W3 and H3, the largest singular value appears outside the study area, There are also singular values with small amplitudes in the study area. From the decomposition results, it can be seen that the non-negative time-frequency amplitude matrix of the data is decomposed by the non-negative matrix decomposition method, and its frequency distribution characteristics and time variation amplitude can be obtained. The resulting global impact components are separated and components more correlated with seismic activity are utilized for seismic anomaly detection.

e. Calculate the energy ratio between the H-matrix row vector seismic study area and the entire orbit, select the component with the largest energy ratio as the local influence component caused by the earthquake, and use the method of overrun rate to judge the abnormal orbit of this component.

The formula for calculating the energy ratio is:

i=1, 2, 3

j=1,2,3,...N

Among them, P _s represents the starting position of the data in the study area, and P _e represents the end position of the data in the study area. Three energy ratio results [η ₁ , η ₂ , η ₃ ] are obtained by calculation, and the H row vector HH corresponding to the maximum energy ratio is selected.

Compute the root mean square value of the row vector HH:

Take k=5, set its threshold as kHH _rms , the points in HH that exceed the threshold can be regarded as abnormal points, and after obtaining the abnormal points, judge whether they only appear in the research area, and the corresponding flag indicates that the abnormality is not due to known Anomalies caused by ionospheric activities. If there are only anomalous points that meet the above conditions, the orbit is an anomalous orbit.

Taking orbit 2 on April 8, 2016 as an example, the decomposition results H are shown in Figure 7, respectively. The energy ratio of H1 is 0.924; the energy ratio of H2 is 0.624; and the energy ratio of H3 is 0.4095. Calculated according to the ratio of the length of the data in the study area to the total data, the energy ratio should be 0.4 when the data is evenly distributed, so H1 is the local component caused by the earthquake, H3 should be the global component caused by the influence of the sun or geomagnetic activity, and H2 is Background components or other distinguishing quantities. According to statistics, the average maximum energy ratio of all orbits is 0.57, an increase of 45% compared with the average value, and the average maximum energy ratio of all abnormal orbits is 0.80; the average value of the second largest energy ratio is 0.39, which is basically the same as the average value of 0.4, and the average value of the minimum energy ratio is 0.39. was 0.23, down 42% from the average. It can be seen that the components of the decomposition result H of the non-negative matrix decomposition can basically reflect the local components caused by the earthquake and the global components caused by the influence of the sun and geomagnetic activities.

Taking orbit 3 on April 9, 2016 and orbit 1 on April 10, as examples, the schematic diagram of abnormal orbit judgment is shown in Figure 9. The solid line in the figure is the HH component data, the dotted line is the threshold value, and the middle of the dotted line is the study area internal data. Figure 9(1) is the orbit 3 on April 9, and no point exceeds the threshold, so it is not judged as an abnormal orbit; Figure 9 (2) is the orbit 1 on April 10, there are points that exceed the threshold, and only in the study area It can be determined to be an anomalous orbit by checking that its flags show no known ionospheric activity.

f. Accumulate the daily number of abnormal orbits, and the result is:

_Nano = [ _Nano (1), _Nano (2), _Nano (3), ... _Nano (d)] (23)

where d is the number of days studied and d=121.

Then accumulate the number of all tracks per day, and the result is:

N _nom = [N _nom (1), N _nom (2), N _nom (3), ... N _nom (d)] (24)

The dispersion standardization is performed on the two curves to achieve normalization, and the dispersion formula is used:

The cumulative normalization result of the daily number of abnormal orbits is calculated as:

MN _ano = [MN _ano (1), MN _ano (2), MN _ano (3), ... MN _ano (d)] (26)

The cumulative normalization result of the number of all tracks per day is:

MN _nom = [MN _nom (1), MN _nom (2), MN _nom (3), ... MN _nom (d)] (27)

Least squares linear fitting is performed on MN _nom , and the fitting line is set as:

y=a ₀ +a ₁ x (28)

The parameters a ₀ and a ₁ are estimated to be solved by:

和

where x is [1, 2, 3...121], y is MN _nom , the parameters of the straight line equation are obtained by calculation

and

和

where x is [1, 2, 3...d], d is the number of days, x _i is the i-th element in x; y is MN _nom , y _i is the i-th element in MN _nom ; x _i y _i is the i-th element in x multiplied by the i-th element in MN _nom ; after calculation, the parameters of the equation of the straight line are obtained

and

The seismic change process is displayed by the deviation of the cumulative curve of the number of abnormal orbits per day from the background fitting straight line. The calculation formula of the deviation from the linear fitting degree is:

i=1, 2, 3, ... 121

where Y _i is the i-th element in MN _ano ;

ρ is an index indicating the degree to which the cumulative curve of the abnormal orbit deviates from the background fitting straight line.

Figure 10 shows the normalized results of the cumulative curve of daily anomalous orbits and the background fitting straight line, where the vertical dotted line is the magnitude 7.8 Ecuador earthquake on April 16. It can be seen from the results that the cumulative curve of the number of abnormal orbits per day began to deviate from the background fitting line abnormally from a period before the earthquake, and reached the maximum deviation near the earthquake. It began to recover a few days after the earthquake, and finally returned to calm. Figure 11 shows the deviation of the daily abnormal track cumulative curve from the background fitting straight line degree index curve, in which the black rectangle is the abnormality detected before the earthquake. The cumulative results of daily anomalous orbits began to appear abnormal about 20 days before the earthquake, began to increase significantly 10 days before the earthquake, reached the maximum two days after the earthquake, and then gradually decreased until it returned to normal levels about 20 days after the earthquake. The results show that the method can effectively detect the anomalies generated by the target earthquake and the changes of anomalies with seismic activity, including the onset of anomalies before the earthquake, the peak of the anomaly two days after the earthquake, and the recovery of the anomaly after the earthquake.

In step g, the curve representing the variation of the index with time indicating the degree of deviation of the daily abnormal track cumulative result from the background fitting straight line is output as the final seismic abnormality detection result.

The above descriptions are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention shall be included in the protection of the present invention. within the range.

Claims (7)

1. A satellite magnetic field data earthquake abnormity detection method based on non-negative matrix factorization is characterized by comprising the following steps:

a. reading magnetic field data of the satellite, and removing invalid data according to the zone bits to be used as original data;

b. subtracting CHAOS-6 magnetic field model data from the original data, and solving a first-order difference of the result to obtain differential data;

c. b, carrying out short-time Fourier transform on the differential data in the step b to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data;

d. decomposing a non-negative time-frequency amplitude matrix of the magnetic field data by using a non-negative matrix decomposition method, and separating local influence components generated by earthquake and global influence components generated by solar activity and geomagnetic activity;

e. selecting local influence components generated by earthquake according to the energy ratio, and judging abnormal orbits of the components by an overrun method;

f. accumulating the number of abnormal orbits every day, and detecting earthquake abnormity according to the deviation degree of the abnormal orbits from the background fitting straight line;

g. and outputting a deviation degree exponential curve result graph.

2. The method of claim 1,

step a, reading the magnetic field data of the satellite, wherein the magnetic field data comprises a north X component B_xEast Y component B_yAnd a vertical Z component B_zData of three components.

3. The method of claim 1,

step b, specifically: subtracting the background field calculated by the CHAOS-6 magnetic field model from the Y component magnetic field data of the original data to obtain residual data,

and then, taking first-order difference to the residual data to remove the residual small background value to obtain the variation of the data of the magnetic field, and calculating by adopting a formula (1):

wherein B is_yAs raw Y-component magnetic field data, B_mBackground field calculated for the CHAOS-6 magnetic field model, t_iIs the time corresponding to the ith data, L is the length of each track data,and solving a preprocessing result after the first-order difference is solved for the subtraction model.

4. The method of claim 1,

step c, carrying out short-time Fourier transform on the differential data to construct a non-negative time-frequency amplitude matrix corresponding to the magnetic field data, wherein the non-negative time-frequency amplitude matrix comprises the following steps:

performing time-frequency transformation on the preprocessed magnetic field data through short-time Fourier transformation to obtain a time-frequency matrix comprising amplitude information and phase information;

and taking a time-frequency amplitude matrix formed by the amplitudes as an input matrix of the subsequent non-negative matrix decomposition.

5. The method of claim 1,

step d, specifically comprising: decomposing any one high-dimensional non-negative time-frequency amplitude matrix into two low-dimensional non-negative matrixes W and H by a non-negative matrix decomposition method, so that the product of W and H is approximate to the original high-dimensional matrix V, wherein W is called a base matrix, H is called a weight matrix, and the amplitude of a time-frequency distribution characteristic matrix and a time-frequency distribution characteristic matrix of data on the time domain is obtained;

the mathematical model of the non-negative matrix factorization method is as follows:

V≈WH (6)

wherein V belongs to a nonnegative matrix with the value of R being more than or equal to 0 and is MxN, W belongs to a nonnegative matrix with the value of R being more than or equal to 0 and is MxR, H belongs to a nonnegative matrix with the value of R being more than or equal to 0 and is RxN, and the selection of the decomposition characteristic number R meets the condition that R is less than MN/(M + N);

formulating an objective function of non-negative matrix factorization based on the K-L divergence, wherein the objective function based on the K-L divergence is as follows:

wherein V_ijRepresents the ith row, j column element of matrix V, [ WH]_ijThe ith row and the j column elements of the matrix are obtained by multiplying the base matrix W and the weight matrix H;

the weight matrix H is constrained by using a minimum determinant criterion, the volume of a vector expansion space of the weight matrix H is constrained to be minimum, and a solving formula of the volume of the vector expansion space of the weight matrix H is as follows:

then the K-L divergence non-negative matrix factorization objective function with the addition of the minimum determinant constraint is:

where det represents the determinant of the matrix, and α is the balance parameter;

optimizing and solving the objective function by adopting a multiplicative iteration method, obtaining an optimal solution by alternately updating variables, and calculating an updating formula by using a gradient descent method:

wherein let mu_wAnd mu_hRespectively as follows:

the final update iteration formula is:

wherein

Represents a dot product of matrix elements;

the normalized formula is:

setting an iteration stop condition that the target function is smaller than a set threshold Th and setting the upper limit of the iteration times to N_iterIf iteration N_iterForcibly stopping iteration if the value of the target function is still not less than the set threshold value after the next time;

and obtaining a base matrix W obtained by separating each track as follows:

weight matrix H

Comprises the following steps:

6. the method of claim 4,

e, selecting local influence components generated by the earthquake according to the energy ratio, and judging abnormal orbits of the components by an overrun method, wherein the method specifically comprises the following steps:

the ratio of the energy in a research area to the energy of the whole orbit in each row vector of a weight matrix H obtained by calculating non-negative time-frequency amplitude matrix decomposition is selected, and a component HH with the largest energy ratio is used as a local influence component generated by an earthquake;

and then selecting the root mean square of the energy average level of the reaction data as a threshold parameter, and judging abnormal points according to whether the value of the maximum component HH is greater than k times of the root mean square of the whole orbit, wherein the abnormal points are considered to be only present in the research area, and the corresponding zone bits show that the abnormal points are not the orbit caused by the known ionospheric activity and are abnormal orbits.

7. The method of claim 1,

f, accumulating the number of abnormal orbits every day, and detecting the earthquake abnormality according to the deviation degree of the abnormal orbits from the background fitting straight line specifically comprises the following steps:

the abnormal orbit every day and the number of all the orbits every day are normalized through dispersion standardization,

the cumulative normalization result of the number of the abnormal tracks per day is obtained by calculation as follows:

MN_ano＝[MN_ano(1)，MN_ano(2)，MN_ano(3)，...MN_ano(d)](26)

performing linear fitting on the accumulated normalization results of all the tracks every day through least squares to obtain a background fitting linear line;

the cumulative normalization result of calculating the number of all the tracks per day is as follows:

MN_nom＝[MN_nom(1)，MN_nom(2)，MN_nom(3)，...MN_nom(d)](27)

for MN_nomPerforming least square linear fitting, wherein a fitting straight line is as follows:

y＝a₀+a₁x (28)

wherein the parameter a₀And a₁The estimation solution formula is:

wherein x is [1, 2, 3.. d ]]D is the number of days, x_iIs the ith element in x; y is MN_nom，y_iIs MN_nomThe ith element in (1); x is the number of_iy_iIs the ith element in x and MN_nomThe ith element of (a); linear equation parameters are obtained through calculation

And

the earthquake change process is displayed through the degree of deviation of the daily abnormal orbit number cumulative curve from a background fitting straight line, and the deviation from linear fitting degree calculation formula is as follows:

wherein Y is_iIs MN_anoThe ith element in (1);

ρ is an index indicating how far the daily abnormal orbit accumulation curve deviates from the background fitting straight line, and a larger value indicates a larger abnormality of the data due to the influence of the earthquake.

Images